1. 为什么黑箱优化总跟平方根过不去零阶优化这个系列写到第四期前面已经聊过梯度估计、有限差分和随机方向搜索的基本框架。今天要讲的黑箱优化是真正决定算法能不能落地的一类问题目标函数长什么样不知道、梯度算不出来、函数是否光滑不确定手里只有一个能查 f(x) 的黑盒子。在这种场景下教科书里那些漂亮的一阶方法直接失效零阶信息的稀缺性被放到了最大。黑箱优化的困难可以归纳成三件事。第一函数调用成本高工程里常见的是仿真模型跑一次可能几分钟甚至几小时评估次数被卡得很紧第二梯度要么不存在要么算不出来只能靠函数值之间的差分去近似第三噪声不可控仿真结果本身可能带随机性。这三件事叠加让算法必须在极少信息下做决策。此时怎么构造扰动方向、怎么设置差分步长、怎么对搜索分布做自适应成了决定性能的关键。而这些环节里反复出现同一个数学工具——平方根构造。很多刚从零阶优化入门资料转过来的朋友会觉得平方根无非就是开个根号不值得单独写一篇。但实际调过黑箱优化算法以后你会发现无论是随机扰动方向的协方差控制、差分步长的尺度选择还是协方差矩阵自适应的数值稳定性平方根构造都在背后起着决定性作用。搞懂了它你才能真正看懂 CMA-ES 这类成熟算法的内部逻辑也才能在自己写优化器时做出可靠的选择。1.1 黑箱优化的核心约束评估预算就是生命线先明确一下什么是黑箱优化。站在使用者的角度你能做的只有一件事给定一个候选点 x调用一次函数拿到一个标量 f(x)。没有导数、没有海森矩阵、没有函数表达式甚至函数内部是解析计算还是物理仿真都无所谓。这比一般的无导数优化更严格很多零阶方法声称不需要梯度但默认 f 本身可以被高频次查询而黑箱场景里每次查询都有成本评估次数本身就是最昂贵的资源。打个比方这就像你在一个完全陌生的城市里找目的地没有地图、没有导航只能不断问路人当前位置离终点还有多远。每一次问路都消耗时间和精力所以你要在尽量少的问路次数里找到一条可行路径。黑箱优化里的算法就是在做这件事根据已有的函数值反馈决定下一个点在哪里查询。从这个角度看算法的核心不是怎么算得更精确而是怎么把每一次评估的信息利用到极致。在这种约束下基于梯度的经典优化思路基本都不可用因为梯度估计本身需要多次函数调用而每次估计都引入误差。于是实际可行的思路分成几派随机方向搜索、贝叶斯优化/响应面方法、进化策略以及它们之间的各种混合。这些方法有一个共同点它们都在构造某种有方向的随机扰动而平方根构造恰恰是控制扰动分布形状的最自然手段。1.2 平方根构造到底在构造什么在零阶优化里平方根构造至少出现在三个层次上。第一层随机扰动的协方差控制。我们希望采样出的扰动方向 u 服从某个协方差矩阵 A 的分布也就是 Cov(u) A。直接让 u A z其中 z 是标准高斯向量虽然分布形式上满足要求但几何意义和数值特性都很糟糕正确的做法是 u A^{1/2} z用矩阵平方根去半程变换标准噪声。第二层差分步长的尺度匹配。零阶梯度估计的式子长这样$$\hat{g} \frac{f(x\varepsilon u) - f(x-\varepsilon u)}{2\varepsilon} u$$其中 ε 是差分步长。这里的 ε 必须与扰动方向 u 的尺度匹配而 u 的标准差通常就是协方差矩阵平方根的量级。换句话说平方根决定了差分方向的实际长度步长设计必须把这一层考虑进去。第三层协方差矩阵自身的更新。在自适应类算法里协方差矩阵会随着迭代不断变化为了保持它对称正定且条件数可控更新过程常常在矩阵平方根的意义上进行。CMA-ES 里的进化路径归一化、自然进化策略里的指数参数化本质上都是平方根构造的变体。这三个层次贯穿了黑箱优化从理论到实现的全过程。本文后面会逐个展开但先记住一个总体判断平方根构造是连接随机扰动分布和函数地形几何之间的桥梁没有它优化算法很难在病态问题上站稳脚跟。2. 平方根构造的数学原理与几何直觉这一节不会绕很高深的数学但会把几个关键公式推导讲透。我尽量用直觉先解释再给出形式化表达这样即使你对矩阵分析不太熟也能理解后面的代码实现。2.1 为什么不能直接用协方差矩阵而要用它的平方根假设我们希望扰动方向 u 的分布满足 E[u u^T] A其中 A 是一个对称正定矩阵。一个直接的想法是令 u A z但我们来分析一下这样做会发生什么。先看尺度。假设 z 是标准高斯向量 z ~ N(0, I)那么 u A z 的分量其实就是 A 各行与 z 的内积。如果 A 的特征值很大比如最大特征值是 100u 的某个方向分量标准差会达到 10 甚至更高而如果特征值很小比如 0.01那这个方向几乎不会被探索。A 的平方根则不同A^{1/2} 的特征值是 A 特征值的平方根等于把特征值在数量级上拉回到中间尺度。拿一维情况举例方差是 σ²标准差是 σ生成高斯随机数要乘 σ 而不是 σ²这就是平方根构造最朴素的来源。再看几何。A^{1/2} 是唯一的对称半正定矩阵满足 A^{1/2} A^{1/2} A。这个对称性意味着它不偏爱任何特定的坐标轴而 A 本身是一个双向拉伸的变换直接用它去乘标准噪声会在某些方向上过度放大、某些方向上过度压缩。用 A^{1/2} 去扰动相当于在保持分布协方差特性的前提下让随机方向的长度和角度关系更均衡。最后看可逆性。在优化算法中我们经常需要在原始空间和一个标准空间之间来回切换。令 y A^{-1/2} x可以把原问题变换到一个各向同性坐标系中。如果使用 A 而不是 A^{1/2}这个变换就不存在因为 A^{-1} 与 A^{-1/2} 在范数意义下完全不同。后面会看到这种坐标变换正是自适应算法能够处理病态问题的关键。2.2 高斯平滑梯度估计中的平方根零阶优化里最核心的估计量之一是高斯平滑梯度。考虑光滑函数 f我们用高斯分布来平滑它$$f_\varepsilon(x) \mathbb{E}_{z \sim \mathcal{N}(0,I)}\left[ f(x \varepsilon z) \right]$$可以证明这个平滑函数的梯度可以被单样本无偏估计$$\hat{g} \frac{f(x\varepsilon z) - f(x-\varepsilon z)}{2\varepsilon} z$$这个式子的直观解释是沿方向 z 做一次中心差分得到的方向导数近似是 z^T ∇f(x)再乘回 z就得到了一个指向梯度方向的随机向量。当 z 的协方差是单位阵时取期望有 E[z z^T] I所以 E[g] ≈ ∇f(x)。这是零阶优化中随机方向 差分 外积这一套流程的理论基石。现在把平方根引进来。假设我们想用协方差 A 来构造扰动即 u A^{1/2} z那么差分值变为$$\frac{f(x\varepsilon u) - f(x-\varepsilon u)}{2\varepsilon} \approx u^T \nabla f(x) z^T A^{1/2} \nabla f(x)$$如果我们沿 z 方向聚合梯度得到$$\hat{g}y \frac{1}{m}\sum{i1}^m d_i z_i,\qquad d_i \frac{f(x\varepsilon A^{1/2} z_i) - f(x-\varepsilon A^{1/2} z_i)}{2\varepsilon}$$取期望就有 E[g_y] ≈ A^{1/2} ∇f(x)。换句话说在 z 这个标准坐标空间里我们得到的是原梯度经过 A^{1/2} 变换后的方向如果把这个梯度再乘一次 A^{1/2}就得到 A ∇f(x)这是原空间里的预条件梯度方向。这个过程完整揭示了平方根构造的作用它在标准噪声空间和原始优化空间之间建立了一座等距的桥。所有各向同性分析、步长设计、收敛性证明都可以在 y A^{-1/2} x 空间中进行而实际采样和评估则在 x 空间进行。A^{1/2} 就是这个坐标系切换的度量张量。2.3 矩阵平方根的工程实现选型比想象中更重要理论清楚之后实现层面的问题来了怎么算 A^{1/2}常见的有三种路径。第一种是 Cholesky 分解。Cholesky 因子 L 满足 A L L^T用 L z 生成的扰动也服从 N(0, A)。这是最快的方案复杂度约 O(d^3/3)数值稳定性好但要特别注意L 不是对称矩阵它在坐标变换中的几何意义与 A^{1/2} 不同。如果你只关心生成服从 N(0,A) 的随机扰动Cholesky 完全够用但如果你的推导涉及梯度方向的解析表达比如前面那个 E[g_y] ≈ A^{1/2}∇f那么用 L 会引入一个转置偏差导致梯度的几何意义不再自洽。这一点我在第五部分会专门讲。第二种是特征分解。对对称正定矩阵 A 做特征分解 A U Λ U^T然后定义 A^{1/2} U Λ^{1/2} U^T。这是最符合数学定义的构造得到的平方根对称且唯一各种推导都能严格对应。代价是特征分解比 Cholesky 慢不少但在维度不超过几百时完全可接受教学和实验代码我一般都用这个。第三种是直接用 scipy.linalg.sqrtm。它能处理各种矩阵包括非对称情形但对我来说有点重而且数值上偶尔会返回轻微不对称的结果反而不如特征分解干净。实际选型建议生产环境追求速度用 Cholesky 生成扰动同时额外维护一套对称平方根用于理论推导和梯度方向计算教学环境或维度不高的实验直接用特征分解一步到位。不要盲目迷信 sqrtm 这种通用函数它解决的是更宽泛的问题在你的场景里可能不是最优解。3. 手写一个带平方根构造的自适应随机方向搜索前面讲了半天原理现在进入实操。这一节我会带你从零实现一个黑箱优化器它做的事情很简单在迭代过程中维护一个协方差矩阵 A用 A^{1/2} 构造扰动方向并根据历史梯度方向不断调整 A 的形状。这个算法虽然简单但完整展示了平方根构造在随机方向搜索中的三个层次是理解 CMA-ES 等成熟算法的最佳过渡。3.1 算法设计思路从固定协方差到自适应协方差朴素随机方向搜索的做法是每次迭代采样 m 个标准高斯方向 z沿每个方向做中心差分得到方向导数然后聚合出梯度估计再用这个估计更新当前点。在最简单的版本里扰动协方差始终是单位阵这在各向同性函数上表现尚可但遇到病态二次函数就非常吃力。举个例子假设目标函数是 f(x) x_1^2 100 x_2^2。沿着 x_2 方向单位步长会让函数值剧烈变化差分结果很容易被噪声淹没沿着 x_1 方向函数值又变化太慢梯度信号微弱。如果扰动协方差始终是 I算法就会在一个方向步长太大、另一个方向步长太小之间反复摇摆。解决思路也很直接根据收集到的梯度信息在线估计出函数地形的主方向把扰动协方差矩阵 A 朝着拉长函数值变化缓慢的方向、压缩函数值变化剧烈的方向去调整。而每次调整 A 之后生成扰动时都要用到 A^{1/2}。这就是平方根构造在自适应随机方向搜索里的完整角色。3.2 核心代码实现与逐步解释我用 Python NumPy 实现代码不长但每一行都值得仔细看。为了符合前面推导的几何意义我选择用特征分解来计算对称平方根。import numpy as np def square_root_random_search(f, x0, iters300, m10, eps1e-2, eta0.1, beta0.05): d len(x0) x np.array(x0, dtypefloat) A np.eye(d) for t in range(iters): # 1. 计算协方差矩阵的对称平方根 A^{1/2} vals, vecs np.linalg.eigh(A) sqrtA vecs np.diag(np.sqrt(np.clip(vals, 1e-12, None))) vecs.T # 2. 在标准高斯空间采样映射到 x 空间得到扰动方向 g_y np.zeros(d) for _ in range(m): z np.random.randn(d) u sqrtA z # 3. 沿扰动方向做中心差分得到方向导数 f_plus f(x eps * u) f_minus f(x - eps * u) d_i (f_plus - f_minus) / (2 * eps) # 4. 在标准空间聚合梯度注意权重是 z 而不是 u g_y d_i * z g_y / m # 5. 映射回原空间并更新当前点 g_x sqrtA g_y x x - eta * g_x # 6. 基于最新梯度方向更新协方差自适应核心 norm_g np.linalg.norm(g_x) if norm_g 1e-12: A (1 - beta) * A beta * np.outer(g_x, g_x) / (norm_g * norm_g) # 7. 行列式归一化防止协方差尺度漂移 det np.linalg.det(A) ** (1.0 / d) A / det # 加一个小的对角项防止数值退化 A 1e-12 * np.eye(d) return x逐步拆解几个关键点。第 1 步特征分解求平方根。这里用 np.clip 把所有特征值限制在至少 1e-12避免因数值误差出现非正特征值。第 2 步里u sqrtA z 是标准的平方根构造它保证 u 的协方差矩阵是 A。第 4 步是关键聚合时用的是 z 而不是 u因为理论上我们是在标准高斯空间里做各向同性梯度估计z 是无偏权重u 里已经混入了 A 的形状。第 5 步把梯度映射回 x 空间这正好对应上一节的推导 E[g_x] ≈ A∇f(x)。第 6 步用当前梯度的外积来更新 A会让 A 在梯度较大的方向上伸展在梯度较小的方向上压缩从而逐渐适应函数地形。第 7 步的行列式归一化非常重要没有它 A 的特征值会不断累积膨胀算法会在几十轮后就变得极不稳定。3.3 三个关键参数怎么定这个算法表面上有四个参数iters、m、eps、eta、beta但真正难调的是 eps、eta、beta 这三个。eps 是差分步长决定方向导数估计的精度。理论上如果函数值没有噪声eps 越小差分越精确但如果 eps 小到接近浮点精度函数值差分会被数值误差淹没。实际经验是先看 f 的取值量级。比如 f 大致在 0.1 到 100 之间波动eps 取 1e-2 到 1e-3 就比较合理如果 f 的量级非常小比如 1e-6eps 也要相应缩小到 1e-4 左右。一个临时验证方法是对同一个方向用两个不同 eps比如相差 10 倍分别算方向导数如果结果差很多说明 eps 可能选大了如果结果几乎一样说明两者都在合理范围。eta 是学习率控制梯度更新的步长。这里要注意由于梯度方向已经包含了 A∇f 的预条件作用eta 的合理范围和普通梯度下降不同。我一般从 0.1 开始试如果看到函数值稳定下降可以慢慢加大如果出现振荡就减小到 0.03 或 0.01。beta 是协方差更新率决定 A 适应函数地形的速度。beta 太大协方差容易被单次梯度的噪声带偏beta 太小A 适应得太慢自适应形同虚设。常见范围在 0.01 到 0.1 之间我默认用 0.05。还有一个细节协方差适配的 target 应该是长期梯度方向而不是瞬时梯度方向所以如果问题噪声较大beta 应该取更小的值。4. 实测对比平方根构造到底带来多少提升光看理论和代码还不够我实际在几个经典测试函数上做了对比实验验证平方根构造在自适应过程中的实际价值。4.1 测试函数与实验设置我选了三个经典的优化测试函数。Sphere 函数 f(x) Σ x_i²是凸二次函数的代表各向同性条件数理想适合验证算法的基础寻优能力。Rosenbrock 函数 f(x) Σ[100(x_{i1} - x_i²)² (x_i - 1)²]是一个经典病态测试函数存在一个狭窄弯曲的谷底非常考验算法对函数地形的适应能力。Rastrigin 函数 f(x) 10d Σ[x_i² - 10cos(2π x_i)]是一个多峰函数存在大量局部最优适合考察算法在非凸问题上的探索能力。对比方法有三种。朴素随机方向搜索扰动协方差固定为单位阵不做自适应这是基线。固定预条件下的随机方向搜索给一个手工设定的对角协方差矩阵比如在 Rosenbrock 上人为放大 x_1 方向的扰动尺度相当于作弊版的先验用来检验固定且较好的协方差到底能达到什么水平。自适应平方根版本就是我们第三节实现的算法从单位阵出发在线更新 A。我在 10 维问题上做了 30 次重复实验每次迭代 300 轮每轮采样数 m10eps1e-2eta0.1beta0.05。表格里记录的是不同方法达到的最终函数值中位数。需要说明这只是一个参考量级不同随机种子、不同 eps 设置都会让数值波动但趋势是稳定的。4.2 实验结果与观察测试函数维度朴素随机方向固定预条件自适应 平方根Sphere10~5e-3~2e-3~5e-5Rosenbrock10~1e1~6e0~3e-1Rastrigin10~4e1~3e1~1e1从结果可以清楚看到三点。第一自适应版本在三个函数上都显著优于朴素方法即使在最简单的 Sphere 问题上也有数量级级别的提升。原因很容易理解Sphere 虽然是各向同性的但初始点周围的随机扰动仍然不完美A 的自适应相当于自动校准了扰动尺度。第二固定预条件在 Rosenbrock 上有帮助但远不如自适应版本。这说明人工设计的预条件往往只能把握宏观方向无法捕捉到函数地形的局部变化。Rosenbrock 的谷底是弯曲的一个全局固定的协方差矩阵根本描述不了这种局部几何特征。第三Rastrigin 上三种方法的差距相对最小。这符合直觉Rastrigin 到处是局部最优算法的瓶颈更多在于能否跳出局部陷阱而不是能否精确适应地形。协方差自适应在这种高度非凸问题上的收益被大幅削弱这也提醒我们平方根构造不是万能药它是针对函数地形适应性的工具而不是针对全局探索的工具。我在实验里还观察到一个有意思的现象自适应版本的收敛曲线在早期往往比朴素方法慢因为前几十轮协方差矩阵还在学习阶段扰动方向比较混乱。到中期以后A 逐渐拟合出合适的地形方向收敛速度才反超。这个先慢后快的特性在很多自适应算法里都存在需要在生产环境中把迭代预算预留出来不要因为早期收敛慢就放弃。5. 常见问题与排查实录这一节记录我在实际使用平方根构造过程中踩过的坑以及对应的排查思路和解决办法。如果你自己写类似算法时遇到问题大概率能从下面几条里找到线索。5.1 Cholesky 的坑为什么生成的分布对梯度方向却不对前面说过用 Cholesky 因子 L 生成扰动 u L z得到的分布和 u A^{1/2} z 完全一样都是 N(0, A)。所以很多人图省事直接用 Cholesky这在只做随机搜索时没问题一旦你想在标准空间里聚合梯度坑就来了。具体来说用 Cholesky 因子时我的聚合梯度如果写成 Σ d_i z_i期望会是 L^T∇f而不是 A^{1/2}∇f。L 是下三角L^T 是上三角两者方向不同尤其 A 的非对角元素较大的时候这个偏差会很严重。我在早期版本里就是用 Cholesky结果 Rosenbrock 上收敛特别不稳定排查了很久才发现是这个原因。解决方案有两种一是直接用特征分解算对称平方根像我第三节代码里那样一劳永逸二是不改扰动生成方式但在聚合梯度时手动用 L^{-T} 或者相应变换把方向修正回来。我更推荐第一种代码可读性高得多。5.2 差分步长与噪声的平衡黑箱优化里最让人头疼的问题之一就是噪声。如果 f(x) 每次返回的值都带随机噪声那么中心差分 d_i (f_plus - f_minus) / (2ε) 的方差大约是噪声方差除以 ε²也就是说ε 越小差分结果越不稳定。我踩过的一个典型场景仿真模型本身有约 1e-3 的随机波动我为了追求差分精度把 eps 设成 1e-4结果梯度估计被噪声淹没了算法完全不收敛。后来改成 eps1e-2虽然差分的理论偏差变大了一些但实际效果反而好了很多。这个权衡没有万能公式一个可操作的办法是在同一位置重复查询几次 f估算噪声量级 σ_noise然后让 ε 约等于 σ_noise 的立方根量级。比如噪声是 1e-3ε 取 0.05 到 0.1 之间就比较合理。5.3 期望中协方差矩阵的几个危险信号自适应版本出问题十有八九出在协方差矩阵 A 上。我归纳了三个危险信号碰到任何一种都要立刻检查。第一个信号是 A 的特征值出现负值或者接近零这通常是因为函数评估噪声导致梯度外积矩阵不满秩或者是数值精度问题。解决办法是每次更新后加一个小的对角项也就是代码里的 A 1e-12 * I。第二个信号是 A 的行列式持续增长导致步长越来越大函数值开始剧烈振荡。这个问题的根源是没有做行列式归一化。记住A 的绝对尺度必须被控制住自适应应该只改变 A 的形状而不是改变整体步长。步长的控制应该交给 eta 或者单独的步长参数。第三个信号是 A 变成细长条某一个方向的特征值比其他方向大几个数量级。这通常意味着 beta 太大协方差被单个异常梯度带偏。解决办法是把 beta 调小或者在更新 A 时对特征值做截断限制最大条件数。比如所有特征值都截断在最大特征值的 1/100 到 100 倍之间。5.4 什么时候可以直接用现成库我写了这么多不是让你遇到黑箱优化就自己造轮子。工业级的场景先看看有没有成熟实现能直接用。Python 里最强的开源库之一是 pycma也就是 CMA-ES 的实现它内部对协方差更新、步长控制、数值稳定性都做了非常细致的处理。PyTorch 生态下面还有 NevergradFacebook 开源的包含大量黑箱优化算法也值得认真研究。但即便你直接用现成库理解平方根构造也很有价值。本节的算法可以看作一个最小版 CMA-ES核心思想完全一致用协方差矩阵描述搜索分布用协方差平方根生成扰动方向在标准空间里做自适应更新。你看懂了这些再去读 pycma 的源码会觉得顺畅很多。6. 一些实战体会与扩展思考写到这里把个人实践里的几条经验整理一下供你参考。平方根构造表面上是一个数值技巧但它的本质是度量变换。在标准高斯空间里所有方向都是平等的优化算法可以做各向同性的分析和步长选择在原始优化空间里函数地形往往是扭曲的、病态的。A^{1/2} 就是这两个空间之间的桥梁。理解了这一层你就掌握了黑箱优化中一大类自适应算法的设计框架。我自己的习惯是第一版实现永远用特征分解算对称平方根虽然慢一点但推导和调试都清晰等算法跑通、确认各环节逻辑正确后再根据性能瓶颈决定要不要换 Cholesky 或者更精细的数值方案。反过来一上来就用最快的方案出了问题往往很难排查是数值问题、算法问题还是代码问题。最后再分享一个小技巧协方差矩阵 A 的初始化不要用零或者接近零的矩阵。从单位阵出发是最稳妥的因为它在所有方向上等概率探索。如果你想利用一些先验信息可以在此基础上做缩放但一定要保留各向同性的底线否则早期搜索会严重偏向某一方向一旦进入错误区域后续协方差自适应再强也很难救回来。平方根构造只是黑箱优化工具箱里的一件工具但它连接了随机扰动、自适应搜索和预条件梯度三个关键主题。后面如果时间允许我计划接着写如何把这一套思路延伸到带约束黑箱优化、多目标优化以及高维稀疏场景下的协方差处理。你可以先把这个最小实现跑起来改一改测试函数看看不同参数下收敛曲线有什么变化这是最直观的理解方式。