这是以 “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、其它涉及昂贵人类输入的场景、机器学习中的自动超参数搜索
参数模型上的贝叶斯优化
所谓参数模型，就是我们可以显式地写出 f 的结构且 f 可以由一组参数 w 决定。
Beta-Bernoulli 先验模型
Beta-Bernoulli 是最简单的参数模型贝叶斯优化问题。假设有 K 种药物，每种的治疗成功率未知，我们希望一边治疗病人，一边找出最有效的药物（找出 f 的最大值）。这是一个多臂老虎机问题（multi-armed bandit problem）。
形式上，多臂老虎机可以描述为输入一个指标，给出一个对应的随机变量。在这个问题中该随机变量是 Bernoulli 变量（成功/失败），设概率为标量 w_a ∈ (0, 1)。那么这个问题中所有不可见参数就是 w ∈ (0, 1)^K，f 就是 a → w_a。
随着实验的进行，我们得到越来越多的 (a_i, y_i) 对，a_i 是药品下标，y_i 是实验结果。我们可以通过观测到的数据计算 w 的后验分布。贝叶斯统计理论中，后验是先验乘以似然（再除以边缘似然），所以我们要先选择先验。为了后续计算的简便，我们选择 Beta 分布，因为它是 Bernoulli 分布的共轭先验 (Conjugate Prior)【附录A】：
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 维特征向量，则我们假设 f 为 f(a) = f(x_a) = x_a^T w。实际上观测得到的是均值为 f(a)，方差为 sigma^2 的高斯分布样本。那么这个问题中所有不可见参数就是 w 和 sigma^2。
（线性模型内在的逻辑是，只要估计出了 w 的分布，就能立刻选出最优化 x_a^T w 的下标 a）
已知高斯分布的共轭先验分布是 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 是固定的已知量；假设先验服从 N(0, V_0)。我们将观测值服从的分布用不含 w 的式子表示：
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