简介这份资源是一套基于MATLAB的BFGS拟牛顿法实现代码面向数值优化学习者、算法研究者以及需要求解无约束优化问题的工程人员。代码共包含5个m文件涉及目标函数定义、梯度计算、线性搜索策略、BFGS迭代更新与结果打印等核心模块压缩包大小仅1KB轻量精简且便于对照学习。目前已有456人浏览学习适合初学者理解拟牛顿法的算法流程也适合进阶者改造复用。BFGS方法通过递推更新正定近似Hessian矩阵避免直接计算二阶导数的昂贵开销资源中配套的线性搜索与收敛判断模块有助于掌握Armijo准则、梯度下降方向选取等关键细节。读者可借助fun.m、dfun.m快速替换自定义目标函数与梯度在二次优化、参数估计等场景中直接运行验证从而深入体会二次收敛性与算法稳定性是一份短小精悍的算法教学与参考工具。1. 为什么高维二次优化场景里BFGS比牛顿法更值得写进代码做参数辨识、最优控制或带罚项的回归估计时第一步往往不是选模型而是选优化器。牛顿法每轮要重新组装并分解真实 Hessian 矩阵维度上千之后O(n^3)的代价非常难看BFGSBroyden-Fletcher-Goldfarb-Shanno用相邻迭代点的梯度差去构造 Hessian 近似既保住牛顿法那种大步流星的方向感又把每步代价压到O(n^2)的矩阵更新和一次矩阵乘向量。真正实现时更常用的还是逆 Hessian 近似矩阵H让搜索方向从“解线性方程组”降级成“一次乘法”代码本身也就短了一大截。这套BFGS.rar里的五个 MATLAB 文件正好把这件事讲透了BFGS.m主循环、fun.m与dfun.m负责目标函数和解析梯度、f_search.m做 Armijo 回溯线搜索、print.m输出迭代日志。适合刚接触拟牛顿法、想在 MATLAB 里把自定义函数塞进求解器的研究者也适合想读一遍源码就复现整套优化流程的工程师。2. 割线方程与正定修正BFGS 的数学结构2.1 为什么拟牛顿法要绕开真实 Hessian牛顿法的迭代步是求解H_k d -g_k其中H_k是目标函数在x_k处的 Hessian 矩阵。对 n 维问题组装 Hessian 需要n^2量级的二阶信息Cholesky 分解又要(1/3)n^3次浮点运算。一旦维度涨到数千这基本就是灾难。拟牛顿族的思路是换一个约束不再碰真正的二阶导数只要求近似矩阵B_k满足割线方程B_{k1} s_k y_k其中s_k x_{k1} - x_ky_k ∇f(x_{k1}) - ∇f(x_k)。这个方程由一维牛顿法中割线的思想推广而来是拟牛顿矩阵唯一必须满足的硬约束。问题在于满足它的对称矩阵非常多选哪一个直接决定了算法的收敛性。经典的对称秩一校正SR1实现最简单但无法保证每一步更新后的矩阵正定BFGS 选择的是对称秩二修正在满足割线方程的同时还能把上一轮矩阵的正定性完整继承下来。下面用一段 MATLAB 代码验证更新后的B_new是否真的满足割线方程% 验证 BFGS 更新后的矩阵是否满足割线方程 B_new * s y s x_new - x; % 位移向量 y grad_new - grad; % 梯度差向量 % B 版本 BFGS 更新公式 B_new B - (B * (s * s) * B) / (s * B * s) (y * y) / (y * s); max(abs(B_new * s - y)) % 应该接近 eps 量级这段代码的验证逻辑很直接第一项(B * s * s * B) / (s * B * s)是一个秩一修正负责把B沿s方向调整到割线方程第二项(y * y) / (y * s)同样是对称秩一矩阵用来补偿y方向上的信息。两者叠加就是完整的秩二更新所以最后max(abs(B_new * s - y))应该落在eps量级。如果计算出的残差明显偏大先检查s和y是不是取自同一对迭代点再检查分母y*s是否接近零。先给一个对比表看 BFGS 在工程选型中的位置方法每步额外存储单步计算量收敛速率适用规模牛顿法n×n真实 HessianO(n^3)分解局部二次收敛中低维度Hessian 可解析BFGSn×n近似矩阵O(n^2)更新局部超线性收敛数千维以内L-BFGSm个向量对O(mn)局部超线性收敛数万到百万维2.2 正定性为什么能一直保持BFGS 更新公式中真正起稳定作用的条件是y * s 0。直观理解在一维情形下这等价于f(x_{k1}) f(x_k)即割线斜率必须为正才能保证二次近似开口向上。推广到多维ys 0是矩阵更新后保持正定的充分条件。对于凸函数只要线性搜索满足 Wolfe 条件的曲率部分这个不等式就自动成立这也是为什么f_search.m不能简单返回一个固定步长它承担着维持矩阵正定的职责。实际工程处理时遇到y*s 0的情况并不少见例如目标函数存在数值噪声或线搜索过于粗糙。稳健做法是跳过本轮更新而不是强行计算if y * s 1e-10 rho 1 / (y * s); H (I - rho * s * y) * H * (I - rho * y * s) rho * (s * s); else % 跳过更新保留上一轮 H end注意跳过更新只是少用一次二阶信息最多让这一步退化为梯度下降强上公式则可能让矩阵失去正定性下一步搜索方向直接变成上山方向。从变分角度看BFGS 更新是所有满足割线方程的正定矩阵中与当前矩阵在加权重范数意义下距离最近的解这也是它比 DFP 在数值表现上更稳的理论原因。2.3 收敛速率超线性还是二次关于 BFGS 的收敛速率有一个需要较真的点很多资料直接写“二次收敛”严格说是不准确的。对一般光滑非线性函数BFGS 在解附近具备的是局部超线性收敛对正定二次函数配合精确一维搜索最多 n 步达到精确解这叫做二次终止性。牛顿法才是真正意义上的局部二次收敛因为它的每一步都基于真实 Hessian。这个区别在实际实验里能直接观察BFGS 进入超线性阶段后线性搜索返回的步长alpha会稳定在 1 附近因为 Hessian 近似已经足够准完整拟牛顿步就是最好的方向。如果跑了十几轮alpha还在 0.1 以下飘忽不定多半是矩阵更新或梯度实现有问题。停机条件同样值得较真。函数值在极小点附近非常平坦用“函数值变化小于阈值”停机容易过早退出更常用的是梯度范数判据norm(g, inf) tol阈值建议取1e-6在双精度浮点下不要低于1e-8否则线性搜索求得的步长会淹没在舍入噪声里。3. BFGS.rar 源码逐文件拆解从主循环到线搜索3.1 五个文件的分工与调用链压缩包里的结构是典型的教学型实现每个文件只干一件事。先看总体分工文件职责常见接口被谁调用BFGS.m主循环初始化、方向计算、矩阵更新、停机判断BFGS(fun, dfun, x0, opts)用户直接调用fun.m目标函数f fun(x)BFGS.m、f_search.mdfun.m解析梯度g dfun(x)BFGS.mf_search.mArmijo 回溯线搜索alpha f_search(fun, dfun, x, d, f0, g0)BFGS.mprint.m迭代日志输出print(k, x, f, ng)BFGS.m调用链是一条直线BFGS.m先调fun.m和dfun.m拿到初始函数值与梯度然后进入循环每次先调f_search.m求步长更新位置后再调fun.m和dfun.m计算新梯度最后用print.m打日志。理解这个结构后换自己的目标函数只需要替换fun.m和dfun.m其余文件基本不用动。3.2 BFGS.m 主循环逆 Hessian 更新为什么是主流教学实现里常见两种矩阵近似 Hessian 的B和近似逆 Hessian 的H。B版本更贴近数学定义但求搜索方向时要解B \ (-g)相当于每轮都做一次线性求解H版本把更新公式翻转过来搜索方向直接是d -H * g省掉了整个求解过程。工程实现几乎清一色选择后者。核心循环长这样% 初始化 x x0(:); f feval(fun, x); g feval(dfun, x); n length(x); H eye(n); % 逆 Hessian 近似初始为单位矩阵 d -H * g; % 搜索方向无需解线性方程组 for k 1:maxiter alpha f_search(fun, dfun, x, d, f, g); x_new x alpha * d; f_new feval(fun, x_new); g_new feval(dfun, x_new); s x_new - x; y g_new - g; yTs y * s; if yTs eps rho 1 / yTs; I eye(n); % BFGS 逆 Hessian 更新公式 H (I - rho * s * y) * H * (I - rho * y * s) rho * (s * s); end x x_new; f f_new; g g_new; d -H * g; if norm(g, inf) tol break; end if mod(k, print_interval) 0 print(k, x, f, norm(g, inf)); end end参数说明s和y都是 n 维列向量rho是标量步长修正系数H始终保持对称正定所以d -H * g一定是下降方向。yTs eps判断保证了更新公式分母安全也维护H的正定性。norm(g, inf)是无穷范数比二范数更容易暴露某个分量没有收敛的情况适合做停机判据。如果拿到的是B版本实现对应替换为B B - (B * (s * s) * B) / (s * B * s) (y * y) / (y * s); d -B \ g;逻辑上两者等价但B \ g的计算代价明显更高。对教学理解B版本更直观对性能敏感场景换成H版本。3.3 f_search.mArmijo 回溯线搜索如何保证充分下降线性搜索是拟牛顿法的“底盘”底盘不稳的话前面的矩阵更新做得再漂亮也会翻车。f_search.m最常见的实现是 Armijo 回溯function alpha f_search(fun, dfun, x, d, f0, g0) c1 1e-4; % Armijo 条件常数经验值 rho 0.5; % 回溯缩放因子 alpha 1; % 优先尝试完整拟牛顿步 while feval(fun, x alpha * d) f0 c1 * alpha * (g0 * d) alpha rho * alpha; if alpha 1e-10 break; % 防止死循环给一个下限 end end end逻辑说明初始步长取 1因为拟牛顿法进入超线性阶段后完整步长通常就是最优选择如果函数没有充分下降就把步长乘以rho0.5退回直到满足充分下降条件。c1经验上取1e-4它控制下降量的严格程度取值太接近 1 会让步长被过度截断取太小则可能接受一个几乎没变化的点。注意到接口里有dfun但函数体内没用它这是刻意保留的扩展位Armijo 条件只需要函数值但换成 Wolfe 条件时就需要梯度信息做曲率检验。教学代码保留这个参数方便后续升级。如果观察alpha序列在迭代后期总是稳定在 1 附近说明矩阵近似已经足够好这是收敛进入快车道的一个信号。3.4 print.m 迭代日志与分析print.m的输出格式决定了你能从一段跑完的实验中读到多少信息。常见实现是function print(k, x, f, ng) fprintf(Iter %4d | f %.8e | ||g||inf %.6e | x1 %.6f\n, ... k, f, ng, x(1)); end字段含义f是目标函数值||g||inf是梯度无穷范数x(1)是第一个分量的当前值。看日志时优先盯两处梯度范数是否单调下降x(1)是否还在明显漂移。如果梯度范数卡在某个量级半天不动先怀疑线搜索容差太大或目标函数存在数值噪声如果函数值先降后升大概率是线搜索接受了不满足足够下降条件的点需要回过去检查c1参数和回溯下限。4. 在 MATLAB 里跑通 BFGS.rar从二维测试到百维扩展4.1 用 Rosenbrock 函数做基准实验校验优化算法最经典的测试函数是 Rosenbrock 函数f(x) 100 * (x(2) - x(1)^2)^2 (1 - x(1))^2它的极小值点在(1, 1)但周围是一条香蕉形峡谷梯度方向与指向极小点的方向并不一致专门考验算法对二阶信息的利用能力。替换fun.m和dfun.mfunction f fun(x) f 100 * (x(2) - x(1)^2)^2 (1 - x(1))^2; end function g dfun(x) g [-400 * x(1) * (x(2) - x(1)^2) - 2 * (1 - x(1)); 200 * (x(2) - x(1)^2)]; end从初始点x0 [-1.2; 1]出发tol 1e-6在常见 64 位 MATLAB 版本下会得到类似下面这种日志节选Iter 0 | f 2.420e01 | ||g||inf 2.153e01 | alpha 1.000 Iter 1 | f 1.820e00 | ||g||inf 5.342e00 | alpha 0.500 Iter 2 | f 4.028e-01 | ||g||inf 2.107e00 | alpha 1.000 Iter 8 | f 3.913e-04 | ||g||inf 6.893e-03 | alpha 1.000 Iter 15 | f 4.581e-11 | ||g||inf 1.729e-05 | alpha 1.000注意一个典型现象第一次迭代步长被削到 0.5因为初始点的二次近似还不够好后面alpha稳定在 1进入超线性收敛阶段。最终大概 12 到 20 步内满足停机条件具体步数与 MATLAB 版本、线搜索参数有关。如果你跑出来的迭代次数明显偏多先看是不是ys eps判断把矩阵更新频繁跳过了。4.2 扩展到 100 维停机条件与内存边界把 Rosenbrock 函数推广到 n 维得到扩展版function [f, g] gen_rosen(x, want_grad) n length(x); f 0; g zeros(n, 1); for i 1:n-1 f f 100 * (x(i1) - x(i)^2)^2 (1 - x(i))^2; if want_grad g(i) g(i) - 400 * x(i) * (x(i1) - x(i)^2) - 2 * (1 - x(i)); g(i1) g(i1) 200 * (x(i1) - x(i)^2); end end end这个实现直接循环累加避免显式构造稀疏 Jacobian逻辑清晰缺点是梯度循环不是向量化写法但作为教学实验完全够用。100 维场景里BFGS 需要维护一个100×100的近似矩阵内存占用约 0.08 MB毫无压力到 1000 维时是 8 MB仍然可以接受5000 维以上n×n矩阵开始触及内存瓶颈就该切换 L-BFGS。停机条件的选择在扩展测试里更值得关注三种常见策略对照如下停机方式判断量风险梯度范数norm(g, inf) 1e-6噪声敏感但最通用函数值变化abs(f_new - f) 1e-10峡谷地形容易早停相对步长norm(x_new - x) 1e-8收敛慢时步长先于精度衰减在扩展 Rosenbrock 这类强峡谷问题上我一般以梯度范数为主判据函数值变化只作参考。否则函数值在峡谷底部变化极其缓慢时会提前触发停机返回的点离真实极小点还很远。4.3 与 fminunc 和 MATLAB 优化工具箱对照MATLAB 优化工具箱自带fminunc在指定解析梯度时走的就是 quasi-newton 路径对应 BFGS 算法。对其做一次对齐测试opts optimoptions(fminunc, ... SpecifyObjectiveGradient, true, ... Display, iter, ... OptimalityTolerance, 1e-6); [xf, ff, exitflag] fminunc(rosen_with_grad, x0, opts);其中rosen_with_grad是同时返回函数值和梯度的函数句柄。不指定Algorithm时fminunc默认走quasi-newton路径内部实现本质上就是 BFGS只是线搜索策略和矩阵更新细节经过工业级打磨比教学版更稳。对照实验通常会得到两个结论提供解析梯度比不提供少约 10% 到 20% 迭代次数因为有限差分梯度在峡谷处会引入数值噪声工具箱实现的线搜索在强非凸问题上更抗造教学版偶尔要手动调c1。提示fminunc在没有提供 Hessian 时会把 quasi-newton 作为默认路径因此它常被当作 MATLAB 里现成的 BFGS 可比对象。工具箱不需要额外安装MATLAB 自带优化工具箱即可。5. 从 BFGS 到 L-BFGS换一种 Hessian 记忆方式5.1 两循环递归去掉 n×n 矩阵BFGS 要维护一个n×n的稠密近似矩阵这是它面对大规模问题时唯一的软肋。L-BFGS 的想法很直接不存矩阵只存最近m组(s_i, y_i)向量对需要计算H_k g时用两循环递归现场算。核心代码是function r lbfgs_direction(H0, g, s_mem, y_mem) q g; m length(s_mem); alpha zeros(m, 1); % 第一次循环从最新到最旧 for i m:-1:1 alpha(i) (s_mem{i} * q) / (y_mem{i} * s_mem{i}); q q - alpha(i) * y_mem{i}; end r H0 * q; % H0 通常取单位矩阵乘尺度系数 % 第二次循环从最旧到最新 for i 1:m beta (y_mem{i} * r) / (y_mem{i} * s_mem{i}); r r (alpha(i) - beta) * s_mem{i}; end endm的经验取值是 5 到 20强凸问题取小值非凸问题取大值。这样内存从O(n^2)降到O(mn)几万维的参数空间也能跑。深度网络参数优化里也经常出现 L-BFGS 的变体虽然大热门是自适应步长类算法但小批量强凸子问题上 L-BFGS 仍有一席之地。5.2 一个收敛精度技巧混合停机条件把工程代码里的停机逻辑从单一梯度范数改成双判据能少踩很多坑rel_step norm(x_new - x_old) / (norm(x_old) 1e-12); if norm(g, inf) tol || (rel_step 1e-8 abs(f_new - f_old) 1e-10) break; end单看梯度范数有个盲区当目标函数由仿真程序返回或者梯度用有限差分近似时数值噪声会让norm(g, inf)永远压不到1e-6以下于是算法一直空转到最大迭代次数。混合判据里相对步长rel_step用来捕捉“位置已经不再移动”的事实函数值变化则补充验证目标确实到了平坦区。对于把 BFGS 或 L-BFGS 集成进实际项目的人来说这一行判断往往比调半天线搜索参数更管用。本文还有配套的精品资源点击获取