1. 从“瞎猜”到“聪明猜”为什么我们需要交叉熵优化算法在解决一个复杂的工程优化问题时比如你要设计一个新型无人机的机翼形状有十几个参数翼展、后掠角、厚度分布等需要调整目标是让它在特定飞行条件下的升阻比最高。你可能会先凭经验或直觉给出一组初始参数然后尝试微调。但很快你会发现参数空间太大了像在茫茫大海里捞针传统的梯度下降法在这里可能因为函数不光滑、不可导而失效或者陷入局部最优解。这时候一种“聪明地瞎猜”的方法就显得尤为重要——这就是交叉熵优化算法。交叉熵优化算法简称CEM听起来名字有点学术但它的核心思想非常直观用概率分布来描述我们对最优解的“猜测”然后通过不断迭代让这个概率分布越来越集中在真正的最优解附近。它属于进化策略和蒙特卡洛方法的一个分支特别擅长处理那些黑箱、非线性、多峰、甚至带有噪声的优化问题。你不需要知道目标函数的梯度只需要能对任何一组输入参数计算出它的“好坏”即目标函数值就行。我第一次接触CEM是在解决一个供应链网络优化的问题上有几十个仓库的选址和库存分配策略需要决策目标是最小化总物流成本。问题维度高约束复杂用常规的线性规划加启发式调整效果总是不理想。尝试了CEM之后虽然前期迭代看起来像在随机采样但很快算法就能“嗅”到优质解区域的大致方向并集中火力在那里探索最终找到了比人工方案成本低15%的配置。这种从“漫无目的”到“精准聚焦”的转变正是CEM的魅力所在。2. CEM的核心原理一场“优胜劣汰”的概率游戏要理解CEM我们可以把它想象成一场不断进化的“选拔赛”。我们不是直接寻找最优解那个“点”而是维护一个描述“哪些区域可能出好解”的概率模型通常是多元高斯分布。然后通过“采样-评估-更新模型”的循环让这个概率模型不断向高分区域收缩。2.1 算法流程拆解一个标准的CEM用于连续函数优化的流程可以分解为以下几步初始化概率模型我们假设解服从一个多元正态分布。一开始我们对最优解在哪里一无所知所以会用一个比较“宽泛”的分布来覆盖整个搜索空间。比如对于每个维度变量我们设定一个均值μ和较大的方差σ²。均值通常可以设为搜索范围的中心点方差则设得足够大以确保初始采样能覆盖整个空间。采样从当前的概率分布N(μ, Σ)其中Σ是协方差矩阵中随机抽取N个候选解。这就像第一轮海选从我们当前认为“可能不错”的区域内随机找一批选手出来。评估用目标函数f(x)计算每一个候选解的性能得分。对于最小化问题得分越低越好对于最大化问题得分越高越好。精英选择从这N个候选解中选出得分最好的前K个K N这K个解被称为“精英样本”。这步是关键我们只相信表现最好的那一小部分样本所透露的信息。K通常取N的10%-20%比如从100个样本中选10-20个最好的。更新概率模型用这K个精英样本来重新估计概率分布的参数均值μ和协方差矩阵Σ。新的μ就是这K个精英样本的均值新的Σ就是它们的协方差。这一步是算法的灵魂我们通过精英样本的分布来更新我们对“好解在哪里”的信念。新的分布会更集中在这K个精英样本的周围。迭代用更新后的概率分布均值μ_new协方差Σ_new替换旧的回到第2步开始新一轮的采样。如此循环直到分布变得非常集中方差很小或达到最大迭代次数。这个过程就像一个不断聚焦的镜头初始分布很模糊覆盖全图每一轮迭代我们都根据当前找到的最清晰的部分精英样本重新调整镜头的焦点和范围让它越来越对准图像中最有价值的细节最优解。2.2 交叉熵的名字从何而来你可能会问这算法为什么叫“交叉熵”这和信息论里的概念有关。在概率分布估计中我们希望找到一个概率分布q使得用它来采样得到精英样本的概率最大。这等价于最小化真实精英样本经验分布与我们所假设分布q之间的交叉熵。在算法步骤中用精英样本的均值和协方差来更新分布正是这种最小化交叉熵思想的一种高效实现。不过对于应用者来说你完全可以暂时不去深究这个数学细节而把它理解为“用表现好的样本来修正我们的猜测模型”。3. 手把手实现用Python为多元函数寻优构建CEM引擎理论说再多不如一行代码。下面我将以一个具体的多元函数——Rastrigin函数——为例展示如何从零实现一个CEM优化器。Rastrigin函数是优化算法领域著名的“考题”它在搜索空间内存在大量的局部极小值点全局最小值位于原点(0,0,...)非常适合测试算法的全局搜索和逃离局部最优的能力。其n维形式的公式为f(x) 10*n Σ_{i1}^{n} [ x_i^2 - 10*cos(2πx_i) ]搜索范围通常设在[-5.12, 5.12]^n。我们的目标是找到最小化该函数的x。3.1 环境准备与核心函数定义首先确保你有基本的科学计算环境。我们将使用numpy进行向量和矩阵运算。import numpy as np import matplotlib.pyplot as plt # 定义Rastrigin函数 (最小化问题) def rastrigin(x): 计算Rastrigin函数值。 参数: x: 一个一维numpy数组代表一个候选解。 返回: 函数值 (float). n len(x) return 10 * n np.sum(x**2 - 10 * np.cos(2 * np.pi * x)) # 为了可视化我们定义一个2维的版本方便画图 def rastrigin_2d(x1, x2): return 20 (x1**2 - 10*np.cos(2*np.pi*x1)) (x2**2 - 10*np.cos(2*np.pi*x2))3.2 CEM优化器类实现我们将CEM封装成一个类这样参数管理和迭代过程会更清晰。class CrossEntropyOptimizer: 用于连续空间优化的交叉熵优化算法实现。 def __init__(self, objective_func, dim, bounds, pop_size100, elite_frac0.2, max_iter100, epsilon1e-3, seedNone): 初始化CEM优化器。 参数: objective_func: 目标函数接受一个一维数组输入返回一个标量最小化。 dim: 问题的维度变量个数。 bounds: 每个变量的上下界list of tuples例如 [(lb1, ub1), (lb2, ub2), ...]。 pop_size: 每代采样的人口大小 (N)。 elite_frac: 精英样本比例用于选择前K个最好解 (K int(pop_size * elite_frac))。 max_iter: 最大迭代次数。 epsilon: 收敛阈值当分布的标准差最大值小于此值时停止。 seed: 随机种子用于复现结果。 self.f objective_func self.dim dim self.bounds np.array(bounds) self.pop_size pop_size self.elite_size max(1, int(pop_size * elite_frac)) # 确保至少有一个精英 self.max_iter max_iter self.epsilon epsilon self.rng np.random.default_rng(seed) # 初始化均值和方差均值设在搜索空间中心方差设为范围的四分之一保证初始覆盖。 self.mu np.mean(self.bounds, axis1) self.sigma (self.bounds[:, 1] - self.bounds[:, 0]) / 4.0 # 初始标准差 # 为了简单我们假设各维度独立协方差矩阵为对角矩阵 diag(sigma^2) self.cov np.diag(self.sigma ** 2) # 记录历史 self.history {mu: [], sigma: [], best_score: [], best_solution: []} def sample_population(self): 从当前多元正态分布 N(mu, cov) 中采样 pop_size 个候选解。 # numpy的multivariate_normal用于从多元正态分布采样 samples self.rng.multivariate_normal(self.mu, self.cov, sizeself.pop_size) # 简单处理边界约束将越界的值裁剪到边界上这是最简单的方法更高级的可以反射或惩罚 samples np.clip(samples, self.bounds[:, 0], self.bounds[:, 1]) return samples def update_distribution(self, elite_samples): 用精英样本更新分布参数 mu 和 cov。 new_mu np.mean(elite_samples, axis0) # 计算样本协方差矩阵ddof0表示除以N是总体协方差的无偏估计在CEM中的常用方式。 # 加入一个很小的正则项防止协方差矩阵退化变成奇异矩阵。 new_cov np.cov(elite_samples, rowvarFalse) 1e-6 * np.eye(self.dim) return new_mu, new_cov def optimize(self): 执行优化主循环。 best_global_score float(inf) best_global_solution None for iteration in range(self.max_iter): # 1. 采样 population self.sample_population() # 2. 评估 scores np.array([self.f(ind) for ind in population]) # 3. 精英选择 elite_indices np.argsort(scores)[:self.elite_size] elite_samples population[elite_indices] elite_scores scores[elite_indices] # 更新全局最优 current_best_idx np.argmin(scores) if scores[current_best_idx] best_global_score: best_global_score scores[current_best_idx] best_global_solution population[current_best_idx].copy() # 4. 更新分布 self.mu, self.cov self.update_distribution(elite_samples) self.sigma np.sqrt(np.diag(self.cov)) # 记录标准差用于收敛判断 # 记录历史 self.history[mu].append(self.mu.copy()) self.history[sigma].append(self.sigma.copy()) self.history[best_score].append(elite_scores[0]) # 记录当代精英中最好的分数 self.history[best_solution].append(best_global_solution.copy()) # 打印进度 if iteration % 10 0: print(fIter {iteration:3d}, Best Score: {best_global_score:.6f}, Mean Sigma: {np.mean(self.sigma):.6f}) # 5. 收敛检查如果所有维度的标准差都足够小则认为收敛 if np.max(self.sigma) self.epsilon: print(fConverged at iteration {iteration}.) break print(f\nOptimization finished.) print(fBest solution found: {best_global_solution}) print(fBest score: {best_global_score}) return best_global_solution, best_global_score, self.history3.3 运行优化与结果分析现在让我们在2维Rastrigin函数上测试我们的CEM优化器。# 定义问题 dim 2 bounds [(-5.12, 5.12)] * dim # 每个变量的边界相同 # 创建优化器实例 cem CrossEntropyOptimizer(objective_funcrastrigin, dimdim, boundsbounds, pop_size50, # 每代采样50个点 elite_frac0.2, # 选择前20%作为精英 (10个) max_iter80, epsilon1e-4, seed42) # 执行优化 best_sol, best_score, history cem.optimize() # 可视化优化过程 # 绘制目标函数曲面2D等高线形式 x np.linspace(-5.12, 5.12, 200) y np.linspace(-5.12, 5.12, 200) X, Y np.meshgrid(x, y) Z rastrigin_2d(X, Y) plt.figure(figsize(15, 5)) # 子图1函数等高线及均值演化路径 plt.subplot(1, 3, 1) plt.contourf(X, Y, Z, levels50, cmapviridis) plt.colorbar(labelRastrigin Function Value) mu_history np.array(history[mu]) plt.plot(mu_history[:, 0], mu_history[:, 1], r-o, markersize4, linewidth1.5, labelMean Evolution) plt.scatter(best_sol[0], best_sol[1], cwhite, edgecolorsred, s200, marker*, labelBest Found) plt.xlabel(x1) plt.ylabel(x2) plt.title(CEM Search Path on Rastrigin Function) plt.legend() plt.grid(True, alpha0.3) # 子图2最佳得分随迭代的变化 plt.subplot(1, 3, 2) plt.plot(history[best_score], b-, linewidth2) plt.xlabel(Iteration) plt.ylabel(Best Score (Elite)) plt.title(Convergence of Best Score) plt.grid(True, alpha0.3) # 子图3标准差随迭代的变化反映分布的收缩 plt.subplot(1, 3, 3) sigma_history np.array(history[sigma]) for d in range(dim): plt.plot(sigma_history[:, d], labelfDim {d1} sigma) plt.xlabel(Iteration) plt.ylabel(Standard Deviation (Sigma)) plt.title(Shrinking of Search Distribution) plt.legend() plt.grid(True, alpha0.3) plt.tight_layout() plt.show()运行这段代码你会看到CEM如何工作。在等高线图上红色的点线显示了概率分布均值μ的移动轨迹。一开始它可能在搜索空间里跳跃但随着迭代它会稳定地朝着全局最优点(0,0)附近移动。同时“最佳得分”曲线会快速下降并趋于平稳“标准差”曲线会逐渐收缩到接近零这表明算法对解的位置越来越确定。4. 关键参数调优与实战避坑指南实现一个能跑的CEM只是第一步让它跑得好、跑得稳才是体现功力的地方。下面这些参数和技巧是我在多个项目实践中总结出来的很多是官方论文里不会细说的“坑”。4.1 核心参数如何设置不是拍脑袋决定的种群大小pop_size与精英比例elite_fracpop_size这是最重要的参数之一。太小采样不足容易陷入局部最优太大计算开销剧增。一个经验法则是pop_size至少是问题维度dim的10倍。对于10维的问题我通常会从100或200开始尝试。对于像Rastrigin这样多峰的函数需要更大的种群来维持多样性。elite_frac通常在0.1到0.3之间。比例越小选择压力越大收敛越快但也更容易早熟陷入局部最优。比例越大保留了更多样性收敛速度慢但探索能力更强。我个人的经验是对于复杂多峰问题可以从0.2开始如果发现早熟可以尝试增大到0.3对于相对平滑的问题可以降低到0.1以加速收敛。初始方差sigma代码中我们用了搜索范围的1/4。这是一个不错的起点。原则是初始分布要能覆盖你感兴趣或你认为最可能包含最优解的整个区域。如果你对最优解的位置一无所知就设大一些。如果你有先验知识比如知道大概在某个区间可以把均值设在那附近方差设小一些能加速收敛。协方差矩阵的处理在我们的简单实现中我们假设各维度独立协方差矩阵是对角阵。这在很多情况下是有效的尤其是当问题变量之间的耦合不强时。但如果变量间存在强相关性比如调整机翼前缘和后缘的参数是联动的使用完整的协方差矩阵能让算法更高效地捕捉这种关系沿着“好解”的等高线方向收缩。我们的实现已经计算了完整协方差这比只使用对角矩阵更强大。收敛阈值epsilon不要设得太小。通常1e-3到1e-5是一个合理范围。有时算法会因为方差收缩得太快而“卡住”此时即使没达到阈值最佳解也不再改善。更稳健的停止条件是同时监控最佳解在连续多次迭代中的改善程度。例如可以增加一个判断如果最近20代的最佳解改善幅度小于某个微小值如1e-6也停止迭代。4.2 实际项目中踩过的“坑”与解决方案坑一方差坍塌与早熟收敛这是CEM最常见的问题。算法过早地聚焦在一个局部最优区域方差迅速变得极小种群失去多样性再也跳不出来了。解决方案方差平滑更新不要直接用精英样本的方差替换旧方差而是采用平滑更新sigma_new alpha * sigma_elite (1-alpha) * sigma_old。其中alpha是一个学习率如0.7。这保留了部分历史信息防止突变。注入噪声在更新后的方差上加上一个小的常数项或按比例缩放确保方差不会低于一个最小值如sigma max(sigma, sigma_min)。这相当于始终保持一点探索的“火种”。重启策略当检测到收敛方差很小但解的质量不满意时保存当前最优解然后重新初始化分布比如把均值设在当前最优解附近但把方差重新放大开始新一轮优化。这相当于“局部搜索”“全局探索”的循环。坑二处理边界约束我们的简单实现用了np.clip进行裁剪。但这会导致一个副作用大量样本被挤压在边界上扭曲了真实的概率分布特别是当最优解就在边界附近时算法性能会下降。更好的解决方案反射边界如果采样点越界将它“反射”回边界内。例如对于变量x在[a, b]内如果采样值x_sampled b则令x 2*b - x_sampled如果x_sampled a则令x 2*a - x_sampled。这能保持样本在边界附近的分布特性。惩罚函数法不修改采样点但在评估目标函数时对越界的点施加一个很大的惩罚值比如加上一个与越界距离成正比的巨大正数。这样越界的点会在精英选择中被自然淘汰分布更新时会自动向可行域内移动。坑三高维灾难当问题维度dim很高时比如上百维CEM会面临挑战。首先采样那么多点pop_size需要更大计算开销大。其次估计一个高维的完整协方差矩阵需要大量的精英样本至少dim1个否则矩阵会是病态的。解决方案使用对角协方差或因子化模型如果假设维度独立可行就只用对角协方差参数从dim*dim降到dim个。或者使用更复杂的因子化分布如Cholesky因子来减少参数。增量更新与正则化使用秩1更新或秩μ更新来迭代更新协方差矩阵而不是每代重新计算。同时必须加入正则化就像我们代码里加的1e-6 * I防止矩阵奇异。5. 超越基础让CEM更强大的进阶技巧掌握了基础CEM后你可以通过以下技巧让它适应更复杂的场景。5.1 混合策略CEM作为局部搜索器CEM的全局探索能力在初期较强但后期精细搜索能力可能不如一些基于梯度的局部方法。一个常见的策略是将CEM与其他算法结合。例如CEM 梯度下降先用CEM进行全局粗搜索找到有希望的盆地区域然后将CEM找到的最优解作为起点用梯度下降法或L-BFGS等进行精细的局部优化。这结合了全局和局部搜索的优点。CEM in CMA-ES Style借鉴CMA-ES协方差矩阵自适应进化策略中的一些思想如进化路径。进化路径记录了均值更新的连续方向利用它来更新协方差矩阵可以更聪明地学习目标函数的拓扑结构加快收敛。这相当于给CEM增加了“记忆”和“动量”。5.2 处理噪声与鲁棒优化在实际工程中目标函数评估可能有噪声比如仿真有随机性或者实验测量有误差。噪声会干扰精英样本的选择可能导致算法不稳定。重采样平滑对同一个候选解进行多次评估取其平均值作为该解的得分。这增加了评估的可靠性但代价是计算量倍增。分位数选择替代Top-K不直接取前K个而是取分数低于某个分位数如75%分位数的所有样本作为“精英集”。这在一定程度上容忍了噪声因为偶尔评估差的“好解”和评估好的“差解”可能都会被纳入但在统计意义上更新方向仍是正确的。5.3 离散与混合优化问题标准的CEM针对连续变量。但很多问题包含离散变量如选择哪台机器、是否开启某个功能。离散化处理对于有序的离散变量如整数1,2,3...可以将其视为连续变量进行CEM优化采样后四舍五入到最近的离散值。更新分布时仍然使用四舍五入前的连续值来计算均值和方差。联合分布对于同时包含连续变量和分类变量的问题可以维护多个概率模型。例如连续部分用高斯分布分类部分用多项分布每个分类有一个概率。采样时分别从各自的分布中采样然后组合成一个完整解。更新时连续部分用高斯分布更新分类部分用精英样本中各类别的出现频率来更新多项分布的概率。6. 多元函数寻优实战一个简单的工程设计案例让我们设想一个更贴近实际的简单案例优化一个箱型梁的截面尺寸在满足强度、刚度要求下最小化重量。假设梁的截面由宽度b、高度h、腹板厚度t_w、翼缘厚度t_f四个变量决定均为连续正数。约束包括最大应力不超过许用应力最大挠度不超过允许值。目标函数是截面面积正比于重量。这是一个带约束的4维优化问题。我们可以用罚函数法将其转化为无约束问题将约束违反量乘以一个大的惩罚系数加到目标函数面积上。def beam_objective(x): 箱型梁优化目标函数罚函数法处理约束。 x [b, h, t_w, t_f] b, h, t_w, t_f x # 计算截面面积 (目标越小越好) area 2*(b*t_f) 2*( (h-2*t_f)*t_w ) # 简化模型 # 计算应力和挠度 (简化公式仅用于示例) # 假设受均布载荷q跨度L材料弹性模量E许用应力sigma_allow许用挠度delta_allow L, q, E, sigma_allow, delta_allow 5.0, 10000, 2.1e11, 200e6, L/500 # 示例参数 I (b*h**3 - (b-2*t_w)*(h-2*t_f)**3)/12.0 # 惯性矩 stress (q * L**2 * h) / (8 * I) # 最大弯曲应力简化 deflection (5 * q * L**4) / (384 * E * I) # 最大挠度简化 # 罚函数 penalty 0.0 if stress sigma_allow: penalty 1e6 * (stress - sigma_allow)**2 # 应力违反惩罚 if deflection delta_allow: penalty 1e6 * (deflection - delta_allow)**2 # 挠度违反惩罚 # 几何约束厚度必须小于宽度/高度的一半等 if t_f h/2 or t_w b/2: penalty 1e9 return area penalty # 定义变量边界 (单位米) bounds [(0.1, 0.5), # b (0.2, 1.0), # h (0.005, 0.02), # t_w (0.008, 0.03)] # t_f # 使用CEM优化 cem_beam CrossEntropyOptimizer(objective_funcbeam_objective, dim4, boundsbounds, pop_size80, elite_frac0.15, max_iter150, epsilon1e-5, seed123) best_design, best_weight, hist cem_beam.optimize() print(\n--- 箱型梁优化结果 ---) print(f最优截面尺寸: 宽度 b {best_design[0]:.4f} m, 高度 h {best_design[1]:.4f} m) print(f 腹板厚 t_w {best_design[2]:.4f} m, 翼缘厚 t_f {best_design[3]:.4f} m) print(f估算的最小截面面积 (重量指标): {best_weight:.6f} m²) # 可以重新计算一下应力和挠度验证是否满足约束通过这个例子你可以看到CEM如何在一个有实际物理意义、带约束的多元函数寻优问题上工作。它不需要你知道应力、挠度公式的梯度只需要能计算出给定尺寸下的“性能”面积惩罚即可。这种黑箱优化的特性使得CEM在仿真驱动设计、控制器参数整定等领域非常有用。最后我个人在多次使用CEM后的体会是把它看作一个“智能的随机搜索框架”。它的优势在于概念简单实现容易并行友好采样和评估可以完全并行并且对目标函数形式要求极低。它的性能很大程度上依赖于参数设置和那些“小技巧”如方差平滑、边界处理。不要期望它在任何问题上都碾压其他算法但在处理中低维度、计算代价昂贵、黑箱或无梯度的复杂优化问题时CEM绝对是一个值得你首先尝试的、强大而实用的工具。当你没有更好选择的时候试试CEM它常常能给你一个不错的起点甚至是一个惊喜。 交叉熵优化算法：从原理到Python实现，解决多元函数寻优难题
1. 从“瞎猜”到“聪明猜”为什么我们需要交叉熵优化算法在解决一个复杂的工程优化问题时比如你要设计一个新型无人机的机翼形状有十几个参数翼展、后掠角、厚度分布等需要调整目标是让它在特定飞行条件下的升阻比最高。你可能会先凭经验或直觉给出一组初始参数然后尝试微调。但很快你会发现参数空间太大了像在茫茫大海里捞针传统的梯度下降法在这里可能因为函数不光滑、不可导而失效或者陷入局部最优解。这时候一种“聪明地瞎猜”的方法就显得尤为重要——这就是交叉熵优化算法。交叉熵优化算法简称CEM听起来名字有点学术但它的核心思想非常直观用概率分布来描述我们对最优解的“猜测”然后通过不断迭代让这个概率分布越来越集中在真正的最优解附近。它属于进化策略和蒙特卡洛方法的一个分支特别擅长处理那些黑箱、非线性、多峰、甚至带有噪声的优化问题。你不需要知道目标函数的梯度只需要能对任何一组输入参数计算出它的“好坏”即目标函数值就行。我第一次接触CEM是在解决一个供应链网络优化的问题上有几十个仓库的选址和库存分配策略需要决策目标是最小化总物流成本。问题维度高约束复杂用常规的线性规划加启发式调整效果总是不理想。尝试了CEM之后虽然前期迭代看起来像在随机采样但很快算法就能“嗅”到优质解区域的大致方向并集中火力在那里探索最终找到了比人工方案成本低15%的配置。这种从“漫无目的”到“精准聚焦”的转变正是CEM的魅力所在。2. CEM的核心原理一场“优胜劣汰”的概率游戏要理解CEM我们可以把它想象成一场不断进化的“选拔赛”。我们不是直接寻找最优解那个“点”而是维护一个描述“哪些区域可能出好解”的概率模型通常是多元高斯分布。然后通过“采样-评估-更新模型”的循环让这个概率模型不断向高分区域收缩。2.1 算法流程拆解一个标准的CEM用于连续函数优化的流程可以分解为以下几步初始化概率模型我们假设解服从一个多元正态分布。一开始我们对最优解在哪里一无所知所以会用一个比较“宽泛”的分布来覆盖整个搜索空间。比如对于每个维度变量我们设定一个均值μ和较大的方差σ²。均值通常可以设为搜索范围的中心点方差则设得足够大以确保初始采样能覆盖整个空间。采样从当前的概率分布N(μ, Σ)其中Σ是协方差矩阵中随机抽取N个候选解。这就像第一轮海选从我们当前认为“可能不错”的区域内随机找一批选手出来。评估用目标函数f(x)计算每一个候选解的性能得分。对于最小化问题得分越低越好对于最大化问题得分越高越好。精英选择从这N个候选解中选出得分最好的前K个K N这K个解被称为“精英样本”。这步是关键我们只相信表现最好的那一小部分样本所透露的信息。K通常取N的10%-20%比如从100个样本中选10-20个最好的。更新概率模型用这K个精英样本来重新估计概率分布的参数均值μ和协方差矩阵Σ。新的μ就是这K个精英样本的均值新的Σ就是它们的协方差。这一步是算法的灵魂我们通过精英样本的分布来更新我们对“好解在哪里”的信念。新的分布会更集中在这K个精英样本的周围。迭代用更新后的概率分布均值μ_new协方差Σ_new替换旧的回到第2步开始新一轮的采样。如此循环直到分布变得非常集中方差很小或达到最大迭代次数。这个过程就像一个不断聚焦的镜头初始分布很模糊覆盖全图每一轮迭代我们都根据当前找到的最清晰的部分精英样本重新调整镜头的焦点和范围让它越来越对准图像中最有价值的细节最优解。2.2 交叉熵的名字从何而来你可能会问这算法为什么叫“交叉熵”这和信息论里的概念有关。在概率分布估计中我们希望找到一个概率分布q使得用它来采样得到精英样本的概率最大。这等价于最小化真实精英样本经验分布与我们所假设分布q之间的交叉熵。在算法步骤中用精英样本的均值和协方差来更新分布正是这种最小化交叉熵思想的一种高效实现。不过对于应用者来说你完全可以暂时不去深究这个数学细节而把它理解为“用表现好的样本来修正我们的猜测模型”。3. 手把手实现用Python为多元函数寻优构建CEM引擎理论说再多不如一行代码。下面我将以一个具体的多元函数——Rastrigin函数——为例展示如何从零实现一个CEM优化器。Rastrigin函数是优化算法领域著名的“考题”它在搜索空间内存在大量的局部极小值点全局最小值位于原点(0,0,...)非常适合测试算法的全局搜索和逃离局部最优的能力。其n维形式的公式为f(x) 10*n Σ_{i1}^{n} [ x_i^2 - 10*cos(2πx_i) ]搜索范围通常设在[-5.12, 5.12]^n。我们的目标是找到最小化该函数的x。3.1 环境准备与核心函数定义首先确保你有基本的科学计算环境。我们将使用numpy进行向量和矩阵运算。import numpy as np import matplotlib.pyplot as plt # 定义Rastrigin函数 (最小化问题) def rastrigin(x): 计算Rastrigin函数值。 参数: x: 一个一维numpy数组代表一个候选解。 返回: 函数值 (float). n len(x) return 10 * n np.sum(x**2 - 10 * np.cos(2 * np.pi * x)) # 为了可视化我们定义一个2维的版本方便画图 def rastrigin_2d(x1, x2): return 20 (x1**2 - 10*np.cos(2*np.pi*x1)) (x2**2 - 10*np.cos(2*np.pi*x2))3.2 CEM优化器类实现我们将CEM封装成一个类这样参数管理和迭代过程会更清晰。class CrossEntropyOptimizer: 用于连续空间优化的交叉熵优化算法实现。 def __init__(self, objective_func, dim, bounds, pop_size100, elite_frac0.2, max_iter100, epsilon1e-3, seedNone): 初始化CEM优化器。 参数: objective_func: 目标函数接受一个一维数组输入返回一个标量最小化。 dim: 问题的维度变量个数。 bounds: 每个变量的上下界list of tuples例如 [(lb1, ub1), (lb2, ub2), ...]。 pop_size: 每代采样的人口大小 (N)。 elite_frac: 精英样本比例用于选择前K个最好解 (K int(pop_size * elite_frac))。 max_iter: 最大迭代次数。 epsilon: 收敛阈值当分布的标准差最大值小于此值时停止。 seed: 随机种子用于复现结果。 self.f objective_func self.dim dim self.bounds np.array(bounds) self.pop_size pop_size self.elite_size max(1, int(pop_size * elite_frac)) # 确保至少有一个精英 self.max_iter max_iter self.epsilon epsilon self.rng np.random.default_rng(seed) # 初始化均值和方差均值设在搜索空间中心方差设为范围的四分之一保证初始覆盖。 self.mu np.mean(self.bounds, axis1) self.sigma (self.bounds[:, 1] - self.bounds[:, 0]) / 4.0 # 初始标准差 # 为了简单我们假设各维度独立协方差矩阵为对角矩阵 diag(sigma^2) self.cov np.diag(self.sigma ** 2) # 记录历史 self.history {mu: [], sigma: [], best_score: [], best_solution: []} def sample_population(self): 从当前多元正态分布 N(mu, cov) 中采样 pop_size 个候选解。 # numpy的multivariate_normal用于从多元正态分布采样 samples self.rng.multivariate_normal(self.mu, self.cov, sizeself.pop_size) # 简单处理边界约束将越界的值裁剪到边界上这是最简单的方法更高级的可以反射或惩罚 samples np.clip(samples, self.bounds[:, 0], self.bounds[:, 1]) return samples def update_distribution(self, elite_samples): 用精英样本更新分布参数 mu 和 cov。 new_mu np.mean(elite_samples, axis0) # 计算样本协方差矩阵ddof0表示除以N是总体协方差的无偏估计在CEM中的常用方式。 # 加入一个很小的正则项防止协方差矩阵退化变成奇异矩阵。 new_cov np.cov(elite_samples, rowvarFalse) 1e-6 * np.eye(self.dim) return new_mu, new_cov def optimize(self): 执行优化主循环。 best_global_score float(inf) best_global_solution None for iteration in range(self.max_iter): # 1. 采样 population self.sample_population() # 2. 评估 scores np.array([self.f(ind) for ind in population]) # 3. 精英选择 elite_indices np.argsort(scores)[:self.elite_size] elite_samples population[elite_indices] elite_scores scores[elite_indices] # 更新全局最优 current_best_idx np.argmin(scores) if scores[current_best_idx] best_global_score: best_global_score scores[current_best_idx] best_global_solution population[current_best_idx].copy() # 4. 更新分布 self.mu, self.cov self.update_distribution(elite_samples) self.sigma np.sqrt(np.diag(self.cov)) # 记录标准差用于收敛判断 # 记录历史 self.history[mu].append(self.mu.copy()) self.history[sigma].append(self.sigma.copy()) self.history[best_score].append(elite_scores[0]) # 记录当代精英中最好的分数 self.history[best_solution].append(best_global_solution.copy()) # 打印进度 if iteration % 10 0: print(fIter {iteration:3d}, Best Score: {best_global_score:.6f}, Mean Sigma: {np.mean(self.sigma):.6f}) # 5. 收敛检查如果所有维度的标准差都足够小则认为收敛 if np.max(self.sigma) self.epsilon: print(fConverged at iteration {iteration}.) break print(f\nOptimization finished.) print(fBest solution found: {best_global_solution}) print(fBest score: {best_global_score}) return best_global_solution, best_global_score, self.history3.3 运行优化与结果分析现在让我们在2维Rastrigin函数上测试我们的CEM优化器。# 定义问题 dim 2 bounds [(-5.12, 5.12)] * dim # 每个变量的边界相同 # 创建优化器实例 cem CrossEntropyOptimizer(objective_funcrastrigin, dimdim, boundsbounds, pop_size50, # 每代采样50个点 elite_frac0.2, # 选择前20%作为精英 (10个) max_iter80, epsilon1e-4, seed42) # 执行优化 best_sol, best_score, history cem.optimize() # 可视化优化过程 # 绘制目标函数曲面2D等高线形式 x np.linspace(-5.12, 5.12, 200) y np.linspace(-5.12, 5.12, 200) X, Y np.meshgrid(x, y) Z rastrigin_2d(X, Y) plt.figure(figsize(15, 5)) # 子图1函数等高线及均值演化路径 plt.subplot(1, 3, 1) plt.contourf(X, Y, Z, levels50, cmapviridis) plt.colorbar(labelRastrigin Function Value) mu_history np.array(history[mu]) plt.plot(mu_history[:, 0], mu_history[:, 1], r-o, markersize4, linewidth1.5, labelMean Evolution) plt.scatter(best_sol[0], best_sol[1], cwhite, edgecolorsred, s200, marker*, labelBest Found) plt.xlabel(x1) plt.ylabel(x2) plt.title(CEM Search Path on Rastrigin Function) plt.legend() plt.grid(True, alpha0.3) # 子图2最佳得分随迭代的变化 plt.subplot(1, 3, 2) plt.plot(history[best_score], b-, linewidth2) plt.xlabel(Iteration) plt.ylabel(Best Score (Elite)) plt.title(Convergence of Best Score) plt.grid(True, alpha0.3) # 子图3标准差随迭代的变化反映分布的收缩 plt.subplot(1, 3, 3) sigma_history np.array(history[sigma]) for d in range(dim): plt.plot(sigma_history[:, d], labelfDim {d1} sigma) plt.xlabel(Iteration) plt.ylabel(Standard Deviation (Sigma)) plt.title(Shrinking of Search Distribution) plt.legend() plt.grid(True, alpha0.3) plt.tight_layout() plt.show()运行这段代码你会看到CEM如何工作。在等高线图上红色的点线显示了概率分布均值μ的移动轨迹。一开始它可能在搜索空间里跳跃但随着迭代它会稳定地朝着全局最优点(0,0)附近移动。同时“最佳得分”曲线会快速下降并趋于平稳“标准差”曲线会逐渐收缩到接近零这表明算法对解的位置越来越确定。4. 关键参数调优与实战避坑指南实现一个能跑的CEM只是第一步让它跑得好、跑得稳才是体现功力的地方。下面这些参数和技巧是我在多个项目实践中总结出来的很多是官方论文里不会细说的“坑”。4.1 核心参数如何设置不是拍脑袋决定的种群大小pop_size与精英比例elite_fracpop_size这是最重要的参数之一。太小采样不足容易陷入局部最优太大计算开销剧增。一个经验法则是pop_size至少是问题维度dim的10倍。对于10维的问题我通常会从100或200开始尝试。对于像Rastrigin这样多峰的函数需要更大的种群来维持多样性。elite_frac通常在0.1到0.3之间。比例越小选择压力越大收敛越快但也更容易早熟陷入局部最优。比例越大保留了更多样性收敛速度慢但探索能力更强。我个人的经验是对于复杂多峰问题可以从0.2开始如果发现早熟可以尝试增大到0.3对于相对平滑的问题可以降低到0.1以加速收敛。初始方差sigma代码中我们用了搜索范围的1/4。这是一个不错的起点。原则是初始分布要能覆盖你感兴趣或你认为最可能包含最优解的整个区域。如果你对最优解的位置一无所知就设大一些。如果你有先验知识比如知道大概在某个区间可以把均值设在那附近方差设小一些能加速收敛。协方差矩阵的处理在我们的简单实现中我们假设各维度独立协方差矩阵是对角阵。这在很多情况下是有效的尤其是当问题变量之间的耦合不强时。但如果变量间存在强相关性比如调整机翼前缘和后缘的参数是联动的使用完整的协方差矩阵能让算法更高效地捕捉这种关系沿着“好解”的等高线方向收缩。我们的实现已经计算了完整协方差这比只使用对角矩阵更强大。收敛阈值epsilon不要设得太小。通常1e-3到1e-5是一个合理范围。有时算法会因为方差收缩得太快而“卡住”此时即使没达到阈值最佳解也不再改善。更稳健的停止条件是同时监控最佳解在连续多次迭代中的改善程度。例如可以增加一个判断如果最近20代的最佳解改善幅度小于某个微小值如1e-6也停止迭代。4.2 实际项目中踩过的“坑”与解决方案坑一方差坍塌与早熟收敛这是CEM最常见的问题。算法过早地聚焦在一个局部最优区域方差迅速变得极小种群失去多样性再也跳不出来了。解决方案方差平滑更新不要直接用精英样本的方差替换旧方差而是采用平滑更新sigma_new alpha * sigma_elite (1-alpha) * sigma_old。其中alpha是一个学习率如0.7。这保留了部分历史信息防止突变。注入噪声在更新后的方差上加上一个小的常数项或按比例缩放确保方差不会低于一个最小值如sigma max(sigma, sigma_min)。这相当于始终保持一点探索的“火种”。重启策略当检测到收敛方差很小但解的质量不满意时保存当前最优解然后重新初始化分布比如把均值设在当前最优解附近但把方差重新放大开始新一轮优化。这相当于“局部搜索”“全局探索”的循环。坑二处理边界约束我们的简单实现用了np.clip进行裁剪。但这会导致一个副作用大量样本被挤压在边界上扭曲了真实的概率分布特别是当最优解就在边界附近时算法性能会下降。更好的解决方案反射边界如果采样点越界将它“反射”回边界内。例如对于变量x在[a, b]内如果采样值x_sampled b则令x 2*b - x_sampled如果x_sampled a则令x 2*a - x_sampled。这能保持样本在边界附近的分布特性。惩罚函数法不修改采样点但在评估目标函数时对越界的点施加一个很大的惩罚值比如加上一个与越界距离成正比的巨大正数。这样越界的点会在精英选择中被自然淘汰分布更新时会自动向可行域内移动。坑三高维灾难当问题维度dim很高时比如上百维CEM会面临挑战。首先采样那么多点pop_size需要更大计算开销大。其次估计一个高维的完整协方差矩阵需要大量的精英样本至少dim1个否则矩阵会是病态的。解决方案使用对角协方差或因子化模型如果假设维度独立可行就只用对角协方差参数从dim*dim降到dim个。或者使用更复杂的因子化分布如Cholesky因子来减少参数。增量更新与正则化使用秩1更新或秩μ更新来迭代更新协方差矩阵而不是每代重新计算。同时必须加入正则化就像我们代码里加的1e-6 * I防止矩阵奇异。5. 超越基础让CEM更强大的进阶技巧掌握了基础CEM后你可以通过以下技巧让它适应更复杂的场景。5.1 混合策略CEM作为局部搜索器CEM的全局探索能力在初期较强但后期精细搜索能力可能不如一些基于梯度的局部方法。一个常见的策略是将CEM与其他算法结合。例如CEM 梯度下降先用CEM进行全局粗搜索找到有希望的盆地区域然后将CEM找到的最优解作为起点用梯度下降法或L-BFGS等进行精细的局部优化。这结合了全局和局部搜索的优点。CEM in CMA-ES Style借鉴CMA-ES协方差矩阵自适应进化策略中的一些思想如进化路径。进化路径记录了均值更新的连续方向利用它来更新协方差矩阵可以更聪明地学习目标函数的拓扑结构加快收敛。这相当于给CEM增加了“记忆”和“动量”。5.2 处理噪声与鲁棒优化在实际工程中目标函数评估可能有噪声比如仿真有随机性或者实验测量有误差。噪声会干扰精英样本的选择可能导致算法不稳定。重采样平滑对同一个候选解进行多次评估取其平均值作为该解的得分。这增加了评估的可靠性但代价是计算量倍增。分位数选择替代Top-K不直接取前K个而是取分数低于某个分位数如75%分位数的所有样本作为“精英集”。这在一定程度上容忍了噪声因为偶尔评估差的“好解”和评估好的“差解”可能都会被纳入但在统计意义上更新方向仍是正确的。5.3 离散与混合优化问题标准的CEM针对连续变量。但很多问题包含离散变量如选择哪台机器、是否开启某个功能。离散化处理对于有序的离散变量如整数1,2,3...可以将其视为连续变量进行CEM优化采样后四舍五入到最近的离散值。更新分布时仍然使用四舍五入前的连续值来计算均值和方差。联合分布对于同时包含连续变量和分类变量的问题可以维护多个概率模型。例如连续部分用高斯分布分类部分用多项分布每个分类有一个概率。采样时分别从各自的分布中采样然后组合成一个完整解。更新时连续部分用高斯分布更新分类部分用精英样本中各类别的出现频率来更新多项分布的概率。6. 多元函数寻优实战一个简单的工程设计案例让我们设想一个更贴近实际的简单案例优化一个箱型梁的截面尺寸在满足强度、刚度要求下最小化重量。假设梁的截面由宽度b、高度h、腹板厚度t_w、翼缘厚度t_f四个变量决定均为连续正数。约束包括最大应力不超过许用应力最大挠度不超过允许值。目标函数是截面面积正比于重量。这是一个带约束的4维优化问题。我们可以用罚函数法将其转化为无约束问题将约束违反量乘以一个大的惩罚系数加到目标函数面积上。def beam_objective(x): 箱型梁优化目标函数罚函数法处理约束。 x [b, h, t_w, t_f] b, h, t_w, t_f x # 计算截面面积 (目标越小越好) area 2*(b*t_f) 2*( (h-2*t_f)*t_w ) # 简化模型 # 计算应力和挠度 (简化公式仅用于示例) # 假设受均布载荷q跨度L材料弹性模量E许用应力sigma_allow许用挠度delta_allow L, q, E, sigma_allow, delta_allow 5.0, 10000, 2.1e11, 200e6, L/500 # 示例参数 I (b*h**3 - (b-2*t_w)*(h-2*t_f)**3)/12.0 # 惯性矩 stress (q * L**2 * h) / (8 * I) # 最大弯曲应力简化 deflection (5 * q * L**4) / (384 * E * I) # 最大挠度简化 # 罚函数 penalty 0.0 if stress sigma_allow: penalty 1e6 * (stress - sigma_allow)**2 # 应力违反惩罚 if deflection delta_allow: penalty 1e6 * (deflection - delta_allow)**2 # 挠度违反惩罚 # 几何约束厚度必须小于宽度/高度的一半等 if t_f h/2 or t_w b/2: penalty 1e9 return area penalty # 定义变量边界 (单位米) bounds [(0.1, 0.5), # b (0.2, 1.0), # h (0.005, 0.02), # t_w (0.008, 0.03)] # t_f # 使用CEM优化 cem_beam CrossEntropyOptimizer(objective_funcbeam_objective, dim4, boundsbounds, pop_size80, elite_frac0.15, max_iter150, epsilon1e-5, seed123) best_design, best_weight, hist cem_beam.optimize() print(\n--- 箱型梁优化结果 ---) print(f最优截面尺寸: 宽度 b {best_design[0]:.4f} m, 高度 h {best_design[1]:.4f} m) print(f 腹板厚 t_w {best_design[2]:.4f} m, 翼缘厚 t_f {best_design[3]:.4f} m) print(f估算的最小截面面积 (重量指标): {best_weight:.6f} m²) # 可以重新计算一下应力和挠度验证是否满足约束通过这个例子你可以看到CEM如何在一个有实际物理意义、带约束的多元函数寻优问题上工作。它不需要你知道应力、挠度公式的梯度只需要能计算出给定尺寸下的“性能”面积惩罚即可。这种黑箱优化的特性使得CEM在仿真驱动设计、控制器参数整定等领域非常有用。最后我个人在多次使用CEM后的体会是把它看作一个“智能的随机搜索框架”。它的优势在于概念简单实现容易并行友好采样和评估可以完全并行并且对目标函数形式要求极低。它的性能很大程度上依赖于参数设置和那些“小技巧”如方差平滑、边界处理。不要期望它在任何问题上都碾压其他算法但在处理中低维度、计算代价昂贵、黑箱或无梯度的复杂优化问题时CEM绝对是一个值得你首先尝试的、强大而实用的工具。当你没有更好选择的时候试试CEM它常常能给你一个不错的起点甚至是一个惊喜。