参考资料：
Expected Improvement formula for Bayesian Optimisation
通俗科普文：贝叶斯优化与SMBO、高斯过程回归、TPE
理解贝叶斯优化
 A Tutorial on Bayesian Optimization

贝叶斯优化是一种求解函数最优值的算法，它最普遍的使用场景是在机器学习过程中对超参数进行调优。贝叶斯优化算法的核心框架是SMBO (Sequential Model-Based Optimization)，而贝叶斯优化（Bayesian Optimization）狭义上特指代理模型为高斯过程回归模型的SMBO。

问题介绍

\[\max_{x\in A}f(x) \]

输入：\(x\in R^d\)，一般\(d\le20\)，即参数在20个以内。
可行集（feasible set）\(A\)： 一般是超多边形（\(x\in R^d: a_i\le x_i \le b_i\)）或d维的simplex（\(x\in R^d: \sum_ix_i=1\)）
在使用贝叶斯优化时，目标函数\(f\)：
- 需要是连续函数
- 评估成本昂贵，能允许执行的次数有限
- 是黑盒，公式/特征未知，例如：不是凹函数、线性函数之类的。
- 因为不对函数做任何假设，所以只能访问f，不能访问其一阶、二阶导数。
- 可以有噪声，但需要假设每个观测值的噪声是独立同分布，均值为0，方差恒定。

SMBO (Sequential Model-Based Optimization)

SMBO是一套优化框架，也是贝叶斯优化所使用的核心框架。它有两个重要组成部分：

一个代理模型（surrogate model），用于对目标函数进行建模。代理模型通常有确定的公式或者能计算梯度，又或者有已知的凹凸性、线性等特性，总之就是更容易用于优化。更泛化地讲，其实它就是一个学习模型，输入是所有观测到的函数值点，训练后可以在给定任意\(x\)的情况下给出对\(f(x)\)的估计。
一个优化策略（optimization strategy），决定下一个采样点的位置，即下一步应在哪个输入\(x\)处观测函数值\(f(x)\)。通常它是通过采集函数（acquisition function） 来实现的：
采集函数通常是一个由代理模型推出的函数，它的输入是可行集（feasible set）\(A\)上的任意值，输出值衡量了每个输入\(x\)有多值得被观测。通常会从以下两方面考虑：
- 有多大的可能性在x处取得最优值
- 评估x是否能减少贝叶斯统计模型的不确定性
  采集函数通常也是容易求最优值的函数（例如：有公式/能算梯度，等），下一个采样点就是可行集上的最大值点，即使采集函数的取最大值的点。

该框架的主要流程是：

代理模型（Surrogate Model）

代理模型可以有很多，高斯过程、随机森林……等等。其中，贝叶斯优化（Bayesian Optimization） 狭义上特指代理模型为高斯过程回归模型的SMBO。

随机过程

随机过程（Stochastic/Random Process）可以理解为一系列随机变量的集合。更具体地说，它是概率空间上的一族随机变量\(\{X(t),t\in T\}\)，其中是t参数，而T又被称作索引集（index set），它决定了构成随机过程的随机变量的个数，大部分情况下，随机变量都是无限个的。

当\(T=\{0,1,2,...\}\)时称之为随机序列或时间序列
参数t经常被解释为时间。
参数空间T是向量集合时，随机过程\(\{X(t),t\in T\}\)称为随机场。
随机过程有时也被称为随机函数（Random Function），因为它也可以被理解为函数值是随机变量的函数。

\(X(t)\)表示系统在时刻t所处的状态。的所有可能状态构成的集合为状态空间，记为S。

随机过程的直观理解

随机过程的直观含义是我们在 t=0 时刻（此时此刻）来考虑未来每个时间点会发生的情况。

例如，当我们把\(t\)解释为时间，而\(X(t)\)解释为粒子位置时，随机过程就可以被直观地理解为粒子的一种随机运动，它在每个时刻下的位置都是随机的，并且不同时刻\(t_1\)和\(t_2\)\(X(t_1)\)\(X(t_2)\)是相关的，它们的相关性由协方差\(cov[X(t_1),X(t_2)]\)定义。

辨析随机变量和随机函数

假设：

\(x\)是一个变量
\(X\)是一个随机变量，它服从一个概率分布\(P(X<t)=p(t)\)
\(f\)是一个函数，例如：\(f(x)=x+1\)
\(g\)是一个随机函数

那么有：

\(f(x)\)代表一个因变量，当\(x\)确定时，也确定。例如：\(f(x)=x+1\)，\(f(1)=2\)
\(f(X)\)代表一个随机变量，它的随机性是随机变量\(X\)带来的，所以它的概率分布也跟的概率分布有关。例如：\(f(x)=x+1\)，\(P(f(X)<t)=P(X+1<t)=P(X<t-1)\)
\(g(x)\)也代表一个随机变量，但它的随机性是随机过程\(g\)带来的。

高斯过程

高斯过程（Gaussian Process）是一类随机过程\(\{F(x),x\in A\}\)，它的任意n维分布\(\{F(x_1),...,F(x_n)\}\)（n也是任意的）都服从多元正态分布，即：对任意有限个\(x_1,...,x_n\in A\)，\(F(x_1),...,F(x_n)\)的任意线性组合\(a_1F(x_1)+...+a_nF(x_n)\)都是一个正态分布。

正如一个正态分布可以通过指定均值和方差来确定，一个高斯过程可以通过指定均值函数 \(m(x)\) 和协方差函数 \(K(x,x')\)唯一确定：

\[\begin{align*} m(x)&=E[F(x)] \\ K(x,x')&=E[(F(x)-m(x))(F(x')-m(x'))] \end{align*} \]

则高斯过程可以表示为：

\[F(x)\sim GP(m(x), K(x, x')) \]

均值函数定义了每个索引\(x\)对应的随机变量（同时也是正态分布变量）\(F(x)\)的均值；而协方差函数不仅定义了每个索引的方差\(K(x,x')\)，还定义了任意两个索引\(x_1\)、\(x_2\)对应的随机变量\(F(x_1 )\)\(F(x_1)\)\(F(x_2)\)之间的相关性\(K(x_1, x_2)\)。

在高斯过程模型里，协方差函数也被称作核函数（Kernel function）。

高斯过程回归

高斯过程回归（Gaussian Process Regression）就是使用高斯过程模型\(F(x)\)去拟合目标函数\(f(x)\)。

让我们先回顾一下，使用正态分布\(N(\mu,\sigma^2)\)去拟合一个随机变量X的步骤：

建模前，我们知道正态分布的均值和方差都是常数，分别假设为\(\mu\) 和\(\sigma\)。

对随机变量X进行t次采样，得到观测值\(x_1,...,x_t\)（简记为\(x_{1:t}\)）

根据观测值计算出最优的\(\mu\) 和\(\sigma\)值，由此确定最终正态分布的样子，从而完成对X的拟合

高斯过程回归也是类似的：

建模前，我们需要事先指定好均值函数 \(m(x)\) 和协方差函数 \(K(x,x')\)，它定义了高斯过程\(F(x)\)的先验分布。
- 有一些超参数会存在于均值函数和协方差函数中。前面提到的正态分布拟合就像是一种简化：\(m(x)=\mu\)，\(K(x,x')=\sigma^2\)，这里\(\mu\) 和\(\sigma\)可以被理解为是两个超参数。
选择t个采样索引\(x_1,...,x_t\)（简记为\(x_{1:t}\)），得到目标函数的观测值\(f(x_1),...,f(x_t)\)（简记为\(f(x_{1:t})\)），它们同样也是高斯过程里对应的随机变量\(F(x_1),...,F(x_t)\) \(F(x_1),...,F(x_t)\)\(F(x_{1:t})\)）的观测值
根据观测值调整均值函数和协方差函数里的超参数，由此确定最终高斯过程的样子，从而完成对函数\(f(x)\)的拟合。

常用均值函数

常数函数

\[m(x)=\mu \]
这里\(\mu\)是一个超参数。即使将均值统一设置为常数，因为有方差的作用，依然能够对数据进行有效建模。
parametric function

\[m(x)=\mu + \sum_{i=1}^P\beta_i\Psi_i(x) \]
这里\(\mu\)，\(\beta_i\)都是超参数。

常用协方差函数

在回归任务里，由于目标函数f是连续函数，因此对协方差函数最基础的要求是有该性质：

x与x’距离越近时，f(x)与f(x’)相关性越大。即：

\[K(x,x_1)>K(x,x_2), \text{ if } ||x-x_1||<||x-x_2|| \]

常用协方差函数有：

Power exponential / Gaussian

\[K(x,x')=\alpha_0 e^{-||x-x'||^2} \]
这里\(||x-x'||^2=\sum_{i=1}^d\alpha_i(x_i-x'_i)^2\)，不是标准的L2范式的写法，d是自变量x的维度数量，即参数的个数。\(\alpha_{0:d}\)是d+1个超参数。
Màtern

\[K(x,x')=\alpha_0 \frac{2^{1-\nu}}{\Gamma(\nu)}(\sqrt{2\nu}||x-x'||)^\nu K_\nu(\sqrt{2\nu}||x-x'||) \]
\(K_\nu\)是modified Bessel functions of the second kind。

如何计算超参数

采用最大后验估计（MAP，Maximum A Posterior estimate），选择在观测值已知的情况下最可能的超参数值。

\[\hat{\eta}=\text{argmax}_{\eta}{P(\eta|F(x_{1:t})=f(x_{1:t}))} \]

\(\eta\)代表超参数集合，\(P(\eta|F(x_{1:t})=f(x_{1:t}))\)代表在得到所有观测值的情况下，超参数的概率分布。

使用贝叶斯公式进行转换：

\[P(\eta|F(x_{1:t})=f(x_{1:t}))=\frac{P(F(x_{1:t})=f(x_{1:t})|\eta)P(\eta)}{P(F(x_{1:t})=f(x_{1:t}))} \]

这也是以高斯过程为代理模型的SMBO又被叫做贝叶斯优化的原因。

由于分母与超参数无关，因此只需要最大化分子：

\[\hat{\eta}=\text{argmax}_{\eta}{P(F(x_{1:t})=f(x_{1:t})|\eta)P(\eta)} \]

\(P(F(x_{1:t})=f(x_{1:t})|\eta)\)是给定参数的情况下，得到所有观测值的概率。

\(P(\eta)\)是超参数取值的先验概率，在某些问题，可以假设某些参数更可能被使用。但更常见的情况是假设其为等概率分布，那么MAP就是退化为最大似然估计（MLE, Maximum Likelihood Estimate），即选择使观测值最可能出现的超参数值：

\[\hat{\eta}=\text{argmax}_{\eta}{P(F(x_{1:t})=f(x_{1:t})|\eta)} \]

由于在高斯过程中，观测点\(x_{1:t}\)对应的随机变量构成的多维随机变量[\(F(x_1),...,F(x_t)\)]（简记作\(F(x_{1:t})\)）服从多元正态分布。其中，

分布的均值是一个t维向量：\(\mathbf{\mu_t}=[m(x_1),...,m(x_t)]\)。其中，\(m(x_i)\)代表第i个随机变量\(F(x_i)\)的期望/均值\(E[F(x_i)]\)。
分布的协方差矩阵是一个\(t\times t\)的矩阵，记作：

\[\mathbf{\Sigma_t} =\begin{bmatrix} K(x_1,x_1)&...&K(x_1,x_t) \\ ...&...&... \\ K(x_t,x_1)&...&K(x_t,x_t) \\ \end{bmatrix} \]

概率密度函数为

\[L(x_{1:t}|\eta)=\frac{1}{\sqrt{(2\pi)^t|\mathbf{\Sigma_t}|}}e^{-\frac{1}{2}(x_{1:t}-\mathbf{\mu_t})^T\mathbf{\Sigma_t}^{-1}(x_{1:t}-\mathbf{\mu_t})} \]

由于观测点\(x_{1:t}\)是已知的，因此，它是一个以超参数为变量的函数，我们可以选择使该概率（密度）最大化的超参数值。

转化为对数似然函数：

\[\ln{L(x_{1:t}|\eta)}=-\frac{1}{2}[t\ln{2\pi}+\ln{|\mathbf{\Sigma_t}|}+(x_{1:t}-\mathbf{\mu_t})^T\mathbf{\Sigma_t}^{-1}(x_{1:t}-\mathbf{\mu_t})] \]

要最大化它，即是最小化：

\[M(x_{1:t}|\eta)=\ln{|\mathbf{\Sigma_t}|}+(x_{1:t}-\mathbf{\mu_t})^T\mathbf{\Sigma_t}^{-1}(x_{1:t}-\mathbf{\mu_t}) \]

它是一个已知公式可以求梯度的函数，因此在编程中可以通过梯度下降等方法去求解其最优值。

预测\(f(x)\)

在SMBO框架下，代理模型预测的不仅仅是函数值\(f(x)\)，而是它的后验分布\(F(x)|F(x_{1:t})=f(x_{1:t})\)，即，在已知t个观测点的值的情况下，\(f(x)\)取不同值的概率。

由于在高斯过程里，也符合多元正态分布，因此就是多元正态分布的条件分布（或边缘分布）。这个分布服从正态分布，其均值和方差如下：

\[\begin{align*} F(x)|F(x_{1:t})=f(x_{1:t})&\sim N(\mu, \sigma^2) \\ \mu(x) &= K(x,x_{1:t})K(x_{1:t},x_{1:t})^{-1}(f(x_{1:t})-m(x_{1:t}))+m(x)\\ \sigma^2(x)&=K(x,x)-K(x,x_{1:t})K (x_{1:t}, x_{1:t})^{-1}K(x_{1:t},x)\end{align*} \]

其中，

\[m(x_{1:t})=[m(x_1),...,m(x_t)] \]

\[K(x_{1:t}, x_{1:t}) =\begin{bmatrix} K(x_1,x_{1:t}) \\ ... \\ K(x_t,x_{1:t}) \end{bmatrix} =\begin{bmatrix} K(x_1,x_1)&...&K(x_1,x_t) \\ ...&...&... \\ K(x_t,x_1)&...&K(x_t,x_t) \\ \end{bmatrix} \]

由于均值函数\(m(x)\)、协方差函数\(K(x,x')\)都是确定公式的函数，公式里的超参数也确定下来了，因此\(\mu (x)\)和\(\sigma^2(x)\)都是确定的函数，对任意给定的\(x\)，都可以计算出值来。

通常，我们使用这个后验分布的期望/均值来作为函数值\(f(x)\)的预测值：

\[f(x)\leftarrow E[F(x)|F(x_{1:t})=f(x_{1:t})]=\mu(x) \]

【命题1】对任意\(A=[a_1,...,a_t]^T\)，\(i=1,...,t\)，都有\(K(x_i,x_{1:t})K(x_{1:t},x_{1:t})^{-1}A=a_i\)。证明见附录。

因此，对于已经观测到的点\(x_i,i=1,...,t\)，我们就有：

\[\begin{align*} \mu(x_i)&=f(x_i)-m(x_i)+m(x_i)=f(x_i)\\ \sigma^2(x_i)&=K(x_i,x_i)-K(x_i,x_i)=0 \\ \end{align*} \]

因此，拟合出来的模型就像这样：

红色阴影部分是预测的95%置信区间。

采集函数（Acquisition Function）

由于代理模型输出了函数f的后验分布\(F(x)|F(x_{1:t})=f(x_{1:t})\)，我们可以利用这个后验分布去评估下一个采样点应该在哪个位置。由于在采集函数阶段我们讨论的都是后验分布，因此后文中将省略条件部分，提到\(F(x)\)时指的都是\(F(x)|F(x_{1:t})=f(x_{1:t})\)。

通常做法是设计一个采集函数\(A(x, F(x)|F(x_{1:t})=f(x_{1:t}))\)，它的输入相当于对每个采样点x进行打分，分数越高的点越值得被采样。

一般来说，采集函数需要满足下面的要求：

在已有的采样点处采集函数的值更小，因为这些点已经被探索过，再在这些点处计算函数值对解决问题没有什么用
在置信区间更宽（方差更大）的点处采集函数的值更大，因为这些点具有更大的不确定性，更值得探索
对最大（小）化问题，在函数均值更大（小）的点处采集函数的值更大，因为均值是对该点处函数值的估计值，这些点更可能在极值点附近。

有非常多种采集函数可供选择，这里简单介绍一些作为例子。

Expected Improvement (EI)

当我们已经采样过t个点之后，总会有一个最优点\(x_m\)，使得：

\[f^*_t=max_{i<t}f(x_i)=f(x_m) \]

假设我们还可以再观测一轮，得到\(F(x)=f(x)\)，最优点将在\(f(x)\)和\(f^*_t\)之间产生。不妨令

\[[F(x)-f^*_t]^+=\max(0, F(x)-f^*_t) \]

由于现在\([F(x)-f^*_t]^+\)是一个随机变量，因此我们可以计算它的期望：

\[\begin{align*} EI_t(x)&=E[[F(x)-f^*_t]^+] \\ &=\sigma(x)\phi(\frac{\mu(x)-f^*_t}{\sigma(x)})+(\mu(x)-f^*_t)\Phi(\frac{\mu(x)-f^*_t}{\sigma(x)}) \\ \end{align*} \]

其中，\(\mu(x)\)和\(\sigma(x)\)是正态分布\(F(x)\)的均值和标准差，即后验均值和标准差\(\sqrt{\sigma^2(x)}\)。

\(\varphi(x)\)为标准正态分布的概率密度函数：

\[\varphi(x)=-\frac{1}{\sqrt{2\pi}}g(x)=\frac{1}{\sqrt{2\pi}}e^{-\frac{x^2}{2}} \]

而\(\Phi(x)\)为标准正态分布的分布函数：

\[\Phi(x)=\frac{1}{\sqrt{2\pi}}\int_{-\infty}^xe^{-\frac{t^2}{2}}dt \]

\(EI_t(x)\)也是一个仅以x为自变量的函数，它的最大值点就是下一个采样点。

\[\hat{x}=\text{argmax}_{x}{EI_t(x)} \]

由于\(EI_t(x)\)有公式，计算不费劲，也可以求梯度，找到它的最大值/极大值有很多种现成的方案可以做到，相比于求原目标函数\(f(x)\)的最值要简单得多。

\(EI_t(x)\)推导过程

由于\(F(x)\)是正态分布，因此可以通过其概率密度计算出期望：

\[\begin{align*} EI_t(x)&=E[[F(x)-f^*_t]^+] \\ &=\frac{1}{\sigma\sqrt{2\pi}}\int_{-\infty}^{\infty}[z-f^*_t]^+e^{-\frac{(z-\mu)^2}{2\sigma^2}}dz \\ &=\frac{1}{\sigma\sqrt{2\pi}}\int_{f_t^*}^{\infty}(z-f^*_t)e^{-\frac{(z-\mu)^2}{2\sigma^2}}dz \end{align*} \]
令\(k=\frac{z-\mu}{\sigma}\)（也是把随机变量\(F(x)\)标准化），通过换元法可以得到：

\[\begin{align*} EI_t(x)&=\frac{1}{\sigma\sqrt{2\pi}}\int_{f_t^*}^{\infty}(z-f^*_t)e^{-\frac{(z-\mu)^2}{2\sigma^2}}dz \\ &=\frac{1}{\sqrt{2\pi}}\int_{\frac{f_t^*-\mu}{\sigma}}^{\infty}(\sigma k+\mu-f^*_t)e^{-\frac{k^2}{2}}dk \\ &=\frac{\sigma}{\sqrt{2\pi}}\int_{\frac{f_t^*-\mu}{\sigma}}^{\infty}ke^{-\frac{k^2}{2}}dk+(\mu-f^*_t)(1-\frac{1}{\sqrt{2\pi}}\int_{-\infty}^{\frac{f_t^*-\mu}{\sigma}}e^{-\frac{k^2}{2}}dk) \end{align*} \]
由于对\(g(x)=-e^{-\frac{x^2}{2}}\)，有\(g'(x)=xe^{-\frac{x^2}{2}}\)。因此

\[\int_{\frac{f_t^*-\mu}{\sigma}}^{\infty}ke^{-\frac{k^2}{2}}dk=g(\infty)-g(\frac{f_t^*-\mu}{\sigma})=-g(\frac{f_t^*-\mu}{\sigma}) \]
继续化简：

\[\begin{align*} EI_t(x)&=E[[F(x)-f^*_t]^+] \\ &=-\frac{\sigma}{\sqrt{2\pi}}g(\frac{f_t^*-\mu}{\sigma})+(\mu-f^*_t)(1-\Phi(\frac{f^*_t-\mu}{\sigma})) \\ &=\sigma\phi(\frac{f^*_t-\mu}{\sigma})+(\mu-f^*_t)[1-\Phi(\frac{f^*_t-\mu}{\sigma})] \\ \end{align*} \]
由于\(\varphi(x)\)是偶函数，且\(\Phi(x)+\Phi(-x)=1\)，则有

\[\begin{align*} EI_t(x)&=E[[F(x)-f^*_t]^+] \\ &=\sigma\phi(\frac{\mu-f^*_t}{\sigma})+(\mu-f^*_t)\Phi(\frac{\mu-f^*_t}{\sigma}) \end{align*} \]

置信区间边界（Confidence Bound）

对于标准正态分布\(X\sim N(0,1)\)，给定任意关于\(x=0\)对称的区间\([-z,z]\)，可以通过对概率密度进行积分计算出一个概率p，表示符合标准正态分布的随机变量的值落在该区间的概率为p。

z和p是一一对应的，可以令\(z=\delta(p)\)，由正态分布的性质可知，显然

\[p=\delta^{-1}(z)=P(X<z)=\int_{-\infty}^z\varphi(t)dt=\frac{1+p}{2},\ z\ge 0 \]

\(\delta^{-1}(z)\)就是\(\delta(p)\)的反函数。

假设对于给定置信水平p，随机变量\(X\)落在置信区间\([-z,z]\)的概率为p，\(z=\delta(p)\)。

由于\(F(x)\sim N(\mu,\sigma)\)，它可以由标准正态分布转化过来：

\[F(x)=\sigma X+\mu \]

那么对于\(F(x)\sim N(\mu,\sigma)\)，当给定置信水平p时，其置信区间就为\([-\mu-z\sigma,\mu+z\sigma]\)。

可以采用该置信区间的边界作为采集函数。例如，在本文中我们讨论最大化f(x)的优化问题，通常就会使用区间的右边界UCB（Upper Confidence Bound）：

\[UCB_t(x)=\mu(x)+z\sigma(x) \]

其中z是超参数，由想要的置信水平p决定：\(z=\delta(p)\)。

UCB越大，说明在该位置有更高的可能性找到更大的值。

同理，在最小化问题中，通常就会使用区间左边界LCB（Low Confidence Bound）：

\[LCB_t(x)=-\mu(x)-z\sigma(x) \]

Appendix

命题1证明

【命题1】对任意\(A=[a_1,...,a_t]^T\)，\(i=1,...,t\)，都有\(K(x_i,x_{1:t})K(x_{1:t},x_{1:t})^{-1}A=a_i\)。

令\(K(x_i,x_{1:t})K(x_{1:t},x_{1:t})^{-1}=[z_1,...,z_t]\)，那么有

\[\begin{align*} K(x_i,x_{1:t})K(x_{1:t},x_{1:t})^{-1}K(x_{1:t},x_{1:t})&=[z_1,...,z_t]K(x_{1:t},x_{1:t}) \\ K(x_i,x_{1:t})&=[\sum_{j=1}^t{z_jK(x_j,x_1)},...,\sum_{j=1}^t{z_jK(x_j,x_t)}] \\ \end{align*} \]

上述式子需对任意\(K(x,x’)\)都成立。列出方程组，得：

\[\left\{\begin{align*} \sum_{j=1}^t{z_jK(x_j,x_1)}&=K(x_i,x_1) \\ ... \\ \sum_{j=1}^t{z_jK(x_j,x_t)}&=K(x_i,x_t) \end{align*}\right. \]

由于\(K(x_{1:t},x_{1:t})\)可以求逆，因此它是满秩的，这意味着上面的方程组一定有唯一解。

经过变换得出：

\[\left\{\begin{align*} \sum_{j\neq i}{z_jK(x_j,x_1)}&=(1-z_i)K(x_i,x_1) \\ ... \\ \sum_{j\neq i}{z_jK(x_j,x_t)}&=(1-z_i)K(x_i,x_t) \end{align*}\right. \]

由于右侧的项都是一样的，而左侧\(K(x_i,x_j),j=1,...t\)各有不同，要想使他们相等，就必须等于0。因此解得：

\[\left\{\begin{align*} z_i&=1& \\ z_j&=0&, j\neq i,j=1,...t \end{align*}\right. \]

故有：对任意\(A=[a_1,...,a_t]^T\)，\(i=1,...,t\)，都有\(K(x_i,x_{1:t})K(x_{1:t},x_{1:t})^{-1}A=a_i\)。