企业官网建设流程全解析

PyMC 协方差函数（Covariance Functions）完全指南：pymc.gp.cov 模块深度解析

【免费下载链接】pymcBayesian Modeling and Probabilistic Programming in Python项目地址: https://gitcode.com/GitHub_Trending/py/pymc

高斯过程（Gaussian Process, GP）的核心在于协方差函数（又称核函数）k(x, x')，它决定了函数先验的平滑性、周期性与各向异性等一切性质。本指南以 PyMC 官方 API 文档 docs/source/api/gp/cov.rst 为骨架，结合 pymc/gp/cov.py 源码与 tests/gp/test_cov.py 测试用例，系统讲解 PyMC 中pymc.gp.cov模块的全部内置协方差函数、input_dim与active_dims的参数语义、核的代数组合（加、乘、标量缩放、幂运算）以及如何在gp.Latent、gp.Marginal等 GP 实现中落地使用。读完本文，你将能够针对不同数据特征选择并组合出合适的协方差函数，构建可解释、可组合的贝叶斯 GP 模型。

协方差函数在 PyMC GP 中的角色

在贝叶斯建模中，有时我们关心的未知量不是标量或定长向量，而是一个连续函数。高斯过程为函数f(x)提供了先验分布：

f(x) ~ GP(m(x), k(x, x'))

其中m(x)是均值函数，k(x, x')是协方差函数。函数值被建模为多元正态分布的一次抽样，其协方差结构完全由k(x, x')决定。PyMC 的 GP 模块正是利用多元正态分布的边缘化与条件化性质来完成推断（边缘分布）与预测（条件分布）的。

PyMC 的 GP 具备"清晰语法 + 高度可组合"两大特点（参见官方指南 docs/source/guides/Gaussian_Processes.rst）：模块内置了大量预定义的协方差函数、均值函数和多种 GP 实现；更重要的是，GP 在 PyMC 中被当作可嵌入更大层次模型中的分布，而非只能独立使用的回归器。

pymc.gp.cov正是全部协方差函数的实现模块。根据 API 文档，该模块通过automodule:: pymc.gp.cov自动生成文档，公开的类包括：

Constant、WhiteNoise、ExpQuad、RatQuad、Exponential、Matern52、Matern32、Linear、Polynomial、Cosine、Periodic、WarpedInput、Gibbs、Coregion、ScaledCov、Kron。此外源码__all__（pymc/gp/cov.py）中还包含Matern12与WrappedPeriodic。

实例化与求值：先"参数化"，后"喂数据"

使用过 GPy 或 GPflow 的用户会对 PyMC 的语法感到熟悉：协方差函数在实例化时只完成参数化，并不接触输入数据；真正求值发生在之后以cov_func(X, Xs)方式调用时。这一"延迟绑定"设计正是为了支持核的组合构造——组合出的新核在被调用前无需关心具体输入。

所有协方差函数都继承自基类BaseCovariance（pymc/gp/cov.py），其统一的调用接口为：

cov_func(X, Xs=None, diag=False)

• X：训练输入（形状(n, input_dim)）；
• Xs：可选的预测输入；若为None则等价于Xs = X；
• diag：为True时只返回协方差矩阵的对角元，用于降低计算量（例如仅需方差时）。

调用内部会分发到两个抽象方法（BaseCovariance）：

• full(X, Xs)：计算完整的n x n（或n x m）协方差矩阵；
• diag(X)：仅计算对角元素。

以最简单的核为例：

import pymc as pm import numpy as np X = np.linspace(0, 1, 10)[:, None] # 10 个一维输入点 cov = pm.gp.cov.ExpQuad(input_dim=1, ls=0.1) K = cov(X).eval() # 10x10 协方差矩阵 Kd = cov(X, diag=True).eval() # 只取对角

input_dim 与 active_dims：告诉核"操作哪些列"

除Constant和WhiteNoise外，绝大多数协方差函数继承自Covariance基类（pymc/gp/cov.py），其构造函数统一接收两个关键参数：

之所以必须显式声明input_dim，是因为协方差函数在构造时尚未看到任何输入数据。active_dims允许同一个核只作用于输入的部分维度，从而可以构建"不同核作用于不同维度、再组合"的模型。

官方指南给出了经典示例（docs/source/guides/Gaussian_Processes.rst）：对包含三个预测变量的矩阵，构造一个只作用于第二、三列的指数二次核，并为每个维度配置独立的长度尺度：

ls = [2, 5] # 第二、三列各自的长度尺度 cov_func = pm.gp.cov.ExpQuad(input_dim=3, ls=ls, active_dims=[1, 2])

源码中_slice方法（pymc/gp/cov.py）会按active_dims对X和Xs做列切片，并在此处给出一个实用警告：若X的实际列数与input_dim不一致，会提示"只有input_dim列被用于计算协方差"，提醒用户确认参数意图。

需要特别说明的参数校验（Covariance.init）：active_dims中的任何值都不能超过input_dim，否则抛出ValueError。

内置协方差函数逐一详解

以下按照 API 文档列出的顺序，结合源码中的数学定义与参数说明逐个展开。

常数与噪声核

Constant（pymc/gp/cov.py）

k(x, x') = c

恒为常数c的协方差函数。它不继承Covariance（无input_dim/active_dims），常作为偏置项参与核组合。求值时通过_alloc生成全c矩阵。

WhiteNoise（pymc/gp/cov.py）

k(x, x') = σ²·I

白噪声协方差函数，sigma为噪声标准差。其对角为σ²；当给定预测输入Xs时，full返回n x m的全零矩阵——这正确反映了不同输入点之间白噪声不相关的性质。常被加在核上建模观测噪声：cov_func = ExpQuad(...) + WhiteNoise(sigma)。

平稳核（Stationary）

所有平稳核继承Stationary基类（pymc/gp/cov.py），其共性是核值只依赖输入点之间的距离，因此diag恒为 1.0。Stationary统一处理长度尺度参数，构造函数签名如下：

Stationary(input_dim, ls=None, ls_inv=None, active_dims=None)

• ls：长度尺度。input_dim > 1时可以是标量列表/数组（每个维度一个），也支持 PyMC 随机变量；input_dim == 1时为标量。
• ls_inv：逆长度尺度，即1 / ls。ls与ls_inv必须且只能提供一个，否则抛出ValueError（pymc/gp/cov.py）。在源码中ls_inv会被转换为ls = 1.0 / ls_inv。

距离计算方面，square_dist先按1/ls缩放输入再计算平方距离，并用pt.clip(sqd, 0.0, np.inf)消除数值误差导致的负值；euclidean_dist在平方距离上开方并加上1e-12数值稳定项（pymc/gp/cov.py）。

ExpQuad（指数二次核，又称平方指数 / RBF 核，pymc/gp/cov.py）

k(x, x') = exp( -(x - x')² / (2ℓ²) )

最常用的光滑核，处处无穷可微。full_from_distance实现为pt.exp(-0.5 * r2)。它还实现了功率谱密度（见后文"谱密度"一节）。

RatQuad（有理二次核，pymc/gp/cov.py）

k(x, x') = ( 1 + (x - x')² / (2αℓ²) )^(-α)

额外参数alpha（α）控制形状：它可以看作无穷多个不同长度尺度的平方指数核的尺度混合，α 越大越接近ExpQuad。源码注释给出了其作为平方指数核关于精度参数 λ ~ Gamma(α, αℓ²) 的混合表示。

Exponential（指数核，pymc/gp/cov.py）

k(x, x') = exp( -||x - x'|| / (2ℓ) )

与 ExpQuad 不同，它基于欧氏距离而非平方距离，对应的是奥恩斯坦-乌伦贝克（OU）过程，样本路径连续但不可导。

Matern52 / Matern32 / Matern12（Matérn 核族）

• Matern52（ν = 5/2，pymc/gp/cov.py）：

k(x, x') = (1 + √5·r/ℓ + 5r²/(3ℓ²)) · exp(-√5·r/ℓ)

• Matern32（ν = 3/2，pymc/gp/cov.py）：

k(x, x') = (1 + √3·r/ℓ) · exp(-√3·r/ℓ)

• Matern12（ν = 1/2，pymc/gp/cov.py）：

k(x, x') = exp(-r/ℓ)

其中r为（经长度尺度缩放后的）欧氏距离。Matérn 族通过 ν 参数刻画函数可微性：ν 越大函数越光滑，Matern52二阶可微、Matern32一阶可微、Matern12（即指数核的变体）不可微。实践中Matern52是平衡光滑性与数值稳定性的常用选择。

Cosine（余弦核，pymc/gp/cov.py）

k(x, x') = cos( 2π·||x - x'|| / ℓ² )

强周期振荡核，适用于周期性但非平滑衰减的信号。

Periodic（周期核，pymc/gp/cov.py）

k(x, x') = exp( -sin²(π|x - x'|/T) / (2ℓ²) )

需要额外参数period（周期T）。源码特别给出了一个易踩坑的说明：PyMC 的系数约定与常见定义不同，指数上是0.5而非常见的2，因此当你想复现标准定义时，初始化时需将长度尺度除以 2（pymc/gp/cov.py）。此外Periodic还实现了power_spectral_density_approx，用于 HSGP 低秩近似的系数计算（基于第一类修正贝塞尔函数I_j）。

非平稳核

Linear（线性核，pymc/gp/cov.py）

k(x, x') = (x - c)(x' - c)

c为偏移中心点。该核产生的函数是输入空间的线性函数，常用于趋势建模。

Polynomial（多项式核，pymc/gp/cov.py）

k(x, x') = [ (x - c)(x' - c) + offset ]^d

在Linear基础上增加次数d与常数项offset，实现方式为对线性核的结果取(linear + offset)^d。

输入变换类核

WarpedInput（pymc/gp/cov.py）

k(x, x') = k_base(w(x), w(x'))

用任意 PyTensor 函数warp_func对输入做非线性扭曲，再交给底层核cov_func计算。参数args用于向warp_func传递额外的标量或 PyMC 变量（内部由handle_args包装，见 pymc/gp/cov.py）。典型应用：将一维输入映射到高维特征空间，从而获得非平稳行为。

Gibbs（pymc/gp/cov.py）

k(x, x') = sqrt( 2ℓ(x)ℓ(x') / (ℓ²(x) + ℓ²(x')) ) · exp( -(x - x')² / (ℓ²(x) + ℓ²(x')) )

使用随输入变化的长度尺度函数lengthscale_func构造非平稳核。源码明确标注"仅在 1 维下测试过"：若active_dims长度大于 1 或input_dim != 1，会直接抛出NotImplementedError（pymc/gp/cov.py）。args同样用于传递额外参数。

ScaledCov（pymc/gp/cov.py）

k(x, x') = φ(x) · k_base(x, x') · φ(x')

用非负的缩放函数φ(x)（scaling_func，PyTensor 可调用对象）对基础核cov_func做点级缩放，实现对函数幅值的非平稳调制。源码中diag实现为cov_diag * φ(X)²，full实现为outer(φ(X), φ(Xs)) * k_base(X, Xs)。

多输出核：Coregion

Coregion（协区域化核，pymc/gp/cov.py）用于内禀/线性协区域化模型（ICM/LCM），其协方差矩阵为：

B = W·Wᵀ + diag(κ)

• W：形状(num_outputs, rank)的低秩矩阵，决定各输出之间的相关性；
• kappa：形状(num_outputs,)的向量，使各输出可独立变化；
• B：形状(num_outputs, num_outputs)的完整矩阵。

约束条件：(W, kappa)与B二者必须恰好提供一个；并且该核要求恰好激活一个维度（active_dims长度必须为 1），否则抛出ValueError（pymc/gp/cov.py）。调用时输入应为整数索引（输出编号），源码将其cast为int32后通过B[index, index2]查表取值。

结构核：Kron

Kron（Kronecker 积核，pymc/gp/cov.py）

与普通乘法（所有核共享同一份输入）不同，Kron 核先把输入按各因子的input_dim切分到各自子空间，再对各子空间分别求核矩阵，最后做 Kronecker 乘积：

Kron([cov1, cov2, ...])

其整体input_dim等于各因子input_dim之和，因子必须是协方差函数或其组合（数组不支持）。它通常配合gp.MarginalKron与gp.LatentKron使用，用于输入为笛卡尔积结构的数据（如网格数据），可显著加速矩阵运算。测试中也使用pymc.math.kronecker来验证其正确性（tests/gp/test_cov.py）。

核的代数组合：可组合性是 PyMC GP 的灵魂

协方差函数在 PyMC 中严格遵循核函数的代数规则——这一点被设计为语言级特性（BaseCovariance 中的运算符重载），用户可以像拼接乐高一样构造复杂核：

1. 两个核相加仍是核（Add，pymc/gp/cov.py）：

cov_func = pm.gp.cov.ExpQuad(input_dim=1, ls=1.0) + pm.gp.cov.Periodic(input_dim=1, period=0.5)

2. 两个核相乘仍是核（Prod，pymc/gp/cov.py）：

cov_func = pm.gp.cov.ExpQuad(input_dim=1, ls=1.0) * pm.gp.cov.Periodic(input_dim=1, period=0.5)

3. 核与标量的乘/加仍是核（标量会自动被包装为Constant）：

cov_func = eta**2 * pm.gp.cov.Matern32(input_dim=1, ls=1.0) # 幅值缩放 cov_func = pm.gp.cov.ExpQuad(input_dim=1, ls=1.0) + 1.0 # 加常数偏置

4. 幂运算（Exponentiated，pymc/gp/cov.py）：cov_func ** p要求底数必须是继承自Covariance的核（Constant/WhiteNoise会报TypeError），且指数必须是标量。

底层实现上，Combination基类（pymc/gp/cov.py）做了两件事：

• 自动推导元信息：组合核的input_dim取各因子之交集（所有因子必须一致，否则报错）；active_dims取各因子激活维度的并集（排序后）。
• 扁平化因子列表：Add套Add、Prod套Prod会被自动展平，避免嵌套过深。

_merge_factors_cov还负责统一处理diag=True时对数组/张量因子取对角（np.diag/pt.diag）的细节，保证对角计算的正确性——这一点被测试用例反复验证（tests/gp/test_cov.py）。

从测试可以看出这些运算的实际数值行为（tests/gp/test_cov.py）：ExpQuad(1, 0.1) + 1的非对角元为1.53940，正是ExpQuad核值0.53940加上常数1；ExpQuad(1, 0.1) + WhiteNoise(sigma=1)的对角元为2（核对角 1 加噪声方差 1）、非对角元保持0.53940。

功率谱密度：从频域理解平稳核

对于平稳核，PyMC 还实现了power_spectral_density(omega)（Stationary），为 HSGP（Hilbert Space Gaussian Process）等谱方法提供理论基础：

• ExpQuad的谱密度为高斯形式（pymc/gp/cov.py）；
• RatQuad的谱密度涉及第二类修正贝塞尔函数K_ν(z)，并处理了z = 0处的奇异性（pymc/gp/cov.py）；
• Matern32/Matern52的谱密度为有理函数形式（pymc/gp/cov.py）。

组合核的谱密度遵循简单规则：Add的谱密度等于各因子谱密度之和；Prod仅当因子中协方差函数不多于一个时才能计算（两个协方差函数相乘的谱没有解析形式，会抛出NotImplementedError，见 pymc/gp/cov.py）。此外，进行谱密度计算时要求求和项中所有核的active_dims完全一致（pymc/gp/cov.py），非平稳核则会收到明确的ValueError提示。

在 GP 模型中使用协方差函数

协方差函数本身只是"零件"，真正发挥威力的是与 GP 实现的配合。PyMC 提供gp.Latent、gp.Marginal等多种实现，统一的使用模式为：先用均值函数与协方差函数实例化 GP 对象，再调用prior、marginal_likelihood或conditional方法构造 PyMC 随机变量。

Latent GP 的完整流程（docs/source/guides/Gaussian_Processes.rst）：

import pymc as pm # 1. 定义核与均值 cov_func = pm.gp.cov.ExpQuad(input_dim=1, ls=1.0) mean_func = pm.gp.mean.Zero() # 2. 实例化 GP gp = pm.gp.Latent(mean_func, cov_func) with pm.Model() as model: # 3. 在观测点 X 上建立函数先验 f f = gp.prior("f", X) # 4. 预测点的条件分布 f* f_star = gp.conditional("f_star", X_star)

prior的第一个参数是随机变量名，第二个是函数输入X（通常来自数据，也可以是 PyMC 随机变量——若输入是张量/随机变量，则必须显式给出shape）。gp.Marginal类没有prior方法，改用marginal_likelihood，需要额外提供观测数据与噪声等参数。

可加 GP（Additive GP）是 PyMC 的另一大特性（docs/source/guides/Gaussian_Processes.rst）：GP 对象之间可以直接相加，从而把复杂函数分解为多个独立成分并分别推断：

with pm.Model() as model: gp1 = pm.gp.Marginal(mean_func1, cov_func1) # 长程趋势 gp2 = pm.gp.Marginal(mean_func2, cov_func2) # 周期成分 gp = gp1 + gp2 # f1 + f2 f = gp.marginal_likelihood("f", X, y, noise) idata = pm.sample(1000)

注意：相加的两个 GP 对象类型必须一致（gp.Marginal不能与gp.Latent相加）。若想得到某个成分（如gp2）的条件分布，需要在conditional中以given={"X": X, "y": y, "noise": noise, "gp": gp}字典形式显式传入缓存的参数（docs/source/guides/Gaussian_Processes.rst）：

with model: f2_star = gp2.conditional("f2_star", X_star, given={"X": X, "y": y, "noise": noise, "gp": gp}) f_star = gp.conditional("f_star", X_star) # 整体条件分布无需 given

整体gp的条件分布无需given，因为参数已在marginal_likelihood调用时被缓存；而gp1、gp2自身从未调用过marginal_likelihood，必须显式提供。

测试驱动的工程细节与注意事项

仓库中的 tests/gp/test_cov.py（共 948 行）是对pymc.gp.cov行为的最权威验证，值得注意的工程细节包括：

• 对角一致性：所有测试都校验cov(X, diag=True)与np.diag(cov(X))严格一致（如 tests/gp/test_cov.py），这是BaseCovariance接口设计正确性的基本保证。
• 标量/矩阵的左右结合：a + cov、cov + a、M + cov、cov + M均被支持（通过__radd__/__rmul__与__array_wrap__），但形状非法的组合（如三维数组）会抛出ValueError: cannot combine...（tests/gp/test_cov.py）。
• 组合核的元信息合并：不同input_dim的核相加会报错；active_dims取并集，这保证组合核始终能正确切片输入。

实战中的几条核心建议：

1. 长度尺度选择：ls可设为标量（各向同性）或按维度设置的数组（自动相关性确定，ARD）；当维度较多时，优先用ls_inv配合正约束先验（如HalfNormal、Gamma）进行推断。
2. 噪声与核分离：观测噪声建议显式使用WhiteNoise或gp.Marginal的noise参数，保持核本身干净，便于解释。
3. 周期数据：优先选用Periodic（注意其 0.5 缩放约定）；若底层核需要更强的可微性，可用WrappedPeriodic将任意Stationary核"卷"成周期核（pymc/gp/cov.py）。
4. 大规模数据：网格化/可分解结构数据用Kron配合gp.MarginalKron/gp.LatentKron；否则可结合Periodic.power_spectral_density_approx与 HSGP 近似（见 pymc/gp/hsgp_approx.py）降低计算负担。

小结

pymc.gp.cov为贝叶斯高斯过程建模提供了完整、可组合、经过充分测试的协方差函数工具箱：从经典的ExpQuad、Matérn 族、周期核，到非平稳的WarpedInput、Gibbs、ScaledCov，再到多输出Coregion与结构化的Kron，全部遵循统一的BaseCovariance调用契约，并支持加、乘、标量缩放与幂运算的代数组合。配合input_dim/active_dims的维度控制与ls/ls_inv的长度尺度参数化，你可以在 pymc/gp/gp.py 提供的Latent、Marginal、MarginalKron等实现中自由搭建从简单回归到可加成分模型的各种 GP 模型。

【免费下载链接】pymcBayesian Modeling and Probabilistic Programming in Python项目地址: https://gitcode.com/GitHub_Trending/py/pymc