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

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

问题介绍

[max_{xin A}f(x)
]

输入：(xin R^d)，一般(dle20)，即参数在20个以内。
可行集（feasible set）(A)： 一般是超多边形（(xin R^d: a_ile x_i le b_i)）或d维的simplex（(xin 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),tin T})，其中是t参数，而T又被称作索引集（index set），它决定了构成随机过程的随机变量的个数，大部分情况下，随机变量都是无限个的。

当(T={0,1,2,…})时称之为随机序列或时间序列
参数t经常被解释为时间。
参数空间T是向量集合时，随机过程({X(t),tin 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
(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)
(g(x))也代表一个随机变量，但它的随机性是随机过程(g)带来的。

高斯过程

高斯过程（Gaussian Process）是一类随机过程({F(x),xin A})，它的任意n维分布({F(x_1),…,F(x_n)})（n也是任意的）都服从多元正态分布，即：对任意有限个(x_1,…,x_nin 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}^Pbeta_iPsi_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||

常用协方差函数有：

Power exponential / Gaussian

[K(x,x’)=alpha_0 e^{-||x-x’||^2}
]

这里(||x-x’||^2=sum_{i=1}^dalpha_i(x_i-x’_i)^2)，不是标准的L2范式的写法，d是自变量x的维度数量，即参数的个数。(alpha_{0:d})是d+1个超参数。
Mtern

[K(x,x’)=alpha_0 frac{2^{1-nu}}{Gamma(nu)}(sqrt{2nu}||x-x’||)^nu K_nu(sqrt{2nu}||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)])。
分布的协方差矩阵是一个(ttimes 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{(2pi)^t|mathbf{Sigma_t}|}}e^{-frac{1}{2}(x_{1:t}-mathbf{mu_t})^Tmathbf{Sigma_t}^{-1}(x_{1:t}-mathbf{mu_t})}
]

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

转化为对数似然函数：

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

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

[M(x_{1:t}|eta)=ln{|mathbf{Sigma_t}|}+(x_{1:t}-mathbf{mu_t})^Tmathbf{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

假设我们还可以再观测一轮，得到(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{2pi}}g(x)=frac{1}{sqrt{2pi}}e^{-frac{x^2}{2}}
]

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

[Phi(x)=frac{1}{sqrt{2pi}}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}{sigmasqrt{2pi}}int_{-infty}^{infty}[z-f^*_t]^+e^{-frac{(z-mu)^2}{2sigma^2}}dz \
&=frac{1}{sigmasqrt{2pi}}int_{f_t^*}^{infty}(z-f^*_t)e^{-frac{(z-mu)^2}{2sigma^2}}dz
end{align*}
]

令(k=frac{z-mu}{sigma})（也是把随机变量(F(x))标准化），通过换元法可以得到：

[begin{align*}
EI_t(x)&=frac{1}{sigmasqrt{2pi}}int_{f_t^*}^{infty}(z-f^*_t)e^{-frac{(z-mu)^2}{2sigma^2}}dz \
&=frac{1}{sqrt{2pi}}int_{frac{f_t^*-mu}{sigma}}^{infty}(sigma k+mu-f^*_t)e^{-frac{k^2}{2}}dk \
&=frac{sigma}{sqrt{2pi}}int_{frac{f_t^*-mu}{sigma}}^{infty}ke^{-frac{k^2}{2}}dk+(mu-f^*_t)(1-frac{1}{sqrt{2pi}}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{2pi}}g(frac{f_t^*-mu}{sigma})+(mu-f^*_t)(1-Phi(frac{f^*_t-mu}{sigma})) \
&=sigmaphi(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]^+] \
&=sigmaphi(frac{mu-f^*_t}{sigma})+(mu-f^*_t)Phi(frac{mu-f^*_t}{sigma})
end{align*}
]

置信区间边界（Confidence Bound）

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

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

[p=delta^{-1}(z)=P(X

(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-zsigma,mu+zsigma])。

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

[UCB_t(x)=mu(x)+zsigma(x)
]

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

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

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

[LCB_t(x)=-mu(x)-zsigma(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_{jneq i}{z_jK(x_j,x_1)}&=(1-z_i)K(x_i,x_1) \
… \
sum_{jneq 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&, jneq 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)。

posted @
2023-10-23 15:38
milliele
阅读(0)
评论(0)
编辑
收藏
举报

var cb_entryId = 17782631, cb_entryCreatedDate = ‘2023-10-23 15:38’, cb_postType = 1, cb_postTitle = ‘一文详解贝叶斯优化（Bayesian Optimization）原理’;
var allowComments = true, cb_blogId = 521207, cb_blogApp = ‘milliele’, cb_blogUserGuid = ’42e647c1-6f64-4fbb-10ba-08d6ef9f5750′;
mermaidRender.render()
markdown_highlight()
zoomManager.apply(“#cnblogs_post_body img:not(.code_img_closed):not(.code_img_opened)”);

会员救园

刷新页面返回顶部

var commentManager = new blogCommentManager();
commentManager.renderComments(0);
fixPostBody();
window.footnoteTipManager.generateFootnoteTips();

window.tocManager.displayDisableTocTips = false;
window.tocManager.generateToc();

setTimeout(function() { countViews(cb_blogId, cb_entryId); }, 50);

deliverT2();
deliverC1C2();
loadNewsAndKb();

LoadPostCategoriesTags(cb_blogId, cb_entryId);

LoadPostInfoBlock(cb_blogId, cb_entryId, cb_blogApp, cb_blogUserGuid);
GetPrevNextPost(cb_entryId, cb_blogId, cb_entryCreatedDate, cb_postType);
loadOptUnderPost();
GetHistoryToday(cb_blogId, cb_blogApp, cb_entryCreatedDate);

服务器托管，北京服务器托管，服务器租用 http://www.fwqtg.net
机房租用，北京机房租用，IDC机房托管， http://www.fwqtg.net