高斯过程与高斯过程回归

本部分内容来自于经典的 Gaussian Processes for machine learning ¹ . 主要说明了高斯过程的建模方法和高斯过程回归是如何做的。

本文当中：

小写粗体字母代表向量 $\boldsymbol{x}\in \mathbb{R}^{n}$ .
大写粗体字母代表矩阵 $\boldsymbol{X}\in \mathbb{R}^{m\times n}$ .

权重空间视角下的高斯过程

权重空间下的高斯过程实质上就是贝叶斯线性模型.

函数空间视角下的高斯过程

高斯过程（随机过程）

函数空间视角下的高斯过程是一类随机过程，定义为：

高斯过程：对随机过程 $f$，如果对于任意 $\boldsymbol{X} = (\boldsymbol{x}_{1}, \boldsymbol{x}_{2},\cdots, \boldsymbol{x}_{n})\in \mathcal{X}^{n}$ ，都有 $(f(\boldsymbol{x}_{1}), \cdots, f(\boldsymbol{x}_{n}))\in \mathbb{R}^{n}$ 是多元高斯随机变量，则 $f$ 称为高斯过程，记为 $f\sim \mathcal{GP}(\mu, \kappa)$. 其中：
$\mu: \mathcal{X}\to \mathbb{R}$ 为均值函数 (mean function)；
$\kappa: \mathcal{X}\times \mathcal{X}\to \mathbb{R}$ 为核函数 (kernel function) 或称为协方差函数 (covariance function).
上述的输入空间通常为 $\mathcal{X} \subset \mathbb{R}^d$，即为实数向量.

这里均值函数和核函数还可以扩展记为

$$ \begin{aligned} \mu(\boldsymbol{X}) & = \mathbb{E}[f(\boldsymbol{X})]\in \mathbb{R}^{n} \\ \kappa(\boldsymbol{X}, \boldsymbol{X}') & = \mathbb{E}[(f(\boldsymbol{X}) - \mu(\boldsymbol{X}))(f(\boldsymbol{X}') - \mu(\boldsymbol{X}'))^{\top}] \in \mathbb{R}^{n\times m} \end{aligned} $$

这里 $\boldsymbol{X}' = (\boldsymbol{x}_{1}', \cdots, \boldsymbol{x}_{m}')\in \mathcal{X}^{m}$，写开就是：

$$ \boldsymbol{\mu}(\boldsymbol{X}) = \begin{pmatrix}\mu(\boldsymbol{x}_{1}) \\ \mu(\boldsymbol{x}_{2}) \\ \vdots \\ \mu(\boldsymbol{x}_{m})\end{pmatrix} $$

以及

$$ \kappa(\boldsymbol{X}, \boldsymbol{X}') = \begin{pmatrix}\kappa(\boldsymbol{x}_{1}, \boldsymbol{x}_{1}') & \kappa(\boldsymbol{x}_{1}, \boldsymbol{x}_{2}') & \cdots & \kappa(\boldsymbol{x}_{1}, \boldsymbol{x}_{m}') \\ \kappa(\boldsymbol{x}_{2}, \boldsymbol{x}_{1}') & \kappa(\boldsymbol{x}_{2}, \boldsymbol{x}_{2}') & \cdots & \kappa(\boldsymbol{x}_{2}, \boldsymbol{x}_{m}') \\ \vdots & \vdots & \ddots & \vdots \\ \kappa(\boldsymbol{x}_{n}, \boldsymbol{x}_{1}') & \kappa(\boldsymbol{x}_{n}, \boldsymbol{x}_{2}') & \cdots & \kappa(\boldsymbol{x}_{n}, \boldsymbol{x}_{m}') \end{pmatrix}\in \mathbb{R}^{n\times m} $$

常见的均值函数

通常而言，均值函数取为 $0$ ，这往往是为了简单起见，方便分析，但是事实上也有更强的均值函数，例如多项式函数 ² ：

$$ \mu(\boldsymbol{x}) = \sum\limits_{i=1}^{p} \beta_{i} f_{i} (\boldsymbol{x}) = [\boldsymbol{f}(\boldsymbol{x})]^{\top} \boldsymbol{\beta}, \quad \boldsymbol{x}\in \mathcal{X} $$

这里的 $\boldsymbol{\beta}\in \mathbb{R}^{p}$ 是对应的系数，假设 $\boldsymbol{x}_{i} = (x_{i,1}, x_{i,2},\cdots, x_{i,n})^{\top}\in \mathcal{X}$ ，$\boldsymbol{X} = (\boldsymbol{x}_{1},\cdots, \boldsymbol{x}_{m})^{\top}\in \mathcal{X}^{m}$，而 $\boldsymbol{f}(\boldsymbol{x}_{i}) = (f_{j}(\boldsymbol{x}_{i}))\in \mathbb{R}^{p}$ .

零阶多项式就是常数均值函数，此时

$$ \mu(\boldsymbol{x}) = c\in \mathbb{R}, \quad \boldsymbol{x}\in \mathcal{X} $$

一阶多项式形式：

$$ f_{1}(\boldsymbol{x}_{i}) = 1, f_{2}(\boldsymbol{x}_{i}) = x_{i,1},\cdots, f_{n+1}(\boldsymbol{x}_{i}) = x_{i, n} $$

那么此时就是一阶多项式的均值函数. 例如

$$ \mu(\boldsymbol{x}_{i}) = 1+ \sum\limits_{j=1}^{n} x_{i, j} $$

二阶多项式形式就相当于：

$$ \begin{cases} f_{1}(\boldsymbol{x}_{i}) = 1 \\ f_{2}(\boldsymbol{x}_{i}) = x_{i,1},\cdots , f_{n+1} (\boldsymbol{x}_{i}) = x_{i,n} \\ f_{n+2} (\boldsymbol{x}_{i}) = x_{i,1}^{2} ,\cdots, f_{2n+1} (\boldsymbol{x}) = x_{i,1} x_{i,n} \\ f_{2n+2} (\boldsymbol{x}_{i}) = \boldsymbol{x}_{i,2}^{2} ,\cdots, f_{3n}(\boldsymbol{x}_{i}) = x_{i,2} x_{i,n} \\ \quad \vdots \end{cases} $$

总共 $\dfrac{(n+1)(n+2)}{2}$ 项.

针对设计矩阵 $\boldsymbol{X}$ ，可以对应有

$$ \boldsymbol{\mu}(\boldsymbol{X}) = \begin{pmatrix}\mu(\boldsymbol{x}_{1}) \\ \mu(\boldsymbol{x}_{2}) \\ \vdots \\ \mu(\boldsymbol{x}_{m})\end{pmatrix} = \boldsymbol{F} \boldsymbol{\beta} $$

这里

$$ \boldsymbol{F} = \begin{pmatrix}f_{1}(\boldsymbol{x}_{1}) & f_{2}(\boldsymbol{x}_{1}) & \cdots & f_{p}(\boldsymbol{x}_{1}) \\ f_{1}(\boldsymbol{x}_{2}) & f_{2}(\boldsymbol{x}_{2}) & \cdots & f_{p}(\boldsymbol{x}_{2}) \\ \vdots & \vdots & \ddots & \vdots \\ f_{1}(\boldsymbol{x}_{m}) & f_{2}(\boldsymbol{x}_{m}) & \cdots & f_{p}(\boldsymbol{x}_{m}) \end{pmatrix}\in \mathbb{R}^{m\times p} $$

这就可以直接进行均值函数的计算了.

[!faq] 例：均值函数计算
设样本为
$$ \boldsymbol{X} = \begin{pmatrix}3 & 2 \\ 5 & 1 \\ 4 & 6\end{pmatrix} $$
若选择一阶多项式作为均值函数，$\boldsymbol{\beta}$ 为全 $1$ 向量，那么请计算 $\boldsymbol{F}$ .

考虑到

$$ \boldsymbol{x}_{1} = \begin{pmatrix}3 \\ 2\end{pmatrix}, \boldsymbol{x}_{2} = \begin{pmatrix}5 \\ 1\end{pmatrix}, \boldsymbol{x}_{3}=\begin{pmatrix}4 \\ 6\end{pmatrix} $$

于是 $\boldsymbol{F}\in \mathbb{R}^{3\times 3}$ ，于是

$$ f_{1}(\boldsymbol{x}_{1}) = 1, f_{2}(\boldsymbol{x}_{1}) = 3, f_{3}(\boldsymbol{x}_{1}) = 2 $$

相应地可以推导出

$$ \boldsymbol{F} = \begin{pmatrix}1 & 3 & 2 \\ 1 & 5 & 1 \\ 1 & 4 & 6\end{pmatrix} $$

这就能得到对应的均值向量为 $(6, 7, 11)$ . $\square$

常见的核函数

核函数的作用主要是表征不同输入点之间的关系,度量样本自身及其相互之间的相似程度,进行未知点处的函数值预测。有效的协方差函数须满足其对应的[[协方差矩阵]]对称、半正定条件.

协方差函数可分为平稳和非平稳协方差函数。平稳协方差函数只与样本点之间的距离 $r = \|\boldsymbol{x}-\boldsymbol{x}'\|_{2}$ 有关,故可将协方差函数记为 $C(r)$，而非平稳协方差与样本点本身有关，如线性核函数.

高斯过程抽样

下面为了对高斯过程的样子有个直观的感受，我们考虑对高斯过程进行一次数值的抽样，来看看会是什么样子.

如果我们有 $n$ 个样本点：$\mathcal{D} = \left\lbrace (\boldsymbol{x}_{i}, y_{i}) \right\rbrace_{i=1}^{n}$ ，设

$$ f\sim \mathcal{GP}(0, \kappa) $$

核函数此时设为 [[RBF 核]]，那么对于 $n$ 个点，此时得到的协方差阵就是：

$$ \boldsymbol{K} = \begin{pmatrix}\kappa(\boldsymbol{x}_{1}, \boldsymbol{x}_{1}) & \cdots & \kappa(\boldsymbol{x}_{1}, \boldsymbol{x}_{n}) \\ \vdots & \ddots & \vdots \\ \kappa(\boldsymbol{x}_{n}, \boldsymbol{x}_{1}) & \cdots & \kappa(\boldsymbol{x}_{n}, \boldsymbol{x}_{n})\end{pmatrix} $$

那么，此时就相当于可以在 $n$ 元正态分布：

$$ \mathcal{N}(\boldsymbol{0}, \boldsymbol{K}) $$

当中进行抽样，一个样本就是 $n$ 元向量，因此只要 $n\to \infty$ ，一个样本就足以绘图.

我们可以绘图查看，对 $n=201$ ，此时对于一个样本，抽样会得到 $201$ 维的向量，针对 $n=1:201$ 进行绘图可以得到：

高斯过程抽样.png

下面给出相应的 MATLAB 代码.

%% 1. 定义核函数（平方指数核）
sf = 1;      % 信号标准差
ell = 0.04;  % 长度尺度
kernel = @(x1, x2) sf^2 * exp(- (x1 - x2).^2 / (2*ell^2));

%% 2. 离散输入点（一维网格）
x = (0:0.001:0.2)';   % 从 0 到 0.2，步长 0.001
n = length(x);

%% 3. 计算协方差矩阵 K (n×n)
K = zeros(n, n);
for i = 1:n
    for j = 1:n
        K(i,j) = kernel(x(i), x(j));
    end
end
% 数值稳定性：添加极小噪声（jitter）
jitter = 1e-8;
K = K + jitter * eye(n);

%% 4. 从多元正态分布中抽取样本
rng(42);               % 固定随机种子，使结果可重复
numSamples = 10;        % 要绘制的样本数量
samples = mvnrnd(zeros(1, n), K, numSamples);  % 每一行是一个样本

%% 5. 绘制样本曲线
figure;
hold on;
for i = 1:numSamples
    % plot(x, samples(i, :), 'LineWidth', 1.2);
    scatter(x, samples(i, :), 3)
end
xlabel('Input x');
ylabel('f(x)');
title(sprintf('GP samples with zero mean, SE kernel (ell=%.2f, sf=%d)', ell, sf));
grid on;
hold off;

高斯过程回归

无噪声情况

假设 $\boldsymbol{X}$ 是 $n$ 个样本点，对应的无噪声观测值 $\boldsymbol{f} = (f_{1}, \cdots, f_{n})^{\top}$ ，即 $f_{i} = f(\boldsymbol{x}_{i}), i=1,2,\cdots,n$，而 $f\sim \mathcal{GP}(\mu, \kappa)$ .

模型输出服从高斯先验分布：

$$ \boldsymbol{f}\sim \mathcal{N}(\mu(\boldsymbol{X}), \kappa(\boldsymbol{X}, \boldsymbol{X})) $$

其密度函数对应记为 $p(\boldsymbol{f}\mid \boldsymbol{X})$，此时，假设有新的待测点 $\boldsymbol{x}_{*}$ 需要知道预测值 $f_{*}$ ，那么此时有联合高斯分布：

$$ (\boldsymbol{f}, f_{*}) \sim \mathcal{N}\left(\begin{pmatrix}\mu(\boldsymbol{X}) \\ \mu(\boldsymbol{x}_{*})\end{pmatrix}, \begin{pmatrix}\kappa(\boldsymbol{X}, \boldsymbol{X}) & \kappa(\boldsymbol{X}, \boldsymbol{x}_{*}) \\ \kappa(\boldsymbol{x}_{*}, \boldsymbol{X}) & \kappa(\boldsymbol{x}_{*}, \boldsymbol{x}_{*})\end{pmatrix}\right) $$

为了表述方便，我们使用 $\boldsymbol{\kappa}_{*}$ 来表达 $\kappa(\boldsymbol{X}, \boldsymbol{x}_{*})$ 列向量，$\kappa_{**}$ 表达数 $\kappa(\boldsymbol{x}_{*}, \boldsymbol{x}_{*})$，用 $\boldsymbol{K}$ 表示 $\kappa(\boldsymbol{X}, \boldsymbol{X})$ .

其联合密度函数描述为 $p(\boldsymbol{f}, f_{*}\mid \boldsymbol{X}, \boldsymbol{x}_{*})$ . 根据：

$$ p(f_{*} \mid \boldsymbol{f}, \boldsymbol{X}, \boldsymbol{x}_{*}) = \dfrac{p(\boldsymbol{f}, f_{*}\mid \boldsymbol{X}, \boldsymbol{x}_{*})}{p(\boldsymbol{f}\mid \boldsymbol{X}, \boldsymbol{x}_{*})} = \dfrac{p(\boldsymbol{f}, f_{*}\mid \boldsymbol{X}, \boldsymbol{x}_{*})}{p(\boldsymbol{f}\mid \boldsymbol{X})} $$

因此后验分布还是高斯分布，推导后可以得到

$$ f_{*} \mid \boldsymbol{X} , \boldsymbol{f}, \boldsymbol{x}_{*} \sim \mathcal{N}(\mu_{*}, \sigma_{*}^{2}) $$

这里

$$ \mu_{*} = \mu(\boldsymbol{x}_{*}) + \boldsymbol{\kappa}_{*}^{\top} \boldsymbol{K}^{-1} (\boldsymbol{f} - \mu(\boldsymbol{X})), \quad \sigma_{*}^{2} = \kappa_{**} - \boldsymbol{\kappa}_{*}^{\top} \boldsymbol{K}^{-1} \boldsymbol{\kappa}_{*} $$

后验分布的计算使得我们不仅有预测值 $\mu_{*}$ ，还有对应的不确定性度量 $\sigma_{*}^{2}$ ，这个方差是很有意义的，在[[贝叶斯优化]]当中，它就是指示下一步采样的利器.

有噪声情形

现在假如原来的函数输出是有随机噪声的，通常表示为：

$$ y = f(\boldsymbol{x}) + \varepsilon, \quad \boldsymbol{x}\in \mathcal{X}, \varepsilon\sim \mathcal{N}(0, \sigma_{\mathrm{n}}^{2}) $$

这里 $\sigma_{\mathrm{n}}^{2}$ 是噪声方差，是一个超参数，是需要调节的，表示对预测目标噪声水平的先验理解 ³.

假设 $\boldsymbol{X}$ 是 $n$ 个样本点，对应的有噪声观测值 $\boldsymbol{y} = (y_{1}, \cdots, y_{n})^{\top}$ ，即 $y_{i} = f(\boldsymbol{x}_{i}) + \varepsilon_{i}, i=1,2,\cdots,n$，而 $f\sim \mathcal{GP}(\mu, \kappa)$ . 模型输出服从高斯先验分布：

$$ \boldsymbol{y}\sim \mathcal{N}(\mu(\boldsymbol{X}), \boldsymbol{K}_{y}) $$

这里我们记 $\boldsymbol{K} = \kappa(\boldsymbol{X}, \boldsymbol{X})$ ，那么

$$ \boldsymbol{K}_{y} = \boldsymbol{K} + \sigma_{\mathrm{n}}^{2}\boldsymbol{I} $$

到后面的推导和无噪声情况是一致的，只是将 $\boldsymbol{K}$ 换成 $\boldsymbol{K}_{y}$ ，因此在有噪声情况下，我们也可以进行不确定性的度量.

高斯过程的调参

有哪些参数？

高斯过程回归的调参一般参数来自多个方面，其中包括：

噪声参数：即对噪声的参数，包括噪声参数；
均值函数参数：均值函数的参数；
核函数参数：核函数的参数，包括尺度参数和信号方差等；

不管有多少参数，我们首先可以都将其直接设为 $\boldsymbol{\theta}$ ，然后再细分参数类型，例如各向同性、带噪声预测以及信号方差的 [[RBF 核]]会有参数：

$$ \boldsymbol{\theta} = \left\lbrace l, \sigma_{f} , \sigma_{\mathrm{n}}\right\rbrace $$

MLE

[[极大似然估计]]法是一个非常常用的高斯过程参数估计方法，它的主要思路就是先建立样本集合的边际似然函数：

$$ p(\boldsymbol{y}\mid \boldsymbol{X}, \boldsymbol{\theta}) \sim \mathcal{N}(\boldsymbol{0}, \boldsymbol{K}+ \sigma_{\mathrm{n}}^{2} \boldsymbol{I}) $$

这里 $\boldsymbol{\theta}$ 就是超参数，默认使用带噪声情形的进行推导，我们建立负对数似然函数：

$$ - \log p(\boldsymbol{y}\mid \boldsymbol{X}, \boldsymbol{\theta}) = \underbrace{\frac{1}{2} \boldsymbol{y}^{\top} \boldsymbol{G}^{-1} \boldsymbol{y}}_{\text{empirical loss}} + \underbrace{\frac{1}{2} \log |\boldsymbol{G}|}_{\text{structural penalty}} + \frac{n}{2} \log 2 \pi $$

式里 $\boldsymbol{G} = \boldsymbol{K}+ \sigma^{2}_{\mathrm{n}} \boldsymbol{I}$ ，这里面的右侧第一项是经验风险，第二项则是对模型复杂度进行惩罚，有了这种做法之后，剩下的就是求导然后求零点，假设对各向异性、带噪声预测以及信号方差的 [[RBF 核]]，其参数为 $\boldsymbol{\theta} = (\boldsymbol{l}, \sigma_{f}, \sigma_{\mathrm{n}})$ ，那么最小化负对数似然就是求下面的零点：

$$ \begin{aligned} \dfrac{\partial \log p(\boldsymbol{y}\mid \boldsymbol{X}, \boldsymbol{\theta})}{\partial l_{k}} & = - \frac{1}{2} \boldsymbol{y}^{\top} \boldsymbol{G}^{-1} \dfrac{\partial \boldsymbol{G}}{\partial l_{k}} \boldsymbol{G}^{-1} \boldsymbol{y} - \frac{1}{2}\mathrm{tr}\left(\boldsymbol{G}^{-1} \dfrac{\partial \boldsymbol{G}}{\partial l_{k}}\right), k=1,2,\cdots, d \\ \dfrac{\partial \log p(\boldsymbol{y}\mid \boldsymbol{X}, \boldsymbol{\theta})}{\partial \sigma_{f}} & = - \frac{1}{2} \boldsymbol{y}^{\top} \boldsymbol{G}^{-1} \dfrac{\partial \boldsymbol{G}}{\partial \sigma_{f}} \boldsymbol{G}^{-1} \boldsymbol{y} - \frac{1}{2}\mathrm{tr}\left(\boldsymbol{G}^{-1} \dfrac{\partial \boldsymbol{G}}{\partial \sigma_{f}}\right) \\ \dfrac{\partial \log p(\boldsymbol{y}\mid \boldsymbol{X}, \boldsymbol{\theta})}{\partial \sigma_{\mathrm{n}}} & = - \frac{1}{2} \boldsymbol{y}^{\top} \boldsymbol{G}^{-1} \dfrac{\partial \boldsymbol{G}}{\partial \sigma_{\mathrm{n}}} \boldsymbol{G}^{-1} \boldsymbol{y} - \frac{1}{2}\mathrm{tr}\left(\boldsymbol{G}^{-1} \dfrac{\partial \boldsymbol{G}}{\partial \sigma_{\mathrm{n}}}\right) \end{aligned} $$

这里 $\boldsymbol{l} = (l_{1},\cdots, l_{d})$，也就是分维度求导，由此可以得到最终的合适的参数 $\widehat{\boldsymbol{\theta}}$ .

MAP

[[最大后验估计]]可以看作引入对 $\boldsymbol{\theta}$ 的先验分布 $p(\boldsymbol{\theta})$，从频率派视角看来就是加入正则化，最大化如下的函数：

$$ \boldsymbol{\theta}^{*} \in \arg\max_{\boldsymbol{\theta}} [\log p(\boldsymbol{y}\mid \boldsymbol{X}, \boldsymbol{\theta}) + \log p(\boldsymbol{\theta})] $$

高斯过程的运算

高斯过程的和可以表达成如下形式，假设

$$ f_{1}\sim \mathcal{GP}(\mu_{1}, \kappa_{1}), \quad f_{2} \sim \mathcal{GP}(\mu_{2}, \kappa_{2}) $$

那么

$$ f := f_{1} + f_{2} \sim \mathcal{GP}(\mu_{1} + \mu_{2} , \kappa_{1} + \kappa_{2}) $$

这是核函数的结构使然. ⁴