这是以 “paper:” 作为标题前缀的论文速读系列。最初的动机方法是参考多份资料学习并写出一份贝叶斯优化入门笔记，结果发现了这篇 paper，尤其适合阅读来入门，于是本学习笔记就以翻译和注解这篇论文的前半部分，也是我能勉强看懂的部分为主了。
Taking the Human Out of the Loop: A Review of Bayesian Optimization
【内容不完整警告】本文目前并不是一个对贝叶斯优化的完整入门笔记
动机和问题描述
希望求出 f(X) 的全局最大值，其中 X 通常不高于20个维度，f 的结构未知，只能求值，而每次求值成本高昂。此外，在任何取值点上，f(X) 还是一个随机变量，每次“求值” 只是从分布中采样一次。
整体求解思路为，首先选择一个 f 的先验概率，然后不断 {决定求值位置、求值、依据观测值得到后验概率 }
形式化地，算法可以描述为
for n=1,2,... do
select new x_{n+1} by optimizing acquisition function alpha: x_{n+1} = argmax_x alpha(x; D_n)
evaluate function at x_{n+1}, value is y_{n+1}
D_{n+1} = D_n.add( (x_{n+1}, y_{n+1}) )
update statistical model
end for常见的实际应用：A/B Testing、其它涉及昂贵人类输入的场景、机器学习中的自动超参数搜索
贝叶斯学派和频率学派
统计学中，贝叶斯学派和频率学派给出了两种截然不同的看待问题的方式。贝叶斯优化作为一个随机性、统计学的框架，继承了贝叶斯学派的视角，因此有必要回顾一下统计问题的“贝叶斯视角”。
频率学派视角相比之下比较直白易懂，所以先聊这个。这个视角下，未知参数通常具有一个确定的数值。我们对这个数值做出断言，然后根据实验结果计算这个断言正确的概率。典型的例子包括：
- 置信区间理论。最朴素的例子，方差已知的随机变量具有一个确定但未知的期望。拿到一组随机变量的观测值，不难断言均值偏离期望一定距离以内的概率大于 95%。这种情况下，我们得到了期望的 95% 水平置信区间。
- 统计假设检验方法。例如“女士品茶”中，提出女士没有能力判断出先加茶还是先加牛奶的原假设，选择“八杯都正确”这一显著性水平为 1.43% 的拒绝域，实验结果落在拒绝域中，因此以 1.43% 的显著性水平拒绝原假设。
贝叶斯视角下，未知参数往往自身就是随机变量。先验、后验分布表示我们对随机变量的信念度，通过贝叶斯公式，可以根据观测结果来更新我们的信念。
参数模型上的贝叶斯优化
所谓参数模型，就是我们可以显式地写出 f 的结构且 f 可以由一组参数 w 决定。此外，在贝叶斯优化的范畴中，f 并不是一个确定性的函数。
概念上的对应
首先，我们将贝叶斯统计学中概念对应到参数化模型的描述中。前面所说的“先验概率”和“后验概率”其实是我们在不同阶段相信参数 w 服从的概率分布。贝叶斯统计学中的“似然函数”是当我们给定一组参数 w 和输入 x 的确切数值时，f(x) 的概率分布函数。
Beta-Bernoulli 先验模型
Beta-Bernoulli 是最简单的参数模型贝叶斯优化问题。假设有 K 种药物，每种的治疗成功率未知，我们希望一边治疗病人，一边找出最有效的药物。这是一个多臂老虎机问题（multi-armed bandit problem）。
形式上，多臂老虎机可以描述为输入一个指标，给出一个对应的随机变量。在这个问题中该随机变量是 Bernoulli 变量（成功/失败），设概率为标量 w_a ∈ (0, 1)。那么这个问题中所有不可见参数就是 w ∈ (0, 1)^K，输入是下标 a。这个问题中，f(a) = Ber(w_a)。
随着实验的进行，我们得到越来越多的 (a_i, y_i) 对，a_i 是药品下标，y_i 是实验结果。我们可以通过观测到的数据计算 w 的后验分布。贝叶斯统计理论中，后验是先验乘以似然（再除以边缘似然）。既然似然是 Bernoulli 分布，我们希望选择相信 w 服从的分布满足好的性质：即先验和后验在同一个分布家族，方便我们在每个观测点更新我们的 belief。这样良好性质的分布称为 Bernoulli 分布的共轭先验 (Conjugate Prior)。在【附录A】中可以证明 Beta 分布是 Bernoulli 分布的共轭先验。到此为止，我们有
P(w | alpha, beta) = Pi_{a=1}^K Beta(w_a | alpha, beta)
where the pdf of Beta is x^{alpha-1} (1-x)^{beta-1} / B(alpha, beta)
通过【附录A】，不难验证后验概率恰好可以解析表示为
P(w | D) = Pi_{a=1}^K Beta(w_a | alpha+n_{a,1}, beta+n{a,0})
where n_{a,1} counts number of successes of index a, n_{a,0} counts number of failures of index a
可以看出，这个情形下先验的 alpha, beta 可以看成虚空成功和失败次数（pseudo-count）。
接下来，我们需要定义采集函数（acquisition function）。例如，我们可以使用 Thompson Sampling，这个策略让样本在后验概率上最优。常常需要使用蒙特卡洛采样来估计采样函数，具体而言，这里我们依据后验概率采样一个 w 出来，然后选择 w 中成功概率最高的药品喂给病人。
总结而言，这个算法和蒙特卡洛树搜索的子节点选择策略大同小异，只不过采集函数不同。
线性模型
Beta-Bernoulli 模型的问题是老虎机的手臂数指数爆炸：手臂数目等于每个特征选项数的连乘。对于有多维特征的 X，我们转而提出线性模型：首先假设有容易正向、反向求值的映射 X → x_a，x_a 是表示 X 的 d 维特征向量，只要估计出了 w 的分布，要能立刻选出最优化 x_a^T w 的下标 a。然后，我们假设 f(a) = N(x_a^T w, sigma^2)。这个问题中所有不可见参数就是 w 和 sigma^2。
根据”概念上的对应”一节的分析，这个问题中的似然函数是高斯分布。已知高斯分布的共轭先验分布是 Normal-inverse-gamma (NIG) 函数。有四个需要输入的先验参数：w_0, V_0, alpha_0, beta_0。它的密度函数定义为：
NIG(w, sigma^2 | w_0, V_0, alpha_0, beta_0) = |2*pi*sigma^2*V_0|^{-1/2} * exp(-\frac{1}{2*sigma^2}(w-w_0)V_0^{-1}(w-w_0)) * \frac{beta_0^alpha_0}{Gamma(alpha_0)(sigma^2)^{alpha_0+1}}* exp(-beta_0/sigma^2)把 n 次观测值的输入记为 n x d 的矩阵 X, 输出记为 n 维向量 y，那么后验概率分布同属 NIG 家族，只不过要将其中的四个参数替换为
w_n = V_n(V_0^{-1}w_0 + X^T y)
V_n = (V_0^{-1} + X^T X)^{-1}
alpha_n = alpha_0 + n/2
beta_n = beta_0 + (1/2)(w_0^T V_0^{-1} w_0 + y^T y - w_n^T V_n^{-1} w_n)仍然用 Thompson Sampling 来采集样本：
a_{n+1} = argmax_a x_a^T \hat{w}, where \hat{w} is sampled from p(w | D_n)
可以发现，这种方法也适用于 X 是连续空间的情形。
还可以通过非线性基函数来扩展线性模型：f(x) = ɸ(x)^T w，比如径向基函数（Radial Basis Function, RBF）: ɸ_j(x) = exp(-(1/2) (x-z_j)^T Gamma (x-z_j))。在计算时，只需将所有的矩阵 X 替换为 Φ(X) 即可。
通用线性模型 Generalized Linear Model, GLM
通用线性模型（GLM）允许实数值以外的观测结果，比如 A/B testing 中的二元选择。方法是引入一个观测空间到实数值的链接映射 (link function) g(·)。为了定义这个函数，通常考虑其反函数 g^{-1}，定义实数值如何决定观测结果。一种常见的选择是 sigmoid 函数：g^{-1}(z) = 1/(1+e^z)。
不幸的是，我们找不到这种观测值对应的共轭先验，因此需要近似。Markov chain Monte Carlo (MCMC) 提供了是一种近似方法。
讨论：如何评价算法的好坏？
上面的问题模型中，尽管目标之一是 “求出 f(X) 的全局最大值”，我们似乎提出了不同的评价算法好坏的方式。最朴素的评价方式似乎是算法输出的值离真正最大值又多近，但它不知为何没有被采用。
目标可能是 观测值和观测结果序列上的效用函数，如 Beta-Bernoulli 模型的累计病人治愈数；更准确地，在多臂老虎机情形中，是指数递减的加权和（Sigma_i gamma^{i-1} y_i）
目标可能是 尽可能确定参数 w 的值，那么形式上就是最小化 w 的后验分布的微分熵 (differential entropy)，这个目标被称为 D-optimality utility。
非参数化模型的贝叶斯优化
首先，我们建立为什么能进行非参数模型贝叶斯优化的基础，其实就是应用 kernel trick 的可能。
假设观测值服从 N(x^T w, sigma^2)，其中 sigma^2 是固定的已知量；假设 w 服从分布 N(0, V_0)。我们将观测值服从的分布用不含 w 的式子表示（注意到变量仅剩下 X）：
p(y | X, sigma^2) = \int N(Xw, sigma^2 I) N(0, V_0)
= N(0, X V_0 X^T + sigma^2 I)
// the convolution of two gaussians (the first is y|w, the second is w) is another gaussian如果对 x 用基函数变为 ɸ(x)，则记 Φ = Φ(X)，上面的式子变为
p(y | X, sigma^2) = N(0,ΦV_0Φ^T + sigma^2 I)其中，ΦV_0Φ^T 是一个元素为 ɸ(x_i)^T ɸ(x_j) = k(x_i, x_j) 的矩阵，由核函数 k 决定。
常见先验：高斯过程
高斯过程是一种随机过程
常见的先验是高斯过程 (Gaussian Process)。高斯过程是一种随机过程（Stochastic Process），随机过程是一个依据指标（时间、空间）返回随机变量的函数。当用任意有限指标集索引得到的随机变量联合概率分布都是多变量高斯分布时，称这个随机过程为高斯过程。
一个示例的高斯过程：X_t = cos(at) E_1 + sin(at) E_2, where E_1, E_2 i.i.d N(0, 1)。
一个重要性质是，期望为 0 的高斯过程可以被协方差函数完全确定下来。具体而言，对于 x, x’ 两个下标，高斯过程给出了 X, X’ 两个随机变量，那么 X, X’ 形成二元高斯分布。由于期望为 0，确定这个分布，仅需要协方差矩阵：
Var[X], Cov(X, X')
Cov(X, X'), Var[Y]由于 Var[X] = Cov(X, X)，我们进而只需要协方差（标量）函数 Cov(X, X’) = E[(X-E[X])(X’-E[X’])]。协方差函数需要满足一些有效条件（Admissibility），保证组合出来的随机变量方差非负。一个示例的协方差函数：Cov(X, X') = X^T X'。协方差函数可以类比 SVM 中的核函数。
高斯过程中的后验分布
我们希望得到任意 x 对应的随机变量 f(x) 的后验分布。为此，利用高斯过程的定义，考察 n 次测试 x_{1:n} 和任意一个求值点 x 的联合分布：
[y,\\ f(x)] = N( [m,\\ u_0(x)], [K+sigma^2 I, k(x_{1:n},x),\\ k(x,x_{1;n}, k(x,x))])协方差矩阵的左上角之所以是 K + sigma^2 I，是先验概率的协方差矩阵加上观测过程造成的方差（回忆 sigma^2 是固定的已知量）。注意这里的协方差矩阵已经被替换为核函数，因为任意一个高斯过程的协方差函数都可以看成一个核函数，而任意选定的有效核函数都可以作为某种高斯过程的协方差函数，两者是同义的。
我们使用多元高斯分布中利用联合分布求条件分布的公式 【附录 B】，得到 f(x) 的期望和方差：
f(x) | y ~ N(mean, variance), where
mean = u_0(x) + k(x,x_{1;n})(K+sigma^2 I)^{-1}(y-m)
variance = k(x,x) - k(x,x_{1;n})(K+sigma^2 I)^{-1}k(x_{1;n},x)这个条件分布，实际上就是后验分布！换言之，当高斯过程的方差固定已知时，高斯分布的共轭先验仍然是高斯分布。
常见的核函数（协方差函数）
常见的核函数通常有一些比较好的性质，比如它导出的高斯过程是平稳过程（Stationary Process，即时间、空间平移后联合分布函数相等）。
Matérn 核函数（亦称为 Matérn 协方差函数）是一种常见的核函数，最简单的情形下可以退化为随半径指数递减函数。
采样函数
TBD
附录
A. Beta 分布是 Bernoulli 分布的共轭先验
我们只考虑第一次实验前后，假设第一次实验选择的下标为 a。容易看出，先验、后验概率分布关于每种药品是独立的，所以可以不做累乘，而且下面的 alpha, beta 可以是对 a 而言的。
如果成功，得到 (a, 1), 则
P(w_a | {(a,1)}) ∝ P({(a,1)} | w_a) P(w_a) ∝ x * x^{alpha-1} (1-x)^{beta-1} = x^{alpha} (1-x)^{beta-1}
=> P(w_a | {(a,1)}) ~ Beta(alpha+1, beta)
如果失败，同理分析。
在这个例子中，似然函数是 Bernoulli 分布，如果先验概率是 Beta 分布，后验概率恰好也是 Beta 分布。这种情况，称 Beta 分布是 Bernoulli 分布的共轭先验（Conjugate Prior）。
B. 多元正态分布的条件分布
假设两个多维随机变量 (x1,\\ x2) 联合分布是期望为 (u1,\\ u2)，协方差矩阵为 (Sigma_11, Sigma_12, \\ Sigma_21, Sigma_22) 的正态分布，则条件分布 x2 | x1 = a 也是正太分布，其期望为 u1 + Sigma_12 Sigma_22^{-1} (a-u2), 协方差矩阵为 Sigma_11 – Sigma_12 Sigma_22^{-1} Sigma_21。
Leave a Reply