贝叶斯优化详解

关于参数搜索

参数搜索是一个开放的问题，假设我们拥有一个函数$f(x;\theta )$，该函数的性质完全由参数$\theta$决定，那么，我们就可以通过对参数进行搜索，得到一个满足我们要求的函数。

这里我所说的满足要求，可以是拟合的误差小于某个阈值，亦可以是其他

当选取一个参数${\theta }_c$之后（这里下标c表示候选，candidate），我们会对这个函数进行一个评价，这里的评价是指，评价这个函数是不是符合我们的要求。

一般来说，这个评价，或者说成评估的过程，往往是代价比较大的，即，计算耗费的时间比较多，不能立即得到结果。所以，在进行参数搜索的时候，我们更希望用一个高效的方法来进行尝试。

在这里，你可能会想到随机梯度下降（Stochastic Gradient Descent），直接求出参数的偏导数即可。但是，我这里所说的参数，有时候不是十分直接地参与函数的计算，计算偏导数可能十分困难，比如神经网络的学习率，每一层神经元的个数等。

对于参数搜索，比较朴素的有，随机搜索以及网格搜索。随机搜索主要是基于一个概率，随机地进行参数的变动。而网格搜索是给定一系列需要尝试的离散值的组合，挨个尝试。

本文介绍的贝叶斯优化，属于一种参数搜索方案，利用先验的知识，帮助我们去对参数进行搜索。

本文的细节比较多，会将其中的推导的详细步骤都展示出来

贝叶斯优化

参数搜索问题的核心在于，如何基于已有的观察，决定下一个需要尝试的参数？

假设我们现在有一个评价函数$g(x)$，通过前$t$次的尝试，我们获得了一系列点 { ( x 1 , g ( x 1 ) ) , ( x 2 , g ( x 2 ) ) , . . . , ( x t , g ( x t ) ) } \{(x_1,g(x_1)),(x_2,g(x_2)),...,(x_t,g(x_t))\} {(x1​,g(x1​)),(x2​,g(x2​)),...,(xt​,g(xt​))}，其中$x_i=(x_{i,1},x_{i,2},...,x_{i,n})\in {\mathbb{R}}^n$

我们的目标是得到一个参数$x^{\ast }=\underset{x}{\mathrm{argmax}}g(x)$，这里不局限于最大，也可以是最小值，取决于你的评价的设计。

在贝叶斯优化之中，有一个重要的假设：
$$x\sim \mathcal{N}({\mu }_x,{\sigma }_x),g(x)\sim \mathcal{N}({\mu }_g,{\sigma }_g)$$

基于这个假设，在已有的$t$个观察点上，我们可以构造出一个针对$g(x)$的高斯分布，参数为
$${\mu }_{1:t}=\frac{1}{t}\sum _{i=1}^{t}g(x_i),{\Sigma }_{1:t}=\left(\begin{array}{cccc}& k(x_1,x_1)& ...& k(x_1,x_t)\\ & k(x_2,x_1)& ...& k(x_2,x_n)\\ & ⋮& \ddots & ⋮\\ & k(x_t,x_1)& \cdots & k(x_t,x_t)\end{array}\right)$$

我们用下标$1:t$表示在从$1$到$t$的观察点上的

其中$k(x_i,x_j)=||x_i-x_j|{|}_p$表示$x_i$和$x_j$通过核函数（kernel function）的内积，用来描述二者相似程度。这样的一个核函数矩阵可以作为协方差矩阵。

范数p可以取1，2等等

于是，关于$g(x_{1:t})$的概率密度函数为$$f_{1:t}(x)=\frac{1}{(2\pi {)}^{t/2}(\det ({\Sigma }_{1:t}){)}^{1/2}}\exp (-\frac{1}{2}(x-{\mu }_{1:t}){\Sigma }_{1:t}^{-1}(x-{\mu }_{1:t}{)}^T)$$

不妨假设，下一个尝试的点为$x_{t+1}$, $x_{1:t+1}$亦是服从正态分布$\mathcal{N}({\mu }_{t+1},{\Sigma }_{t+1})$，于是，$x_{1:t}$和$x_{t+1}$的联合分布的参数应该为：
$${\mu }_{1:t+1}=\left(\begin{array}{c}{\mu }_{1:t}\\ {\mu }_{t+1}\end{array}\right),{\Sigma }_{1:t+1}=\left(\begin{array}{cc}{\Sigma }_{1:t}& \Sigma \\ {\Sigma }^T& {\Sigma }_{t+1}\end{array}\right)$$

其中$\Sigma =[k(x_1,x_{t+1}),k(x_2,x_{t+1}),...,k(x_t,x_{t+1}){]}^T$以及${\Sigma }_{t+1}=k(x_{t+1},x_{t+1})$

我们的目的其实是去在$x_{1:t}$的条件下估计$g(x_{t+1})$

$g(x_{t+1})$的概率密度函数为
$$f_{t+1}(x)=\frac{1}{(2\pi {)}^{1/2}(det({\Sigma }_{t+1}){)}^{1/2}}\exp (-\frac{1}{2}(x-{\mu }_{t+1}){\Sigma }_{t+1}^{-1}(x-{\mu }_{t+1}{)}^T)$$

$g(x_{1:t})$和$g(x_{t+1})$的联合概率密度函数为
$$f_{1:t,t+1}(x)=\frac{1}{(2\pi {)}^{(t+1)/2}(det({\Sigma }_{1:t+1}){)}^{1/2}}\exp (-\frac{1}{2}(x-{\mu }_{1:t+1}){\Sigma }_{1:t+1}^{-1}(x-{\mu }_{1:t+1}{)}^T)$$

进一步地，我们可以得到一个条件概率密度函数
$$f_{t+1|1:t}(x)=\frac{f_{1:t,t+1}(x)}{f_{1:t}(x)}$$

我们将结果分成两个部分：系数部分和指数部分

系数部分有：
$$\frac{(2\pi {)}^{t/2}(\det ({\Sigma }_{1:t}){)}^{1/2}}{(2\pi {)}^{(t+1)/2}(\det ({\Sigma }_{1:t+1}){)}^{1/2}}=\frac{1}{(2\pi {)}^{1/2}(\det ({\Sigma }_{1:t+1})/\det ({\Sigma }_{1:t}){)}^{1/2}}$$

注意到分块矩阵的一个性质：$\det ({\Sigma }_{1:t+1})=\det ({\Sigma }_{1:t})\det ({\Sigma }_{t+1}-{\Sigma }^T{\Sigma }_{1:t}^{-1}\Sigma )$

$\left(\begin{array}{cc}A& B\\ C& D\end{array}\right)\left(\begin{array}{cc}I& -A^{-1}B\\ 0& I\end{array}\right)=\left(\begin{array}{cc}A& 0\\ C& D-C{A}^{-1}B\end{array}\right)$ 两边取行列式得到$\left|\begin{array}{cc}A& B\\ C& D\end{array}\right|=|A||D-C{A}^{-1}B|$

从而系数部分结果为：
$$\frac{1}{(2\pi {)}^{1/2}(\det ({\Sigma }_{t+1}-{\Sigma }^T{\Sigma }_{1:t}^{-1}\Sigma ){)}^{1/2}}$$

从这个结果可以看出来方差为${\Sigma }_{t+1|1:t}={\Sigma }_{t+1}-{\Sigma }^T{\Sigma }_{1:t}^{-1}\Sigma$

对于指数部分
$$E=\exp (-\frac{1}{2}((x-{\mu }_{1:t+1}){\Sigma }_{1:t+1}^{-1}(x-{\mu }_{1:t+1}{)}^T-(x-{\mu }_{1:t}){\Sigma }_{1:t}^{-1}(x-{\mu }_{1:t}{)}^T\left)\right)$$

为了记录方便，我们这里使用一些记号：$x-{\mu }_{1:t}={\hat{x}}_{1:t},x-{\mu }_{t+1}={\hat{x}}_{t+1}$

从而
$$x-{\mu }_{1:t+1}=\left(\begin{array}{cc}{\hat{x}}_{1:t}& {\hat{x}}_{t+1}\end{array}\right)$$

这里需要用到分块矩阵的逆
${\left(\begin{array}{cc}A& B\\ C& D\end{array}\right)}^{-1}=\left(\begin{array}{cc}H& -HB{D}^{-1}\\ -D^{-1}CH& D^{-1}(I+CHB{D}^{-1})\end{array}\right)$ 其中$H=(A-B{D}^{-1}C{)}^{-1}$

分块矩阵的逆怎么来的？
若存在一个逆矩阵，不妨记作$\left(\begin{array}{cc}E& F\\ G& H\end{array}\right)$
由于$\left(\begin{array}{cc}A& B\\ C& D\end{array}\right)\left(\begin{array}{cc}E& F\\ G& H\end{array}\right)=\left(\begin{array}{cc}I& 0\\ 0& I\end{array}\right)$
便有方程组$\begin{array}{c}AE+BG=I\\ AF+BH=0\\ CE+DG=0\\ CF+DH=I\end{array}$
$(AE+BG)-B{D}^{-1}(CE+DG)=I=AE-B{D}^{-1}CE=(A-B{D}^{-1}C)E$
进而$E=A-B{D}^{-1}C$，类似地可以得到其他值

基于分块矩阵的逆，对于${\Sigma }_{1:t+1}$，我们可以得到
$${\Sigma }_{1:t+1}^{-1}=\left(\begin{array}{cc}H& -H\Sigma {\Sigma }_{t+1}^{-1}\\ -{\Sigma }_{t+1}^{-1}{\Sigma }^{T}H& {\Sigma }_{t+1}^{-1}(I+\Sigma H{\Sigma }^T{\Sigma }_{t+1}^{-1})\end{array}\right)$$
其中$H=({\Sigma }_{1:t}-{\Sigma }^T{\Sigma }_{t+1}^{-1}\Sigma {)}^{-1}$

对称地，可以得到
$${\Sigma }_{1:t+1}^{-1}=\left(\begin{array}{cc}{\Sigma }_{1:t}^{-1}(I+{\Sigma }^{T}H\Sigma {\Sigma }_{1:t}^{-1})& -{\Sigma }_{1:t}^{-1}{\Sigma }^{T}H\\ -H\Sigma {\Sigma }_{1:t}^{-1}& H\end{array}\right)$$
其中$H=({\Sigma }_{t+1}-\Sigma {\Sigma }_{1:t}^{-1}{\Sigma }^T{)}^{-1}$

将指数部分拆解合并得到
$$\begin{gathered} E=\exp (-\frac{1}{2}({\hat{x}}_{1:t}({\Sigma }_{1:t}^{-1}({\Sigma }^{T}H\Sigma {\Sigma }_{1:t}^{-1})){\hat{x}}_{1:t}^T+{\hat{x}}_{t+1}(-H\Sigma {\Sigma }_{1:t}^{-1}){\hat{x}}_{1:t}^T \\ +{\hat{x}}_{1:t}(-{\Sigma }_{1:t}^{-1}{\Sigma }^{T}H){\hat{x}}_{t+1}^T+{\hat{x}}_{t+1}H{\hat{x}}_{t+1}^T)) \end{gathered}$$

利用之前我们得到的方差${\Sigma }_{t+1|1:t}={\Sigma }_{t+1}-{\Sigma }^T{\Sigma }_{1:t}^{-1}\Sigma$

猜测均值为${\mu }^{\ast }$,有$({\hat{x}}_{t+1}-{\mu }^{\ast }){\Sigma }_{t+1|1:t}^{-1}({\hat{x}}_{t+1}-{\mu }^{\ast }{)}^T$

展开之后对齐系数项可以得到
$${\mu }^{\ast }={\hat{x}}_{1:t}{\Sigma }_{1:t}^{-1}{\Sigma }^T$$

可以得到最终的均值估计为
$${\mu }_{t+1|1:t}={\mu }_{t+1}+{\hat{x}}_{1:t}{\Sigma }_{1:t}^{-1}{\Sigma }^T$$

我们便得到了两个重要的参数估计
$$\begin{gathered} {\Sigma }_{t+1|1:t}={\Sigma }_{t+1}-{\Sigma }^T{\Sigma }_{1:t}^{-1}\Sigma \\ {\mu }_{t+1|1:t}={\mu }_{t+1}+{\hat{x}}_{1:t}{\Sigma }_{1:t}^{-1}{\Sigma }^T \end{gathered}$$

其中方差${\Sigma }_{t+1|1:t}$可以通过$x_{t+1}$和$x_{1:t}$通过核函数计算得到
均值${\mu }_{t+1|1:t}$的估计中的${\mu }_{t+1}$可以使用${\mu }_{1:t}$近似计算

然后我们可以计算出每个未知点的概率为0.95的置信区间

通常为$[\mu -1.96\frac{\sigma }{\sqrt{n}},\mu +1.96\frac{\sigma }{\sqrt{n}}]$

到了这里，我们基本上获得了每一个未知点的贝叶斯后验估计，在这之上需要评估到底该选择哪些未知点进行尝试

这里介绍两种，分别是Probability of Improvement以及Expected Improvement
简记为PI方法以及EI方法

对于PI方法来说，我们的目标为
$$x^{\ast }=\underset{x}{\mathrm{argmax}}P(g(x)\ge g(x^{+}))$$
其中$x^{+}=\underset{x\in x_{1:t}}{\mathrm{argmax}}g(x)$表示着在已有观测值里面，对应评价最高的参数。

$P(X\le x)={\int }_{-\infty }^x\varphi (x)dx$, 其中$\varphi (x)$是概率密度函数

我们之前对$g(x)$进行高斯分布的建模，于是，上述的概率可以写成分布函数的形式
$$\begin{gathered} P(g(x)\ge g(x^{+}))=1-P(g(x)\le g(x^{+})) \\ =1-{\int }_{-\infty }^{g(x^{+})}\frac{1}{\sigma \sqrt{(2\pi )}}\exp (-\frac{1}{2{\sigma }^2}(x-\mu {)}^2)dx \\ =1-\Phi (\frac{\mu -g(x^{+})}{\sigma }) \end{gathered}$$

$\Phi (x)=\frac{1}{\sqrt{2\pi }}{\int }_{-\infty }^x\exp (-\frac{u^2}{2})du$ 是标准的正态累积分布函数

在上述式子中，变换到标准的正态累积分布函数的过程并不是很直接，这里给出中间的变换步骤
$$\Phi (\frac{\mu -g(x^{+})}{\sigma })=\frac{1}{\sqrt{2\pi }}{\int }_{-\infty }^{\frac{\mu -g(x^{+})}{\sigma }}\exp (-\frac{x^2}{2})dx$$
令$x=\frac{x}{\sigma }$, 便有
$$\Phi (\frac{\mu -g(x^{+})}{\sigma })=\frac{1}{\sqrt{2\pi }}{\int }_{-\infty }^{\mu -g(x^{+})}\exp (-\frac{x^2}{2{\sigma }^2})d(\frac{x}{\sigma })=\frac{1}{\sigma \sqrt{2\pi }}{\int }_{-\infty }^{\mu -g(x^{+})}\exp (-\frac{x^2}{2{\sigma }^2})dx$$
再令$x=x-\mu$，便有
$$\Phi (\frac{\mu -g(x^{+})}{\sigma })=\frac{1}{\sigma \sqrt{2\pi }}{\int }_{-\infty }^{g(x^{+})}\exp (-\frac{(x-\mu {)}^2}{2{\sigma }^2})d(x-\mu )=\frac{1}{\sigma \sqrt{2\pi }}{\int }_{-\infty }^{g(x^{+})}\exp (-\frac{(x-\mu {)}^2}{2{\sigma }^2})dx$$

我们会去计算所有未知的待搜索点的这个概率，进而选择一个最高概率的值进行探索。

PI方法被认为搜索的范围较为小，只去搜索局部的区域，因此有了搜索区域更广的EI方法

首先定义一个获益函数
I ( x ) = max ⁡ { 0 , g ( x ) − g ( x + ) } I(x) = \max\{0, g(x) - g(x^+)\} I(x)=max{0,g(x)−g(x+)}
代表着新尝试的值对于已有的最佳值的提升，注意到该函数的值服从正态分布，其概率密度函数为
$$\frac{1}{\sigma \sqrt{2\pi }}\exp (-\frac{(I(x)-(\mu -g(x^{+})){)}^2}{2{\sigma }^2})$$

我们去计算这个函数的期望

$$\begin{gathered} \mathbb{E}(I(x))={\int }_0^{+\infty }I(x)\frac{1}{\sigma \sqrt{2\pi }}\exp (-\frac{(I(x)-(\mu -g(x^{+})){)}^2}{2{\sigma }^2})d(I(x)) \\ ={\int }_0^{+\infty }(I(x)-\mu +g(x^{+})+\mu -g(x^{+}))\frac{1}{\sigma \sqrt{2\pi }}\exp (-\frac{(I(x)-(\mu -g(x^{+})){)}^2}{2{\sigma }^2})d(I(x)) \\ ={\int }_0^{+\infty }(I(x)-\mu +g(x^{+}))\frac{1}{\sigma \sqrt{2\pi }}\exp (-\frac{(I(x)-(\mu -g(x^{+})){)}^2}{2{\sigma }^2})d(I(x)) \\ +{\int }_0^{+\infty }(\mu -g(x^{+}))\frac{1}{\sigma \sqrt{2\pi }}\exp (-\frac{(I(x)-(\mu -g(x^{+})){)}^2}{2{\sigma }^2})d(I(x)) \end{gathered}$$

我们分为两个部分进行解析，首先是第一个部分
$$\begin{gathered} {\int }_0^{+\infty }(I(x)-\mu +g(x^{+}))\frac{1}{\sigma \sqrt{2\pi }}\exp (-\frac{(I(x)-(\mu -g(x^{+})){)}^2}{2{\sigma }^2})d(I(x)) \\ ={\int }_{-\mu +g(x^{+})}^{+\infty }u\frac{1}{\sigma \sqrt{2\pi }}\exp (-\frac{u^2}{2{\sigma }^2})du \\ =\frac{1}{\sigma \sqrt{2\pi }}{\int }_{-\mu +g(x^{+})}^{+\infty }u\exp (-\frac{u^2}{2{\sigma }^2})du \\ =\frac{1}{\sigma \sqrt{2\pi }}{\int }_{(-\mu +g(x^{+}){)}^2}^{+\infty }\frac{1}{2}\exp (-\frac{u^2}{2{\sigma }^2})d{u}^2 \\ =\frac{1}{2\sigma \sqrt{2\pi }}{\int }_{(-\mu +g(x^{+}){)}^2}^{+\infty }\exp (-\frac{t}{2{\sigma }^2})dt \\ =\frac{1}{2\sigma \sqrt{2\pi }}[-2{\sigma }^2\exp (-\frac{t}{2{\sigma }^2}){]}_{(-\mu +g(x^{+}){)}^2}^{+\infty } \\ =\frac{1}{2\sigma \sqrt{2\pi }}(0-(-2{\sigma }^2\exp (-\frac{(-\mu +g(x^{+}){)}^2}{2{\sigma }^2}))) \\ =\frac{\sigma }{\sqrt{2\pi }}\exp (\frac{(-\mu +g(x^{+}){)}^2}{2{\sigma }^2}) \\ =\sigma \phi (\frac{\mu -g(x^{+})}{\sigma }) \end{gathered}$$

其中$\phi (x)=\frac{1}{\sqrt{2\pi }}\exp (-\frac{x^2}{2})$是标准正态分布的概率密度函数

第二个部分为
$$\begin{gathered} {\int }_0^{+\infty }(\mu -g(x^{+}))\frac{1}{\sigma \sqrt{2\pi }}\exp (-\frac{(I(x)-(\mu -g(x^{+})){)}^2}{2{\sigma }^2})d(I(x)) \\ =(\mu -g(x^{+})){\int }_0^{+\infty }\frac{1}{\sigma \sqrt{2\pi }}\exp (-\frac{(I(x)-(\mu -g(x^{+})){)}^2}{2{\sigma }^2})d(I(x)) \\ =(\mu -g(x^{+}))(1-{\int }_{-\infty }^0\frac{1}{\sigma \sqrt{2\pi }}\exp (-\frac{(I(x)-(\mu -g(x^{+})){)}^2}{2{\sigma }^2})d(I(x))) \end{gathered}$$
令$I(x)=I(x)+g(x^{+})$, 则有
$$\begin{gathered} (\mu -g(x^{+}))(1-{\int }_{-\infty }^0\frac{1}{\sigma \sqrt{2\pi }}\exp (-\frac{(I(x)-(\mu -g(x^{+})){)}^2}{2{\sigma }^2})d(I(x))) \\ =(\mu -g(x^{+}))(1-{\int }_{-\infty }^{g(x^{+})}\frac{1}{\sigma \sqrt{2\pi }}\exp (-\frac{(I(x)-\mu {)}^2}{2{\sigma }^2})d(I(x))) \\ =(\mu -g(x^{+}))(1-\Phi (\frac{\mu -g(x^{+})}{\sigma })) \end{gathered}$$

所以期望的结果为
$$\mathbb{E}(I(x))=(\mu -g(x^{+}))(1-\Phi (\frac{\mu -g(x^{+})}{\sigma }))+\sigma \phi (\frac{\mu -g(x^{+})}{\sigma })$$

同样的，我们会对所有未知点计算出期望，然后选择期望最高的未知点进行探索。