引言
在机器学习实践中，模型超参数的选择直接决定了最终性能。网格搜索与随机搜索虽然简单，但在高维空间中效率极低。贝叶斯优化（Bayesian Optimization）是一种高效的全局优化方法，通过构建目标函数的概率代理模型，在有限的评估次数内找到近似最优解。它已被广泛应用于神经网络架构搜索、超参数调优、材料设计等领域。
本文将用Java实现完整的贝叶斯优化框架，核心包括高斯过程回归建模目标函数的不确定性，以及采集函数（UCB/EI）平衡探索与利用的采样策略。
核心原理
1. 问题定义
寻找黑箱函数的最小值：
$$x^* = \arg\min_{x \in \mathcal{X}} f(x)$$
其中 $f(x)$ 计算代价高昂，无闭式表达式，仅能通过采样获取函数值（可能含噪声）。
2. 高斯过程回归（Gaussian Process Regression）
高斯过程是定义在函数空间上的概率分布，任意有限个采样点的函数值服从联合高斯分布：
$$f(x) \sim GP(m(x), k(x, x’))$$
其中：
– $m(x)$：均值函数，通常设为0
– $k(x, x’)$：核函数（协方差函数），衡量两点间的相似度
常用的径向基核函数（RBF）：
$$k(x, x’) = \sigma^2 \exp\left(-\frac{|x – x’|^2}{2l^2}\right)$$
给定观测数据 $D = {(x_i, y_i)}{i=1}^n$，对新点 $x$ 的预测分布为：
$$\mu(x_) = k_^T (K + \sigma_n^2 I)^{-1} y$$
$$\sigma^2(x_) = k(x_, x_) – k_^T (K + \sigma_n^2 I)^{-1} k_$$
其中 $K_{ij} = k(x_i, x_j)$，$k_ = [k(x_, x_1), …, k(x_*, x_n)]^T$。
3. 采集函数（Acquisition Function）
采集函数指导下一个采样点的选择，在探索（高不确定性区域）与利用（高期望收益区域）之间取得平衡。
上置信界（UCB）
$$UCB(x) = \mu(x) + \kappa \cdot \sigma(x)$$
其中 $\kappa$ 控制探索程度，值越大越偏向探索未知区域。
期望改进（Expected Improvement, EI）
$$EI(x) = E[\max(0, f_{best} – f(x))]$$
对于高斯分布，EI有闭式解：
$$EI(x) = (f_{best} – \mu(x)) \Phi(Z) + \sigma(x) \phi(Z)$$
其中 $Z = \frac{f_{best} – \mu(x)}{\sigma(x)}$，$\Phi$ 和 $\phi$ 分别为标准正态分布的CDF和PDF。
Java实现
import java.util.*;
/**
* 贝叶斯优化框架
* 高斯过程回归 + UCB/EI采集函数
*/
public class BayesianOptimizer {
private final int dim; // 参数维度
private final double[] lowerBound;
private final double[] upperBound;
private final Kernel kernel;
private final AcquisitionFunction acquisition;
private final double noiseSigma; // 观测噪声
// 观测数据
private final List<double[]> X = new ArrayList<>();
private final List<Double> y = new ArrayList<>();
// 高斯过程缓存
private double[][] K; // 核矩阵
private double[][] KInv; // 核矩阵逆
private double[] alpha; // (K + sigma^2 I)^{-1} y
public BayesianOptimizer(int dim, double[] lowerBound, double[] upperBound,
String kernelType, String acqType, double noiseSigma) {
this.dim = dim;
this.lowerBound = lowerBound.clone();
this.upperBound = upperBound.clone();
this.noiseSigma = noiseSigma;
// 初始化核函数
if ("rbf".equalsIgnoreCase(kernelType)) {
this.kernel = new RBFKernel(1.0, 1.0); // 默认长度尺度和信号方差
} else {
this.kernel = new RBFKernel(1.0, 1.0);
}
// 初始化采集函数
if ("ucb".equalsIgnoreCase(acqType)) {
this.acquisition = new UCB(2.0); // kappa = 2.0
} else if ("ei".equalsIgnoreCase(acqType)) {
this.acquisition = new ExpectedImprovement();
} else {
this.acquisition = new UCB(2.0);
}
}
/**
* 核函数接口
*/
interface Kernel {
double eval(double[] x1, double[] x2);
}
/**
* RBF径向基核函数
*/
static class RBFKernel implements Kernel {
private final double lengthScale;
private final double signalVariance;
RBFKernel(double lengthScale, double signalVariance) {
this.lengthScale = lengthScale;
this.signalVariance = signalVariance;
}
@Override
public double eval(double[] x1, double[] x2) {
double distSq = 0;
for (int i = 0; i < x1.length; i++) {
double diff = x1[i] - x2[i];
distSq += diff * diff;
}
return signalVariance * Math.exp(-distSq / (2 * lengthScale * lengthScale));
}
}
/**
* 采集函数接口
*/
interface AcquisitionFunction {
double eval(double mean, double std, double bestY);
}
/**
* 上置信界（UCB）采集函数
*/
static class UCB implements AcquisitionFunction {
private final double kappa;
UCB(double kappa) {
this.kappa = kappa;
}
@Override
public double eval(double mean, double std, double bestY) {
// 注意：这里最小化问题，所以用LCB
return mean - kappa * std;
}
}
/**
* 期望改进（EI）采集函数
*/
static class ExpectedImprovement implements AcquisitionFunction {
@Override
public double eval(double mean, double std, double bestY) {
if (std < 1e-10) return 0;
double z = (bestY - mean) / std;
double ei = (bestY - mean) * normalCDF(z) + std * normalPDF(z);
return ei;
}
}
/**
* 标准正态分布PDF
*/
static double normalPDF(double x) {
return Math.exp(-0.5 * x * x) / Math.sqrt(2 * Math.PI);
}
/**
* 标准正态分布CDF（使用误差函数近似）
*/
static double normalCDF(double x) {
return 0.5 * (1 + erf(x / Math.sqrt(2)));
}
/**
* 误差函数近似（Abramowitz and Stegun）
*/
static double erf(double x) {
double sign = x >= 0 ? 1 : -1;
x = Math.abs(x);
double a1 = 0.254829592, a2 = -0.284496736;
double a3 = 1.421413741, a4 = -1.453152027, a5 = 1.061405429;
double p = 0.3275911;
double t = 1 / (1 + p * x);
double y = 1 - (((((a5 * t + a4) * t) + a3) * t + a2) * t + a1) * t * Math.exp(-x * x);
return sign * y;
}
/**
* 添加观测数据点
*/
public void observe(double[] x, double y) {
X.add(x.clone());
this.y.add(y);
updateGP();
}
/**
* 更新高斯过程模型
*/
private void updateGP() {
int n = X.size();
K = new double[n][n];
for (int i = 0; i < n; i++) {
for (int j = 0; j < n; j++) {
K[i][j] = kernel.eval(X.get(i), X.get(j));
if (i == j) K[i][j] += noiseSigma * noiseSigma;
}
}
// 计算 K^{-1} y
KInv = invertMatrix(K);
alpha = new double[n];
for (int i = 0; i < n; i++) {
for (int j = 0; j < n; j++) {
alpha[i] += KInv[i][j] * y.get(j);
}
}
}
/**
* 高斯过程预测：均值与标准差
*/
public double[] predict(double[] x) {
int n = X.size();
double[] kStar = new double[n];
for (int i = 0; i < n; i++) {
kStar[i] = kernel.eval(x, X.get(i));
}
double mean = 0;
for (int i = 0; i < n; i++) {
mean += kStar[i] * alpha[i];
}
double kStarStar = kernel.eval(x, x);
double v = 0;
for (int i = 0; i < n; i++) {
double sum = 0;
for (int j = 0; j < n; j++) {
sum += kStar[j] * KInv[j][i];
}
v += sum * kStar[i];
}
double variance = kStarStar - v;
double std = Math.sqrt(Math.max(0, variance));
return new double[]{mean, std};
}
/**
* 选择下一个采样点（最大化采集函数）
*/
public double[] suggestNext() {
double bestY = Double.POSITIVE_INFINITY;
for (double val : y) {
if (val < bestY) bestY = val;
}
// 使用随机采样 + 局部优化寻找采集函数最大值
double[] bestX = null;
double bestAcq = Double.NEGATIVE_INFINITY;
Random rand = new Random();
int nCandidates = 1000;
for (int i = 0; i < nCandidates; i++) {
double[] candidate = new double[dim];
for (int d = 0; d < dim; d++) {
candidate[d] = lowerBound[d] + rand.nextDouble() * (upperBound[d] - lowerBound[d]);
}
double[] pred = predict(candidate);
double acq = acquisition.eval(pred[0], pred[1], bestY);
if (acq > bestAcq) {
bestAcq = acq;
bestX = candidate.clone();
}
}
// 简单局部优化：在最佳候选点附近微调
if (bestX != null) {
for (int step = 0; step < 100; step++) {
double[] neighbor = bestX.clone();
int d = rand.nextInt(dim);
double delta = (rand.nextDouble() - 0.5) * (upperBound[d] - lowerBound[d]) * 0.1;
neighbor[d] += delta;
neighbor[d] = Math.max(lowerBound[d], Math.min(upperBound[d], neighbor[d]));
double[] pred = predict(neighbor);
double acq = acquisition.eval(pred[0], pred[1], bestY);
if (acq > bestAcq) {
bestAcq = acq;
bestX = neighbor;
}
}
}
return bestX;
}
/**
* 获取当前最优解
*/
public double[] getBest() {
int bestIdx = 0;
for (int i = 1; i < y.size(); i++) {
if (y.get(i) < y.get(bestIdx)) bestIdx = i;
}
return X.get(bestIdx).clone();
}
public double getBestValue() {
double best = Double.POSITIVE_INFINITY;
for (double val : y) {
if (val < best) best = val;
}
return best;
}
/**
* 矩阵求逆（高斯-约旦消元法）
*/
private double[][] invertMatrix(double[][] matrix) {
int n = matrix.length;
double[][] aug = new double[n][2 * n];
for (int i = 0; i < n; i++) {
System.arraycopy(matrix[i], 0, aug[i], 0, n);
aug[i][n + i] = 1;
}
for (int i = 0; i < n; i++) {
double pivot = aug[i][i];
for (int j = 0; j < 2 * n; j++) aug[i][j] /= pivot;
for (int k = 0; k < n; k++) {
if (k == i) continue;
double factor = aug[k][i];
for (int j = 0; j < 2 * n; j++) {
aug[k][j] -= factor * aug[i][j];
}
}
}
double[][] inv = new double[n][n];
for (int i = 0; i < n; i++) {
System.arraycopy(aug[i], n, inv[i], 0, n);
}
return inv;
}
/**
* 主优化流程
*/
public double[] optimize(ObjectiveFunction obj, int nIterations, int nInit) {
Random rand = new Random();
// 1. 随机初始化
for (int i = 0; i < nInit; i++) {
double[] x = new double[dim];
for (int d = 0; d < dim; d++) {
x[d] = lowerBound[d] + rand.nextDouble() * (upperBound[d] - lowerBound[d]);
}
observe(x, obj.evaluate(x));
}
// 2. 迭代优化
for (int iter = 0; iter < nIterations; iter++) {
double[] nextX = suggestNext();
observe(nextX, obj.evaluate(nextX));
}
return getBest();
}
/**
* 目标函数接口
*/
interface ObjectiveFunction {
double evaluate(double[] x);
}
/**
* 测试：优化经典测试函数
*/
public static void main(String[] args) {
// 测试1：Branin函数（2维，3个全局最小值）
BayesianOptimizer bo = new BayesianOptimizer(
2,
new double[]{-5, 0},
new double[]{10, 15},
"rbf", "ei", 1e-5
);
ObjectiveFunction branin = x -> {
double a = x[1] - 5.1/(4*Math.PI*Math.PI)*x[0]*x[0] + 5/Math.PI*x[0] - 6;
return a*a + 10*(1-1/(8*Math.PI))*Math.cos(x[0]) + 10;
};
double[] best = bo.optimize(branin, 30, 5);
System.out.println("Branin函数最优值: " + bo.getBestValue());
System.out.println("最优参数: [" + best[0] + ", " + best[1] + "]");
System.out.println("观测次数: " + bo.X.size());
// 测试2：Rosenbrock函数（高维山谷）
BayesianOptimizer bo2 = new BayesianOptimizer(
3,
new double[]{-2, -2, -2},
new double[]{2, 2, 2},
"rbf", "ucb", 1e-5
);
ObjectiveFunction rosenbrock = x -> {
double sum = 0;
for (int i = 0; i < x.length - 1; i++) {
sum += 100 * Math.pow(x[i+1] - x[i]*x[i], 2) + Math.pow(1 - x[i], 2);
}
return sum;
};
double[] best2 = bo2.optimize(rosenbrock, 50, 10);
System.out.println("\nRosenbrock函数最优值: " + bo2.getBestValue());
System.out.println("最优参数: [" + best2[0] + ", " + best2[1] + ", " + best2[2] + "]");
}
}
算法复杂度分析
时间复杂度
- 高斯过程预测：$O(n^2)$，其中 $n$ 为观测点数（核向量与逆矩阵的乘积）
- 核矩阵求逆：$O(n^3)$，使用高斯-约旦消元法。当 $n$ 较大时可改用Cholesky分解或稀疏近似
- 单次优化迭代：$O(C \cdot n^2)$，$C$ 为候选点数量
- 总优化时间：$O(N \cdot n^3)$，$N$ 为总迭代次数
空间复杂度
- 核矩阵及逆矩阵：$O(n^2)$
- 观测数据：$O(n \cdot d)$，$d$ 为维度
- 总空间：$O(n^2 + n \cdot d)$
扩展与改进
1. 核函数自适应
通过最大化边际似然（Marginal Likelihood）自动优化核函数的超参数（长度尺度 $l$ 和信号方差 $\sigma^2$）：
$$\log p(y|X) = -\frac{1}{2}y^T K^{-1} y – \frac{1}{2}\log|K| – \frac{n}{2}\log(2\pi)$$
2. 稀疏高斯过程
当观测点超过1000时，使用诱导点（Inducing Points）或随机特征（Random Features）近似核矩阵，将复杂度降至 $O(m^2 n)$ 或 $O(n)$。
3. 并行采集函数
在分布式计算场景下，使用qEI（批量期望改进）或局部惩罚方法同时推荐多个评估点，加速并行优化。
4. 约束处理
对于带约束的优化问题，使用约束采集函数如CEI（Constrained Expected Improvement），同时建模目标函数和约束违反概率。
总结
贝叶斯优化通过概率建模将”黑箱”目标函数转化为可分析的不确定性分布，以采集函数智能引导采样，在样本效率上远超传统网格/随机搜索。本文实现的Java框架包含完整的高斯过程回归、UCB/EI采集函数及优化流程，可直接用于超参数调优、实验设计等场景。其核心思想——在不确定性中做决策——也是强化学习、主动学习等领域的共同基础。
