贝叶斯优化

贝叶斯优化（Bayesian optimization）为评价昂贵的未知目标建立概率代理模型，再用采集函数选择下一次评价位置。概率描述的是模型对尚未评价函数值的不确定性；它不是凭空获得了真实目标的导数，也不保证每次新评价都改善当前最好结果。

一次观测后的代理预测

设候选集合是教学用的 $\{0,1,2\}$，真实目标只在评价时可见，本例为 $f(x)=(x-1{)}^2$，要最小化它。先观察到 $f(0)=1$。选择零均值高斯过程先验，协方差核为 $k(x,x^{\prime })=\exp [-(x-x^{\prime }{)}^2/2]$，观测无噪声。条件正态公式给出 $$\mu (x)=k(x,0),s^2(x)=1-k(x,0{)}^2.$$ 于是候选一的均值约0.606531、标准差约0.795060；候选二的均值约0.135335、标准差约0.990800。观测点零的后验方差为零。图中的带宽来自声明的核模型，并不是由这一个观测验证出的真实误差界。

虚线为一个后验标准差带，表示模型内不确定性。

一般已观察 $y=(f(x_1),\dots ,f(x_t))$，先验均值为零，核矩阵 $K_{ij}=k(x_i,x_j)$ 可逆，向量 $k_x=(k(x,x_i){)}_i$。无噪声后验为 $\mu (x)=k_x^{𝖳}{K}^{-1}y$、$s^2(x)=k(x,x)-k_x^{𝖳}{K}^{-1}{k}_x$。实际计算通过线性系统求解，不应显式求逆；有噪声时需要在观测协方差中加入噪声模型。

期望改进的完整积分

记目前最小观测值为b。在某候选点，模型给出 $Y\sim N(\mu ,s^2)$，改进量为 $(b-Y{)}_{+}=max(b-Y,0)$。当 $s>0$，设 $z=(b-\mu )/s$，标准化后 $$\mathrm{EI}={\int }_{-\mathrm{\infty }}^z(b-\mu -su)\phi (u)\mathrm{d}u=(b-\mu )\Phi (z)+s\phi (z).$$ 最后一步使用 ${\phi }^{\prime }(u)=-u\phi (u)$，所以 ${\int }_{-\mathrm{\infty }}^{z}u\phi (u)\mathrm{d}u=-\phi (z)$。当s为零时，改进量确定，EI为 $max(b-\mu ,0)$。EI与目标同单位，不能解释为一个概率。

有限候选一、二中，二的EI更大；连线展示同一核模型下的采集函数。

主例取b为一，计算两点EI并选择较大的候选二。实际评价却得到 $f(2)=1$，并没有改善最好值；它仍提供了新信息。再条件化于零和二的观测后，只剩候选一未评价，取值为零。这个有限集合的最优性最终来自三个候选都已评价，而不是EI值本身给出的全局证书。

探索变式与模型边界

变式比较两个候选：A的均值0.8且标准差零，EI为0.2；B的均值一、标准差一，EI为 $1/\sqrt{2\pi }\approx 0.398942$。模型均值较差的B反而更值得探索，因为其下尾存在较大改进机会。反过来，过于自信或不合适的核也可能低估某区域的不确定性，导致有限预算内错过真实谷底。

代理模型、核参数估计、采集函数和采集函数的数值优化共同决定算法。噪声观测下“最小观测值”可能只是偶然偏低，不能无说明地沿用无噪声EI。连续域中的有限观测不排除尚未看见的窄谷；模型内的后验区间也不自动成为对任意目标成立的频率覆盖保证。

三个候选的可运行演示

代码只用标准库；输入为互异的三个有限标量候选，固定从第一个开始评价。返回观测与每次后验/采集记录，不改写输入。小线性系统用带主元的消元法，失败时显式报错。这份朴素实现每个候选各解系统，t个观测、M个候选时约为 $O(M{t}^3)$；正式实现复用分解可减少重复成本。内存主要为 $O(t^2+M)$，不含返回历史。

Python 3

import math

def solve_system(A,b):
    n = len(b)
    M = [list(A[i])+[float(b[i])] for i in range(n)]
    for j in range(n):
        pivot = max(range(j,n),key=lambda i:abs(M[i][j]))
        M[j],M[pivot] = M[pivot],M[j]
        q = M[j][j]
        if abs(q) < 1e-12:
            raise ArithmeticError("singular covariance")
        M[j] = [v/q for v in M[j]]
        for i in range(n):
            if i != j:
                q = M[i][j]
                M[i] = [v-q*w for v,w in zip(M[i],M[j])]
    return [row[-1] for row in M]

def bayes_three(candidates=(0.0,1.0,2.0)):
    points = tuple(map(float,candidates))
    if len(points) != 3 or len(set(points)) != 3:
        raise ValueError("three distinct candidates required")
    kernel = lambda x,y:math.exp(-(x-y)**2/2)
    objective = lambda x:(x-1)**2
    observed = [(points[0],objective(points[0]))]
    trace = []
    while len(observed) < 3:
        K = [[kernel(x,y) for y,_ in observed] for x,_ in observed]
        weights = solve_system(K,[y for _,y in observed])
        best = min(y for _,y in observed)
        trials = []
        for x in points:
            if any(x == old for old,_ in observed):
                continue
            k = [kernel(x,old) for old,_ in observed]
            mu = sum(a*b for a,b in zip(k,weights))
            variance = max(0.0,1-sum(a*b for a,b in zip(k,solve_system(K,k))))
            sd = math.sqrt(variance)
            if sd < 1e-12:
                ei = max(best-mu,0.0)
            else:
                z = (best-mu)/sd
                ei = (best-mu)*(1+math.erf(z/math.sqrt(2)))/2+sd*math.exp(-z*z/2)/math.sqrt(2*math.pi)
            trials.append((x,mu,sd,ei))
        x = max(trials,key=lambda row:row[3])[0]
        trace.append(trials)
        observed.append((x,objective(x)))
    return observed,trace

observations,trace = bayes_three()
for row in trace:
    print("candidates",row)
print("observed",observations)
assert observations[-1] == (1.0,0.0)