简介：面向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-6；maxit：最大迭代次数，防止死循环，默认 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 的默认值，目标函数值域很小时，h=1e-6 会引入大舍入误差，可以把 h 调到 1e-4 左右再观察结果。
2.4 非二次函数与方向重置
共轭梯度在严格二次函数上有有限步收敛保证，工程问题大多不是二次函数，迭代若干轮后共轭关系被破坏。标准做法是每 n 步把搜索方向重置为负梯度。部分源码实现里没有写重置逻辑，如果发现iter达到maxit仍不收敛，不要无脑增大最大迭代次数，先在方向更新处加一句：
if mod(k, numel(x0)) == 0 beta = 0; end这里k是当前迭代步数，numel(x0)是变量维度。beta=0表示放弃之前的共轭修正，让搜索方向回到最速下降方向。这个操作对 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_{k+1} = 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 求梯度后等价于解正规方程X'X w = X'y。因此调用方式可以写成：
% 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 / GMRES | Jacobi 必须对角占优 |
从工程角度看，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)^2，x 是列向量没问题，但梯度函数返回行向量时，与 x 做运算会触发维度不一致。统一把变量写成列向量，梯度函数返回[...; ...]而不是[... , ...]。第二，内点法初始点不满足严格可行约束，log 函数直接产生 NaN。解决方法是先画约束函数，找一个离边界距离较大的点作为 x0，或者对边界点加一个小扰动。第三，外点法惩罚因子初值设得太大，第一轮迭代就让增广函数病态。从 rho0=0.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 里分节执行，每次迭代的中间变量都会保留在工作区，调试时直接看变量变化比命令窗口连续打印更直观。
本文还有配套的精品资源，点击获取
