高斯过程回归的权空间视点推导及代码实现

1.高斯过程简介

1.1定义

高斯过程是随机变量的集合，其中任意有限数量的随机变量具有联合高斯分布。

在函数空间(function-space view)的观点中，高斯过程可以看作是一个定义在函数上的分布，并且直接在函数空间中进行inference。

与之等价的观点是权空间观点(weight-space view)，在权空间中进行推断,对权向量的不同抽样将产生不同的函数，这与函数空间观点是一致的,但是函数空间的观点更为直接和抽象。

注：为方便起见，本文没有刻意区分概率密度和概率质量函数。向量用等小写字母表示，矩阵用等大写字母表示，标量会专门说明。

2.部分基础知识（已具备的直接跳至第3节）

2.1 部分矩阵计算基础

2.1.1 分块矩阵求逆

有兴趣可以看推导过程，否则直接看最终结论。

设置矩阵

是可逆矩阵，求矩阵的逆矩阵如下（推导是在逆矩阵存在的假设下进行的）。

认为
,
可用的

从、购买
将带入(1)并移项整理可得，

把带入整理得到，

令，其实就是关于的舒尔补(The Shur complements)。

单独的订单

其中是单位矩阵。
可以得到原始分块矩阵的逆矩阵

2.1.2 矩阵求逆引理

其中矩阵是可逆的。
证明：

，其中是单位矩阵，那么我们可以得到

可以从得到，带入整理

回来得到

2.2 多元高斯分布

2.2.1 联合分布

设为维向量，则其概率密度函数为

其中，和分别为随机向量的协方差矩阵和均值向量。

多维高斯分布具有非常好的特性：

1. 边际分布满足高斯分布。
   
2. 条件分布满足高斯分布。
   
3. 分量的线性组合也是一个高斯随机变量。
   

2.2.2 条件概率分布

令随机向量符合多维高斯分布，将其分成不相交的两个部分
，

的平均向量

协方差矩阵是

精度矩阵是

其中协方差矩阵是正定的，因为它的对称性，，。

注意这是一个分块矩阵，我们不能简单地将每个矩阵分块求逆，我们使用公式，我们可以得到

接下来，我们找到给定的的条件概率分布。注意，高斯分布的形式主要取决于指数项，所以我们使用匹配的方法找到对应的均值和协方差矩阵，而不考虑归一化系数，得到条件分布的形式。

索引项是

观察式(9)中指数项的形式，可发现精度矩阵出现在的二次项中，精度矩阵和均值的乘积出现在的线性项中，因此我们可得

由式(10)(11)可得

这样我们得到的分布，我们发现它的协方差与无关，而均值是的线性函数，实际上是线性高斯模型的一个例子。

2.2.3 简单的线性高斯模型及贝叶斯定理

贝叶斯公式：

我们设置

其中，是的精度矩阵，的均值是的线性函数。
接下来，我们想知道的联合分布。

依然使用配方的方法，关注于指数项的系数。根据式(15)可得，

所以观察

整理次要项后，有

因此，可以得到精度矩阵

根据公式(5)可得协方差矩阵，

再次查看项目

由此并根据式(13’)(16)可得均值

这个结果也符合我们的直觉，可以得到的边际分布为

3.高斯过程回归的权空间观点推导

先回忆一下一般的线性回归模型，我们就不介绍基函数了，

其中是代表实际数据值的一维变量，代表高斯噪声，我们假设

因此，由于训练样本是一个定量，那么

下面规定是由实际数据值组成的向量，是由输入组成的矩阵。这里我们异常规定的每一列都是输入，样本是

我们首先对 做一个先验分布假设，设置

Y的似然函数为

根据贝叶斯定理以及式(19)(20)可得的后验分布

得到的后验分布后，我们需要进行预测，即得到预测分布。给定测试样本，我们仍然考虑测试样本点有高斯噪声的情况。从根本上说，我们最终想要的不是带有噪声的样本值，而是生成这些数据的函数，这也符合定义：高斯过程是定义在函数上的分布。假设预测函数为，
但


这同样是一个线性高斯模型，我们需要求解边缘概率分布，由式(19)(20)可得

在

接下来，我们引入基函数，用基函数将样本输入映射到高维特征空间，用来表示映射后的特征向量(feature vector，与eigenvector区分)，用表示特征向量组成的矩阵。
但

在

但是显示表示一个合适的基函数(basis function)不是一件容易的事情，更别说加上一个先验的协方差矩阵，因此我们隐式的引入这一式子。

我们设，我们先对式(24)进行处理。
等式两边同时乘以左边的和右边的，整理并带入得到

将(25)式代入(23)的均值部分可得，

利用矩阵求逆引理(6)可得并带入(23)的协方差部分，可得

最后我们引入核函数的概念，设置核函数，
以此类推，，

我们这里介绍常用的高斯核函数，其中是高斯过程的length-scale，这里不做过多解释。

最终形式是


至此，权重空间视角的推导过程结束。

下面是python实现代码：

gaussian_process_regression.py

import numpy as np


class GaussianProcessRegressor:
    """
    kernel: RBF(sigma_overall, l_scale)
    alpha: noise, 1-D array or scaler
    """
    def __init__(self, kernel, sigma_overall, l_scale, alpha=0.):
        self.kernel = kernel(sigma_overall, l_scale)
        self.alpha = alpha

    def fit(self, X, y):
        X = np.asarray(X)
        y = np.asarray(y)
        self.train_x_ = X
        self.train_y_ = y

    def predict(self, X, return_cov=True, return_std=False):
        if return_cov and return_std:
            raise RuntimeError("return_cov, return_std can't be True in the same time")
        if not hasattr(self, 'train_x_'):
            y_mean = np.zeros(X.shape[0])
            if return_cov:
                y_cov = self.kernel(X, X)
                return y_mean, y_cov
            elif return_std:
                y_cov = self.kernel(X, X)
                return y_mean, np.sqrt(np.diag(y_cov))
            else:
                return y_mean
        K = self.kernel(self.train_x_, self.train_x_)
        L = np.linalg.cholesky(K + self.alpha * np.eye(self.train_x_.shape[0]))
        alpha = np.linalg.solve(L, self.train_y_)
        alpha = np.linalg.solve(L.T, alpha)
        y_mean = self.kernel(self.train_x_, X).T @ alpha
        v = np.linalg.solve(L, self.kernel(self.train_x_, X))
        y_cov = self.kernel(X, X) - v.T @ v + self.alpha * np.eye(X.shape[0])
        if return_cov:
            return y_mean, y_cov
        elif return_std:
            return y_mean, np.sqrt(np.diag(y_cov))
        else:
            return y_mean

    def sample_func(self, X, n_samples=1):
        y_mean, y_cov = self.predict(X, return_cov=True, return_std=False)
        sampled_y = np.random.multivariate_normal(y_mean, y_cov, size=n_samples)
        return sampled_y

kernel.py

import numpy as np


class RBFKernel:
    def __init__(self, sigma, scale):
        self.sigma = sigma
        self.scale = scale

    def __call__(self, x1: np.ndarray, x2: np.ndarray):
        m, n = x1.shape[0], x2.shape[0]
        K_matrix = np.zeros((m, n), dtype=float)
        for i in range(m):
            for j in range(n):
                K_matrix[i, j] = self.sigma * np.exp(-0.5 * np.sum((x1[i] - x2[j]) ** 2) / self.scale)
        return K_matrix

测试代码：

from gaussian_process import RBFKernel, GaussianProcessRegressor
import numpy as np
import matplotlib.pyplot as plt


def get_y(x, alpha):
    return np.cos(x)*0.3 + np.random.normal(0, alpha, size=x.shape)

observation_size = 6
gpr = GaussianProcessRegressor(RBFKernel, sigma_overall=0.04, l_scale=0.5, alpha=1e-4)
sample_size = 3

test_x = np.linspace(0, 10, 100)
prior_mean, prior_cov = gpr.predict(test_x, return_cov=True)
sample_ys = gpr.sample_func(test_x, n_samples=sample_size)
uncertainty = 1.96 * np.sqrt(np.diag(prior_cov))

plt.plot(test_x, prior_mean, label='mean')
plt.fill_between(test_x, prior_mean-uncertainty, prior_mean+uncertainty, alpha=0.1)
for i in range(sample_size):
    plt.plot(test_x, sample_ys[i], linestyle='--')
plt.legend()
plt.title('prior')
plt.show()  # 至此绘制先验分布

train_x = np.array([3, 1, 4, 5, 7, 9])
train_y = get_y(train_x, alpha=1e-4)

gpr.fit(train_x, train_y)
y_mean, y_cov = gpr.predict(test_x, return_cov=True)
sample_ys = gpr.sample_func(test_x, n_samples=sample_size)
uncertainty = 1.96 * np.sqrt(np.diag(y_cov))

plt.plot(test_x, y_mean, label='mean', linewidth=3)
plt.fill_between(test_x, y_mean-uncertainty, y_mean+uncertainty, alpha=0.1)
for i in range(sample_size):
    plt.plot(test_x, sample_ys[i], linestyle='--')
plt.scatter(train_x, train_y, c='red', marker='x', label='observation', linewidths=5)
plt.legend()
plt.title('posterior')
plt.show()  # 绘制后验分布

运行效果如下

下次我将添加功能空间的观点。

我是大二学生，水平有限。不严谨请见谅。

参考

[1]C。E。Rasmussen & C。K。I。Williams.(2006).Gaussian Processes for Machine Learning.7-29.

[2]Christopher M. Bishop.(2007).Pattern Recognition and Machine Learning. 78-93.