贝叶斯优化 (Bayesian Optimization)

范叶亮 / 2020-06-06

分类: 机器学习, 最优化 / 标签: 高斯分布, Gaussian Distribution, 正态分布, Normal Distribution, 边缘化, Marginalization, 条件化, Conditioning, 高斯过程, Gaussian Processes, 高斯过程回归, Gaussian Processes Regression, 主动学习, Active Learning, 代理模型, Surrogate Model, 贝叶斯优化, Bayesian Optimization, 采集函数, Acquisition Functions / 字数: 5043

本文内容主要参考自：
1. 从高斯分布到高斯过程、高斯过程回归、贝叶斯优化
2. A Visual Exploration of Gaussian Processes
3. Gaussian Process Regression
4. Exploring Bayesian Optimization

一元高斯分布

若随机变量 X 服从一个均值为 μ，方差为 σ2 的高斯分布，则记为：

(1)X∼N(μ,σ2)

其概率密度函数为：

(2)f(x)=1σ2πe−(x−μ)22σ2

图片来源：https://zh.wikipedia.org/wiki/正态分布

二元高斯分布

若随机变量 X,Y 服从均值为 μ=(μX,μY)⊤，方差为 μ=(σX,σY)⊤ 的高斯分布，则记为：

(3)(X,Y)∼N(μ,σ)

其概率密度函数为：

(4)f(x,y)=12πσXσY1−ρ2e−12(1−ρ2)[(x−μX)2σX2+(y−μY)2σY2−2ρ(x−μX)(y−μX)σXσY]

其中，ρ 是 X 和 Y 之间的相关系数，σX>0 且 σY>0。

图片来源：Bayesian tracking of multiple point targets using expectation maximization

多元高斯分布

若 K 维随机向量 X=[X1,⋯,XK]⊤ 服从多元高斯分布，则必须满足如下三个等价条件：

1. 任何线性组合 Y=a1X1+⋯aKXK 均服从高斯分布。
2. 存在随机向量 Z=[Z1,⋯,ZL]⊤（每个元素服从独立标准高斯分布），向量 μ=[μ1,⋯,μK]⊤ 以及 K×L 的矩阵 A，满足 X=AZ+μ。
3. 存在 μ 和一个对称半正定矩阵 Σ 满足 X 的特征函数 ϕX(u;μ,Σ)=exp⁡(iμ⊤u−12u⊤Σu)

如果 Σ 是非奇异的，则概率密度函数为：

(5)f(x1,⋯,xk)=1(2π)k|Σ|e−12(x−μ)⊤Σ−1(x−μ)

其中 |Σ| 表示协方差矩阵的行列式。

边缘化和条件化

高斯分布具有一个优秀的代数性质，即在边缘化和条件化下是闭合的，也就是说从这些操作中获取的结果分布也是高斯的。边缘化（Marginalization）和条件化（Conditioning）都作用于原始分布的子集上：

(6)PX,Y=[XY]∼N(μ,Σ)=N([μXμY],[ΣXXΣXYΣYXΣYY])

其中，X 和 Y 表示原始随机变量的子集。

对于随机向量 X 和 Y 的高斯概率分布 P(X,Y)，其边缘概率分布为：

(7)X∼N(μX,ΣXX)Y∼N(μY,ΣYY)

X 和 Y 两个子集各自只依赖于 μ 和 Σ 中它们对应的值。因此从高斯分布中边缘化一个随机变量仅需从 μ 和 Σ 中舍弃相应的变量即可：

(8)pX(x)=∫ypX,Y(x,y)dy=∫ypX|Y(x|y)pY(y)dy

条件化可以用于得到一个变量在另一个变量条件下的概率分布：

(9)X|Y∼N(μX+ΣXYΣYY−1(Y−μY),ΣXX−ΣXYΣYY−1ΣYX)Y|X∼N(μY+ΣYXΣXX−1(X−μX),ΣYY−ΣYXΣXX−1ΣXY)

需要注意新的均值仅依赖于作为条件的变量，协方差矩阵和这个变量无关。

边缘化可以理解为在高斯分布的一个维度上的累加，条件化可以理解为在多元分布上切一刀从而获得一个维数更少的高斯分布，如下图所示：

高斯过程（Gaussian Process）是观测值出现在一个连续域（例如时间或空间）的随机过程。在高斯过程中，连续输入空间中每个点都是与一个正态分布的随机变量相关联。此外，这些随机变量的每个有限集合都有一个多元正态分布，换句话说它们的任意有限线性组合是一个正态分布。高斯过程的分布是所有那些（无限多个）随机变量的联合分布，正因如此，它是连续域（例如时间或空间）上函数的分布。

简单而言，高斯过程即为一系列随机变量，这些随机变量的任意有限集合均为一个多元高斯分布。从一元高斯分布到多元高斯分布相当于增加了空间维度，从高斯分布到高斯过程相当于引入了时间维度。一个高斯过程可以被均值函数 m(x) 和协方差函数 K(x,x′) 共同唯一确定：

(10)m(x)=E[f(x)]K(x,x′)=E[(f(x)−m(x))(f(x′)−m(x′))]

则高斯过程可以表示为：

(11)f(x)∼GP(m(x),K(x,x′))

均值函数决定了样本出现的整体位置，如果为零则表示以 y=0 为基准线。协方差函数描述了不同点之间的关系，从而可以利用输入的训练数据预测未知点的值。常用的协方差函数有：

• 常数：Kc(x,x′)=C
• 线性：KL(x,x′)=x⊤x′
• 高斯噪声：KGN(x,x′)=σ2δx,x′
• 指数平方：KSE(x,x′)=exp⁡(−|d|22ℓ2)
• Ornstein-Uhlenbeck：KOU(x,x′)=exp⁡(−|d|ℓ)
• Matérn：KMatern (x,x′)=21−νΓ(ν)(2ν|d|ℓ)νKν(2ν|d|ℓ)
• 周期：KP(x,x′)=exp⁡(−2sin2⁡(d2)ℓ2)
• 有理平方：KRQ(x,x′)=(1+|d|2)−α,α≥0

高斯过程回归

回归任务的目标是给定一个输入变量 x∈RD 预测一个或多个连续目标变量 y 的值。更确切的说，给定一个包含 N 个观测值的训练集 X={xn}1N 和对应的目标值 Y={yn}1N，回归的目标是对于一个新的 x 预测对应的 y。目标值和观测值之间通过一个映射进行关联：

(12)f:X→Y

在贝叶斯模型中，我们通过观测数据 D={(xn,yn)}n=1N 更新先验分布 P(Θ)。通过贝叶斯公式我们可以利用先验概率 P(Θ) 和似然函数 P(D|Θ) 推导出后验概率：

(13)p(Θ|D)=p(D|Θ)p(Θ)p(D)

其中 p(D) 为边际似然。在贝叶斯回归中我们不仅希望获得未知输入对应的预测值 y∗ ，还希望知道预测的不确定性。因此我们需要利用联合分布和边缘化模型参数 Θ 来构造预测分布：

(14)p(y∗|x∗,D)=∫p(y∗,Θ|x∗,D)dΘ=∫p(y∗|x∗,Θ,D)p(Θ|D)dΘ

通常情况下，由于积分形式 p(Θ|D) 不具有解析可解性（Analytically Tractable）：

(15)p(D)=∫p(D|Θ)p(Θ)dΘ

但在高斯似然和高斯过程先验的前提下，后验采用函数的高斯过程的形式，同时是解析可解的。

对于高斯过程回归，我们构建一个贝叶斯模型，首先定义函数输出的先验为一个高斯过程：

(16)p(f|X,θ)=N(0,K(X,X))

其中 K(⋅,⋅) 为协方差函数，θ 为过程的超参数。假设数据已经变换为零均值，因此我们不需要在先验中设置均值函数，则令似然形式如下：

(17)p(Y|f)∼N(f,σn2I)

假设观测值为独立同分布的高斯噪音的累加，则整个模型的联合分布为：

(18)p(Y,f|X,θ)=p(Y|f)p(f|X,θ)

虽然我们并不关心变量 f，但由于我们需要对不确定性进行建模，我们仍需考虑 Y 和 f 以及 f 和 X 之间的关系。高斯过程作为一个非参数模型，其先验分布构建于映射 f 之上，f 仅依赖于核函数的超参数 θ，且这些超参数可以通过数据进行估计。我们可以将超参数作为先验，即：

(19)p(Y,f|X,θ)=p(Y|f)p(f|X,θ)p(θ)

然后进行贝叶斯推断和模型选择，但是通常情况下这是不可解的。David MacKay 引入了一个利用最优化边际似然来估计贝叶斯平均的框架，即计算如下积分：

(20)p(Y|X,θ)=∫p(Y|f)p(f|X,θ)df

其中，高斯似然 p(Y|f) 表示模型拟合数据的程度，p(f|X,θ) 为高斯过程先验。经过边缘化后，Y 不在依赖于 f 而仅依赖于 θ。

假设采用零均值函数，对于一个高斯过程先验，我们仅需指定一个协方差函数。以指数平方协方差函数为例，选择一系列测试输入点 X∗，利用协方差矩阵和测试输入点可以生成一个高斯向量：

(21)f∗∼N(0,K(X∗,X∗))

从高斯先验中进行采样，我们首先需要利用标准正态来表示多元正态：

(22)f∗∼μ+BN(0,I)

其中，BB⊤=K(X∗,X∗)，B 本质上是协方差矩阵的平方根，可以通过 Cholesky 分解获得。

上图（左）为从高斯先验中采样的 10 个序列，上图（右）为先验的协方差。如果输入点 xn 和 xm 接近，则对应的 f(xn) 和 f(xm) 相比于不接近的点是强相关的。

我们关注的并不是这些随机的函数，而是如何将训练数据中的信息同先验进行合并。假设观测数据为 {(xi,fi)|i=1,…,n}，则训练目标 f 和测试目标 f∗ 之间的联合分布为：

(23)[ff∗]∼N(0,[K(X,X)K(X,X∗)K(X∗,X)K(X∗,X∗)])

根据观测值对联合高斯先验分布进行条件化处理可以得到高斯过程回归的关键预测方程：

(24)f∗|X,X∗,f∼N(f―∗,cov⁡(f∗))

(25)f―∗≜E[f∗|X,X∗,f]=K(X∗,X)K(X,X)−1fcov⁡(f∗)=K(X∗,X∗)−K(X∗,X)K(X,X)−1K(X,X∗)

函数值可以通过对联合后验分布采样获得。

我们以三角函数作为给定的函数，并随机采样一些训练数据 {(xi,fi)|i=1,…,n}，如下图所示：

我们希望将训练数据和高斯过程先验进行合并得到联合后验分布，我们可以通过在观测值上条件化联合高斯先验分布，预测的均值和协方差为：

(26)f―∗=K(X∗,X)K(X,X)−1fcov⁡(f∗)=K(X∗,X∗)−K(X∗,X)K(X,X)−1K(X,X∗)

Rasmussen 和 Williams 给出了一个实现高斯过程回归的实用方法：

\begin{algorithm}
\caption{高斯过程回归算法}
\begin{algorithmic}
\REQUIRE \\
    输入 $\mathbf{X}$ \\
    目标 $\mathbf{y}$ \\
    协方差函数 $k$ \\
    噪音水平 $\sigma^2_n$ \\
    测试输入 $\mathbf{x}_*$
\ENSURE \\
    均值 $\bar{f}_*$ \\
    方差 $\mathbb{V}\left[f_{*}\right]$
\FUNCTION{GaussianProcessRegression}{$\mathbf{X}, \mathbf{y}, k, \sigma^2_n, \mathbf{x}_*$}
\STATE $L \gets \text{cholesky} \left(K + \sigma^2_n I\right)$
\STATE $\alpha \gets L^{\top} \setminus \left(L \setminus \mathbf{y}\right)$
\STATE $\bar{f}_* \gets \mathbf{k}^{\top}_* \alpha$
\STATE $\mathbf{v} \gets L \setminus \mathbf{k}_*$
\STATE $\mathbb{V}\left[f_{*}\right] \gets k \left(\mathbf{x}_*, \mathbf{x}_*\right) - \mathbf{v}^{\top} \mathbf{v}$
\RETURN $\bar{f}_*, \mathbb{V}\left[f_{*}\right]$
\ENDFUNCTION
\end{algorithmic}
\end{algorithm}

高斯过程后验和采样的序列如下图所示：

先验的协方差矩阵和后验的协方差矩阵可视化如下图所示：

本小结代码请参见这里。

贝叶斯优化

在很多机器学习问题中，数据标注往往需要耗费很大成本。主动学习（Active Learning）在最大化模型准确率时最小化标注成本，例如对不确定性最高的数据进行标注。由于我们仅知道少量数据点，因此我们需要一个代理模型（Surrogate Model）来建模真正的模型。高斯过程因其灵活性和具有估计不确定性估计的特性不失为一个常用的代理模型。

在估计 f(x) 的过程中，我们希望最小化评估的次数，因此我们可以通过主动学习来“智能”地选择下一个评估的数据点。通过不断的选择具有最高不确定性的数据点来获得 f(x) 更准确的估计，直至收敛或达到停止条件。下图展示了利用主动学习估计真实数据分布的过程：

贝叶斯优化问题

贝叶斯优化的核心问题是：基于现有的已知情况，如果选择下一步评估的数据点？在主动学习中我们选择不确定性最大的点，但在贝叶斯优化中我们需要在探索不确定性区域（探索）和关注已知具有较优目标值的区域之间进行权衡（开发）。这种评价的依据称之为采集函数（Acquisition Functions），采集函数通过当前模型启发式的评估是否选择一个数据点。

贝叶斯优化的目标是找到一个函数 f:Rd↦R 最大值（或最小值）对应的位置 x∈Rd。为了解决这个问题，我们遵循如下算法：

1. 选择一个代理模型用于建模真实函数 f 和定义其先验。
2. 给定观测集合，利用贝叶斯公式获取后验。
3. 利用采集函数 α(x) 确性下一个采样点 xt=arg⁡maxxα(x)。
4. 将采样的点加入观测集合，重复步骤 2 直至收敛或达到停止条件。

• Probability of Improvement (PI)

Probability of Improvement (PI) 采集函数会选择具有最大可能性提高当前最大的 f(x+) 值的点作为下一个查询点，即：

(27)xt+1=arg⁡max(αPI(x))=arg⁡max(P(f(x))≥(f(x+)+ϵ))

其中，P(⋅) 表示概率，ϵ 为一个较小的正数，x+=arg⁡maxxi∈x1:tf(xi)，xi 为第 i 步查询点的位置。如果采用高斯过程作为代理模型，上式则转变为：

(28)xt+1=arg⁡maxxΦ(μt(x)−f(x+)−ϵσt(x))

其中，Φ(⋅) 表示标准正态分布累积分布函数。PI 利用 ϵ 来权衡探索和开发，增加 ϵ 的值会更加倾向进行探索。

• Expected Improvement (EI)

PI 仅关注了有多大的可能性能够提高，而没有关注能够提高多少。Expected Improvement (EI) 则会选择具有最大期望提高的点作为下一个查询点，即：

(29)xt+1=arg⁡minxE(‖ht+1(x)−f(x∗)‖|Dt)

其中，f 为真实函数，ht+1 为代理模型在 t+1 步的后验均值，Dt={(xi,f(xi))},∀x∈x1:t 为训练数据，x∗ 为 f 取得最大值的真实位置。

上式中我们希望选择能够最小化与最大目标值之间距离的点，由于我们并不知道真实函数 f，Mockus 1 提出了一种解决办法：

(30)xt+1=arg⁡maxxE(max{0,ht+1(x)−f(x+)}|Dt)

其中，f(x+) 为到目前为止遇见的最大函数值，如果采用高斯过程作为代理模型，上式则转变为：

(31)EI(x)={(μt(x)−f(x+)−ϵ)Φ(Z)+σt(x)ϕ(Z), if σt(x)>00 if σt(x)=0Z=μt(x)−f(x+)−ϵσt(x)

其中 Φ(⋅) 表示标准正态分布累积分布函数，ϕ(⋅) 表示标准正态分布概率密度函数。类似 PI，EI 也可以利用 ϵ 来权衡探索和开发，增加 ϵ 的值会更加倾向进行探索。

• 对比和其他采集函数

上图展示了在仅包含一个训练观测数据 (0.5,f(0.5)) 情况下不同点的采集函数值。可以看出 αEI 和 αPI 的最大值分别为 0.3 和 0.47。选择一个具有较小的 αPI 和一个较大的 αEI 的点可以理解为一个高的风险和高的回报。因此，当多个点具有相同的 αEI 时，我们应该优先选择具有较小风险（高 αPI）的点，类似的，当多个点具有相同的 αPI 时，我们应该优先选择具有较大回报（高 αEI）的点。

其他采集函数还有 Thompson Sampling 2，Upper Confidence Bound (UCB)，Gaussian Process Upper Confidence Bound (GP-UCB) 3，Entropy Search 4，Predictive Entropy Search 5 等，细节请参见原始论文或 A Tutorial on Bayesian Optimization 6。

← 马尔可夫决策过程 (Markov Decision Process) 利用动态规划求解马尔可夫决策过程 (Planning by Dynamic Programming) →