高斯过程回归（Gaussian Process Regression, GPR）是使用高斯过程（Gaussian Process, GP）先验对数据进行回归分析的非参数模型（non-parameteric model）。

与GPR有关的早期研究大致可分为3部分，在时间序列分析中，安德雷·柯尔莫哥洛夫（Андрей Николаевич Колмогоров）和罗伯特·维纳（Robert Wiener）在二十世纪40年代提出了估计0均值平稳高斯过程信号的滤波技术。随后出现的卡尔曼滤波（Kalman filter）提供了估计高斯隐变量（latent variable）的有效方法。在地统计学领域，1963年提出的克里金法（Kriging）首次实现了平稳随机场的非参数回归，其后在二十世纪90年代出现的贝叶斯克里金提供了与GPR接近或等价的高斯随机场估计。在机器学习方面，包含无限节点的贝叶斯神经网络（Bayesian neural network）和各类核学习（kernel learning）方法例如核岭回归（kernel ridge regression）为GPR的核函数超参数求解带来了启发。

GPR作为机器学习一般方法的提出来自剑桥大学（University of Cambridge）学者Carl E. Rasmussen和爱丁堡大学（University of Edinburgh）学者Christopher K. I. Williams，二者在1996年发表的研究按函数空间观点给出了GPR的求解系统并讨论了GPR超参数的极大似然估计和蒙特卡罗方法求解。

这里给出GPR在权重空间（weight-space）观点和函数空间（function-space）观点下的模型推导，前者由贝叶斯线性回归（Bayesian linear regression）出发，通过特征空间的映射得到GPR的预测形式，后者由高斯过程出发，由正态分布的边缘分布性质得到与前者等价的结果。

权重空间观点

参见：贝叶斯线性回归

在权重空间（weight-space）观点下，GPR可以由正态假设的贝叶斯线性回归导出。具体地，给定相互独立的N组学习样本： ，贝叶斯线性回归是如下形式的多元线性回归模型：

式中 为权重系数， 为残差或噪声。贝叶斯线性回归假设残差服从独立同分布（independent and identically distributed, iid）的0均值正态分布： ，由此可得贝叶斯线性回归的似然：

在此基础上，为模型权重赋予0均值正态先验： 。由贝叶斯定理（Bayes' theorem）可知，模型权重的后验正比于似然和先验的乘积： 。由正态分布的共轭性（conjugacy）可知，对正态分布的似然和方差已知的正态分布先验，其后验也为正态分布，因此带入似然和先验的解析形式并按正态分布整理可得：

给定测试样本 ，贝叶斯线性回归通过对可通过边缘化模型权重，即按其后验积分得到测试结果 的概率分布：

贝叶斯线性回归是一个线性参数模型，为使其表示样本间的非线性关系，可以使用给定的函数将 映射至高维空间：

由于映射函数 是固定的，即与模型权重无关，因此可直接带入贝叶斯线性回归的结果得到：

使用核方法（kernel method），即定义核函数 可改写上式：

上述改写结果即是GPR对 的均值和协方差进行预测的形式，式中 表示由学习样本得到的格拉姆矩阵（Gram matrix），即核矩阵， 表示带入测试样本后得到的行向量和列向量：

函数空间观点

参见：高斯过程

对回归模型 ，若函数 的形式不是固定的，则其为潜函数（latent function）。潜函数的每个取值都是函数空间（function-space）的一个测度。GPR取该函数空间的先验为高斯过程（Gaussian Process, GP），不失一般性，这里表示为0均值高斯过程：

式中 为学习样本，其在高斯过程中的测度是高斯过程的有限维分布（finite-dimensional distribution），由定义可知，该有限维分布是联合正态分布：

式中 为核函数，0均值高斯过程由其核函数完全决定：

给定N组学习样本 ，作为对算法的推导，假设回归残差服从iid正态分布：，则GPR在高斯过程先验和正态分布似然下求解回归模型的后验：，并对测试样本的测试结果进行估计（非正态似然的情形参见算法部分）。具体地，由回归模型和高斯过程的定义， 和 的概率分布为：

因此二者的联合概率分布是如下形式的联合正态分布：

对上述联合分布取 的边缘分布，由联合正态分布的边缘分布性质（marginalization property）可得：

上式是GPR的预测形式，也是函数空间后验对测试样本的有限维分布。式中协方差的展开使用了分块矩阵求逆公式。注意到，上式与权重空间观点下GPR的预测形式相同，即假设先验和iid残差为0均值正态分布的贝叶斯线性回归在引入核函数后，等价于使用相同核函数和0均值高斯过程先验的GPR。

由上述预测形式的推导可知，GPR在假设0均值高斯过程先验时，解得的后验往往不是0均值的，因此0均值先验具有泛用性。但作为说明，GPR也可使用均值不为0的高斯过程先验，此时常见的做法是将潜函数表示为一组显式基函数（explicit basis function）： 。在给定基函数的正态分布先验 时， 是具有如下形式的高斯过程：

参照先前的推导，此时GPR的预测形式为：

式中 ，具体步骤和讨论可参见Rasmussen and Williams (2006), pp.27-30。

参见：核函数

GPR使用高斯过程作为先验，即假设了学习样本是高斯过程的采样，因此其估计结果与核函数有密切联系。GPR中核函数的实际意义为协方差函数（covariance function），描述了学习样本间的相关性，因此其不是通过核方法简化计算的手段，而是模型假设的一部分。这里对GPR常见的核函数进行简单介绍，并引入核函数计算的加速方法。

核函数的选择

核函数的选择要求满足Mercer定理（Mercer's theorem），即核函数在样本空间内的任意格拉姆矩阵（核矩阵）为半正定矩阵（semi-positive definite）。

若GPR的先验为平移不变（transformation invariant）的平稳高斯过程，可用的核函数包括径向基函数核（RBF kernel）、马顿核（Matérn kernel）、指数函数核（exponential kernel）、二次有理函数核（rational quadratic kernel, RQ kernel）等，以RBF核为例，其解析形式如下：

式中 为RBF核的超参数，表示其带宽（width）和特征长度尺度（characteristic length-scale）， 是表征各向异性的矩阵函数，式中的3种形式分别对应各向同性和两类各向异性。常见地，对各向同性的情形，RBF核表示为：

GPR也可使用非平稳高斯过程，此时常见的核函数选择为周期核（periodic kernel）与多项式函数核（polynominal kernel），二者分别赋予高斯过程周期性和旋转不变性（rotation invariant）。

核矩阵求逆的计算简化

由模型推导部分可知，GPR要求计算核矩阵的逆矩阵： ，若使用Cholesky分解（Cholesky decomposition）求逆，则对N组学习样本，其计算复杂度为 。因此随着学习样本的增加，GPR的计算开销呈非线性增长。这里简单介绍核矩阵求逆的简化方法以应对上述问题，一般地，这些方法除GPR外也适用于其它核学习（kernel learning）方法，例如核岭回归、支持向量机（Support Vector Machine, SVM）等。

1. 选取数据子集（subset of datapoints, SD）

由于GPR可以在少量学习样本的情形下给出回归问题的可靠估计，因此从学习样本中合理选择一个子集可以降低核矩阵的大小，但不对GPR的表现造成明显影响。选取SD的常见方法是通过寻找微分熵（differential entropy）最大的学习样本作为支持向量定义SD，选取过程的计算复杂度为 ，在选取较小子集的情形下简化了计算，使用SD优化的高斯过程建模被称为稀疏高斯过程（sparse Gaussian process），或信息向量机（Informativevector machine, IVM）。

2. 低秩近似（low-rank approximation）

低秩近似将核矩阵近似为一系列低秩矩阵的乘积以简化求逆计算，具体地，N×N的核矩阵可以近似表示为： ，其中 的大小分别为N×m、m×m、m×N，此时由Sherman-Morrison-Woodbury定理可得：

式中求逆部分的矩阵大小为m×m，因此计算复杂度减少至 ，在m取小值时简化了计算。

确定低秩矩阵的常见方法包括Nyström法（Nyström method）、子集回归法（subset ofregressors, SR）、映射过程近似（projected process approximation, PP）、BCM（Bayesian committee machine）等，这里对前两者进行介绍。Nyström法将核矩阵分解为4个分块矩阵并在构建低秩近似后直接带入GPR的预测形式： SR按与Nyström法相同的方式构建分块矩阵，但其首先将GPR的预测形式改写为核矩阵子集的回归并带入低秩矩阵，例如对GPR的均值部分有： 。

GPR的求解也被称为超参学习（hyper-parameter learning），是按贝叶斯方法通过学习样本确定核函数中超参数的过程。由贝叶斯定理（Bayes' theorem）可知，GPR的超参数后验有如下表示：

式中 表示GPR的超参数，包括核函数的超参数和残差的方差 。 为似然（liklihood），作为非参数模型，这里的似然是对GPR的输出边缘化得到的边缘似然（marginal liklihood）： 由模型推导可知，在残差服从iid正态分布时，GPR的似然服从正态分布： ，此时其常见的求解方法是包含非线性优化的极大似然估计。

极大似然估计（Maximum Likelihood Estimation, MLE）

MLE通过令GPR的似然取极大值实现对超参数的单点估计。不失一般性，对0均值先验的GPR，带入联合正态分布的解析形式并求其负自然对数可得：

式中第一项包含学习样本，是数据拟合项，表示模型的经验风险（empirical risk）；式中第二项仅与回归模型有关，回归模型的核矩阵越复杂其取值越高，表示模型的结构风险（structural risk），因此按统计学习理论，GPR是一个包含正则化的求解系统，但二者的比例，即残差的方差由MLE给出。求解GPR似然极大值等价于求解负自然对数的极小值，因此这里给出其对 的偏导数：

非线性优化算法可以求解GPR的MLE问题，常见的选择为共轭梯度法（conjugate gradient method）和拟牛顿法（quasi-Newton method）。有研究使用随机优化方法，例如随机梯度下降（stochastic gradient descent）、遗传算法（genetic algorithm）和粒子群算法（particle swarm optimization）求解GPR。

GPR的对数似然不是凸函数，且其优化复杂度随学习样本的增加而增大，在学习样本较多的情形下，可能会发现多个局部最优，且局部最优解的差异很大的情形，考虑解出的核函数超参数通常表示高斯过程的特征长度尺度，多个局部最优意味着学习样本可以按多个尺度进行回归，其中小尺度的模型复杂度高但考虑更多局部特征，即具有高方差（high variance），大尺度的模型反之，具有高偏差（high bias）。

除MLE外，GPR也可通过极大后验估计（Maximum A Posteriori estimation, MAP）求解。MAP是与MLE相近的估计方法，其优化目标是似然和先验的乘积。因为MAP与MLE在正态后验时得到的结果等价，在其他情形下的结果也相近，却引入了先验和额外的计算，因此在GPR研究中被较少提及。

极大伪似然估计（maximum pseudo-likelihood estimator, MPLE）

使用MPLE求解GPR的过程也被称为交叉验证（cross-validation），相比于MLE，其改变是对学习样本使用留一法（Leave One Out, LOO）计算其伪似然（pseudo-likelihood）的负自然对数：

式中 表示除 以外的所有样本， 表示选取该矩阵的第i行第j列元素，由 的形式可知其为留一数据由GPR估计的均值和方差，式中第2行为留一数据的似然，GPR的伪似然为所有留一数据似然的和。MPLE求解GPR在计算开销上与MLE相当，因为两者共同的计算效率瓶颈是核矩阵求逆步骤（参见理论部分）。在求解方法上，MPLE独立第考察了每个学习样本，因此其本质上是频率估计（frequentist inference），有研究表明MPLE在GPR模型选择不恰当（miss-specified）时得到的结果更好，MLE在模型选择恰当（well-specified）时得到的结果更好。

编程实现

这里提供一个在Python 3环境下使用GPy进行GPR建模的例子：

在似然不服从正态分布，例如处理异方差噪声数据时，GPR对测试数据的后验没有解析形式，此时可以使用的解析近似（analytical approximation）方法，例如拉普拉斯近似（Laplace Approximation, LA）和期望传播（Expectation Propagation, EP）将非正态后验近似表示为正态后验。LA和EP的推导可参见Rasmussen and Williams (2006), pp.41-60。

数值方法可用在非正态似然，或核函数超参数过多，不利于MLE问题优化的情形下求解GPR。选择包括蒙特卡罗方法（Monte Carlo）和变分贝叶斯估计（Variational Bayesian Inference, VBI）。VBI使用Jensen不等式得到GPR对数似然的凸函数下界作为代理损失并使用梯度算法求解，其计算复杂度接近，在学习样本较多时节约了计算成本。蒙特卡罗方法被认为是在GPR数值解法中拥有较高的准确度和泛用性的方法，其原理是通过计算开销迭代生成随机数逼近GPR后验。Hamiltonian蒙特卡罗（Hamiltonian Monte Carlo, HMC）或混合蒙特卡罗（Hybrid Monte Carlo）是GPR的常见解法，其步骤接近于Metropolis-Hastings算法，但通过“动量”项限制了随机游走的取值。HMC求解GPR的例子可参见Rasmussen (1996)。

这里给出与高斯过程回归有关的扩展算法，但更一般地，这些方法是高斯过程模型的扩展，除回归问题外也可被应用于其它机器学习主题。

扭曲高斯过程（Warped Gaussian Process, WGP）

WGP通过对GPR进行非线性变换将其拓展至非高斯过程的情形，具体地，WGP首先在特征空间选择一个单调扭曲函数（monotonic warping function），将学习样本变换至封装函数指定的潜空间后再求解包含封装函数超参数的极大似然。WGP的中封装函数的选择通常为双曲正切函数（hyperbolic tangent）线性组合，这里给出WGP对应于GPR的模型和算法：

式中 为预先给定的WGP封装次数，在求解超参数后，类比GPR可得到WGP的预测形式。除上述形式的封装函数外，使用更复杂的扭曲函数，例如 可进一步提升WPG的灵活性。在一些测试数据中，WGP的表现优于GPR。

半参数高斯过程（Semi-parametric Gaussian Processes, SGP）

SGP是在回归问题中将高斯过程模型与参数模型线性组合以获得两者优点——GPR的灵活性和参数模型在高维问题中利用数据的有效性——的算法。具体地，SGPR构建了如下求解系统：

式中 为参数回归部分，常见形式为线性回归，第二项为0均值高斯过程先验的GPR，最后一项为残差。SGPR的学习分为两步，先对学习样本使用参数模型进行回归，随后使用GPR拟合参数模型的残差。SGP被使用的常见情形是在已知学习数据和学习目标的确定性关系时考虑噪声带来的随机效应，例如自动控制系统的预测问题。

深度高斯过程（deep Gaussian Process, DGP）

由性质可知，高斯过程可以视为一个拥有单隐含层和无限个隐含层节点的多层感知器（Multi-Layer Perceptron, MLP），DGP是MLP由单隐含层推广至多隐含层时得到的高斯过程，即拥有多个隐含层和无限宽度的MLP。按定义，DGP是一组从高斯过程先验中选择的函数的分布：

式中 表示第l层的输出。DGP在结构上是由非参数的高斯过程基函数表示隐含层的神经网络，并通过样本求解核函数超参数。DGP的另一种等价观点/结构是控制MLP中的多层无限节点保持不变并在每个隐含层后加入一组加权求和，由样本求解求和部分，即 的权重。

可加高斯过程（Additive Gaussian Process, AGP）

GPR的常见核函数，例如RBF核、马顿核等，被认为在高维特征空间的学习问题中泛化能力较差，而AGP即是为解决GPR高维学习问题而提出的扩展方法。AGP利用核函数性质，在构建高斯过程先验时对核函数进行了组合和结构优化。介绍性地，AGP在不同维度将不同结构的核函数相加作为先验，并通过贝叶斯阶层核学习（Bayesian Hierarchical Kernel Learning, HKL）更新超参数以去除对学习样本解释不利的结构。以RBF核为例，其1阶可加高斯过程先验可表示为：

式中 表示数据维度。作为核函数结构的优化方法，AGP的特点在于其构建的核函数具有可解析的参数化形式，容易对模型进行解释。在回归问题中，AGP被认为具有与其它现代机器学习算法相当的预测能力。

高斯过程分类（Gaussian Process Classification, GPC）

GPC是使用高斯过程作为先验的分类器，在二元分类中，GPC在对潜函数进行估计后将其作为Sigmoid函数的输入得到分类概率；在多元分类中，则将Sigmoid函数替换为归一化指数函数。以二元分类为例，因为标签数据是取值的二项分布随机变量，所以对应的GPC似然是潜函数对学习样本的因子乘积：

GPC通过最大化似然求解超参数并得到后验，带入Sigmoid函数的表达式可知，上式不是正态分布，因此不同于GPR，GPC的后验没有解析形式，要求使用非正态似然的求解方法，例如解析近似和蒙特卡罗方法。

克里金法（Kriging）

主词条：克里金法

克里金法是与GPR相近的非参数模型，常见于随机场的插值问题。若协方差函数的形式等价，简单克里金（simple Kriging）和普通克里金（ordinary Kriging）的输出与GPR在正态似然下输出的数学期望相同，若克里金法使用高斯随机场假设，则给出的置信区间也与GPR相同。克里金法与GPR的不同点在于，前者假设随机场为固有平稳过程（intrinsically stationary process）并给出其对测试样本的最优无偏估计（Best Linear Unbiased Prediction, BLUP）；后者假设随机场为高斯过程并给出其对测试样本的完整后验。此外，考虑在地统计学中应用的常见情形，克里金法通常不加入噪声，其估计结果在学习样本处与学习目标完全相同。

有一些克里金法的扩展方法与GPR等价，例如Handcock and Stein (1993)和Handcock and Wallis (1994)提出的包含贝叶斯推断的克里金法使用了各向同性的马顿核高斯过程先验，并在正态似然假设下按MLE求解超参数。

在使用特定核函数时，GPR是一个通用近似（universal approximator），具体地，若实数域内有紧致空间 和赋范向量空间（normed vector space） ，则可以证明，对任意函数 ，等价核（Equivalent Kernel, EK）构成的再生核希尔伯特空间（Reproducing Kernel Hilbert Space, RKHS）内总是存在 能够按任意精度逼近原函数，即RKHS内总是存在线性组合： 。等价核的例子包括RBF核、二次有理函数核、指数函数核等。

作为非参数高斯过程模型的性质：GPR的复杂程度取决于学习样本，因此天然避免了过拟合（overfitting）问题。当测试样本与学习样本在特征空间的距离趋于无穷时，GPR的估计结果趋于其高斯过程先验，例如对0均值的各向同性高斯过程先验有： 。由于高斯过程可以通过协方差函数模拟随机变量间的相关性，因此GPR可应用于平移不变、旋转不变或周期性数据的回归问题，但同时，由于核函数/非参数方法的局限，GPR在高维数据的学习中表现不佳，该现象有时被称为“维数诅咒（curse of dimensionality）”。此外由理论部分可知，GPR的计算复杂度较高。一些研究给出了计算量更低的改进算法并取得了效果，但总体而言，GPR是一个为小样本设计的机器学习算法。

作为贝叶斯方法的性质：GPR是一个包含全贝叶斯特性（full Bayesian）的回归模型，可以给出测试数据的完整后验。此外，在似然为正态分布时，GPR给出的后验是具有闭合形式的联合正态分布，即GPR具有可解析性，该性质在非参数模型中并不多见。

与核学习方法的比较：GPR与适用于回归问题的核学习（kernel learning）方法，例如核岭回归、支持向量回归（support vector regression）等都使用了核函数与核方法，且可以通过相同的方式简化核矩阵求逆计算，但对后者，其核函数是输入空间向高维特征空间的映射，而GPR中的核函数是高斯过程协方差函数的模型，表示随机变量间的相关性。上述不同在求解中的体现是，GPR的“学习”是通过数据确定核函数超参数的过程，而核方法的“学习”，是在给定核函数超参数后求解模型权重或寻找支持向量的过程。

与人工神经网络（Artificial Neural Network, ANN）的关系：ANN是参数模型，在理论和结构上与GPR有较大差异，但二者在应用中有重叠，例子包括时间序列分析和计算机视觉问题。若学习样本充足，ANN由于使用随机梯度下降方法进行学习，在求解效率上优势明显；但在要求给出概率输出时，GPR是更合适的方法。ANN中的多层感知器（Multi-Layer Perceptron, MLP）与高斯过程存在联系，高斯过程可以由拥有单隐含层和无限个隐含层节点的MLP导出。具体地，对如下预测形式的MLP： 若其权重为iid随机变量且均值为0，方差为有限值，则由中心极限定理（central limit theorem）可推出，当隐含层节点数n趋于无穷时，模型输出的联合分布趋于0均值的联合正态分布： 该结论不要求iid权重服从正态分布，但当权重服从iid正态分布时，有限个隐含层节点的MLP也可表示为高斯过程。另一方面，MLP也可以由高斯过程导出，由Mercer定理可知，核函数可以表示为映射函数的内积： ，因此可视为以 为隐含层特征（输出）的MLP。

GPR可用于一般意义上的低维回归问题，尤其是时间序列数据的预测。例子包括太阳辐射的有关变量、风速、土壤温度、冒纳罗亚观测站（Mauna Loa Observatory，MLO）的全球对流层平均二氧化碳浓度等。在图像处理（image processing）方面，GPR被用于图像去噪和生成超分辨率图像。在自动控制方面，GPR被用于机械臂（robotic arm）数据的实时学习，也有研究开发了GPR的机器人学习（robot learning）系统。

包含GPR的编程模块

基于Python开发的机器学习模块scikit-learn提供封装的GPR工具GaussianProcessRegressor，其求解方法为BFGS算法（Broyden-Fletcher-Goldfarb-Shanno）的MLE，拥有自定义的求解器函数接口。在核函数方面，scikit-learn包含RBF核、马顿核、二次有理核、指数-周期核、多项式核以及自定义核函数的API工具。

由谢菲尔德机器学习团队（Sheffield machine learning group）开发的Python高斯过程建模工具GPy提供了GPR的完整求解系统和核函数，包括异方差噪声GPR、隐变量模型，解析近似、VBI和蒙特卡罗求解、稀疏高斯过程和低秩近似。此外有基于GPy的模型超参数优化工具GPyOpt可用。

pyGPs是基于Python开发的面向对象高斯过程建模应用，包含可视化学习界面，在功能上包括了常见核函数、稀疏高斯过程和解析近似。