贝叶斯优化
贝叶斯优化（Bayesian Optimization, BO）是一类面向中低维度昂贵黑盒函数优化问题的序贯决策方法。它所解决的问题看似简单：
\[x^\star = \arg\min_{x\in\mathcal X} f(x)\]真正困难之处在于，我们通常既不知道 $f(x)$ 的解析表达式，也无法获得其梯度，甚至每计算一次 $f(x)$ 都意味着一次昂贵的真实实验。
例如，$x$ 可以是一组温度、压力、浓度和时间参数，而 $f(x)$ 是实验结束后才能测得的材料性能；也可以是一组机器学习超参数，而 $f(x)$ 是完成一次训练之后才能获得的验证集误差。
真实函数评估次数非常昂贵，因此必须尽可能让每一次实验都提供最大的信息价值。
贝叶斯优化正是围绕这一目标构建的。
它并不直接优化未知的真实函数 $f(x)$，而是建立一个计算成本很低的概率代理模型（surrogate model），利用已经完成的实验推断未知区域的函数值及其不确定性，再通过采集函数（acquisition function）决定下一次最值得进行的实验。
可以将整个过程概括为：
\[\text{已有实验数据} \rightarrow \text{概率代理模型} \rightarrow \text{采集函数} \rightarrow \text{下一实验点} \rightarrow \text{真实实验} \rightarrow \text{更新代理模型}\]背景
黑盒优化问题
考虑一个 $d$ 维输入空间：
\[x=(x_1,x_2,\ldots,x_d)^\mathsf T\in\mathcal X\]真实目标函数为：
\[y=f(x)\]如果实验中存在随机误差，则更合理的模型是：
\[y=f(x)+\varepsilon\]其中噪声通常近似表示为：
\[\varepsilon\sim\mathcal N(0,\sigma_n^2)\]在普通数值优化中，我们可能通过梯度下降不断计算 $\nabla f(x)$。
但对于黑盒实验，$f$ 本身未知，更不用说梯度。我们真正拥有的是有限个已经测量过的数据：
\[\mathcal D_n=\{(x_i,y_i)\}_{i=1}^{n}\]因此贝叶斯优化面对的核心问题不是：
已知 $f(x)$，怎样最快求出最优解？
而是：
只知道有限个函数观测值时，下一次应该在哪里测量，才能在尽可能少的后续实验次数的情况下接近全局最优解？
网格搜索
假设每个变量仅离散为 $m$ 个可能值，而问题有 $d$ 个维度，那么完整网格需要 $N=m^d$ 次实验。
例如 $d=8$、每个变量仅测试 $10$ 个水平，就已经需要 $10^8$ 个组合。因此，高成本实验几乎不可能通过穷举解决。
随机搜索虽然避免了规则网格的一些问题，但随机搜索并不会利用已经获得的实验结果。第 $100$ 次实验的选择方式和第 $1$ 次实验没有采样上的区别。
贝叶斯优化则不同：
\[\mathcal D_n \rightarrow p(f\mid\mathcal D_n) \rightarrow x_{n+1}\]每完成一次实验，都会改变模型对未知函数的认识，从而改变下一次实验的位置。
概率视角
贝叶斯推断
贝叶斯推断的一般形式可以写成：
\[p(\theta\mid\mathcal D) = \frac{p(\mathcal D\mid\theta)p(\theta)} {p(\mathcal D)}\]其中 $p(\theta)$ 是先验，$p(\mathcal D\mid\theta)$ 是似然，而 $p(\theta\mid\mathcal D)$ 是观察数据之后得到的后验。
贝叶斯优化中的特殊之处在于，我们感兴趣的未知对象并不只是一个有限维参数 $\theta$，而是未知函数 $f$。
因此可以抽象地写成：
\[p(f\mid\mathcal D) \propto p(\mathcal D\mid f)p(f)\]即：
- 首先对 函数族 function family 建立先验；
- 完成若干真实实验；
- 使用实验结果更新函数分布；
- 根据后验分布决定下一次实验点选择。
Gaussian Process 是最常被选用的函数族。
Gaussian Process
从随机变量扩展到随机函数
普通概率模型描述的通常是有限个随机变量，而高斯过程（Gaussian Process, GP）描述的是一组由输入 $x$ 索引的随机变量：
\[\{f(x)\mid x\in\mathcal X\}\]若随机函数 $f$ 服从均值函数 $m$ 和协方差函数 $k$ 所定义的高斯过程，通常记为：
\[f\sim\mathcal{GP}(m,k)\]其中均值函数定义为：
\[m(x)=\mathbb E[f(x)]\]而协方差函数，也称核函数（kernel function），定义为：
\[k(x,x') = \operatorname{Cov}\left(f(x),f(x')\right)\]这里 $x$ 和 $x’$ 表示输入空间 $\mathcal X$ 中任意两个输入点。$k(x,x’)$描述的是 $f(x)$ 与 $f(x’)$ 之间的统计相关性。
高斯过程的严格定义是：对于输入空间中任意有限的一组输入点
\[x_1,x_2,\ldots,x_n\in\mathcal X\]对应的函数值组成随机向量：
\[\mathbf f= \begin{bmatrix} f(x_1)\\ f(x_2)\\ \vdots\\ f(x_n) \end{bmatrix} \in\mathbb R^n\]该随机向量服从一个多元高斯分布：
\[\mathbf f \sim \mathcal N(\mathbf m,K)\]其中均值向量为：
\[\mathbf m= \begin{bmatrix} m(x_1)\\ m(x_2)\\ \vdots\\ m(x_n) \end{bmatrix}\]协方差矩阵为：
\[K= \begin{bmatrix} k(x_1,x_1)&k(x_1,x_2)&\cdots&k(x_1,x_n)\\ k(x_2,x_1)&k(x_2,x_2)&\cdots&k(x_2,x_n)\\ \vdots&\vdots&\ddots&\vdots\\ k(x_n,x_1)&k(x_n,x_2)&\cdots&k(x_n,x_n) \end{bmatrix}\]即：
\[K_{ij}=k(x_i,x_j)\]对于 $d$ 维输入，每个实验点本身是一个 $d$ 维向量：
\[x_i= \begin{bmatrix} x_{i1}\\ x_{i2}\\ \vdots\\ x_{id} \end{bmatrix} \in\mathbb R^d\]将 $n$ 个实验点排列起来，可形成设计矩阵：
\[X= \begin{bmatrix} x_1^\mathsf T\\ x_2^\mathsf T\\ \vdots\\ x_n^\mathsf T \end{bmatrix} \in\mathbb R^{n\times d}\]这里 $n$ 表示已经完成的实验数量，$d$ 表示每个实验包含的参数数量。
需要特别注意，虽然：
\[X\in\mathbb R^{n\times d}\]但 GP 的协方差矩阵仍然满足：
\[K\in\mathbb R^{n\times n}\]因为矩阵 $K$ 描述的是 $n$ 个实验点对应的 $n$ 个函数值之间的两两协方差，而不是 $d$ 个输入变量之间的协方差。
一个高斯过程由均值函数 $m(x)$ 和合法的协方差函数 $k(x,x’)$共同确定，其中 $k$ 必须产生对任意有限输入集合都为半正定的协方差矩阵。
在实际贝叶斯优化中，经常在经过中心化处理之后的prior mean采用零均值假设：$m(x)=0$
此时 GP 可以简写为：
\[f\sim\mathcal{GP}(0,k)\]于是未知函数的结构假设主要由核函数 $k$ 决定。
RBF 核
该黑盒优化器使用的是典型的 RBF 核，即：
\[k(x,x') = \sigma_f^2 \exp\left( -\frac{\lVert x-x'\rVert^2}{2\ell^2} \right)\]其中：
- $\sigma_f^2$ 是信号方差；
- $\ell$ 是长度尺度（length scale）；
- $\lVert x-x’\rVert$ 是两个输入点之间的欧氏距离。
当两个点非常接近时：
\[\lVert x-x'\rVert\rightarrow0\]于是：
\[k(x,x')\rightarrow\sigma_f^2\]说明模型认为两点的函数值高度相关。
反之，当两点相距很远：
\[\lVert x-x'\rVert\gg\ell\]则：
\[k(x,x')\rightarrow0\]意味着模型不再认为两处函数值具有明显相关性。
因此，长度尺度 $\ell$ 实际控制的是模型对函数“变化速度”的判断。
较大的 $\ell$ 对应较平滑、缓慢变化的函数；较小的 $\ell$ 则允许函数在较短距离内发生明显变化。
输入归一化
实际实验变量通常具有完全不同的单位和数量级。
例如温度变量：
\[T\in[100,300]\]而另一个变量可能满足：
\[c\in[1,10]\]如果直接计算欧氏距离，第一个变量的数值变化天然比第二个变量大得多，于是 RBF 核会错误地认为第一个变量更加重要。
因此实现中首先对每一个维度进行：
\[z_j= \frac{x_j-l_j}{u_j-l_j}\]其中 $l_j$ 和 $u_j$ 分别是第 $j$ 个变量的上下界。
这样所有变量都被映射到：
\[z_j\in[0,1]\]这样可以避免不同量纲互相压制。
如果 GP 在归一化空间使用统一长度尺度 $\ell$，那么它对应到第 $j$ 个原始变量上的有效尺度为：
\[\ell_j^{\mathrm{raw}} = \ell(u_j-l_j)\]因此 $\ell=0.25$ 意味着
输入变化达到该维度完整搜索区间的约 $25\%$ 时，函数相关性开始明显下降。
这使长度尺度具有统一的无量纲解释。
GP 后验分布推导
已知实验点与未知预测点的联合分布
假设已经有 $n$ 个实验：
\[X= \begin{bmatrix} x_1\\ \vdots\\ x_n \end{bmatrix}\]观测结果为：
\[\mathbf y= \begin{bmatrix} y_1\\ \vdots\\ y_n \end{bmatrix}\]现在考虑一个新的候选点 $x_\ast$。
根据 GP 定义，训练点函数值与未知点函数值具有联合高斯分布：
\[\begin{bmatrix} \mathbf f\\ f_\ast \end{bmatrix} \sim \mathcal N \left( \begin{bmatrix} \mathbf 0\\ 0 \end{bmatrix}, \begin{bmatrix} K & k_\ast\\ k_\ast^\mathsf T & k_{\ast\ast} \end{bmatrix} \right)\]其中：
\[k_\ast= \begin{bmatrix} k(x_1,x_\ast)\\ \vdots\\ k(x_n,x_\ast) \end{bmatrix}\]并且：
\[k_{\ast\ast}=k(x_\ast,x_\ast)\]如果实验观测包含独立高斯噪声，那么训练数据的协方差变为：
\[K_y=K+\sigma_n^2I\]于是联合分布写为：
\[\begin{bmatrix} \mathbf y\\ f_\ast \end{bmatrix} \sim \mathcal N \left( \mathbf 0, \begin{bmatrix} K_y & k_\ast\\ k_\ast^\mathsf T & k_{\ast\ast} \end{bmatrix} \right)\]条件高斯分布
对于一般联合高斯变量：
\[\begin{bmatrix} a\\ b \end{bmatrix} \sim \mathcal N \left( \begin{bmatrix} \mu_a\\ \mu_b \end{bmatrix}, \begin{bmatrix} A&C\\ C^\mathsf T&B \end{bmatrix} \right)\]条件分布满足：
\[b\mid a \sim \mathcal N \left( \mu_b+C^\mathsf TA^{-1}(a-\mu_a), B-C^\mathsf TA^{-1}C \right)\]将 GP 中的对应项代入，可以得到：
\[f_\ast\mid X,\mathbf y,x_\ast \sim \mathcal N(\mu_\ast,\sigma_\ast^2)\]其中后验均值为：
\[\mu_\ast = k_\ast^\mathsf T K_y^{-1} \mathbf y\]后验方差为：
\[\sigma_\ast^2 = k_{\ast\ast} - k_\ast^\mathsf T K_y^{-1} k_\ast\]这两个公式构成整个 GP 贝叶斯优化的数学基础。
后验均值
后验均值：
\[\mu_\ast = k_\ast^\mathsf T K_y^{-1} \mathbf y\]可以理解成对已有实验结果的核加权推断。
如果 $x_\ast$ 接近某些已经实验过的点，那么对应的核值较大，这些实验结果会对预测产生较强影响。
如果 $x_\ast$ 位于一个完全没有探索过的区域，则 $k_\ast$ 整体较小，模型就不会过度相信已有观测能够准确预测该区域。
后验方差
后验方差：
\[\sigma_\ast^2 = k_{\ast\ast} - k_\ast^\mathsf T K_y^{-1} k_\ast\]表达的不是实验本身的目标值，而是模型对于预测有多不确定。
在已经充分采样的区域，候选点和已有实验点高度相关，因此第二项较大，$\sigma_\ast^2$ 较小。
在几乎没有实验数据的区域，第二项较小，于是 $\sigma_\ast^2$ 较大。
因此 GP 同时给出了：
\[\text{预测值} + \text{预测不确定性}\]“哪里还不知道”与“哪里看起来最好”同样重要。
实现中使用 Cholesky 分解
理论公式中频繁出现：$K_y^{-1}\mathbf y$ 但数值计算时通常不应该显式求逆矩阵。
实现首先对协方差矩阵进行 Cholesky 分解：$K_y=LL^\mathsf T$
其中 $L$ 是下三角矩阵。随后求解：$Lz=\mathbf y$ 以及：$L^\mathsf T\alpha=z$
于是：$\alpha=K_y^{-1}\mathbf y$ GP 后验均值即可写成：$\mu_\ast=k_\ast^\mathsf T\alpha$
对于方差，同样先求解：$Lv=k_\ast$ 则：$\sigma_\ast^2=k_{\ast\ast}-v^\mathsf Tv$
使用 Cholesky 分解不仅是性能优化，更重要的是数值稳定性。
标准化
实现并没有直接使用原始 $y$ 建立 GP，而是首先计算：
\[\bar y = \frac{1}{n} \sum_{i=1}^{n}y_i\]以及样本标准差 $s_y$，然后构造：
\[z_i= \frac{y_i-\bar y}{s_y}\]GP 实际拟合的是标准化后的 $z_i$。
预测完成后再恢复到原始尺度：
\[\mu_y = \bar y+s_y\mu_z\]以及：
\[\sigma_y=s_y\sigma_z\]实现中的预测函数正执行了这一逆变换。
这样做有两个明显好处。
第一，GP 的信号方差和噪声参数不需要随着 $y$ 的数量级剧烈变化。
第二，Expected Improvement 中的探索参数 $\xi$ 可以具有相对稳定的无量纲含义。
如果一个问题中的 $y$ 大约为 $0.01$，另一个问题中的 $y$ 大约为 $10^6$，直接使用同一个 $\xi$ 显然没有可比性。标准化以后，这一问题大幅减弱。
超参数
RBF GP 至少存在三个重要参数：
\[\theta= (\ell,\sigma_f^2,\sigma_n^2)\]其中 $\ell$ 决定函数变化尺度，$\sigma_f^2$ 决定潜在函数变化幅度，$\sigma_n^2$ 则描述观测噪声。
对数边际似然
GP 可以通过边际似然评价一组超参数是否合理：
\[p(\mathbf y\mid X,\theta) = \mathcal N(\mathbf y\mid0,K_\theta)\]其中：
\[K_\theta = K(X,X;\theta)+\sigma_n^2I\]展开多元高斯密度：
\[p(\mathbf y\mid X,\theta) = \frac{1} {(2\pi)^{n/2}|K_\theta|^{1/2}} \exp \left( -\frac{1}{2} \mathbf y^\mathsf T K_\theta^{-1} \mathbf y \right)\]取对数得到：
\[\log p(\mathbf y\mid X,\theta) = -\frac{1}{2} \mathbf y^\mathsf T K_\theta^{-1} \mathbf y -\frac{1}{2} \log|K_\theta| -\frac{n}{2}\log(2\pi)\]第一项 $-\frac{1}{2} \mathbf y^\mathsf T K_\theta^{-1} \mathbf y $ 衡量模型对已有数据的拟合程度。
第二项 $ -\frac{1}{2}\log|K_\theta| $ 可以理解为对模型复杂度的惩罚（正则化）。
第三项对于固定数据集而言只是常数。
因此最大化边际似然是在 数据拟合 与 模型复杂度 之间自动进行权衡。
当前实现如何选择超参数
该实现没有运行连续的梯度优化器，而是定义长度尺度、信号方差和噪声方差的离散候选集合，然后计算每个组合的对数边际似然，选择分数最高的一组参数。
例如长度尺度候选覆盖从很短到大于整个归一化搜索区间的多个尺度，而噪声也覆盖近似无噪声到较明显噪声的情况。
这是一种工程上相对稳健的策略。它牺牲了一部分超参数估计精度，却避免了在浏览器端再实现一个复杂的连续优化器。
因此这里应该区分两个层面的优化：
\[\text{外层：寻找最优实验参数 }x\]以及：
\[\text{内层：寻找合理 GP 超参数 }\theta\]两者并不是同一个问题。
贝叶斯优化的权衡：Exploration 与 Exploitation
假设 GP 已经完成拟合。对于每个候选点 $x$，我们现在拥有：
\[f(x)\mid\mathcal D \sim \mathcal N(\mu(x),\sigma^2(x))\]接下来应该在哪里做实验？一种策略是直接选择预测均值最优的位置。例如最小化问题可以使用：
\[x_{n+1} = \arg\min_x\mu(x)\]但这实际上非常危险。因为代理模型早期只有少量数据，预测均值本身可能存在很大误差。如果始终选择当前预测最优区域，就可能永远不会探索那些“模型暂时不知道，但实际上可能更好”的区域。
因此贝叶斯优化必须同时考虑两个因素：
Exploitation： $\mu(x)$ 候选点看起来是否足够好。
Exploration： $\sigma(x)$ 候选点是否存在足够大的未知性。
贝叶斯优化的采集函数就是把这两种需求转换成一个可优化的标量。
实验价值
假设有两个候选点。
点 A：
\[\mu_A=0.20,\qquad\sigma_A=0.01\]点 B：
\[\mu_B=0.24,\qquad\sigma_B=0.15\]如果目标是最小化，只看预测均值显然会选择 A。
但 A 的不确定性极小，意味着模型几乎已经知道那里会发生什么。
B 虽然平均预测略差，却存在很大的未知性，因此真实结果可能远好于当前最优值。
贝叶斯优化实际上在不断计算：
\[\text{已知的好} \quad\text{与}\quad \text{未知的可能更好}\]之间的价值。
综合评估这种价值的函数是采集函数，常见的采集函数有：
- EI (Expected Improvement, 期望提升)
- UCB (Upper Confidence Bound, 上置信界)
- TS (Thompson Sampling, 汤普森采样)
- Entropy Search (熵搜索) MES (Max-value Entropy Search) PES (Predictive Entropy Search)
- PI (Probability of Improvement, 改进概率)
Expected Improvement 的完整推导
当前实现使用 Expected Improvement（EI）作为采集函数。
EI 的思想非常直接：
如果在某一点进行实验，相对于当前最优结果，它平均能够带来多少改进？
Improvement 随机变量
先考虑最小化问题。设当前最好观测值为 $f_{\min}$。为了要求新实验至少产生一定程度的改进，引入探索参数 $\xi\ge0$。
定义 improvement：$I(x)=\max\left(f_{\min}-f(x)-\xi,0\right)$
由于 GP 认为 $f(x)$ 是随机变量：$f(x)\sim\mathcal N(\mu(x),\sigma^2(x))$
所以 $I(x)$ 也是随机变量。
Expected Improvement 定义为：$EI(x)=\mathbb E[I(x)]$
即：
\[EI(x) = \mathbb E \left[ \max \left( f_{\min}-f(x)-\xi, 0 \right) \right]\]只有当 $f(x)<f_{\min}-\xi$ 时 improvement 才大于零。
因此：
\[EI(x) = \int_{-\infty}^{f_{\min}-\xi} \left( f_{\min}-f-\xi \right) p(f)\,df\]令：
\[a=f_{\min}-\xi\]则：
\[EI(x) = \int_{-\infty}^{a} (a-f) p(f)\,df\]由于：
\[f\sim\mathcal N(\mu,\sigma^2)\]进行标准化：
\[z= \frac{f-\mu}{\sigma}\]以及：
\[Z= \frac{a-\mu}{\sigma} = \frac{f_{\min}-\mu-\xi}{\sigma}\]于是：
\[f=\mu+\sigma z\]代入积分：
\[EI(x) = \int_{-\infty}^{Z} (a-\mu-\sigma z) \phi(z)\,dz\]拆开得到：
\[EI(x) = (a-\mu) \int_{-\infty}^{Z} \phi(z)\,dz - \sigma \int_{-\infty}^{Z} z\phi(z)\,dz\]根据标准正态分布性质：
\[\int_{-\infty}^{Z}\phi(z)\,dz = \Phi(Z)\]标准正态密度满足：
\[\phi'(\varepsilon) = -\varepsilon\phi(\varepsilon)\]因此：
\[\int_{-\infty}^{Z} \varepsilon\phi(\varepsilon)\,d\varepsilon = -\phi(Z)\]因此：
\[EI(x) = (a-\mu)\Phi(Z) + \sigma\phi(Z)\]代回 $a=f_{\min}-\xi$：
\[EI(x) = (f_{\min}-\mu-\xi)\Phi(Z) + \sigma\phi(Z)\]其中：
\[Z= \frac{f_{\min}-\mu-\xi}{\sigma}\] \[EI(x) = \sigma \left[ Z\Phi(Z)+\phi(Z) \right]\]这就是最小化问题的 Expected Improvement。
为什么 EI 自动实现探索与利用
\[EI = \sigma \left[ Z\Phi(Z)+\phi(Z) \right]\]可以看到 EI 实际由两个因素决定。
一个是不确定性的绝对尺度：$\sigma$
另一个是标准化后的改善优势：$Z=\frac{f_{\min}-\mu-\xi}{\sigma}$
$Z$ 可以理解成：
\[\frac{\partial EI}{\partial Z} = \sigma \cdot \frac{d}{dZ}[Z\Phi(Z) + \phi(Z)] = \sigma[\underbrace{\Phi(Z) + Z\phi(Z)}_{\frac{d}{dZ}[Z\Phi(Z)]} + \underbrace{(-Z\phi(Z))}_{\frac{d}{dZ}\phi(Z)}] = \sigma \Phi(Z) > 0\]当前候选点的预测均值相对于 improvement threshold 有多大的优势，并且这个优势是多少个 posterior standard deviation。
EI 关于 $\sigma$ 和 $Z$ 都是增函数，这意味着 EI 平衡了探索（Exploration）与利用（Exploitation）。
$\xi$ 的作用：虽然 EI 自动平衡探索与利用，但超参数 $\xi$ (exploration parameter) 仍然可以微调这种平衡：
- $\xi = 0$：纯利用，只关注期望改进。
- $\xi > 0$：增加探索，随着 $\xi$ 增大，一个候选点必须具有更明显的潜在改善，才能仅依靠预测均值获得较大的 EI。因此模型会相对更加重视不确定性 $\sigma$，搜索行为趋向探索。
该网页实现 adaptiveXi 并没有让 $\xi$ 永远固定，而是根据实验次数逐渐衰减。基础探索强度分为低、中、高三个等级，然后根据实验数量 $n$ 使用衰减因子：
并令：
\[\xi_n=\xi_0d(n)\]实验早期数据不足，应该允许更多探索；随着数据积累，模型逐渐可信，可以提高 exploitation 的比例。
候选点采样 EI 最大值
理论上，下一实验点应该满足：
\[x_{n+1} = \arg\max_{x\in\mathcal X}EI(x)\]问题在于，EI 本身是由 GP 后验均值和方差构成的非线性函数。在高维空间中，它通常是非凸的，并可能存在多个局部极值。
怎样找到采集函数的最大值？
成熟 BO 框架通常可以使用多起点梯度优化、DIRECT、CMA-ES、L-BFGS 或其他全局优化方法。
当前实现采用了一种更容易在浏览器中简单稳定执行的方法：候选点采样。
生成大量候选点：${x^{(1)},x^{(2)},\ldots,x^{(M)}}$
计算：$EI_i=EI(x^{(i)})$
选择：$x_{n+1}=\arg\max_i EI_i$
如果所有候选点都均匀随机生成，那么算法具有不错的全局覆盖能力，但在最优区域附近的搜索精度可能不足。如果所有候选点都围绕当前最优点生成，则算法容易陷入局部最优。因此实现采用混合策略。
大约 $70\%$ 的候选点在整个参数空间中随机产生，而约 $30\%$ 的候选点围绕当前最优实验点进行局部扰动。
1
2
3
4
5
6
7
8
9
10
11
12
13
参数空间
+ 表示全局候选点
· 表示局部候选点
+ + + + ← global
+ + +
· ·
· · · ·
· · ★ · · ★ = 当前真实实验最优点
· · ·
·
+ + + ← global
+ +
局部扰动表示为：
\[x'_j = x_j^{\mathrm{best}} + \epsilon_j s_n (u_j-l_j)\]其中：
\[\epsilon_j\sim\mathcal N(0,1)\]而 $s_n$ 是随实验数量缓慢变化的步长比例。$(u_j-l_j)$ 是每个参数自身的区间宽度缩放。
在 $d$ 维单位超立方体中，固定数量的随机点随着 $d$ 增大将迅速变得稀疏。这就是高维空间中的“维度灾难”。实现采用经验公式：$M(d) \approx 1800d^{1.45} $ 同时将候选点数量限制在一定范围内。否则在高维空间中，仅使用几千个随机候选点可能根本无法充分覆盖高 EI 区域。但另一方面，这也暴露出候选采样法的局限：维度继续升高以后，单纯增加候选数量并不能真正克服维度灾难。
初始实验使用 Latin Hypercube Sampling
贝叶斯优化不能从完全没有数据的状态直接得到有意义的 GP 后验。因此通常需要首先执行若干初始实验。最简单的方法是完全随机采样，但随机点可能发生明显聚集，导致有限实验预算被浪费。
Latin Hypercube Sampling（LHS）则把每个维度划分为 $n$ 个区间，并保证每一个区间恰好被使用一次。
对于第 $j$ 个维度，可以将第 $i$ 个样本写成：
\[z_{ij} = \frac{\pi_j(i)+u_{ij}}{n}\]其中 $\pi_j$ 是该维度上的随机排列，而：
\[u_{ij}\sim U(0,1)\]因此：
\[z_{ij}\in[0,1]\]最后再映射回真实参数范围：
\[x_{ij} = l_j+z_{ij}(u_j-l_j)\]与简单随机采样相比，LHS 能确保每个单独维度都具有较好的边际覆盖。
当前实现正是按照这一方法生成初始实验点。对于昂贵实验，这一点尤其重要，因为最初几个实验决定了 GP 对整个参数空间的第一轮认识。
算法流程：从初始化到贝叶斯优化
第一步：搜索空间
指定：
\[x_j\in[l_j,u_j]\]并将每个维度归一化为：
\[z_j\in[0,1]\]第二步：生成初始设计
使用 LHS 获得：
\[X_0=\{x_1,\ldots,x_{n_0}\}\]完成真实实验：
\[y_i=f(x_i)+\varepsilon_i\]得到初始数据集：
\[\mathcal D_{n_0} = \{(x_i,y_i)\}_{i=1}^{n_0}\]第三步：标准化输出
计算：
\[z_i= \frac{y_i-\bar y}{s_y}\]以提高 GP 与 EI 参数的尺度稳定性。
第四步：估计 GP 超参数
寻找：
\[\theta^\star = \arg\max_\theta \log p(\mathbf y\mid X,\theta)\]得到长度尺度、信号方差与噪声水平。
第五步：计算 GP 后验
对于任意候选点 $x$：
\[f(x)\mid\mathcal D_n \sim \mathcal N(\mu_n(x),\sigma_n^2(x))\]第六步：通过采集函数EI选择下一实验点
对于所有候选点计算EI值
选择：
\[x_{n+1} = \arg\max_x EI_n(x)\]第七步：执行真实实验
获得：
\[y_{n+1} = f(x_{n+1})+\varepsilon_{n+1}\]并更新：
\[\mathcal D_{n+1} = \mathcal D_n \cup \{(x_{n+1},y_{n+1})\}\]然后重新拟合 GP。
整个过程不断循环：
\[\mathcal D_n \rightarrow GP_n \rightarrow EI_n \rightarrow x_{n+1} \rightarrow y_{n+1} \rightarrow \mathcal D_{n+1}\]原实现的交互方式实际上严格遵循这一闭环：系统给出一个待实验点，用户完成真实实验并填写 $y$，数据被加入历史记录，随后 GP 参数和后续推荐根据新数据重新计算。
贝叶斯优化适用场景
单次实验成本高
这是最典型的适用条件。
例如：
- 材料配方实验；
- 化学反应条件优化；
- 制造工艺参数优化；
- 生物实验条件优化；
- 机器人真实环境控制参数；
- CFD 或有限元仿真参数；
- 机器学习超参数优化；
- 需要数分钟、数小时甚至数天才能完成一次评估的问题。
如果总共只能进行几十次或几百次实验，那么每一个实验点都需要尽可能有价值。这正是 BO 相对于随机搜索最明显的优势。
如果一次函数计算只需要几微秒，那么 BO 中 GP 拟合、采集函数计算等额外开销可能反而没有意义。
黑盒：没有解析表达式或梯度
贝叶斯优化只需要：$x\rightarrow y$
而不要求知道：$\nabla f(x)$
因此它天然适用于黑盒系统。
输入维度中低
经典 GP 贝叶斯优化通常最适合低维到中等维度问题。（一般是小于20维）
没有一个绝对维度界限，但随着 $d$ 增大，两个问题会越来越明显：$\text{空间覆盖难度}\uparrow$ $\text{采集函数优化难度}\uparrow$
因此几十维甚至上百维问题通常需要更专门的方法，例如稀疏 GP、随机嵌入、trust-region BO 或结构化核。
对于当前这种依靠随机候选集寻找最大 EI 的轻量级实现，高维问题尤其需要谨慎。
哪些问题不适合使用这种 BO
目标函数极不连续
RBF 核隐含了较强的平滑性假设。
如果目标函数存在大量突变：
\[x\rightarrow x+\delta\]却可能导致：
\[|f(x+\delta)-f(x)|\]突然变得非常大，那么 RBF GP 可能产生错误的平滑外推。
大量离散或类别变量
例如：
- 算法 A、B、C；
- 材料类型；
- 开关变量；
- 树状条件参数。
这些变量不存在普通欧氏距离意义上的“接近”。
直接把类别编码为 $0,1,2$ 会隐含错误的距离关系。
这类问题需要 categorical kernel、mixed-variable BO 或其他专门方法。
超高维搜索空间
当 $d$ 很大时，有限候选点很难覆盖空间。
即使：
\[M(d)\propto d^{1.45}\]也无法从根本上消除维度灾难。
每次实验几乎没有成本
如果 $f(x)$ 可以高速计算，那么直接使用成熟的全局优化方法、进化算法甚至大规模随机搜索可能更加合适。
贝叶斯优化的优势建立在：
\[C_{\mathrm{experiment}} \gg C_{\mathrm{BO}}\]这一成本关系上。
当前实现中的其他工程化设计
避免重复实验
在昂贵实验中，推荐一个几乎已经测试过的参数组合通常没有意义。因此实现计算候选点与已有实验点在归一化空间中的最近距离，并过滤极度接近已有数据的候选点。
\[d_{\min}(x) = \min_i \lVert z(x)-z(x_i)\rVert_2\]若 $ d_{\min}(x)<\varepsilon $ 则跳过该候选点。
GP 协方差矩阵 数值jitter
实际计算协方差矩阵时，即使理论上矩阵应该是正定的，也可能由于浮点误差变得接近奇异。
因此工程实现使用：
\[K \leftarrow K+ (\sigma_n^2+\epsilon)I\]其中 $\epsilon$ 是一个非常小的正数。
它与实验噪声 $\sigma_n^2$ 的统计意义不同，主要目的是保证 Cholesky 分解稳定。
一维情况下可视化
对于一维问题，系统绘制：$\mu(x)$
以及：$\mu(x)\pm1.96\sigma(x)$
作为近似 $95\%$ 的后验区间。
这个可视化过程实际上很好地展示了 BO 的核心：
- 实验点附近，不确定性收缩；
- 未探索区域，不确定性扩大；
- 新实验会重新塑造整条后验均值曲线；
- 采集函数EI 会综合后验均值与不确定性决定下一点。
当前轻量实现局限性
当前实现属于一个相对轻量的 BO 系统。但如果用于更加复杂的生产问题，还存在进一步扩展空间。
单一长度尺度
当前 RBF 核使用统一的 $\ell$：
\[k(x,x') = \sigma_f^2 \exp \left( -\frac{\lVert x-x'\rVert^2}{2\ell^2} \right)\]这意味着所有归一化维度具有相同的相关长度。
更一般的 Automatic Relevance Determination（ARD）核可以写成：
\[k(x,x') = \sigma_f^2 \exp \left( -\frac{1}{2} \sum_{j=1}^{d} \frac{(x_j-x'_j)^2}{\ell_j^2} \right)\]每个变量都有自己的 $\ell_j$。
如果某个 $\ell_j$ 很大，意味着目标函数对第 $j$ 个变量不敏感；如果 $\ell_j$ 很小，则说明该变量可能引起快速变化。
网格搜索超参数
当前实现：
\[\theta^\star = \arg\max_{\theta\in\Theta_{\mathrm{grid}}} \log p(\mathbf y\mid X,\theta)\]工业级实现则常直接求解：
\[\theta^\star = \arg\max_{\theta} \log p(\mathbf y\mid X,\theta)\]并使用梯度优化或多起点优化。
网格搜索更简单稳定，连续优化则具有更高精度和更大的参数空间。
采集函数优化
当前方法从有限候选集选择：
\[x^\star = \arg\max_{x\in\mathcal C}EI(x)\]而更完整的方法会尝试求：
\[x^\star = \arg\max_{x\in\mathcal X}EI(x)\]这两者的差异在高维空间尤其明显。
单点 BO 与批量 BO
当前流程一次推荐一个实验：
\[x_{n+1}\]实验完成以后才推荐下一个。
如果实验室可以同时运行 $q$ 个实验，更适合 batch Bayesian Optimization：
\[X_{\mathrm{batch}} = \{x_{n+1},\ldots,x_{n+q}\}\]此时需要处理不同候选实验之间的信息冗余，常使用 qEI 等方法，而不能简单地取 EI 最大的前 $q$ 个点。
计算复杂度
对于标准 GP，训练阶段需要对 $n\times n$ 协方差矩阵进行 Cholesky 分解。其时间复杂度大约为：$O(n^3)$ 存储复杂度约为：$O(n^2)$
而对一个新的候选点进行预测，在已有 Cholesky 分解的情况下通常需要约：$O(n^2)$ 的计算。
因此 GP 贝叶斯优化非常适合：
\[n\text{ 较小，但每次实验非常昂贵}\]的场景。
这恰好与 BO 的设计目标一致。
总结
贝叶斯优化可以最终浓缩为三个数学对象。
第一个是未知目标函数：
\[f(x)\]第二个是基于已有数据得到的概率代理模型：
\[f(x)\mid\mathcal D_n \sim \mathcal N(\mu_n(x),\sigma_n^2(x))\]第三个是将“预测值”和“不确定性”转化为实验价值的采集函数：
\[a_n(x) = EI_n(x)\]然后选择：
\[x_{n+1} = \arg\max_x a_n(x)\]完成真实实验以后：
\[\mathcal D_{n+1} = \mathcal D_n \cup \{(x_{n+1},y_{n+1})\}\]再重新开始。
这也是贝叶斯优化能够在昂贵黑盒问题中，用有限实验预算逐渐逼近最优实验超参数设计的根本原因。
