去年我接过一个数据拟合的活训练点一共只有13个输入输出之间明显是非线性关系。用多项式回归跑了几轮曲线倒是贴着数据走可对方紧接着问了句这个区间里的预测值误差带是多少我当场卡壳——传统回归给的是点估计预测区间得额外假设误差分布才能算出来而且误差方差在所有位置都一样明显不合理。后来把工具换成高斯过程才真正把“曲线拟合”这件事做完整了它不只是给出一条拟合线连每个位置的置信区间都是一起带出来的。这篇文章我就把高斯过程、高斯过程采样、高斯过程回归一次讲清楚附上可以直接跑的MATLAB代码讲明白每一步背后的数学逻辑和我在实际项目中踩过的坑。1. 为什么函数拟合需要一套“概率世界观”如果你只用最小二乘拟合隐含的假设是存在一个确定性的函数 f(x)我们目测或根据经验给它定一个参数形式然后调整参数让它尽量贴合数据。这个思路在数据量大、函数形态明确时没问题但在小样本、高噪声的场景下会很别扭。1.1 传统点估计的根本局限最小二乘给的是 f(x) 的一个点估计比如二次多项式拟合出来一堆系数代入 x 就得到 y 的预测值。但真实世界不是这样运作的。你在某个位置测到 y3.2另一天再测可能变成 3.5噪声本身就是随机的。更麻烦的是拟合出的预测值在不同位置的可靠程度是不一样的有数据点的地方预测可信度高两个数据点中间是“盲区”数据点之外几乎是外推。传统回归对这些信息几乎没有表达能力。做过实验设计的人都懂采样点怎么布、传感器布在哪都要考虑预测精度的空间分布。如果模型能告诉我们“哪个区域不确定性高”就可以引导下一步往哪里补数据点。这就是主动学习和贝叶斯优化的核心思想而高斯过程天然适合干这件事因为它把函数本身当作一个随机过程来对待。1.2 把函数当作随机变量的集合高斯过程严格的定义是一个随机过程其中任意有限个输入点上的函数值都服从一个多元高斯分布。听起来有点绕拆开说就是——假设我们有一组输入点 x₁, x₂, ..., xₙ它们对应的函数值 f(x₁), f(x₂), ..., f(xₙ) 不是一系列独立的随机数而是联合服从一个 n 维高斯分布。这个 n 维高斯分布的均值向量描述了函数在每一点的“中心趋势”协方差矩阵描述了不同点之间函数值的统计相关性。这种“概率世界观”最大的好处是预测的时候不光给出均值还给出方差方差大小直接反映模型在那个位置的自信程度。高斯过程回归预测结果天然是一个高斯分布后验均值和后验方差可以写成一个闭环公式不需要做采样逼近。这也是它在小样本机器学习任务里始终没被淘汰的原因。2. 从多元高斯到无限维函数先验分布的构建要用高斯过程描述函数必须完整定义它的先验。这一步对初学者来说最容易懵因为一下子从“多元高斯分布”跳到了“无限维”的抽象概念。我建议先回顾多元高斯的关键性质再慢慢把维度推高。2.1 由均值函数和核函数确定的高斯过程多元高斯分布由两个要素决定均值向量 μ 和协方差矩阵 Σ。推广到函数空间后高斯过程也由两个要素决定均值函数 m(x) 和协方差函数也叫核函数k(x, x)。对于任意一组点 x₁...xₙ函数值的先验是p(f) N(f | m, K)其中 f [f(x₁), ..., f(xₙ)]ᵀm 的每个分量是 m(xᵢ)K 的每个元素 Kᵢⱼ k(xᵢ, xⱼ)。多数时候均值函数取常数 0因为核函数已经能表达足够的结构均值函数强行搞得很复杂反而容易引入错误先验。核函数才是重头戏。2.2 核矩阵如何编码相似性核函数的直观含义是两个输入点距离越近对应的函数值越相关。一个典型的平方指数核RBF核也叫SE核写出来是k(xᵢ, xⱼ) σ_f² · exp(−‖xᵢ − xⱼ‖² / (2l²))这里的 σ_f² 叫信号方差控制函数的整体波动幅度l 叫长度尺度控制函数在输入空间上的变化快慢。l 越小函数扭曲越厉害两点但凡离得稍远相关性就迅速衰减l 越大函数越平缓两条曲线在长距离上仍然像。拿温度场拟合来说两个观测点距离 100 米温度变化通常不会太大长度尺度 l 就对应了这种物理上的“空间相关长度”。这不是玄学而是核函数把领域知识编码进了统计模型的结构里。在程序里核函数算出来的 K 是一个 n×n 对称正定矩阵。正定性是高斯过程一切推导的基础因为协方差矩阵必须合法Cholesky 分解才能稳定执行。这个矩阵也是后面所有计算的核心后验均值、方差、边际似然全都建立在它的逆或分解之上。3. 核函数决定函数的“性格”常见核与选择思路核函数是高斯过程模型里最重要的设计决策。我见过不少工程上跑出来的预测结果很离谱最后排查发现不是代码写错而是核函数跟数据本身的规律不匹配。3.1 平方指数核与 Matérn 核的真实差别平方指数核RBF是高斯完整且无限可微的意味着它建模出来的函数非常“光滑”在实际数据里偶尔会过于理想——很多真实物理过程的函数值变化并没有那么平缓而是带有某种粗糙度。Matérn 核通过参数 ν 控制函数的光滑程度。ν 3/2 时表达式是k(xᵢ, xⱼ) σ_f² · (1 √3·‖xᵢ−xⱼ‖/l) · exp(−√3·‖xᵢ−xⱼ‖/l)ν 5/2 时是k(xᵢ, xⱼ) σ_f² · (1 √5·‖xᵢ−xⱼ‖/l 5·‖xᵢ−xⱼ‖²/(3l²)) · exp(−√5·‖xᵢ−xⱼ‖/l)它的好处是对函数光滑程度的假设没有那么强处理工程数据时经常比 RBF 更稳。如果你不确定数据到底有多光滑优先试 Matérn-5/2很多项目中它的默认表现就不错。3.2 周期核、线性核与核函数的加减组合数据有明显周期性时比如昼夜温度变化、设备运行周期振动就应该用周期核k(xᵢ, xⱼ) σ_f² · exp(−2·sin²(π·‖xᵢ−xⱼ‖/p) / l²)其中 p 是周期长度。核函数还可以相加、相乘形成新的核。相加表示不同效应叠加比如长周期趋势 短周期波动相乘表示两个因素同时满足时才相关比如空间位置和时间的交互效应。不要一开始就把核函数组合得很复杂工程上先跑出基线结果再逐步增加结构否则容易过度参数化。3.3 选择建议先看数据形状再定核我给一个很个人的经验拿到数据后先画散点图。如果函数曲线肉眼可见地平滑用 RBF如果表面看起来有细碎的波动优先 Matérn-3/2如果存在明显重复趋势用周期核如果整体是线性的线性核或者“线性核 RBF”组合会更可靠。核函数不是越复杂越好它本质上是先验假设错误的假设会带来不靠谱的预测区间。4. 先验采样怎么做Cholesky分解与MATLAB代码高斯过程先验本身已经定义了函数分布但光看公式很难建立直觉。最好的理解方式是做一件事从先验里“抽函数”出来画在图上。这就是高斯过程采样。4.1 为什么用 Cholesky 分解而不是直接求逆要从高斯分布 N(m, K) 里抽一个向量标准做法是先做 Cholesky 分解 K L Lᵀ再生成一组标准正态随机向量 z然后算 f m L z。为什么不用矩阵开方或者直接对 K 求逆因为 Cholesky 分解是三角分解数值上比求逆稳定得多速度也快。K 是对称正定矩阵这个分解一定能做成功前提是不出现病态。如果 K 里有数值误差导致微小负特征值Cholesky 通常会报错这反而是在提醒你数值出了问题比求逆后悄悄给一个错误结果好得多。4.2 MATLAB采样代码与样本可视化在 MATLAB 里跑一次先验采样非常直接% 采样点 X linspace(0, 1, 30); % 超参数 sigma_f 1.0; l 0.2; % 计算核矩阵 K rbf_kernel(X, X, sigma_f, l); % Cholesky分解 L chol(K, lower); % 抽取一组函数值 z randn(30, 1); f L * z; % 画图 plot(X, f, LineWidth, 1.5);其中核矩阵函数function K rbf_kernel(X1, X2, sigma_f, l) n1 size(X1, 1); n2 size(X2, 1); X1sq sum(X1.^2, 2); X2sq sum(X2.^2, 2); sqdist repmat(X1sq, 1, n2) repmat(X2sq, n1, 1) - 2 * X1 * X2; K sigma_f^2 * exp(-sqdist / (2 * l^2)); end注意这里我用了 repmat 加平移扩展的技巧避免双重循环。数据量几百个点以内时这种写法足够高效代码还短。4.3 多画几次你就会理解超参数多运行几次上面的代码把每次的曲线叠加到同一张图上你就能直观看到先验分布描述的函数集合。长度尺度 l 0.2 时曲线在 0 到 1 区间上上下下摆动很多次改成 l 0.8 时曲线会变得非常温柔两条邻近点的函数值相差不多。先验采样看起来只是可视化手段但它有实际用途在贝叶斯优化里需要从函数先验里生成初始候选点在蒙特卡洛积分时从高斯过程先验抽取函数再计算函数期望有些GPU计算场景下还会用固定数量的随机特征去近似核函数本质上也是先验采样思想。5. 含噪观测下的回归后验均值和方差推导高斯过程回归的目标是在给定观测数据后计算任意预测点上的函数值分布。这一步是整个方法的核心。5.1 观测模型的建立假设观测是潜在函数值加高斯噪声yᵢ f(xᵢ) ε, ε ~ N(0, σₙ²)观测点的联合分布是p(y) N(y | 0, K(X, X) σₙ²I)这里 K(X, X) 是训练点之间的核矩阵再加上噪声方差构成观测协方差。5.2 联合分布与条件分布的推导对于预测点 X*m 个点训练输出 y 和预测函数值 f* 的联合分布是[y, f*]ᵀ ~ N(0, [[K σₙ²I, K*], [K*ᵀ, K**]])其中 K* K(X, X*) 大小为 n×mK** K(X*, X*) 大小为 m×m。利用多元高斯条件分布的公式条件后验均值是μ* K*ᵀ · (K σₙ²I)⁻¹ · y条件后验协方差是Σ* K** − Kᵀ · (K σₙ²I)⁻¹ · K这个公式不需要采样就能得到解析解这就是高斯过程回归比贝叶斯神经网络更让人舒服的地方。5.3 从公式看模型行为后验均值 μ* 有一种很直观的解释它在每个预测点上计算训练数据所有观测值的加权线性组合权重由核函数 k(xᵢ, x*) 决定。离预测点越近的训练点权重大离得远的权重小。这和RBF插值的加权平均思想是一脉相承的但高斯过程在此基础上又给出了方差。后验方差 Σ* 等于预测点的先验方差 K** 减去观测数据带来的信息量。越靠近训练点的区域后验方差越小离所有训练点都很远的区域后验方差会退回接近先验方差。这是高斯过程最迷人的地方——预测不确定性与空间位置直接挂钩。抽样工程里布置传感器可以先算一下哪些区域的预测方差最大然后专门去那里补点这就是主动学习的标准玩法。6. 超参数优化的实质在边际似然面上寻优高斯过程不是“无参模型”核函数里藏着超参数σ_f、l 以及观测噪声 σₙ。这些参数不能靠拍脑袋盯着一张曲线调标准做法是最大化对数边际似然。6.1 负对数边际似然的分解与含义观测数据 y 的边际似然是对潜在函数积分后的结果log p(y | X, θ) −½ yᵀ · (K σₙ²I)⁻¹ · y − ½ log|K σₙ²I| − n/2 · log(2π)右边第一项是“数据拟合项”衡量模型解释数据的程度数据拟合得越好这一项越小。第二项是“复杂度惩罚项”行列式的对数会随核矩阵特征值变大而变大相当于对模型复杂度做自然惩罚。第三项是归一化常数与超参数无关。高斯过程不会被复杂度惩罚项一巴掌拍死而是自动在拟合度和复杂度之间找均衡。我在调参时经常说别只盯着预测误差调边际似然里已经内置了奥卡姆剃刀。6.2 优化技巧对数域参数与多起点超参数必须为正所以实际优化的是 log(σ_f)、log(l)、log(σₙ)保证指数回退后永远是正数。目标函数定义为负对数边际似然然后用 fminsearch 或者 fminunc 求解。我一开始犯过错误在原始空间里优化迭代几次长度尺度变成负数核矩阵彻底崩掉。换成对数域以后再也没出过这种问题。另一个容易踩的坑是局部最优。边际似然面在长度尺度小的地方经常有高方差区域搜索很容易陷进去。稳妥起见用多个随机起点分别优化比较最终负对数边际似然值选最小的那组超参数。MATLAB 里负对数边际似然可以这样写function nll neg_log_marginal_likelihood(theta, X, y) sigma_f exp(theta(1)); l exp(theta(2)); sigma_n exp(theta(3)); K rbf_kernel(X, X, sigma_f, l) sigma_n^2 * eye(size(X,1)); L chol(K, lower); alpha L \ (L \ y); nll 0.5 * y * alpha sum(log(diag(L))) 0.5 * size(X,1) * log(2*pi); end用 chol 返回 L 之后sum(log(diag(L))) 等价于 0.5·log|K|都是数值稳定的写法。注意不要直接用 log(det(K))矩阵接近奇异时 det 会下溢或上溢完全不可靠。7. 完整的MATLAB实现与结果解读把前面所有环节串起来给一个小而完整的回归示例。假设我们的潜在函数是 f(x) sin(2πx)训练点 10 个噪声标准差 0.1。7.1 主程序与预测% 生成带噪声的观测数据 rng(42); X linspace(0, 1, 10); y sin(2 * pi * X) 0.1 * randn(10, 1); % 预测点 Xs linspace(-0.5, 1.5, 100); % 初始化超参数对数域 theta_init log([1.0, 0.2, 0.1]); % 优化超参数 theta_opt fminsearch((t) neg_log_marginal_likelihood(t, X, y), theta_init); % 用优化后的超参数计算后验 sigma_f exp(theta_opt(1)); l exp(theta_opt(2)); sigma_n exp(theta_opt(3)); K rbf_kernel(X, X, sigma_f, l) sigma_n^2 * eye(size(X,1)); L chol(K, lower); alpha L \ (L \ y); K_star rbf_kernel(X, Xs, sigma_f, l); mu_star K_star * alpha; % 计算后验方差 v L \ K_star; var_star diag(rbf_kernel(Xs, Xs, sigma_f, l)) - sum(v.^2, 1); std_star sqrt(var_star); % 画图 figure; hold on; fill([Xs; flipud(Xs)], [mu_star - 1.96*std_star; flipud(mu_star 1.96*std_star)], ... [0.8, 0.9, 1.0], EdgeColor, none); plot(Xs, mu_star, b-, LineWidth, 2); plot(X, y, ro, MarkerSize, 8, MarkerFaceColor, r); xlabel(x); ylabel(y); legend(95%置信区间, 后验均值, 观测数据, Location, best);7.2 看结果时要注意的几个数字跑完之后你会看到优化得到的长度尺度 l 大约在 0.15 到 0.25 之间对应数据波动的空间频率。后验均值在 x0.5 附近能较准确地描述正弦波峰但在两个边界 x0 和 x1 的区域置信区间快速张大因为那里没有观测数据模型在向外推时表现出了诚实的不确定性。95%置信区间宽度由 1.96·std_star 决定整个区间带在训练数据点附近明显收紧在两个相邻训练点中间略微变宽。这种区间宽度的空间变化是高斯过程相对普通回归的卖点。7.3 训练点数量对计算量的影响必须强调核矩阵求逆或 Cholesky 分解的计算复杂度是 O(n³)n 是训练样本数。300 个点以内基本秒开3000 个点就开始吃力30000 个点不可行。如果数据量上来了标准高斯过程要上稀疏近似或者随机特征扩展。你如果只是做小样本拟合、贝叶斯优化代理模型、实验设计这个复杂度完全够用。8. 实际使用中的数值稳定性坑与优化建议高斯过程看起来公式干净但工程实现有大量细节我把真正踩过的坑和习惯做法一次性列出来。8.1 矩阵病态与Jitter处理核矩阵对角线加上噪声方差 σₙ²之后理论上正定但实际中两个非常接近的观测点会让矩阵几乎奇异。此时 chol 函数直接报错。常规做法是在 K 的对角线上加一个很小的常数比如 1e-10甚至 1e-6这个操作叫加 jitter。加了之后对预测均值几乎没有影响但能让矩阵稳定分解。注意不要为了省事把 jitter 设成过大值比如和信号方差同量级那样会人为增大噪声压低后验方差的表达能力。工程经验是先设为 1e-8如果 chol 报错再逐步加大。8.2 特征值分解判断问题位置如果加了 jitter 还是报错先检查核函数的距离计算是否有问题。我之前写过一版 rbf_kernel距离矩阵里有很多 NaN结果核矩阵对角线全是 NaNCholesky 直接崩溃。打印 K 的前几行总能第一时间发现问题。如果问题确实来自病态用 eig 快速查看最小特征值有多大。最小特征值接近机器精度时意味着某些超参数下核函数把两个点定义得过于接近可以考虑修改长度尺度初始化值或者重新检查输入数据。8.3 数据归一化几乎必须做高斯过程对输入输出尺度很敏感。如果 x 范围在 0 到 1e6长度尺度 l 会被迫学习到百万量级如果 y 的量纲差距太大σ_f 的优化也会很不稳。实际项目里我会把输入和输出都归一化到零均值、单位方差训练完成后再反变换回原始量纲。有人会问归一化是不是改变了核函数的解释性确实会但你换来的是优化稳定性。只要把长度尺度等参数归一化后的值记录下来反归一化时可以换算出原始尺度下的相关长度。8.4 复杂度惩罚与超参数初始化的关系负对数边际似然函数非常容易在长度尺度很小区间产生多个局部极小点。别只用一组初始值我会随机抽 5 到 10 组初始值分别优化保留最优的一组。这样做的好处不仅仅是找到更好的参数还能顺带观察多组解之间后验曲线的差异如果差异很大说明数据对超参数不够敏感这时候对预测区间的解读要格外谨慎。8.5 预测点数量与内存控制预测时如果用一个特别密的网格K** 会变得巨大。设计实验时我会先只计算方差的对角线而不算了整个 K**因为大多数预测场景只需要单点方差。在上面的代码里用到的是 var_star diag(...) - sum(v.^2, 1)这已经是省内存的写法。如果连预测点也上百万就分段做预测避免一次性把 K_star 矩阵全部放进内存。最后聊聊我的实际体会高斯过程这东西最初学的时候容易被数学绕晕但真正上手发现核心就三件事选核函数、算核矩阵、解一个线性系统。在实际工程里它最出彩的地方不是给你一条好看的拟合曲线而是告诉你哪里预测可信、哪里不可信。做设备标定的时候我用它来确定下一组实验点在哪个位置最划算做仿真代理模型的时候我用它来代替每轮都跑重型仿真。每次跑通都有一种“这钱省得真值”的感觉。如果在你的项目里数据量不大但预测置信区间非常重要高斯过程值得作为第一选择。
去年我接过一个数据拟合的活训练点一共只有13个输入输出之间明显是非线性关系。用多项式回归跑了几轮曲线倒是贴着数据走可对方紧接着问了句这个区间里的预测值误差带是多少我当场卡壳——传统回归给的是点估计预测区间得额外假设误差分布才能算出来而且误差方差在所有位置都一样明显不合理。后来把工具换成高斯过程才真正把“曲线拟合”这件事做完整了它不只是给出一条拟合线连每个位置的置信区间都是一起带出来的。这篇文章我就把高斯过程、高斯过程采样、高斯过程回归一次讲清楚附上可以直接跑的MATLAB代码讲明白每一步背后的数学逻辑和我在实际项目中踩过的坑。1. 为什么函数拟合需要一套“概率世界观”如果你只用最小二乘拟合隐含的假设是存在一个确定性的函数 f(x)我们目测或根据经验给它定一个参数形式然后调整参数让它尽量贴合数据。这个思路在数据量大、函数形态明确时没问题但在小样本、高噪声的场景下会很别扭。1.1 传统点估计的根本局限最小二乘给的是 f(x) 的一个点估计比如二次多项式拟合出来一堆系数代入 x 就得到 y 的预测值。但真实世界不是这样运作的。你在某个位置测到 y3.2另一天再测可能变成 3.5噪声本身就是随机的。更麻烦的是拟合出的预测值在不同位置的可靠程度是不一样的有数据点的地方预测可信度高两个数据点中间是“盲区”数据点之外几乎是外推。传统回归对这些信息几乎没有表达能力。做过实验设计的人都懂采样点怎么布、传感器布在哪都要考虑预测精度的空间分布。如果模型能告诉我们“哪个区域不确定性高”就可以引导下一步往哪里补数据点。这就是主动学习和贝叶斯优化的核心思想而高斯过程天然适合干这件事因为它把函数本身当作一个随机过程来对待。1.2 把函数当作随机变量的集合高斯过程严格的定义是一个随机过程其中任意有限个输入点上的函数值都服从一个多元高斯分布。听起来有点绕拆开说就是——假设我们有一组输入点 x₁, x₂, ..., xₙ它们对应的函数值 f(x₁), f(x₂), ..., f(xₙ) 不是一系列独立的随机数而是联合服从一个 n 维高斯分布。这个 n 维高斯分布的均值向量描述了函数在每一点的“中心趋势”协方差矩阵描述了不同点之间函数值的统计相关性。这种“概率世界观”最大的好处是预测的时候不光给出均值还给出方差方差大小直接反映模型在那个位置的自信程度。高斯过程回归预测结果天然是一个高斯分布后验均值和后验方差可以写成一个闭环公式不需要做采样逼近。这也是它在小样本机器学习任务里始终没被淘汰的原因。2. 从多元高斯到无限维函数先验分布的构建要用高斯过程描述函数必须完整定义它的先验。这一步对初学者来说最容易懵因为一下子从“多元高斯分布”跳到了“无限维”的抽象概念。我建议先回顾多元高斯的关键性质再慢慢把维度推高。2.1 由均值函数和核函数确定的高斯过程多元高斯分布由两个要素决定均值向量 μ 和协方差矩阵 Σ。推广到函数空间后高斯过程也由两个要素决定均值函数 m(x) 和协方差函数也叫核函数k(x, x)。对于任意一组点 x₁...xₙ函数值的先验是p(f) N(f | m, K)其中 f [f(x₁), ..., f(xₙ)]ᵀm 的每个分量是 m(xᵢ)K 的每个元素 Kᵢⱼ k(xᵢ, xⱼ)。多数时候均值函数取常数 0因为核函数已经能表达足够的结构均值函数强行搞得很复杂反而容易引入错误先验。核函数才是重头戏。2.2 核矩阵如何编码相似性核函数的直观含义是两个输入点距离越近对应的函数值越相关。一个典型的平方指数核RBF核也叫SE核写出来是k(xᵢ, xⱼ) σ_f² · exp(−‖xᵢ − xⱼ‖² / (2l²))这里的 σ_f² 叫信号方差控制函数的整体波动幅度l 叫长度尺度控制函数在输入空间上的变化快慢。l 越小函数扭曲越厉害两点但凡离得稍远相关性就迅速衰减l 越大函数越平缓两条曲线在长距离上仍然像。拿温度场拟合来说两个观测点距离 100 米温度变化通常不会太大长度尺度 l 就对应了这种物理上的“空间相关长度”。这不是玄学而是核函数把领域知识编码进了统计模型的结构里。在程序里核函数算出来的 K 是一个 n×n 对称正定矩阵。正定性是高斯过程一切推导的基础因为协方差矩阵必须合法Cholesky 分解才能稳定执行。这个矩阵也是后面所有计算的核心后验均值、方差、边际似然全都建立在它的逆或分解之上。3. 核函数决定函数的“性格”常见核与选择思路核函数是高斯过程模型里最重要的设计决策。我见过不少工程上跑出来的预测结果很离谱最后排查发现不是代码写错而是核函数跟数据本身的规律不匹配。3.1 平方指数核与 Matérn 核的真实差别平方指数核RBF是高斯完整且无限可微的意味着它建模出来的函数非常“光滑”在实际数据里偶尔会过于理想——很多真实物理过程的函数值变化并没有那么平缓而是带有某种粗糙度。Matérn 核通过参数 ν 控制函数的光滑程度。ν 3/2 时表达式是k(xᵢ, xⱼ) σ_f² · (1 √3·‖xᵢ−xⱼ‖/l) · exp(−√3·‖xᵢ−xⱼ‖/l)ν 5/2 时是k(xᵢ, xⱼ) σ_f² · (1 √5·‖xᵢ−xⱼ‖/l 5·‖xᵢ−xⱼ‖²/(3l²)) · exp(−√5·‖xᵢ−xⱼ‖/l)它的好处是对函数光滑程度的假设没有那么强处理工程数据时经常比 RBF 更稳。如果你不确定数据到底有多光滑优先试 Matérn-5/2很多项目中它的默认表现就不错。3.2 周期核、线性核与核函数的加减组合数据有明显周期性时比如昼夜温度变化、设备运行周期振动就应该用周期核k(xᵢ, xⱼ) σ_f² · exp(−2·sin²(π·‖xᵢ−xⱼ‖/p) / l²)其中 p 是周期长度。核函数还可以相加、相乘形成新的核。相加表示不同效应叠加比如长周期趋势 短周期波动相乘表示两个因素同时满足时才相关比如空间位置和时间的交互效应。不要一开始就把核函数组合得很复杂工程上先跑出基线结果再逐步增加结构否则容易过度参数化。3.3 选择建议先看数据形状再定核我给一个很个人的经验拿到数据后先画散点图。如果函数曲线肉眼可见地平滑用 RBF如果表面看起来有细碎的波动优先 Matérn-3/2如果存在明显重复趋势用周期核如果整体是线性的线性核或者“线性核 RBF”组合会更可靠。核函数不是越复杂越好它本质上是先验假设错误的假设会带来不靠谱的预测区间。4. 先验采样怎么做Cholesky分解与MATLAB代码高斯过程先验本身已经定义了函数分布但光看公式很难建立直觉。最好的理解方式是做一件事从先验里“抽函数”出来画在图上。这就是高斯过程采样。4.1 为什么用 Cholesky 分解而不是直接求逆要从高斯分布 N(m, K) 里抽一个向量标准做法是先做 Cholesky 分解 K L Lᵀ再生成一组标准正态随机向量 z然后算 f m L z。为什么不用矩阵开方或者直接对 K 求逆因为 Cholesky 分解是三角分解数值上比求逆稳定得多速度也快。K 是对称正定矩阵这个分解一定能做成功前提是不出现病态。如果 K 里有数值误差导致微小负特征值Cholesky 通常会报错这反而是在提醒你数值出了问题比求逆后悄悄给一个错误结果好得多。4.2 MATLAB采样代码与样本可视化在 MATLAB 里跑一次先验采样非常直接% 采样点 X linspace(0, 1, 30); % 超参数 sigma_f 1.0; l 0.2; % 计算核矩阵 K rbf_kernel(X, X, sigma_f, l); % Cholesky分解 L chol(K, lower); % 抽取一组函数值 z randn(30, 1); f L * z; % 画图 plot(X, f, LineWidth, 1.5);其中核矩阵函数function K rbf_kernel(X1, X2, sigma_f, l) n1 size(X1, 1); n2 size(X2, 1); X1sq sum(X1.^2, 2); X2sq sum(X2.^2, 2); sqdist repmat(X1sq, 1, n2) repmat(X2sq, n1, 1) - 2 * X1 * X2; K sigma_f^2 * exp(-sqdist / (2 * l^2)); end注意这里我用了 repmat 加平移扩展的技巧避免双重循环。数据量几百个点以内时这种写法足够高效代码还短。4.3 多画几次你就会理解超参数多运行几次上面的代码把每次的曲线叠加到同一张图上你就能直观看到先验分布描述的函数集合。长度尺度 l 0.2 时曲线在 0 到 1 区间上上下下摆动很多次改成 l 0.8 时曲线会变得非常温柔两条邻近点的函数值相差不多。先验采样看起来只是可视化手段但它有实际用途在贝叶斯优化里需要从函数先验里生成初始候选点在蒙特卡洛积分时从高斯过程先验抽取函数再计算函数期望有些GPU计算场景下还会用固定数量的随机特征去近似核函数本质上也是先验采样思想。5. 含噪观测下的回归后验均值和方差推导高斯过程回归的目标是在给定观测数据后计算任意预测点上的函数值分布。这一步是整个方法的核心。5.1 观测模型的建立假设观测是潜在函数值加高斯噪声yᵢ f(xᵢ) ε, ε ~ N(0, σₙ²)观测点的联合分布是p(y) N(y | 0, K(X, X) σₙ²I)这里 K(X, X) 是训练点之间的核矩阵再加上噪声方差构成观测协方差。5.2 联合分布与条件分布的推导对于预测点 X*m 个点训练输出 y 和预测函数值 f* 的联合分布是[y, f*]ᵀ ~ N(0, [[K σₙ²I, K*], [K*ᵀ, K**]])其中 K* K(X, X*) 大小为 n×mK** K(X*, X*) 大小为 m×m。利用多元高斯条件分布的公式条件后验均值是μ* K*ᵀ · (K σₙ²I)⁻¹ · y条件后验协方差是Σ* K** − Kᵀ · (K σₙ²I)⁻¹ · K这个公式不需要采样就能得到解析解这就是高斯过程回归比贝叶斯神经网络更让人舒服的地方。5.3 从公式看模型行为后验均值 μ* 有一种很直观的解释它在每个预测点上计算训练数据所有观测值的加权线性组合权重由核函数 k(xᵢ, x*) 决定。离预测点越近的训练点权重大离得远的权重小。这和RBF插值的加权平均思想是一脉相承的但高斯过程在此基础上又给出了方差。后验方差 Σ* 等于预测点的先验方差 K** 减去观测数据带来的信息量。越靠近训练点的区域后验方差越小离所有训练点都很远的区域后验方差会退回接近先验方差。这是高斯过程最迷人的地方——预测不确定性与空间位置直接挂钩。抽样工程里布置传感器可以先算一下哪些区域的预测方差最大然后专门去那里补点这就是主动学习的标准玩法。6. 超参数优化的实质在边际似然面上寻优高斯过程不是“无参模型”核函数里藏着超参数σ_f、l 以及观测噪声 σₙ。这些参数不能靠拍脑袋盯着一张曲线调标准做法是最大化对数边际似然。6.1 负对数边际似然的分解与含义观测数据 y 的边际似然是对潜在函数积分后的结果log p(y | X, θ) −½ yᵀ · (K σₙ²I)⁻¹ · y − ½ log|K σₙ²I| − n/2 · log(2π)右边第一项是“数据拟合项”衡量模型解释数据的程度数据拟合得越好这一项越小。第二项是“复杂度惩罚项”行列式的对数会随核矩阵特征值变大而变大相当于对模型复杂度做自然惩罚。第三项是归一化常数与超参数无关。高斯过程不会被复杂度惩罚项一巴掌拍死而是自动在拟合度和复杂度之间找均衡。我在调参时经常说别只盯着预测误差调边际似然里已经内置了奥卡姆剃刀。6.2 优化技巧对数域参数与多起点超参数必须为正所以实际优化的是 log(σ_f)、log(l)、log(σₙ)保证指数回退后永远是正数。目标函数定义为负对数边际似然然后用 fminsearch 或者 fminunc 求解。我一开始犯过错误在原始空间里优化迭代几次长度尺度变成负数核矩阵彻底崩掉。换成对数域以后再也没出过这种问题。另一个容易踩的坑是局部最优。边际似然面在长度尺度小的地方经常有高方差区域搜索很容易陷进去。稳妥起见用多个随机起点分别优化比较最终负对数边际似然值选最小的那组超参数。MATLAB 里负对数边际似然可以这样写function nll neg_log_marginal_likelihood(theta, X, y) sigma_f exp(theta(1)); l exp(theta(2)); sigma_n exp(theta(3)); K rbf_kernel(X, X, sigma_f, l) sigma_n^2 * eye(size(X,1)); L chol(K, lower); alpha L \ (L \ y); nll 0.5 * y * alpha sum(log(diag(L))) 0.5 * size(X,1) * log(2*pi); end用 chol 返回 L 之后sum(log(diag(L))) 等价于 0.5·log|K|都是数值稳定的写法。注意不要直接用 log(det(K))矩阵接近奇异时 det 会下溢或上溢完全不可靠。7. 完整的MATLAB实现与结果解读把前面所有环节串起来给一个小而完整的回归示例。假设我们的潜在函数是 f(x) sin(2πx)训练点 10 个噪声标准差 0.1。7.1 主程序与预测% 生成带噪声的观测数据 rng(42); X linspace(0, 1, 10); y sin(2 * pi * X) 0.1 * randn(10, 1); % 预测点 Xs linspace(-0.5, 1.5, 100); % 初始化超参数对数域 theta_init log([1.0, 0.2, 0.1]); % 优化超参数 theta_opt fminsearch((t) neg_log_marginal_likelihood(t, X, y), theta_init); % 用优化后的超参数计算后验 sigma_f exp(theta_opt(1)); l exp(theta_opt(2)); sigma_n exp(theta_opt(3)); K rbf_kernel(X, X, sigma_f, l) sigma_n^2 * eye(size(X,1)); L chol(K, lower); alpha L \ (L \ y); K_star rbf_kernel(X, Xs, sigma_f, l); mu_star K_star * alpha; % 计算后验方差 v L \ K_star; var_star diag(rbf_kernel(Xs, Xs, sigma_f, l)) - sum(v.^2, 1); std_star sqrt(var_star); % 画图 figure; hold on; fill([Xs; flipud(Xs)], [mu_star - 1.96*std_star; flipud(mu_star 1.96*std_star)], ... [0.8, 0.9, 1.0], EdgeColor, none); plot(Xs, mu_star, b-, LineWidth, 2); plot(X, y, ro, MarkerSize, 8, MarkerFaceColor, r); xlabel(x); ylabel(y); legend(95%置信区间, 后验均值, 观测数据, Location, best);7.2 看结果时要注意的几个数字跑完之后你会看到优化得到的长度尺度 l 大约在 0.15 到 0.25 之间对应数据波动的空间频率。后验均值在 x0.5 附近能较准确地描述正弦波峰但在两个边界 x0 和 x1 的区域置信区间快速张大因为那里没有观测数据模型在向外推时表现出了诚实的不确定性。95%置信区间宽度由 1.96·std_star 决定整个区间带在训练数据点附近明显收紧在两个相邻训练点中间略微变宽。这种区间宽度的空间变化是高斯过程相对普通回归的卖点。7.3 训练点数量对计算量的影响必须强调核矩阵求逆或 Cholesky 分解的计算复杂度是 O(n³)n 是训练样本数。300 个点以内基本秒开3000 个点就开始吃力30000 个点不可行。如果数据量上来了标准高斯过程要上稀疏近似或者随机特征扩展。你如果只是做小样本拟合、贝叶斯优化代理模型、实验设计这个复杂度完全够用。8. 实际使用中的数值稳定性坑与优化建议高斯过程看起来公式干净但工程实现有大量细节我把真正踩过的坑和习惯做法一次性列出来。8.1 矩阵病态与Jitter处理核矩阵对角线加上噪声方差 σₙ²之后理论上正定但实际中两个非常接近的观测点会让矩阵几乎奇异。此时 chol 函数直接报错。常规做法是在 K 的对角线上加一个很小的常数比如 1e-10甚至 1e-6这个操作叫加 jitter。加了之后对预测均值几乎没有影响但能让矩阵稳定分解。注意不要为了省事把 jitter 设成过大值比如和信号方差同量级那样会人为增大噪声压低后验方差的表达能力。工程经验是先设为 1e-8如果 chol 报错再逐步加大。8.2 特征值分解判断问题位置如果加了 jitter 还是报错先检查核函数的距离计算是否有问题。我之前写过一版 rbf_kernel距离矩阵里有很多 NaN结果核矩阵对角线全是 NaNCholesky 直接崩溃。打印 K 的前几行总能第一时间发现问题。如果问题确实来自病态用 eig 快速查看最小特征值有多大。最小特征值接近机器精度时意味着某些超参数下核函数把两个点定义得过于接近可以考虑修改长度尺度初始化值或者重新检查输入数据。8.3 数据归一化几乎必须做高斯过程对输入输出尺度很敏感。如果 x 范围在 0 到 1e6长度尺度 l 会被迫学习到百万量级如果 y 的量纲差距太大σ_f 的优化也会很不稳。实际项目里我会把输入和输出都归一化到零均值、单位方差训练完成后再反变换回原始量纲。有人会问归一化是不是改变了核函数的解释性确实会但你换来的是优化稳定性。只要把长度尺度等参数归一化后的值记录下来反归一化时可以换算出原始尺度下的相关长度。8.4 复杂度惩罚与超参数初始化的关系负对数边际似然函数非常容易在长度尺度很小区间产生多个局部极小点。别只用一组初始值我会随机抽 5 到 10 组初始值分别优化保留最优的一组。这样做的好处不仅仅是找到更好的参数还能顺带观察多组解之间后验曲线的差异如果差异很大说明数据对超参数不够敏感这时候对预测区间的解读要格外谨慎。8.5 预测点数量与内存控制预测时如果用一个特别密的网格K** 会变得巨大。设计实验时我会先只计算方差的对角线而不算了整个 K**因为大多数预测场景只需要单点方差。在上面的代码里用到的是 var_star diag(...) - sum(v.^2, 1)这已经是省内存的写法。如果连预测点也上百万就分段做预测避免一次性把 K_star 矩阵全部放进内存。最后聊聊我的实际体会高斯过程这东西最初学的时候容易被数学绕晕但真正上手发现核心就三件事选核函数、算核矩阵、解一个线性系统。在实际工程里它最出彩的地方不是给你一条好看的拟合曲线而是告诉你哪里预测可信、哪里不可信。做设备标定的时候我用它来确定下一组实验点在哪个位置最划算做仿真代理模型的时候我用它来代替每轮都跑重型仿真。每次跑通都有一种“这钱省得真值”的感觉。如果在你的项目里数据量不大但预测置信区间非常重要高斯过程值得作为第一选择。
平方复杂度无法被常规优化消除｜TaoToken 统一 Key 通道实测)
