MATLAB数值优化源码解析：共轭梯度法与罚函数实现

简介面向MATLAB优化算法学习与研究者的实用源码包覆盖梯度法、内点法、外点法、罚函数及线性梯度法等经典约束与无约束优化方法。全部程序为可直接运行的.m脚本用户只需在命令窗口按提示输入参数即可得到结果免去重复编写调试的麻烦适合算法入门、课程实验及论文复现使用。资源共6个文件均为MATLAB源程序压缩包仅4KB轻量易用当前已有1734人学习下载广受认可代码实现包含共轭梯度迭代、内点惩罚函数、外点惩罚函数等典型策略并配有二维示例用于可视化对比不同算法的收敛路径。通过运行这些程序可直观体会梯度法沿负梯度方向迭代、内点法从可行域内部逼近、外点法从外部逐步修正约束等核心思想理解罚函数如何通过惩罚项将约束问题转化为无约束问题各脚本逻辑清晰、注释简洁便于结合理论推导进行单步调试深入掌握优化算法的实现细节与适用场景是提升MATLAB编程能力和解决实际优化问题的得力工具。1. 从梯度法到罚函数这套MATLAB源码包解决什么问题如果你正在做机械优化设计或复现一本数值优化教材里的算例大概率遇到过这种情况fmincon 能给出结果却看不到惩罚因子怎么变、内点迭代路径如何贴着约束边界走论文需要的过程曲线一张都画不出来。这套源码包把六种基础数值方法拆成独立 .m 文件包括共轭梯度法二维与通用版、内点罚函数、外点罚函数、线性系数回归和 Jacobi 迭代。解压后按提示输入目标函数、约束函数和初始点每一步迭代的目标值、惩罚因子和迭代点都能被拿到。它不依赖优化工具箱基础 MATLAB 环境就能跑适合课程设计、教材复现和想手写优化器的人边读边改。2. 梯度法与共轭梯度两个共轭梯度文件的调用与验证2.1 从最速下降到共轭方向收敛性差在哪最速下降法沿着负梯度方向搜索在二次函数等值线为长椭圆时会出现明显的锯齿效应。原因很简单相邻两次迭代的梯度方向不满足正交性搜索路径反复纠正步长越来越小收敛极慢。共轭梯度法Conjugate Gradient构造一组关于 Hessian 矩阵共轭的方向在 n 维二次函数上理论上最多 n 步收敛而且不需要存储 Hessian 矩阵只需要每次迭代做一次矩阵向量乘因此成为中等规模无约束优化问题的常用选择。压缩包里有Conjugate_grad_2d.m和Conjugate_grads_method.m两个文件。前者的定位是二维演示方便把迭代点画在等高线上后者是通用 n 维实现名字里带复数形式内部大概率实现了 Fletcher-Reeves 或 Polak-Ribiere 公式并且包含一维搜索。理解这两个文件的区别比直接拿过来跑更重要因为很多教材里共轭梯度法的收敛性证明基于二次函数而实际工程目标函数很少是二次的。2.2 打开文件先看签名别急着运行拿到.m文件后先在编辑器中双击打开看第一行function的定义。这决定了传入参数顺序和返回值个数。常见签名大概是function [x_opt, iter, hist] Conjugate_grads_method(fun, grad, x0, tol, maxit)如果实际文件里变量名不同以头部定义为准。参数含义如下fun目标函数句柄必须写成(x) ...形式x 按列向量传入grad梯度函数句柄返回与 x 同维的列向量如果该文件支持数值差分可以不传x0初始点列向量长度与变量维度一致tol终止阈值常指相邻两次迭代点差值的二范数默认可取 1e-6maxit最大迭代次数防止死循环默认 100 到 500。调用前先在命令窗口确认which Conjugate_grads_method.m能返回路径否则后面所有脚本都会报 Undefined function。2.3 一个可直接复现的二维无约束算例为了验证共轭梯度法的收敛性我习惯用椭圆等值线的二次函数测试% demo_cg_2d.m f (x) (x(1) - 2).^2 3*(x(2) 1).^2; grad (x) [2*(x(1) - 2); 6*(x(2) 1)]; x0 [-3; 2]; tol 1e-8; maxit 50; [x_opt, iter, hist] Conjugate_grads_method(f, grad, x0, tol, maxit); fprintf(最优解: [%.6f, %.6f]\n, x_opt(1), x_opt(2)); fprintf(迭代次数: %d\n, iter);这段代码里3*(x(2)1)^2人为拉长了等值线让问题从圆形变成椭圆。共轭梯度法应当一步或两步就收敛到[2, -1]而最速下降法需要几十步甚至上百步。运行后如果迭代次数超过 10常见原因是梯度函数写错比如二次函数求导后忘记乘系数。把hist打印出来看目标值序列是否单调下降如果不单调说明一维搜索的步长没有精确实现。对于Conjugate_grad_2d.m它通常只接受目标函数和初始点内部用有限差分近似梯度省去手写梯度函数。这种情况下需要注意差分步长 h 的默认值目标函数值域很小时h1e-6 会引入大舍入误差可以把 h 调到 1e-4 左右再观察结果。2.4 非二次函数与方向重置共轭梯度在严格二次函数上有有限步收敛保证工程问题大多不是二次函数迭代若干轮后共轭关系被破坏。标准做法是每 n 步把搜索方向重置为负梯度。部分源码实现里没有写重置逻辑如果发现iter达到maxit仍不收敛不要无脑增大最大迭代次数先在方向更新处加一句if mod(k, numel(x0)) 0 beta 0; end这里k是当前迭代步数numel(x0)是变量维度。beta0表示放弃之前的共轭修正让搜索方向回到最速下降方向。这个操作对 Rosenbrock 这类非二次函数特别有效也是教材里很少强调但实际调试必用的技巧。3. 内点法与外点法两个罚函数程序的原理和参数设置3.1 罚函数法的统一形式内点法InteriorPenaltyFunctionMethod.m和外点法ExteriorPenaltyFunctionMethod.m都属于序列无约束极小化技术。核心思想是把约束条件以惩罚项形式叠加到目标函数上构造一个新函数min f(x) rho_k * P(x)其中 rho_k 是惩罚因子随外层迭代逐步增大。内点法要求迭代点始终在可行域内部P(x) 在边界处趋于无穷大常见形式是-log(g(x))或1/g(x)外点法则不限制迭代点位置允许在可行域外计算P(x) 通常取max(0, g(x))^2或等式约束的h(x)^2。两者的区别直接决定了初始点选择和参数调整方式。3.2 内点罚函数严格可行初始点是硬前提内点法的对数障碍函数在g(x) 0时未定义因此初始点必须严格可行这一点经常被忽略。调用示例如下% demo_interior.m % 约束 g(x) x1 x2 - 1 0 f (x) x(1).^2 x(2).^2; constraint (x) x(1) x(2) - 1; x0 [1; 2]; % 严格可行点 rho0 1; c 5; tol 1e-6; [x_opt, iter] InteriorPenaltyFunctionMethod(f, constraint, x0, rho0, c, tol);参数说明rho0是初始惩罚因子控制迭代点与边界的距离取值太小会导致第一次无约束优化时目标函数主导迭代点几乎贴着边界c是外层循环中rho c * rho的放大系数一般取 5 到 10。tol是外层收敛阈值通常检查约束残差和目标函数变化量。如果运行返回 NaN先检查x0代入constraint是否严格大于 0因为浮点在边界附近算 log 会溢出。3.3 外点罚函数对初值宽容但参数敏感外点法不要求初始点可行甚至可以从远离可行域的点开始这是它相对内点法最大的工程优势。调用方式类似% demo_exterior.m f (x) x(1).^2 x(2).^2; h (x) x(1) - 3; % 等式约束 x1 3 x0 [0; 0]; rho0 1; c 8; tol 1e-6; [x_opt, iter] ExteriorPenaltyFunctionMethod(f, h, x0, rho0, c, tol);这里h是等式约束函数外点法用rho/2 * h(x)^2作为惩罚项。原问题最优解是[3; 0]但外点法得到的结果通常会在 3 附近震荡rho 越大越逼近精确解却又越容易让增广函数病态。判断收敛不能只看目标函数值还应该检查abs(h(x_opt))是否小于 tol。如果源码内部的无约束优化用的是最速下降步长需要设置为 0.01 量级否则外层 rho 跳动太大时容易发散。3.4 内点法与外表法的适用性对比对比项内点法外点法初始点位置必须严格可行任意点惩罚项形式对数或倒数障碍函数二次惩罚项迭代点轨迹始终在可行域内部可能在可行域外等式约束处理不方便直接处理直接处理终止条件障碍项趋于 0约束残差趋于 0典型问题不等式约束为主混合约束或初值难找从源码结构上讲两个文件的内层循环基本一致区别只在于增广函数的构造方式。如果想把内点法改成处理等式约束需要额外引入等式障碍项这已经不是简单的参数调整而是算法实现层面的改动。提示内点法初始点必须是严格内部点边界点或不可行点会让 log 障碍函数直接返回 Inf先画出可行域或者随机抽样一个可行点比反复改 rho 更有效。3.5 rho 序列怎么调才不炸罚函数法最常踩的坑是 rho 初值和放大倍率不匹配。我一般建议初始rho0取 1 以下尤其是目标函数值域在 100 以上时先对目标函数做归一化否则第一轮无约束优化就被惩罚项主导。放大倍率c不要超过 10比较稳妥的是 5 到 8。停止条件最好同时检查约束残差和梯度范数if abs(constraint(x)) 1e-4 norm(grad(x)) 1e-4 break; end这种双条件判断比只看步长或者只看目标函数更可靠也是我处理惩罚函数发散问题时最先加的判断。4. 线性梯度法与 Jacobi 迭代配套的线性代数工具4.1 优化和线性方程组的关联压缩包里的LinearCofficientMethod.m和Jacobi_iterative_method.m表面上与前几章的优化算法无关实际上它们解决的是同一类问题。求解线性方程组Ax b可以等价为极小化二次函数0.5 * x A x - b x当 A 对称正定时这个二次函数的极小点就是方程组的解。这意味着梯度法、共轭梯度法都可以直接用来解线性系统反过来Jacobi 迭代又是最基础的无约束迭代方法把两者放在一起既能用来做数值实验也能在罚函数内部实现一维搜索时提供参照。4.2 Jacobi 迭代的收敛前提对角占优Jacobi 迭代把矩阵 A 分解为对角阵 D 和剩余部分 R迭代公式为x_{k1} D^{-1} (b - R x_k)当 A 严格对角占优时收敛保证较强。调用示例% demo_jacobi.m A [5, 1, 0; 1, 4, 1; 0, 1, 3]; b [6; 6; 4]; x0 zeros(3, 1); tol 1e-10; maxit 200; [x, iter, residual] Jacobi_iterative_method(A, b, x0, tol, maxit);这个矩阵每一行对角线元素的绝对值都大于该行其他元素绝对值之和是验证 Jacobi 收敛的标准用例。运行后残差norm(A*x - b)应持续下降迭代次数通常在 20 次以内。如果换成非对角占优矩阵残差会先下降再上升最后变成 NaN。这说明问题本身不适合 Jacobi应该换 Gauss-Seidel 或 SOR不能靠增加maxit解决。4.3 用线性梯度法求回归系数LinearCofficientMethod.m我通常把它理解为求解线性最小二乘问题也就是线性模型的系数估计。目标是最小化||Xw - y||^2对 w 求梯度后等价于解正规方程XX w Xy。因此调用方式可以写成% demo_lineargrad.m X randn(50, 3); % 50 个样本3 个特征 y X * [1; -2; 0.5] 0.01 * randn(50, 1); A X * X; g X * y; % 如果函数接受矩阵和右端项 beta LinearCofficientMethod(A, g);参数说明A是 Gram 矩阵g是右端项函数内部用迭代法求解A * beta g。需要注意数据尺度问题如果两个特征数量级差很多Gram 矩阵条件数会很大梯度法收敛很慢。先对 X 做 zscore 标准化算完系数后还原比直接改迭代次数更有效。如果源码里的函数要求传入X和y而不是A和g内部会自动构造正规方程。4.4 选择迭代法的判断标准场景推荐方法原因教学演示、低维问题Jacobi实现简单路径直观大规模稀疏线性系统共轭梯度无需显式存储分解因子对称正定矩阵共轭梯度理论有限步收敛非对称或对角占优Jacobi / GMRESJacobi 必须对角占优从工程角度看Jacobi 的实际性能远不如 Krylov 子空间方法但它的价值在于代码简单可以当作验证矩阵性质的工具。我写新的优化算法时会先用 Jacobi 确认矩阵和右端项没有拼错再切换到共轭梯度。5. 让这些程序跑起来路径设置、初值选择与三个易错点5.1 解压、加路径与签名确认把压缩包解压到不含中文和空格的目录例如D:\optimization_src然后在 MATLAB 命令窗口执行addpath(genpath(D:\optimization_src));之后用which Conjugate_grads_method.m验证路径是否生效。每个.m文件在运行前都要检查头部 function 行确认返回值是多个变量还是一个结构体这决定了后续如何取出迭代次数和优化结果。5.2 三个高频报错点第一匿名函数使用错误。目标函数写成(x) x(1)^2 x(2)^2x 是列向量没问题但梯度函数返回行向量时与 x 做运算会触发维度不一致。统一把变量写成列向量梯度函数返回[...; ...]而不是[... , ...]。第二内点法初始点不满足严格可行约束log 函数直接产生 NaN。解决方法是先画约束函数找一个离边界距离较大的点作为 x0或者对边界点加一个小扰动。第三外点法惩罚因子初值设得太大第一轮迭代就让增广函数病态。从 rho00.1 开始观察约束残差变化残差下降缓慢再增大 c。5.3 用 Rosenbrock 函数做收尾验证用 Rosenbrock 函数验证整套流程是稳妥的收尾方式f (x) 100*(x(2) - x(1)^2)^2 (1 - x(1))^2; grad (x) [-400*x(1)*(x(2) - x(1)^2) - 2*(1 - x(1)); 200*(x(2) - x(1)^2)]; x0 [-1.2; 1]; [x_opt, iter, hist] Conjugate_grads_method(f, grad, x0, 1e-6, 100);Rosenbrock 函数的最优解是[1; 1]但谷底是一条弯曲的弧形最速下降会走 Z 字形共轭梯度配合方向重置能在几十步内收敛。如果结果偏差大于 1e-4在迭代循环里临时加一行fprintf(%d: %.6f\n, k, norm(grad(x)))观察梯度范数是否下降。出现方向不下降时按 2.4 节把 beta 重置为 0 再跑。这套源码不依赖优化工具箱基础 MATLAB 环境直接运行比安装完整工具箱要省事得多。如果使用的是新版 MATLAB建议把脚本复制到 Live Editor 里分节执行每次迭代的中间变量都会保留在工作区调试时直接看变量变化比命令窗口连续打印更直观。本文还有配套的精品资源点击获取