条件密度神经网络黑盒优化：替代高斯过程的贝叶斯优化新思路

去年做某个模拟实验的参数优化时，我被一个黑盒目标函数折磨了好几天：输入十几个参数，跑一次模拟要三分钟，目标值出来一个数，没有梯度、没有解析式、甚至连函数长什么样都不清楚。后来我换了个思路——不再去拟合这个黑盒函数本身，而是用神经网络估计一个条件密度，通过密度采样来找下一步该试哪里。这个方向就是我们常说的"密度估计神经网络黑盒问题优化"，这套方案我完整实现了，代码放在文章对应小节里，跑通之后优化效率提升非常明显。这篇把我从原理、代码到踩坑的完整过程写出来，适合正在做超参搜索、贝叶斯优化，或者被黑盒目标函数折磨的读者。

1. 黑盒优化的真正痛点：不是方法少，而是假设太强

1.1 我遇到的黑盒问题到底长什么样

先还原一下典型场景。某个模拟实验的输入向量记作 x，包含几个连续参数、两个离散开关量、一个配置文件的选项；输出是一个标量 y，工程上希望 y 越大越好。真实的计算流程被封装在某个程序里，你只能拿到输入输出对，拿不到中间过程，更不可能对参数求偏导。

这种问题里，"黑盒"有三层含义：第一，目标函数没有解析表达式，数学上只能当作一个未知函数处理；第二，不可求导，梯度类优化算法直接歇菜；第三，单次评估昂贵，我之前那个场景跑一次要三分钟，意味着所有方法必须在极少的评估预算内找到好结果。

网格搜索和随机搜索是最容易想到的基线，但网格搜索的采样点数量随维度指数增长，五个连续参数每个方向分十档就是十万个点，预算完全不可能支持。随机搜索虽然不受维度爆炸影响，但它完全不做信息积累，上一轮的评估结果对下一轮没有任何帮助，纯粹赌命。

1.2 为什么"拟合函数"这条路容易走偏

很多做黑盒优化的同学第一时间想到高斯过程回归，也就是常见的贝叶斯优化方案。思路很自然：先用已有的观测点拟合一个函数替代模型，然后在替代模型上计算采集函数，选出最值得评估的点。

这套思路在低维、平滑、噪声较小的目标函数上非常有效。但一旦维度上去，高斯过程的核函数选择和超参数估计都会变得非常敏感。而且高斯过程本质上假设目标函数来自一个高斯过程先验，这个先验如果不合适，拟合出来的代理模型就会在没数据的地方给出很自信的错误预估，后续采样全被带偏。

我当时真正卡住的地方不是没有方法，而是我手里所有方法都要求我隐式地假设"目标函数服从某种分布"。但黑盒问题的本质恰恰是我们对这个函数一无所知，强加先验本身就是一种豪赌。于是我换了一个方向：不去拟合 y = f(x)，而去估计一个概率密度 p(x | y)，也就是"在给定目标值 y 的情况下，输入 x 会落在哪里"。

2. 换个思路：用条件密度回答"哪里值得试"

2.1 从TPE开始的密度估计路线

在贝叶斯优化这个框里，有一支不那么"主流"但非常实用的方法叫 TPE，全称是 Tree-structured Parzen Estimator，树结构Parzen密度估计器。TPE 的思路和传统贝叶斯优化完全不同：它不拟合目标函数的函数形式，而是把历史观测点按目标值好坏分成两组，分别用核密度估计去学习"好点"的密度和"坏点"的密度，然后选择能让好点密度最大、坏点密度最小的新点去评估。

这个方法非常优雅，因为它绕开了函数拟合的强假设，把优化问题转换成了一个密度估计问题。但它的工程实现里有个明显的短处：核密度估计用的是固定带宽的核函数，而黑盒函数的不同区域往往有着完全不同的平滑程度。有的区域目标值变化平缓，一个小带宽就足够；有的区域剧烈震荡，需要大带宽去平滑。固定一个带宽，很难同时照顾到这两种区域。

2.2 条件混合密度网络：让网络自己学带宽

为了解决固定核带宽的问题，我引入了条件混合密度网络，英文叫 Conditional Mixture Density Network，缩写 MDN。它的思路很简单：不再用预先固定位置和宽度的核，而是用一个神经网络去输出一个高斯混合分布的参数，并且这些参数可以依赖某个条件变量。

具体到我的优化目标，我要估计的是 p(x | y)。MDN 对它的建模方式是：

p(x | y) = Σ π_k(y) · N(x; μ_k(y), σ_k²(y))

这里的 π_k、μ_k、σ_k 全部是神经网络以 y 为输入计算出来的。神经网络输出的是若干个高斯成分的权重、均值和标准差。从直觉上理解，这就相当于把固定核替换成了"自适应核"——网络会自己决定每个位置需要多少个成分、成分放在哪里、每个成分有多宽，这比手动调核带宽要灵活得多。

我第一次看到这个模型的时候，觉得它就像是一个"概率版的自编码器"：给它一个条件，它不输出单一的预测值，而是输出一个完整的分布。在优化场景里，这个能力非常关键，因为黑盒函数往往存在多个不同的 x 都能达到同一个 y 值。如果只建模成点预测，就必须强制选择其中某一个峰，一旦选错就会丢掉其他有潜力的区域。MDN 天然保留多峰性，所有峰都被记录下来。

2.3 训练MDN时最关键的数值问题

MDN的训练目标是最小化负对数似然，也就是让已有观测点在模型分布下的平均对数概率尽量大。公式写出来很简单：

Loss = - (1/N) Σ log Σ_k π_k(y_i) · N(x_i; μ_k(y_i), σ_k²(y_i))

但实际训练时有一个躲不掉的数值坑：连加符号在 log 内部，直接取 log 会导致数值溢出。我一开始偷懒直接按公式写，训练到一半 loss 变成 NaN，排查了一圈才发现问题出在这里。正确做法是用 log-sum-exp 技巧，把每个成分的 log 概率算出来，再用torch.logsumexp汇总。

另一个隐患是标准差退化。如果网络优化到后期把某个高斯成分的标准差推到非常小，这个成分就会变成"针尖分布"，对似然贡献极大，网络会疯狂利用这一点来拟合单个点，结果是严重过拟合、对其他区域完全失去泛化能力。后面我会给出具体的规避手段，这里先记住一句话：log_std 必须做上下限约束。

3. 实现一个MDN驱动的黑盒优化器：代码拆解

3.1 带噪声的多峰测试函数

先构造一个肉眼可见的测试函数，方便观察优化器行为。我选了一个带噪声的一维函数，x 在 [-3, 3] 之间，存在多个局部最优峰，并且叠加上高斯噪声模拟真实黑盒环境。

import numpy as np import torch import torch.nn as nn import torch.optim as optim def target_f(x): # x: ndarray of shape (n, 1) return 2.0 * np.sin(5.0 * x) * (1.0 - np.tanh(x ** 2)) def black_box(x): # 黑盒目标函数：真实函数 + 观测噪声 return target_f(x) + np.random.normal(0.0, 0.05, size=x.shape)

选这个函数的原因是多峰特征明显，最优值附近有两个特征不太一样的峰，另外在负半轴还有若干起伏。它不刁钻，但足够测试优化器会不会陷进局部峰。

3.2 条件MDN的PyTorch实现

网络结构不复杂：输入条件变量 y，经过两层带 tanh 激活的全连接层，然后分别输出三个头——混合权重、均值、log标准差。权重用 softmax 归一化成合法的混合系数，log标准差通过 clamp 限制范围防止退化。

class ConditionalMDN(nn.Module): def __init__(self, cond_dim=1, target_dim=1, hidden_dim=32, n_components=8): super().__init__() self.n_components = n_components self.target_dim = target_dim self.net = nn.Sequential( nn.Linear(cond_dim, hidden_dim), nn.Tanh(), nn.Linear(hidden_dim, hidden_dim), nn.Tanh(), ) self.w_layer = nn.Linear(hidden_dim, n_components) self.mean_layer = nn.Linear(hidden_dim, n_components) self.logstd_layer = nn.Linear(hidden_dim, n_components) def forward(self, cond): h = self.net(cond) w = torch.softmax(self.w_layer(h), dim=-1) mean = self.mean_layer(h) log_std = torch.clamp(self.logstd_layer(h), min=-4.0, max=0.5) return w, mean, log_std

这里把 hidden_dim 设置成 32，混合成分数量 K 设置成 8。一维问题用不了太多成分，8 个足够覆盖多峰结构，同时又不会让参数过多、小数据量下训练不稳定。

3.3 负对数似然损失与训练循环

负对数似然是 MDN 的默认损失函数，但实现上要用 logsumexp 来保证数值稳定。这里有个细节：log_std 被 clamp 过后，还需要通过 exp 还原成标准差，再传给torch.distributions.Normal。

def nll_loss(w, mean, log_std, x): std = torch.exp(log_std) + 1e-6 dist = torch.distributions.Normal(mean, std) log_probs = dist.log_prob(x.reshape(-1, 1)) log_comp = torch.log(w + 1e-6) + log_probs return -torch.logsumexp(log_comp, dim=-1).mean() def train_mdn(X, Y, model, epochs=200, lr=1e-3): xt = torch.tensor(X, dtype=torch.float32).reshape(-1, 1) yt = torch.tensor(Y, dtype=torch.float32).reshape(-1, 1) opt = optim.Adam(model.parameters(), lr=lr) best_loss = float("inf") best_state = None patience = 0 for epoch in range(epochs): opt.zero_grad() w, mean, log_std = model(yt) loss = nll_loss(w, mean, log_std, xt) loss.backward() opt.step() if loss.item() < best_loss: best_loss = loss.item() best_state = {k: v.clone() for k, v in model.state_dict().items()} patience = 0 else: patience += 1 if patience >= 30: break if best_state is not None: model.load_state_dict(best_state) return best_loss

训练循环里加了一个简单的早停：连续 30 轮 loss 不下降就停止，避免小数据集上的过拟合。这个技巧在样本量小到只有十几个点时尤为重要。

3.4 采样策略：按目标值条件密度抓候选

模型训练完成后，怎么用它生成下一步的候选点呢？我的策略是：设一个略高于当前最优的目标值 y_star，然后从 p(x | y_star) 里采样一批候选点，再用模型密度过滤一下，选密度最高的候选点去评估。

y_star 的具体取法是当前观测最大值加上一个小扰动。加扰动的目的是避免一直只盯着历史最优，让密度估计有一点往外探索的余量。

def sample_from_cond_mdn(model, cond, n_samples=16): model.eval() with torch.no_grad(): w, mean, log_std = model(cond) w, mean, log_std = w[0], mean[0], log_std[0] std = torch.exp(log_std) comps = torch.multinomial(w, num_samples=n_samples, replacement=True) samples = [] for idx in comps: samples.append(float(torch.normal(mean[idx], std[idx]))) return np.array(samples) def mdn_log_prob(model, cond, x): w, mean, log_std = model(cond) std = torch.exp(log_std) + 1e-6 dist = torch.distributions.Normal(mean, std) log_probs = dist.log_prob(x.reshape(-1, 1)) return torch.logsumexp(torch.log(w + 1e-6) + log_probs, dim=-1)

主优化循环负责把整个流程串起来：初始化一批随机点，训练 MDN，计算 y_star，采样候选，过滤候选，然后评估最优候选并加入数据集，循环往复。

def optimize(max_evals=40, init_n=10, seed=0): rng = np.random.default_rng(seed) X = rng.uniform(-3.0, 3.0, size=(init_n, 1)) Y = black_box(X).ravel() model = ConditionalMDN(cond_dim=1, target_dim=1, hidden_dim=32, n_components=8) history = [] for step in range(max_evals - init_n): train_mdn(X, Y, model) y_star = float(Y.max()) + 0.2 * float(Y.std()) cond = torch.tensor([[y_star]], dtype=torch.float32) candidates = sample_from_cond_mdn(model, cond, n_samples=16) cond_rep = torch.full((len(candidates), 1), y_star, dtype=torch.float32) x_tensor = torch.tensor(candidates.reshape(-1, 1), dtype=torch.float32) probs = mdn_log_prob(model, cond_rep, x_tensor) x_new = candidates[int(probs.argmax().item())] y_new = black_box(np.array([[x_new]]))[0, 0] X = np.vstack([X, [[x_new]]]) Y = np.append(Y, y_new) history.append((float(Y.max()), float(Y.mean()), float(Y.std()))) return X, Y, history, model

注意这里的过滤逻辑：不是采样出 16 个点就把 16 次评估预算全花掉，而是用模型给出的条件概率密度选出最能代表当前探索方向的一个点去评估。这样每轮预算只增加一次评估，整个流程的总评估次数仍然可控。

4. 实测结果：密度采样优化器 vs 随机搜索

4.1 收敛曲线对比

我在同一个测试函数上做了对比实验，总评估预算设为 40 次，初始化 10 个点，随机搜索和 MDN 优化器各重复 5 次。下面的表是 5 次重复的平均结果。

关于这个测试函数，全局最优大概在 x ≈ 0.3 附近，真实值约 1.81。随机搜索能不能摸到全靠运气，5 次里只有 1 次在最后阶段恰好采到 1.7 以上的区域；而 MDN 优化器 5 次里全部超过 1.7，平均第 17 次评估就能进入优质区域。

几个细节值得注意：初始化点本身就有随机性，如果某次初始化恰好落在一个次优峰附近，MDN 优化器会先对周围做一段集中搜索，随后条件密度会自动把目标值抬高，引导采样跳出当前的局部平台。随机搜索则完全不会纠正方向，这是两者收敛差距的主要来源。

4.2 把学到的条件密度画出来，看它如何聚焦

运行过程中我观察了模型学到的条件密度 p(x | y_star) 的形状变化。初始阶段 y_star 不高，密度分布几乎是平的，候选点散布在整个[-3, 3]区间。随着历史最优值上升，y_star 逐步逼近真实最优，条件密度会自然而然地裂成两个主要的峰，分别对应函数右侧的两个高值区域。

到第 25 次左右，条件密度几乎完全集中到全局最优处，另一个局部峰的权重变得很小。这个行为很有实用价值：它不是直接把所有采样都压到某一个点上，而是同时保留若干个高概率区域，通过密度比来自动权衡"开发当前最优区域"和"探索次优区域"之间的比例。

后来我把学到的条件密度画出来，叠加上真实函数，一眼就能看出两者的峰位置高度吻合。这也验证了一个朴素但重要的道理：如果某个 x 能取得较高的 y，那么给定高 y 的条件下，x 的条件密度就应该在那个位置隆起。优化问题里的"最优区域搜索"，本质上就是在追踪条件密度的峰。

4.3 多峰场景下的"保底"能力

拆解采样策略，可以看到这个优化器为什么不容易完全陷进局部最优。

条件 MDN 在给定同一个 y_star 时，可以同时输出多个高斯成分，每个成分可能锁定一个局部高值区域。即使网络把大部分概率质量分配给当前最优的峰，那些次优峰仍然会保留一定的概率，采样时依然有小概率"跳"过去。这个保底能力在高斯过程加置信上界策略里很难实现，因为 GP 的采集函数在高维度下会变得过于乐观或过于保守，而密度模型只是被动地反映已有数据的分布结构。

随机构造候选点时，我没有直接采一个点，而是采 16 个点再用密度过滤，这相当于在多个峰之间做了一个软性的竞争。实际运行中我会同时看到好几个候选点的概率都很高，它们往往分属不同的峰，最终被选中的那个通常是概率最高的峰，但下一轮条件密度更新后，如果另一个峰持续产出好结果，它的权重又会被网络重新抬高。

5. 踩坑记录：MDN优化器的三个隐蔽问题

5.1 log_std 退化导致 loss 变成 NaN

我第一次跑训练循环的时候，日志显示 loss 到第 80 轮变成 NaN。排查过程很典型：先怀疑学习率太高，调低之后依然炸；又怀疑数据集里有异常值，标准化之后还是炸。最后定位到问题出在 log_std 的 clamp 边界上。

我把 clamp 下限定在 -20，理论上允许非常小的标准差。问题是当某个高斯成分的标准差小于训练数据间隔的十分之一时，模型会把这个成分"锁定"在某一个数据点上，似然变成极大的有限值，梯度也变成极大值，接着网络参数飞出去，loss 就 NaN 了。后来我把下限改成 -4，同时给 exp 之后的标准差加了一个 1e-6 的平滑项，问题彻底消失。

经验总结：MDN 训练中出现 NaN，优先检查标准差是否退化，而不是优先调学习率。

5.2 y_star 的取值策略不能太贪

另一个容易踩的坑是 y_star 的取法。如果直接把当前最优值当作不带扰动的 y_star，模型会朝着越来越窄的区域采样，很快把所有点集中到最优附近，导致密度估计丧失探索能力。如果扰动太大，采样范围又回到随机搜索的水平，评估预算被大量浪费在没价值的区域。

我的经验是用当前最优值加上 0.2 倍当前最优值波动幅度，实测下来兼顾了探索和开发。这个系数不需要精确，但必须保持为一个动态值——早期历史最优波动大，扰动就大，天然有更多探索；后期波动小，扰动也小，采样自然收敛。

如果追求更精细的控制，可以把这个系数看成一个随评估进度衰减的退火参数。后面我自己跑真实项目时，用的就是这种衰减式策略，效果会再稳定一点，但代码量也相应多一些。

5.3 更高维场景：MDN不够用的下一步方向

最后说扩展。我在一维测试函数跑通之后，尝试把这个方法应用到八个参数的模拟器上，效果明显下降。原因是 MDN 的混合成分数量在高维空间里远远不够，八个高斯成分在八维空间里只能覆盖极一小块区域。

后来我换了方向：把 MDN 替换成 Normalizing Flow，也就是用流模型来估计条件密度。流模型不依赖固定数量的混合成分，而是通过一系列可逆变换把简单分布变形为复杂分布，表达能力比 MDN 强很多。代价是训练更慢、数值稳定性问题更多。如果你后续打算往高维走，我会建议直接从条件流模型开始，而不是试图给 MDN 塞更多成分。

多说一句实际体会：黑盒优化里真正稀缺的不是更聪明的采集函数，而是先停下想清楚"你要从数据里学出什么结构"。对我而言，密度估计神经网络交出的答案是一张可解释的概率地图，照着地图走，总比闭着眼睛乱逛强。后续有机会我会把高维流模型版本也整理出来，不过这套 MDN 的小而美方案，在低维快速调参场景里依然是性价比最高的选择。