随机森林贝叶斯优化 贝叶斯参数优化

关于参数搜索

参数搜索是一个开放的问题，假设我们拥有一个函数$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_{i} = (x_{i,1},x_{i,2}, ..., x_{i,n}) \in \mathbb{R}^n$

我们的目标是得到一个参数$x^* = \argmax_{x}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}=\begin{pmatrix}&k(x_1,x_1) &... &k(x_1,x_t)\\ &k(x_2,x_1) &... &k(x_2,x_n)\\&\vdots &\ddots &\vdots\\ &k(x_t,x_1) &\cdots &k(x_t,x_t)\end{pmatrix}$

我们用下标$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}_{1:t}(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} = \begin{pmatrix} \mu_{1:t}\\ \mu_{t+1} \end{pmatrix}, \Sigma_{1:t+1} = \begin{pmatrix} \Sigma_{1:t} &\Sigma \\ \Sigma^T&\Sigma_{t+1} \end{pmatrix}$

其中$\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^{-1}_{t+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}_{1:t+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)$

$\begin{pmatrix} A&B\\C&D \end{pmatrix}\begin{pmatrix} I&-A^{-1}B\\0&I \end{pmatrix} = \begin{pmatrix} A&0 \\ C&D-CA^{-1}B \end{pmatrix}$ 两边取行列式得到$\begin{vmatrix} A&B\\C&D \end{vmatrix} = |A| |D-CA^{-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\Big(-\frac{1}{2}\big((x-\mu_{1:t+1})\Sigma^{-1}_{1:t+1}(x-\mu_{1:t+1})^T - (x-\mu_{1:t})\Sigma^{-1}_{1:t}(x-\mu_{1:t})^T\big)\Big)$

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

从而
$x-\mu_{1:t+1} = \begin{pmatrix} \hat{x}_{1:t} & \hat{x}_{t+1} \end{pmatrix}$

这里需要用到分块矩阵的逆
$\begin{pmatrix} A&B\\C&D \end{pmatrix}^{-1} = \begin{pmatrix} H &-HBD^{-1}\\-D^{-1}CH&D^{-1}(I+CHBD^{-1}) \end{pmatrix}$ 其中$H = (A-BD^{-1}C)^{-1}$

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

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

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

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

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

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

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

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

我们便得到了两个重要的参数估计
$\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}_{1:t}\Sigma^T$

其中方差$\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^* = \argmax_{x} P(g(x) \geq g(x^+))$
其中$x^+ = \argmax_{x\in x_{1:t}}g(x)$表示着在已有观测值里面，对应评价最高的参数。

$P(X\leq x) = \int_{-\infty}^{x}\phi(x)dx$, 其中$\phi(x)$是概率密度函数

我们之前对$g(x)$进行高斯分布的建模，于是，上述的概率可以写成分布函数的形式
$P(g(x)\geq g(x^+)) = 1 - P(g(x)\leq 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})$

$\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^+)\}$
代表着新尝试的值对于已有的最佳值的提升，注意到该函数的值服从正态分布，其概率密度函数为
$\frac{1}{\sigma\sqrt{2\pi}}\exp(-\frac{(I(x)-(\mu-g(x^+)))^2}{2\sigma^2})$

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

$\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))$

我们分为两个部分进行解析，首先是第一个部分
$\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})du^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\varphi(\frac{\mu-g(x^+)}{\sigma})$

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

第二个部分为
$\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^+))\Big(1 - \int_{-\infty}^{0}\frac{1}{\sigma\sqrt{2\pi}}\exp(-\frac{(I(x)-(\mu-g(x^+)))^2}{2\sigma^2})d(I(x))\Big)$
令$I(x) = I(x) + g(x^+)$, 则有
$(\mu-g(x^+))\Big(1 - \int_{-\infty}^{0}\frac{1}{\sigma\sqrt{2\pi}}\exp(-\frac{(I(x)-(\mu-g(x^+)))^2}{2\sigma^2})d(I(x))\Big)\\=(\mu-g(x^+))\Big(1 - \int_{-\infty}^{g(x^+)}\frac{1}{\sigma\sqrt{2\pi}}\exp(-\frac{(I(x)-\mu)^2}{2\sigma^2})d(I(x))\Big)\\=(\mu-g(x^+))\Big(1 - \Phi(\frac{\mu-g(x^+)}{\sigma})\Big)$

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

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