通常我们会遇到模型的拟合问题,即给定一系列数据 $(x_i,y_i)$,根据预先设定的模型 $$y_i=f(x_i,\theta)+\epsilon$$ 对模型参数进行最小二乘拟合,其中 $\theta$ 为模型参数,$\epsilon$ 为满足正态分布的噪声,即 $\epsilon \sim \mathcal N(0,\sigma^2_n)$。以线性模型拟合为例,$$y_i= a_0+a_1x_i+\epsilon$$,其中 $a_0$,$a_1$ 为线性模型的参数。然而在绝大多数情况下我们并不知道预先设定的模型是什么样子,因此模型参数也就无从定义,进而无法根据模型来进行预测。针对上述问题,高斯过程回归(Gaussian Process Regression, GPR) 可在无模型参数的情况下对未知量进行预测。下文主要结合高斯过程回归的 Python 程序包 GPy(The Gaussian processes framework in Python) 来阐述高斯过程回归的基本原理。

高斯过程回归的基本原理

已知一系列数据 $(x_i,y_i)$ 和 $x^*$,如何求 $y^*$ 及其置信区间?针对该问题的求解,我们引入高斯过程回归的基本原理。

核函数

假定存在一个最佳的映射 f,使得 $x_i \rightarrow f(x_i)$,那么根据 $x_i$,可求得 $f(x_i)$,其中 $i=1…n$。记 $\mathbf x = [x_1 \ldots x_n]^T$,$\mathbf f(\mathbf x) = [f(x_1) \ldots f(x_n)]^T$,我们规定 $\mathbf f(\mathbf x)$ 满足 n 维正态分布 $\mathcal{N}(\mathbf 0,\Sigma)$,其中协方差矩阵 $$\Sigma = K_{xx} + \sigma_n^2 I = \left [
\begin{matrix}
k(x_1,x_1) & k(x_1,x_2) & \ldots & k(x_1,x_n) \\
k(x_2,x_1) & k(x_2,x_2) & \ldots & k(x_2,x_n) \\
\vdots & \vdots & \vdots & \vdots \\
k(x_n,x_1) & k(x_n,x_2) & \ldots & k(x_n,x_n)
\end{matrix}
\right ] + \sigma_n^2 I$$,其中 $k(x_i,x_j)$ 为核函数,表示任意两个数据之间的相关关系。核函数满足如下性质:

  • 对称性,即 $k(x_i,x_j) = k(x_j,x_i)$;
  • 正定性,即核函数需要使 $K_{xx}$ 为正定矩阵;

核函数多种多样,例如

  • Squared Exponential Kernel:$k_{SE}(x_i,x_j) = \sigma^2 e^{-\dfrac{(x_i-x_j)^2}{2l^2}}$,其中 $\sigma^2$ 为振幅因子,$l$ 为长度尺度因子,都是用来调整协方差矩阵的超参数(Hyperparameters),记为 $\gamma$;
  • Linear Kernel:$k_{Lin}(x_i,x_j) = \sigma^2 x_i^T x_j$;
  • Periodic Kernel:$k_{Per}(x_i,x_j) = \sigma^2e^{−\dfrac{2\sin^2(\nu\pi|x_i−x_j|)}{l^2}}$,其中 $\nu$ 为频率;

一些复杂的核函数可通过两个核函数相加而得到,例如线性核函数与周期核函数的组合。

高斯过程

高斯过程回归从字面上理解为构建高斯分布并进行回归分析。构建高斯分布指的是用数据 $(x_i,y_i)$ 构建一个 n 维高斯分布,数据量越大,高斯分布的维度越高。回归分析则是通过构建的高斯分布对未知量进行预测。

构建高斯分布主要取决于均值和协方差矩阵。通常我们取均值为零即可,协方差矩阵则是由核函数决定,而核函数则是由数据的先验信息而定。例如数据同时表现为趋势性和周期性,那么核函数为 $$k(x_i,x_j) = \sigma_1^2e^{-\dfrac{(x_i-x_j)^2}{2l_1^2}} + \sigma_2^2e^{−\dfrac{2\sin^2(\nu\pi|x_i−x_j|)}{l_2^2}}$$。一旦确定了核函数,那么也就确定了映射 $f \sim \mathcal{GP}(\mathbf 0, k(x_i,x_j))$,整个过程即为高斯过程(GP)。已知 $\{ x_i \}_{i=1 \ldots 50}$ 序列,我们可以根据高斯过程来生成无穷多个样本(事实上无穷多个样本在实际中无法实现,但是生成尽量多的样本则没有问题,例如一万个样本),步骤如下:

  • 根据序列 $\{x_i\}_{i=1 \ldots 50}$ 计算核函数 $k(x_i,x_j)$;
  • 将核函数 $k(x_i,x_j)$ 组合成协方差矩阵 $K_{xx}$ 的形式;
  • 根据正态分布 $\mathcal{N}(\mathbf 0,K_{xx})$ 生成一万个样本;

因为每一个样本均服从正态分布 $\mathcal{N}(\mathbf 0,K_{xx})$,因此对这一万个样本取平均值,我们期望它们的平均值为零,事实上当取无穷多个样本时,它们的平均值确实为零。

高斯过程回归

经过整个高斯过程后,我们需要对未知量进行回归分析。既然 $\mathbf f(\mathbf x)\sim\mathcal{N}(\mathbf 0,K_{xx})$,那么加入未知量 $f(x^*)$ 后,我们得到一个 n+1 维的变量 $$\left [
\begin{matrix}
y_1 \\
\vdots \\
y_n \\
f(x^*)
\end{matrix}
\right ]
\sim
\mathcal N \left ( \left [
\begin{matrix}
\mathbf 0 \\
0
\end{matrix}
\right ],
\left [
\begin{matrix}
K_{xx} + \sigma_n^2 I & K_{xx^*}^T\\
K_{xx^*} & k(x^*,x^*)
\end{matrix}
\right ]
\right )
$$,其中 $$K_{xx^*} = \left [
\begin{matrix}
k(x_1,x^*)
\ldots
k(x_n,x^*)
\end{matrix}
\right ]
$$。在已知一系列数据 $(x_i,y_i)$ 和 $x^*$的情况下估算 $y^*$ 也就转化为条件概率的问题,即求解 $f(x^*)|\mathbf y(\mathbf x)$,经过运算可得 $$f(x^*)|\mathbf y(\mathbf x) \sim \mathcal N (K_{xx^*} M^{−1}\mathbf y(\mathbf x), k(x^*, x^*) – K_{xx^*} M^{−1} K_{x^*x}^T)
$$,其中 $M = K_{xx}+\sigma_n^2I$,因此 $y^*$ 的期望值 $$E(y^*) = K_{xx^*} M^{−1}\mathbf y(\mathbf x)$$,方差 $$Var(y^*) = k(x^*, x^*) – K_{xx^*} M^{−1} K_{x^*x}^T$$。

超参数的选取

如何选取核函数中超参数 $\gamma$ 的值?在噪声方差未知的情况下如何选取噪声方差 $\sigma_n^2$ 的值?为了确定最似然的超参数和噪声方差(在噪声分布已经明确的情况下只需要确定最似然的超参数即可)的值,我们需要使 $-\log P(\mathbf y(\mathbf x)|\gamma,\sigma_n)$ 为最小值,即 $$ \underset{\gamma,\sigma_n}{Min} [-\log P(\mathbf y(\mathbf x)|\gamma,\sigma_n)]= \underset{\gamma,\sigma_n} {Min} [ \frac{1}{2} \mathbf y(\mathbf x)^T M^{-1} \mathbf y(\mathbf x) + \frac{1}{2}\log|M| + \frac{n}{2} \log(2\pi)]$$。

用 GPy 实现高斯过程回归

我们先随机生成 [-3,3] 的 20 个横坐标 X,然后根据公式 $Y = \sin X + \epsilon$ 生成纵坐标 Y,其中 $\epsilon$ 为位于 [-0.05,0.05] 的观测噪声。

选取 Squared Exponential Kernel 作为本例的核函数。在 GPy 中默认数据的维度、核函数的振幅因子 $\sigma^2$ 和长度尺度因子 $l$ 都为 1,此处我们将长度尺度因子设置为 0.2。

输出的核函数超参数为

rbf.valueconstraintspriors
variance1.0+ve
lengthscale0.2+ve

从上表中可知核函数的振幅因子为 1.0,长度尺度因子为 0.2;这两个超参数都被限定为正值。输出的核函数图像为

Squared Exponential Kernel

图中画出了核函数 $k(x,0) = e^{-\dfrac{x^2}{2l^2}}$ 的函数图像,其中 $l$ 为 0.2。下面根据序列 $\{0,0.02 \ldots 1.0\}$ 共 51 个点来生成 10 个样本。

输出的图像为

核函数样本

图中共显示了 10 个样本,理论上这 10 个样本的均值应该为零。在明确核函数后,我们用之前的 X、Y 数据来构建一个高斯过程。

输出结果为

Model: GP regression
Objective: 25.344523930584618
Number of Parameters: 3
Number of Optimization Parameters: 3
Updates: True

GP_regression.valueconstraintspriors
rbf.variance1.0+ve
rbf.lengthscale0.2+ve
Gaussian_noise.variance1.0+ve

该高斯过程回归共有 3 个参数(2 个超参数、1 个观测噪声方差),GPy 默认观测噪声的方差为 1.0。Objective 表示接近最大似然的程度,值越小越接近。高斯过程回归的图像为

超参数没有最佳化前的高斯过程回归

从图中可看到拟合的曲线与 $Y = \sin X$ 相差甚远,这是因为超参数与观测噪声方差没有选取合适的值导致的。下面我们根据上文描述的最大似然方法来最佳化超参数与观测噪声方差的值,并重新构建高斯过程回归。

输出结果为

Model: GP regression
Objective: -19.5130492961455
Number of Parameters: 3
Number of Optimization Parameters: 3
Updates: True

GP_regression.valueconstraintspriors
rbf.variance0.5438471489551423+ve
rbf.lengthscale1.370942670247255+ve
Gaussian_noise.variance0.0012174968855384238+ve

与之前输出的 3 个参数相比,最佳化的超参数和观测噪声方差的值有明显的变化,并且 Objective 明显变小,这说明超参数和观测噪声方差的值更加接近最大似然。最佳化超参数和观测噪声方差后的高斯过程回归图像为

超参数最佳化后的高斯过程回归

在选取合适的超参数与观测噪声方差后,从图中可看到拟合的曲线与 $Y = \sin X$ 非常接近;图中浅蓝色的区域代表 95% 的置信区间,并且置信区间在两端明显变宽,这是由于两端没有观测数据导致的。

任意在 [-3,3] 中给定一个或多个点 x,计算这些点的 y 值及其方差。例如在样本点处计算均值和方差

输出结果为

由于样本点共有 51 个,因此上面的计算结果共有 51 个 y 值和方差值。

高维的高斯过程回归

高维高斯过程回归指的是数据不再是 $(x_i,y_i)$ 的形式,而是 $(\mathbf x_i,y_i)$ 的形式,其基本原理与一维高斯过程回归的基本原理类似。以二维高斯过程回归为例,记 $\mathbf X = [\mathbf x_1 \ldots \mathbf x_n]^T$,$\mathbf f(\mathbf X) = [f(\mathbf x_1) \ldots f(\mathbf x_n)]^T$,其中 $\mathbf x_i=[x^1_i,x^2_i]^T$,我们规定 $\mathbf f(\mathbf X)$ 满足 n 维正态分布 $\mathcal{N}(\mathbf 0,\Sigma)$,核函数变为 $k(\mathbf x_i,\mathbf x_j)$。同样在已知一系列数据 $(\mathbf x_i,y_i)$ 和 $\mathbf x^*$的情况下估算 $y^*$ 也就转化为条件概率的问题,即求解 $f(\mathbf x^*)|\mathbf y(\mathbf X)$,经过运算可得 $$f(\mathbf x^*)|\mathbf y(\mathbf X) \sim \mathcal N (K_{xx^*} (K_{xx}+\sigma_n^2I)^{−1}\mathbf y(\mathbf X), k(\mathbf x^*, \mathbf x^*) – K_{xx^*} (K_{xx}+\sigma_n^2I)^{−1} K_{x^*x}^T)
$$ ,因此 $y^*$ 的期望值和方差分别为 $$E(y^*) = K_{xx^*} (K_{xx}+\sigma_n^2I)^{−1}\mathbf y(\mathbf X) \\ Var(y^*) = k(\mathbf x^*, \mathbf x^*) – K_{xx^*} (K_{xx}+\sigma_n^2I)^{−1} K_{x^*x}^T$$。


参考

0 0 vote
Article Rating
Subscribe
Notify of
guest
5 Comments
oldest
newest most voted
Inline Feedbacks
View all comments
乔木
乔木
August 2, 2018 12:57

看着这个运算过程想到了数学王子高斯

xiaogua
xiaogua
November 5, 2019 16:37

你好,本人刚刚接触Python,请问为何我的程序提示
“NameError: name ‘display’ is not defined”

恳请大大告知~

Muratatsu
Muratatsu
January 19, 2020 16:21

你好,在进行多维输入的高斯过程回归时,程序的运行上是否只需替换输入数组即可?