从数学到实现，全面回顾高斯过程当中的函数最优化

高斯过程能够被认为是一种机器学习算法，它利用点与点之间同质性的度量做为核函数，以从输入的训练数据预测未知点的值。本文从理论推导和实现详细地介绍了高斯过程，并在后面提供了用它来近似求未知函数最优解的方法。文章选自efavdb，做者： Jonathan Landy，机器之心编译。

咱们回顾了高斯过程（GP）拟合数据所需的数学和代码，最后得出一个经常使用应用的 demo——经过高斯过程搜索法快速实现函数最小化。下面的动图演示了这种方法的动态过程，其中红色的点是从红色曲线采样的样本。使用这些样本，咱们试图利用 GP 尽快找到曲线的最小值。

附录包含（i）高斯回归后验推导； (ii)SKLearn 的 GP 实现；(iii) GP 分类器的快速回顾。

前言

高斯过程（GP）是处理如下通常问题的一个工具：函数 f(x) 在 n 个点采样，并获得一组有噪声 [1] 的函数度量值 {f(xi)=y_i ± σ_i,i=1,…,n}。那么若给定这些可用的样本，且 f hat 为一些候选函数，咱们是否就能估计出 f =f hat 的几率?

为了逐步明确上述问题，咱们首先应用贝叶斯法则，

上式左边的数值是所求几率的简写，即给定样本函数值 {y} 的条件下 f=f hat 的几率。在上式的右边，分子中的第一项须要咱们对测量过程当中的偏差来源作一些假设，分子中的第二项是先验几率，在这里咱们必须采用最合理的假设。例如，咱们将在下面看到先验几率有效地决定了 f 函数在给定平滑度的几率。

在 GP 方法中，右式中两个分子都服从多元正态/高斯分布。模型能够选择高斯的具体参数来达到良好的拟合度，但特征的正态族假设对于解决数学问题是必不可少的。采用这种方法，咱们能够经过分析写出后验几率，而后用在一些应用和计算中。例如，咱们使用这种方法来得到文中最开始处图片上的曲线，即经过拟合 GP 后验几率中随机抽样而获得曲线估计，在两个收缩点处被固定为相等的测量值。后验样本对于可视化和求蒙特卡洛的平均值都颇有用。

在本文中，咱们作的工做有：

（i）回顾计算上述后验几率所需的数学运算；

（ii）讨论数值评估，并使用 GP 来拟合一些实例数据；

（iii）回顾拟合的 GP 如何快速最小化成本函数，例如机器学习中的交叉验证分。

附录包括高斯过程回归推导，SKLearn 的 GP 实现和 GP 分类器的快速回顾。

咱们的 GitHub 提供了简单的高斯过程示例：github.com/EFavDB/gaus…

注意：为了理解这篇文章中的数学细节，应该熟悉多元正态分布。但若是主要对应用感兴趣，能够忽略这些细节。

后验分析评估

为了计算 (1) 式左边的值，咱们要先计算右边的值。由于分母不依赖 f hat，咱们只须要考虑分子中的项。这意味着分母必须是全部候选函数共有的归一化因子。在本节中，咱们先将分子两项的估计公式写出来，而后考虑后验几率。

咱们要作的第一个假设是，假如实际函数是 f hat，那么咱们的测量值 y 关于 f hat 是独立的而且服从高斯分布。这个假设意味着方程 (1) 右边的第一项是：

上式中的 y_i 是咱们样本点的实际测量值，σ_i 方是它们的方差不肯定度。

第二个假设是，假设先验几率 p（f hat）的公式。咱们将注意力集中在一组数据点点 {x_i : i=1,…,N} 上，其中前 n 个点是已被抽样的点，剩下的（N-n）个点是其余位置的测试点，即咱们用来估计 f 函数联合统计数据的点。接下来，简单地假设这些点服从 f 的多元正态分布，由协方差矩阵Σ来控制，由此获得

这里，咱们简写 f_i≡f（x_i）。请注意，咱们已经默认上面正态分布的均值是零。这是为了简便起见：若是非零均值合适，那么能够与均值相加分析，或者从原始 f 中减去非零均值使新的分布均值为零。

Σ的特殊形式是在 GP 建模时最须要建模者观察力和首创性的地方。对研究主题很是熟悉的研究者能够构建很是好且很是复杂的先验几率，而这种先验一般是各项求和的形式，其中每一项都在所讨论问题的数据基础上加入了一些物理相关的贡献。在这篇文章中，咱们假设一个简单的公式，

注意，这个假设下，若是 x_i 和 x_j 很接近，那么指数近似等于 1。这确保了附近的点高度相关，从而使全部高几率函数变得平滑。当两测试点远离时，式 (4) 中的衰减速率由长度参数 l 控制。若是 l 很大（小），曲线将在一个很长（短）的距离上平滑。咱们将在下一节中说明这些问题，并在下下节中解释如何从已有的样本数据中推断合适的长度参数。

如今，若是咱们把式 (2) 和式 (3) 代入式 (1)，将获得后验几率 p(f1|{y}) 的表达式。这是一个指数函数，它的独立变量是 f_i 函数中的二次项。也能够说，与前验几率同样，后验几率服从多变量正态分布。经过简单计算，就能够获得这个分布均值和协方差的明确表达式：使用块（符号），其中 0 对应于采样点，1 对应于测试点，测试点的边缘分布是

y 是测量值的向量，长度为 n,

方程 (5) 是高斯过程回归的一个主要结果——有了这个结果，咱们就能够评估后验几率了。注意，在采样值 y 中全部点的均值是线性的，而且在测量值附近每一个点处的方差减少。若是你有兴趣仔细推导这个结果，能够参考咱们的附录，在那里有两个推导。可是，在接下来的正文中，咱们仅简单地探讨这个公式的应用。

后验几率的数值计算

在本节中，咱们将介绍式 (5) 的两个典型应用：（i）在测试点 x 处评估后验分布的均值和标准差，（ii）从后验几率中直接采样函数 f_hat。前者能够得到 f 函数在全部位置的置信区间，然后者能够用来实现可视化和从后验几率中得到通常的蒙特卡洛平均值。这两个概念都在这篇文章的标题图片中进行了说明：图片中，咱们将 GP 拟合一个已有两个测量点的一维函数。蓝色阴影区域表示每一个位置函数值的一个σ置信区间，彩色曲线是后验样本。

咱们的 SimpleGP fitter 类的代码能够在 GitHub 上找到。咱们将在下文中解释他是如何工做的，但若是对细节感兴趣应该仔细查阅代码。

区间

下面的代码对咱们的 SimpleGP 类进行了初始化，定义了一些样本位置、样本值和不肯定性，而后评估了一组测试点后验几率的均值和标准差。简而言之，这个过程以下：经过拟合评估了出如今式（5）中的逆矩阵

，并保留结果供之后使用，这能够避免在每一个测试点中从新评估这个逆矩阵。接下来，经过调用区间，针对每一个测试点再次评估式 (5)。

# Initialize fitter -- set covariance parameters
WIDTH_SCALE = 1.0
LENGTH_SCALE = 1.0
model = SimpleGP(WIDTH_SCALE, LENGTH_SCALE, noise=0) 

# Insert observed sample data here, fit 
sample_x = [-0.5, 2.5]
sample_y = [.5, 0]
sample_s = [0.01, 0.25]
model.fit(sample_x, sample_y, sample_s) 

# Get the mean and std at each point in x_test
test_x = np.arange(-5, 5, .05)
means, stds = model.interval(test_x)复制代码

在以上代码块中，WIDTH_SCALE 和 LENGTH_SCALE 用来指定协方差矩阵式（4）。前者对应于等式中的σ，后者对应于 l。增长 WIDTH_SCALE 对应于未知函数大小不肯定度的增长，增长 LENGTH_SCALE 对应于增长咱们指望的函数平滑程度。下图说明这些问题：这里，经过设置 WIDTH_SCALE = LENGTH_SCALE = 1 得到蓝色区间，经过设置 WIDTH_SCALE = 0.5 和 LENGTH_SCALE = 2 得到橙色区间。结果是橙色相对蓝色后验估计更加紧密平滑。在这两幅图中，实曲线表示后验分布均值，竖线表示一个σ置信区间。

后验采样

为了从后验几率中采样实际函数，咱们将再次简单地评估式 (5) 中的均值和协方差矩阵，此次是对咱们所求采样函数的多个测试点进行。一旦咱们有了这些测试点后验几率的均值和协方差矩阵，咱们可使用多元正态采样的外部库从 (5) 中抽取样本——为此，咱们使用了 python 中 numpy。下面代码的最后一步执行这些步骤。

# Insert observed sample data here. 
sample_x = [-1.5, -0.5, 0.7, 1.4, 2.5, 3.0]
sample_y = [1, 2, 2, .5, 0, 0.5]
sample_s = [0.01, 0.25, 0.5, 0.01, 0.3, 0.01]

# Initialize fitter -- set covariance parameters
WIDTH_SCALE = 1.0
LENGTH_SCALE = 1.0
model = SimpleGP(WIDTH_SCALE, LENGTH_SCALE, noise=0)
model.fit(sample_x, sample_y, sample_s)

# Get the mean and std at each point in test_x
test_x = np.arange(-5, 5, .05)
means, stds = model.interval(test_x)

# Sample here
SAMPLES = 10
samples = model.sample(test_x, SAMPLES)
复制代码

注意在第 2-4 行中，咱们添加了一些附加的函数样本点（为了好玩）。最后的区间和后验样本以下图所示。注意到在采样点附近，后验结果很是理想。然而，在图的左侧，一旦咱们移动的距离≥1（即协方差矩阵 (4) 中的长度参数），后验就趋近于先验。

选择协方差超参数

以前，咱们证实了咱们的协方差的长度参数显着地影响后验几率的区间形状以及其中的样本。适当设置这些参数是使用 GP 的一个广泛难点。在这里，咱们描述两种方法，能够巧妙地设置超参数，并给出一些采样数据。

交叉验证

交叉验证是设置超参数的标准方法。这须要将可用的样本数据分为训练集和验证集。训练集经过一系列超参数进行 GP 拟合，而后在已有的验证集上评估模型的准确性。而后，经过选择不一样的超参数重复这个过程，选择可使验证集表现最优的一组。

边缘似然最大化

一般状况下，人们倾向于将 GP 应用于评估样本有较高成本的状况。这意味着人们一般在只有少数样本可用的状况下使用 GP。这种状况下，随着训练点数量的增长，最优超参数能够快速变化，意味着经过交叉验证获得的最优选择可能远不如训练一个完整样本集获得的最优集合 [3]。

设置超参数的另外一种常见方法是使边缘似然最大化。这就是说，咱们试图最大化已关察样本的可能性，于是在可用超参数的基础上进行优化。具体来讲，就是经过对未知的 f hat 进行积分来评估边缘似然 [4]。

咱们能够像附录中评估后验分布那样直接进行积分。但更快的方法是注意到 f 积分后，y 值服从以下的正态分布

其中σ^2 * I_00 在式（6）中定义，由此得出，

上面两项是矛盾的：第二项经过找出使指数最大的协方差矩阵来减少，最大化指数使数据过分拟合。然而，第一项与第二项相反，第一项是高斯积分的归一化因子，它随着衰减长度变短和对角线误差下降而变大，至关于抑制复杂度的正则项。

实际中，为了使式 (10) 最大，一般利用梯度的解析表达式和梯度降低法，这是 SKLearn 采起的方法。模型的一个优势是可以优化 GP 的超参数。然而，式（10）不必定是凸的，一般存在多个局部最小值。为了得到一个好的最小值，能够尝试在一些优秀的初始点上进行初始化。或者能够在随机点重复初始化梯度降低，最后选择最优解。

函数最小搜索和机器学习

如今咱们将介绍 GP 的一个经常使用的应用：快速地搜索函数最小值。在这个问题中，咱们能够迭代得到函数的噪声样本，从而尽快识别函数的全局最小值。梯度降低能够应用于这种状况，可是若是函数不具有凸性，一般须要重复采样。为了减小所需的步骤/样本的数量，能够尝试应用更通常的探索式策略，即平衡「优化当前已知最小值的目标」与「寻找可能更小的新局部最小值的目标」。GP 后验为开发这样的策略的提供了一个自然的起点。

GP 搜索法的想法是在 GP 后验的基础上得到一个得分函数。这个得分函数用来对搜索给定点的信息进行编码，它能够对探索（explore）和利用（exploit）造成一种权衡。一旦每一个点都进行评分，那么具备最大（或最小，最合适的）分数的点将会被采样。而后迭代重复该过程直到找到一个符合要求的解为止。咱们将在下面讨论四种可能的选择，并给出一个例子。

高斯置信下界（GLCB）

GLCB 在每点的评分方式为

这里，μ和σ是函数在 x 处的均值和标准差的 GP 后验估计值，κ是控制参数。请注意，第一项μ（x）鼓励利用最可靠的局部最小值，并在它的周围执行搜索。相似地，第二项κσ鼓励在当前 GP 最不肯定真实函数值的点上进行探索。

改进的高斯几率（GPI）

若是目前为止所看到的最小值是 y，则能够利用该点处的真实函数值小于 y 的几率来给每一个点评分。也就是说，咱们能够写为

高斯预期改进（EI）

上式常见的变形叫作预期改进，定义为

这个得分函数倾向于鼓励更多地去探索而不是改善几率，由于它更重视不肯定性。

几率最小值

要获得的最终得分函数是问题中最小值的几率。得到这个分数的一个方法是进行屡次后验采样。对于每一个样本，首先标记它的全局最小值，而后采起多数投票方法来决定接下来的样本。

本文最开始处的动图展现了一个实际的 GP 搜索，使用 skopt[5] 在 python 中执行。左边的红色曲线是正在寻找全局最小值的（隐藏）曲线 f。红点是目前已经得到的样本，绿色阴影曲线是每一个点的 GP 后验置信区间，该置信区间会随着更多样本的得到而逐渐改进。右边是经过在 GP 后验基础上分析获得的每点的预期改进（EI）得分函数——该例中用于指导搜索的得分函数。该过程用五个随机样本进行初始化，而后进行引导搜索。请注意，随着过程的演变，前几个样本集中利用已知的局部最小值。但通过几回迭代后，继续对这些位置采样的收益减少，探索中间点的需求占了上风——中间点是发现的实际全局最小值的点。

讨论

在这篇文章中，咱们归纳了大部分 GP 的数学运算：获得后验几率所需的数学，如何进行后验采样，最后讨论了如何实际应用后验几率。

总的来讲，GP 表明了一个能够拟合任何函数的强大工具。实际中，使用这个工具的挑战主要在于合适超参数的选择，寻找合适的参数常常被困在局部最小，使拟合失效。不过，若是选择得当，GP 的应用能够提供一些有价值的性能提高。

附录中讨论了关于 GP 的其余话题。若是对更多的细节感兴趣，咱们能够推荐 Rasmussen 和 Williams 的免费在线文本 [6]。

附录 A：后验推导

本附录中，咱们提出后验推导（5）的两种方法。

方法 1

先平方，结合式（2）和式（3），简单地计算得出

这里，(1/σ^2) * I 在式（6）中被定义，但在样本集外的全部行中都为零。为了获得式（5），咱们必须统一为正文中逆矩阵的分块结构。

首先，咱们能够得出

这里咱们使用了分块表示法。为了计算上述的逆矩阵，咱们将利用分块矩阵求逆公式，

矩阵（A2）中块 C = 0、D = 1，这大大简化了上述过程。代入后获得

利用这个结果和（A1），咱们能够获得测试集的平均值

其中第二行的分子分母同时乘了 (1/σ^2) * I_00 的逆。相似地，测试集的协方差由（A3）的右下块给出。获得，

由（A4）和（A5）得出式（5）。

方法 2

在第二种方法中，咱们考虑一组测试点 f_1 和一组观测样本 f_0 的联合分布，咱们再次假设密度函数的均值为零。那么二者的联合几率密度是

如今，咱们利用结果

式（5）的得出须要利用上述两个表达式。主要的困难在于平方，相似于以前的推导，这能够经过分块矩阵求逆公式来完成。

附录 B：SKLearn 实现和其余内核

SKLearn 提供了包含 GaussianProcessRegressor 的类。使得能够在任何维度上进行拟合和采样——即它比咱们的最小类更通常，由于它能够在多维上拟合特征向量。另外，SKLearn 类的拟合方法尝试为一组给定的数据找到一组最佳的超参数。如上所述，这是经过边缘似然的最大化来完成的。这里，咱们提供一些关于这个类的基本注释和能够用来定义（3）中协方差矩阵Σ的内核函数，以及一段说明调用的简单代码。

预约义的核函数

• 径向基函数（RBF）：这是默认值，至关于式（4）。RBF 由一个尺度参数 l 表征，多维状况下，能够是一个容许各向异性相关长度的向量。

• White kernel：The White Kernel 用于噪声估计——文档建议用于估计全局噪声水平，但不是逐点。

• Matern：这是一个广义的指数衰减，其中指数是分离距离的幂律。特殊限制包括 RBF 和绝对距离指数衰减。

• 有理二次方程：（1+（d / l）2）α。

• Exp-Sine-Squared：它容许模拟周期性函数。相似于 RBF，但其距离是实际距离的正弦。存在周期性参数和「方差」——高斯抑制（Gaussian suppression）的尺度。

• 点积核函数：格式为 1 +xi⋅xj。它不是稳定的，就是说若是加入一个常量的平移，结果就会改变。若是把 N（0,1）的先验值放在系数上，将获得线性回归分析的结果。

• 核函数做为对象：能够支持核函数之间的二进制操做以建立更复杂的核函数，例如加法、乘法和指数（后者只是将初始核函数提高为幂）。你能够经过一些辅助函数来访问核函数中的全部参数，例如 kernel.get_params().kernel.hyperparameters 是全部超参数的列表。
  

参数

• n_restarts_optimizer：从新拟合的次数，用于探索多个局部最小值，默认值是零。
  

• alpha：这个可选参数容许每一个测量都传递不肯定性。
  

• normalize_y：用来表示咱们正在寻找的 y 的平均值不必定是零。
  

调用示例

下面的代码进行了一次简单拟合，结果是本文最开始展现的图片。

from sklearn.gaussian_process.kernels import RBF, ConstantKernel as C
from sklearn.gaussian_process import GaussianProcessRegressor
import numpy as np

# Build a model
kernel = C(1.0, (1e-3, 1e3)) * RBF(10, (0.5, 2))
gp = GaussianProcessRegressor(kernel=kernel, n_restarts_optimizer=9)

# Some data
xobs = np.array([[1], [1.5], [-3]])
yobs = np.array([3, 0, 1])

# Fit the model to the data (optimize hyper parameters)
gp.fit(xobs, yobs)

# Plot points and predictions
x_set = np.arange(-6, 6, 0.1)
x_set = np.array([[i] for i in x_set])
means, sigmas = gp.predict(x_set, return_std=True)

plt.figure(figsize=(8, 5))
plt.errorbar(x_set, means, yerr=sigmas, alpha=0.5)
plt.plot(x_set, means, 'g', linewidth=4)
colors = ['g', 'r', 'b', 'k']
for c in colors:
y_set = gp.sample_y(x_set, random_state=np.random.randint(1000))
 plt.plot(x_set, y_set, c + '--', alpha=0.5)
复制代码

关于 sklearn 实现的更多细节能够在这里找到：scikit-learn.org/stable/modu…

附录 C：GP 分类器

这里，咱们说明一般 GP 如何被用来拟合二进制类型数据，响应变量 y 能够取 0 或 1 的数据。GP 分类器的数学运算不像 GP 回归那样清楚。由于 0/1 响应不是高斯分布的，意味着后验几率也不是。为了利用该程序，能够经过拉普拉斯（Laplace）近似正常地对后验几率近似。

首先写这样一个公式，表示在 x 处看到一个给定的 y 值的几率。具体以下，

这个公式是 logistic 回归的一个正常非线性泛化。此外，f 的先验几率再次获得等式（3）。使用此式和（A8），咱们能够获得 f 的后验几率

利用这个公式，能够很容易地从近似后验中得到置信区间和样本，相似于回归。

原文连接：efavdb.com/gaussian-pr…