BFGS拟牛顿法结合Armijo线搜索：MATLAB非线性优化算法实现

前阵子做机器人标定参数反演，目标函数有三十多个参数，每次求值都要跑一整遍运动学正解，梯度存在但带有数值噪声。我一开始图省事直接甩给fminunc，结果某些初始点收敛精度上不去；换最速下降法倒是稳定了，可收敛速度慢得让人怀疑人生。最后自己用 BFGS 拟牛顿方向加上 Armijo 回溯线搜索实现了一套优化器，前后跑了几轮实验，才把问题彻底搞定。

这篇文章就把这套非线性优化算法在 MATLAB 里的实现完整拆开讲一遍。从“为什么要用 BFGS + Armijo 的组合”，到数学原理、再到可直接复制运行的完整代码，以及我实际跑实验时踩过的坑和调参记录，全部写出来。适合两类人看：一类是想搞懂拟牛顿法和线搜索原理的初学者，另一类是手里有真实优化问题、不想被工具箱黑盒限制的工程向同学。

1. 为什么是 BFGS + Armijo，而不是别的组合

很多初学优化的人一上来就写最速下降法，因为代码实在太简单了，方向取负梯度，步长用固定值或手动调。但真用起来，最速下降法面对稍微复杂一点的目标函数就会让人崩溃。

1.1 最速下降法的尴尬：稳定但慢得离谱

最速下降法的迭代格式是：

x_{k+1} = x_k + α_k * (-∇f(x_k))

每一步只用了当前点的一阶信息，也就是梯度方向。当目标函数的等高线是椭球形状、条件数比较大时，最速下降法会走出非常明显的之字形路线，来回震荡，收敛速度是线性的，而且线性收敛的常数往往非常接近 1。

我用 Rosenbrock 函数做过对比实验，这个函数长这样：

f(x) = 100 * (x₂ - x₁²)² + (1 - x₁)²

从 (-1.2, 1) 出发，最速下降法动辄几百上千次迭代才能逼近极小点。如果你把中间路径画出来，会发现它像在一条狭长山谷里反复横跳，效率极低。在实际工程里，目标函数每求值一次可能要跑几分钟的仿真，这种收敛速度完全不可接受。

1.2 牛顿法的代价：精度高，但二阶信息算不起

牛顿法的想法很直接：既然最速下降法只用了局部的一阶信息，那我把二阶曲率信息也用上，不就能更准地判断下山方向了吗？

x_{k+1} = x_k - [∇²f(x_k)]⁻¹ * ∇f(x_k)

其中 ∇²f(x_k) 是 Hessian 矩阵。理论上牛顿法在极小点附近能实现二次收敛，收敛速度快得惊人，迭代几步就到了。但工程里用纯牛顿法有三个很现实的问题：

• 目标函数可能是仿真函数、黑盒函数，根本没有解析二阶导。
• 即使能算二阶导，Hessian 矩阵也未必正定。如果 Hessian 有负特征值，牛顿方向甚至可能不是下降方向。
• 每次迭代都要计算并存储 n×n 的 Hessian，当参数有几千个时，这个计算量和存储量都是大问题。

所以工程上更多使用“拟牛顿法”这个折中路线。

1.3 BFGS 的思路：用梯度变化去近似 Hessian

BFGS 是拟牛顿法里最经典、最稳健的一种。它的核心思想是：不直接算 Hessian，而是利用相邻两次迭代的梯度变化量，迭代出一个近似 Hessian 矩阵 B。这个 B 会在迭代过程中越来越接近真实的 Hessian，从而让搜索方向逐渐逼近牛顿方向，同时避开了二阶导计算。

它和 Armijo 线搜索是天生一对。因为 BFGS 给出的只是一个“近似”牛顿方向，刚开始迭代时 Hessen 近似可能很不准，如果固定步长，很容易出现函数值不降反升的情况。Armijo 回溯线搜索能保证每一步都满足充分下降条件，把整个算法从“可能发散”变成“一定下降”，最终收敛到稳定点。

1.4 为什么不用 fminunc 一把梭

MATLAB 优化工具箱里的fminunc确实内置了 BFGS 实现，用起来也方便。我自己的经验是：工具箱适合快速拿结果，但一旦涉及以下场景，自己实现反而更合适。

• 需要自定义线搜索准则或者变体策略。
• 需要每次迭代输出中间变量、保存历史、做可视化。
• 需要把优化器嵌入到更大框架里，比如参数辨识循环、模型预测控制等。
• 自己的目标函数有很多特殊情况，比如梯度存在噪声、函数在部分区域不可导等。

而且从学习角度，亲手实现一遍 BFGS + Armijo，你会真正明白拟牛顿法每一步在做什么。工具箱封装得太好，反而容易让人雾里看花。

2. BFGS 数学细节：割线方程、更新公式与正定性

要在 MATLAB 里写出正确的代码，光看懂伪代码不够，还得理解公式从哪来的。否则遇到数值问题，你根本不知道怎么诊断。

2.1 从牛顿法到割线方程

我们记目标函数为 f(x)，梯度为 g(x) = ∇f(x)，Hessian 为 H(x) = ∇²f(x)。牛顿法的迭代格式是：

x_{k+1} = x_k - H(x_k)⁻¹ * g(x_k)

BFGS 的思路是用一个对称正定矩阵 B_k 来近似 H(x_k)。那么 B 应该满足什么条件呢？把梯度在 x_{k+1} 处做一阶泰勒展开：

g(x) ≈ g(x_{k+1}) + H(x_{k+1}) * (x - x_{k+1})

令 s_k = x_{k+1} - x_k，y_k = g(x_{k+1}) - g(x_k)，把上面的式子整理一下，得到：

B_{k+1} * s_k = y_k

这是拟牛顿法的“割线方程”，也常被称为拟牛顿方程。它说的是：近似 Hessian 矩阵作用在位移向量上，应该等于梯度变化量。这个方程是用一阶信息约束二阶近似的基本手段。

2.2 BFGS 更新公式的推导直觉

满足割线方程的 B_{k+1} 有很多个，BFGS 的想法是：在所有满足条件的对称正定矩阵中，选一个与当前 B_k 距离最近的矩阵。这个问题最终可以解出显式公式。

BFGS 更新公式有两个常用形式。一个是直接更新 B，也就是 Hessian 近似：

B_{k+1} = B_k + y_k * y_kᵀ / (y_kᵀ s_k) - B_k * s_k * s_kᵀ * B_kᵀ / (s_kᵀ * B_k * s_k)

另一个是更新 Hessian 逆矩阵的近似：

H_{k+1} = (I - ρ * s_k * y_kᵀ) * H_k * (I - ρ * y_k * s_kᵀ) + ρ * s_k * s_kᵀ

其中 ρ = 1 / (y_kᵀ s_k)。

这两个公式本质等价。第一个公式在理解和教学上更直观，因为它的形式直接对应“加一项、减一项”的秩二修正；第二个公式在实现时更便于直接计算搜索方向，因为搜索方向 p_k = -H_k * g_k 不需要求解线性方程组。

2.3 正定性：这是 BFGS 能稳定的关键

BFGS 有一个优秀的性质：只要 B_k 对称正定，并且满足曲率条件 y_kᵀ s_k > 0，那么更新后的 B_{k+1} 仍然对称正定。这一步可以从数学上严格证明，直觉上也很容易理解。

正定性为什么那么重要？因为搜索方向是 p_k = -B_k⁻¹ * g_k。当 B_k 正定时，p_k 与负梯度的夹角一定小于 90 度，也就是说这是一个严格的下降方向。如果 B 失去了正定性，搜索方向可能直接朝山上走，整个算法就彻底崩了。

这里我要特别提醒一点：虽然 Armijo 条件在理论上能保证收敛到满足曲率条件的区域，但在数值实现中，由于浮点误差或者目标函数本身不够光滑，y_kᵀ s_k 有时候确实会小于等于 0。所以代码里必须处理这种情况，后面我会详细讲。

2.4 搜索方向计算的正确姿势

公式上，BFGS 的搜索方向是：

p_k = -B_k⁻¹ * g_k

但实际写代码时，千万别用inv(B) * g这种写法。数值分析领域有个共识：永远不要显式求矩阵逆。正确的做法是解线性方程组：

p = -B \ g;

在 MATLAB 里，反斜杠运算符会自动选择合适的求解器，数值稳定性更好，速度也更快。对于中等规模的问题，一次求解的代价可以接受。如果你的问题规模特别大，那就是后面要说的 L-BFGS 的舞台了。

3. Armijo 回溯线搜索：让每一步都真正“下山”

有了搜索方向 p_k，下一步的问题是：沿着这个方向走多远？定的步长太长，可能一步跨过山谷甚至发散；太短，则迭代缓慢。线搜索就是解决这个问题的。

3.1 为什么需要线搜索

有人可能会问：BFGS 方向的长度不是由 B 的尺度决定吗？是的，理论上当 B 逼近真实 Hessian 时，最优步长接近 1。但在迭代初期，B 还很粗糙，如果直接按步长 1 走，可能走出一个函数值上升的坏点。

线搜索的作用就是在给定方向上找到一个合适的步长，使得函数值有足够的下降。它把“方向好不好”和“走多远”两个问题分开处理，让算法更稳健。

3.2 Armijo 条件的直观理解

Armijo 条件有时翻译为“充分下降条件”，表达式是：

f(x_k + α * p_k) ≤ f(x_k) + c * α * ∇f(x_k)ᵀ * p_k

其中 c 是一个常数，通常取 1e-4 到 1e-3。右边的项 f(x_k) + c * α * g_kᵀ p_k 是一条斜率为 c * g_kᵀ p_k 的直线。因为 p_k 是下降方向，g_kᵀ p_k < 0，所以这条直线的斜率是负的。条件的意思是：函数真实下降的量，至少要达到这条直线所代表的下降量的一定比例。

这个条件保证了步长不能太大。如果 α 太大，函数值通常无法满足这个不等式。

3.3 回溯策略：从 1 开始，不行就缩

Armijo 回溯线搜索的执行逻辑非常简单直接：

1. 初始化 α = 1。
2. 计算试探点 x_trial = x + α * p。
3. 检查 Armijo 条件是否满足。
4. 如果满足，接受这个步长；如果不满足，令 α = ρ * α，其中 ρ 通常取 0.5 到 0.8。
5. 重复步骤 2-4，直到满足条件或者 α 小于某个下限。

为什么要从 α = 1 开始？因为在 BFGS 算法收敛后期，B 已经非常接近真实 Hessian，此时最优步长就应该接近 1。从 1 开始可以最大限度地保留 BFGS 的快速收敛性。如果你一开始就给一个很小的初始步长，那 BFGS 会退化成一阶方法，收敛速度大打折扣。

3.4 参数选择的经验值

我在实验中的默认设置是：

ρ 的取值会影响回溯速度。ρ 越小，步长收缩越快，但可能错过合理步长；ρ 太接近 1，回溯次数会变多。0.5 是个非常经典的选择，绝大多数情况都不用改。

还有一个容易被忽略的点：线搜索时目标函数每次都要重新求值。如果你的目标函数很贵，最好在函数内部做一点缓存，或者在线搜索函数里把已经计算过的值传出去，避免重复计算。这个在后面代码里也会体现。

4. MATLAB 代码逐块实现

理论知识差不多了，下面是核心内容：完整的 MATLAB 实现。我把它拆成几个子函数，这样逻辑清晰，也方便你单独调试。

4.1 主函数框架

主函数接收目标函数句柄、梯度句柄、初始点和配置参数，输出最优解、最优值和迭代信息。我习惯把迭代历史也保存下来，方便后续画收敛曲线。

function [x_opt, f_opt, info] = bfgs_armijo(fun, gfun, x0, opts) % BFGS + Armijo backtracking line search for unconstrained optimization. % % Input: % fun : R^n -> R, 目标函数句柄 % gfun : R^n -> R^n, 梯度函数句柄 % x0 : n x 1 初始列向量 % opts : 结构体，可选字段包括: % tol_g 梯度停机阈值，默认 1e-6 % max_iter 最大迭代次数，默认 500 % c Armijo 条件常数，默认 1e-4 % rho 回溯缩小因子，默认 0.5 % max_backtrack 最大回溯次数，默认 30 % % Output: % x_opt : 近似最优解 % f_opt : 最优目标函数值 % info : 结构体，包含迭代次数、历史、退出标记等 % --- 处理默认参数 --- if nargin < 4 opts = struct(); end if ~isfield(opts, 'tol_g'), opts.tol_g = 1e-6; end if ~isfield(opts, 'max_iter'), opts.max_iter = 500; end if ~isfield(opts, 'c'), opts.c = 1e-4; end if ~isfield(opts, 'rho'), opts.rho = 0.5; end if ~isfield(opts, 'max_backtrack'), opts.max_backtrack = 30; end % --- 初始状态 --- x = x0(:); % 统一成列向量 n = length(x); f = fun(x); g = gfun(x); B = eye(n); % 初始 Hessian 近似取单位矩阵 f_hist = zeros(opts.max_iter, 1); g_norm_hist = zeros(opts.max_iter, 1); f_hist(1) = f; g_norm_hist(1) = norm(g, inf); exit_flag = 0; for k = 1:opts.max_iter % --- 停机判断：梯度无穷范数 --- if norm(g, inf) <= opts.tol_g exit_flag = 1; % 成功收敛 break; end % --- 计算搜索方向 --- p = -B \ g; % --- Armijo 回溯线搜索 --- alpha = armijo_backtracking(fun, x, p, g, f, opts.c, opts.rho, opts.max_backtrack); % --- 更新迭代点 --- x_new = x + alpha * p; f_new = fun(x_new); g_new = gfun(x_new); % --- 构造 BFGS 更新数据 --- s = x_new - x; y = g_new - g; % --- 更新 Hessian 近似 --- B = bfgs_update(B, s, y); % --- 更新状态 --- x = x_new; g = g_new; f = f_new; f_hist(k+1) = f; g_norm_hist(k+1) = norm(g, inf); if norm(s, inf) <= 1e-12 && norm(g, inf) <= opts.tol_g * 10 exit_flag = 2; % 位移极小，视为收敛 break; end end % --- 收敛后记录迭代次数 --- if exit_flag == 0 k = opts.max_iter; end info = struct(); info.iter = k; info.exit_flag = exit_flag; info.f_hist = f_hist(1:k+1); info.g_norm_hist = g_norm_hist(1:k+1); x_opt = x; f_opt = f; end

关于停机准则，上面我用了两个判据：一是梯度的无穷范数小于阈值，这是最经典的；二是位移很小且梯度也接近收敛，用来处理一些函数在极小点附近梯度下降缓慢的情况。这个双保险在工程里很有用。

4.2 BFGS 更新子函数

更新公式按前面写的来，但加了一个曲率条件的保护判断：

function B = bfgs_update(B, s, y) % BFGS update with curvature condition protection. % % 如果 y' * s 太小或为负，说明割线条件不稳定，此时跳过更新， % 保持原先的 Hessian 近似不变。 ys = y' * s; if ys <= 1e-12 * norm(y) * norm(s) return; end Bs = B * s; B = B + (y * y') / ys - (Bs * Bs') / (s' * Bs); end

这里有个细节要注意：Bs * Bs'是外积，维度是 n×n；s' * Bs是一个标量。MATLAB 对矩阵乘法很敏感，如果变量维度没对上，运行时会直接报错。这个在代码里我已经写清楚了。

4.3 Armijo 回溯线搜索子函数

子函数单独抽出来，好处是逻辑独立，以后想换成 Wolfe 条件或者 Goldstein 条件，动这一块就行。

function alpha = armijo_backtracking(fun, x, p, g, f, c, rho, max_backtrack) % Armijo backtracking line search. % % 从 alpha = 1 开始测试，不满足充分下降条件就按 rho 缩小。 % 如果到达最大回溯次数，返回最后一次尝试的步长。 alpha = 1.0; for i = 1:max_backtrack x_trial = x + alpha * p; f_trial = fun(x_trial); if f_trial <= f + c * alpha * (g' * p) return; end alpha = rho * alpha; end % 兜底：到达最大回溯次数但未满足条件，仍返回当前步长 % 此时调用方需要小心，可能目标函数在方向上有数值问题 end

如果线搜索一直不满足 Armijo 条件，说明可能出现两种情况：一是搜索方向不是有效的下降方向，这可能意味着 Hessian 近似出了问题；二是目标函数本身有噪声或不可导点。兜底返回一个很小的步长总比死循环强。

4.4 测试代码：Rosenbrock 函数

下面给一个完整的可运行测试脚本，里面包含目标函数、解析梯度、数值梯度验证和调用主函数。

% test_bfgs_armijo.m % Rosenbrock 函数测试 % % f(x) = 100*(x2 - x1^2)^2 + (1 - x1)^2 clear; clc; % 目标函数 fun = @(x) 100 * (x(2) - x(1).^2).^2 + (1 - x(1)).^2; % 解析梯度 gfun = @(x) [ -400 * x(1) * (x(2) - x(1).^2) - 2 * (1 - x(1)); 200 * (x(2) - x(1).^2) ]; % 数值梯度验证（可选） % g_num = numerical_gradient(fun, [-1.2; 1]); % g_ana = gfun([-1.2; 1]); % fprintf('解析梯度: [%f, %f]\n', g_ana(1), g_ana(2)); % fprintf('数值梯度: [%f, %f]\n', g_num(1), g_num(2)); % 初始点 x0 = [-1.2; 1]; % 调用 BFGS [x_opt, f_opt, info] = bfgs_armijo(fun, gfun, x0); fprintf('最优解: [%.10f, %.10f]\n', x_opt(1), x_opt(2)); fprintf('最优值: %.10e\n', f_opt); fprintf('迭代次数: %d\n', info.iter); fprintf('退出标记: %d\n', info.exit_flag); % 收敛曲线 figure; semilogy(abs(info.f_hist - f_opt) + 1e-16); xlabel('迭代次数'); ylabel('目标函数误差 (log scale)'); title('BFGS + Armijo 收敛曲线'); grid on;

4.5 数值梯度验证函数

在写优化代码时，我强烈建议先做一次梯度验证。解析梯度手写很容易出错，第 i 个分量的中心差分公式是：

g_num(i) = (f(x + h*e_i) - f(x - h*e_i)) / (2h)

MATLAB 里可以这么写：

function g = numerical_gradient(fun, x) n = length(x); g = zeros(n, 1); h = 1e-6; for i = 1:n x_plus = x; x_minus = x; x_plus(i) = x_plus(i) + h; x_minus(i) = x_minus(i) - h; g(i) = (fun(x_plus) - fun(x_minus)) / (2 * h); end end

注意：如果目标函数值很小，或者在不同尺度上差异很大，固定步长 h 的中心差分可能不够准。更稳妥的做法是用相对步长，比如h = 1e-6 * max(1, abs(x(i)))。但对于大多数测试场景，上面的版本已经够用。

5. 数值实验：Rosenbrock 函数上的收敛行为

代码写完了，跑实验验证是必须的。我的测试条件如下。

5.1 实验设置

• 目标函数：Rosenbrock
• 初始点：(-1.2, 1)
• 梯度阈值 tol_g：1e-6
• Armijo 参数 c = 1e-4，ρ = 0.5
• 最大迭代次数：500

5.2 一次典型运行的结果

在我的机器上（MATLAB R2023a，Windows 11），运行结果大致是：

最优解: [1.0000000000, 1.0000000000] 最优值: 0.0000000000e+00 迭代次数: 34 退出标记: 1

从 (-1.2, 1) 出发，经历 34 次迭代就达到了 1e-6 的梯度阈值。如果换最速下降法，同样的阈值下通常要跑到 500 次以上，有时甚至卡在狭长山谷里很久出不来。

收敛曲线画出来是平滑下降的，前期下降速度很快，后期接近最优解时函数值呈指数衰减。这正好体现了 BFGS 的拟牛顿性质：随着迭代进行，B 越来越接近真实的 Hessian，搜索方向越来越接近牛顿方向。

5.3 不同初始点的表现

我把初始点换成了几组不同位置，结果如下：

这说明在无约束光滑问题上，BFGS + Armijo 对初始点的依赖相对较小。当然，这不是绝对的。如果目标函数有很多局部极小点，或者非常非凸，初始点的选择仍然至关重要。

5.4 与最速下降法的直观对比

我做过一个简单对照实验，同样从 (-1.2, 1) 出发，同样梯度阈值，最速下降法实现里线搜索可以用精确线搜索（比如直接用黄金分割）来公平对比。结果是：

注意这里最速下降法用的是精确线搜索，每一步都找到当前方向的最优步长，仍然远远比不过 BFGS。这说明问题不在线搜索，而在于方向本身。最速下降法每次只按负梯度方向走，在条件数大的问题里本质上就有天花板。

6. 工程实践中的坑与调参建议

前面代码已经能跑通，但实际工程里不会总是这么顺利。下面这些坑是我自己在不同项目里踩过的，逐个列出来供你排查。

6.1 梯度写错：最隐蔽的坑

这是所有优化代码里最容易出问题的点。解析梯度手写一长串，特别容易少一个负号、写错一个系数。更麻烦的是，如果梯度是错的，优化器可能依然能够收敛，只是收敛到一个错误的位置。

我的建议是：跑任何优化之前，先做梯度验证。写一个脚本，随机采样几个点，对比解析梯度和数值梯度。如果相对误差在 1e-6 量级，说明解析梯度大概率没错。如果误差在两位数，那百分之百写错了。

有一个更隐蔽的情况：你的目标函数是很多项求和得到的复合函数，单点梯度验证通过了，但程序里某处修改了函数定义没同步修改梯度函数，也会出问题。所以建议每次修改目标函数后，都重新做一次梯度验证。

6.2 曲率条件被破坏导致 BFGS 更新失败

前面提到 BFGS 保持正定性的前提是 yᵀ s > 0。对凸函数来说，这个条件在 Armijo 线搜索下基本能保证。但对非凸函数，尤其是有振荡的目标函数，yᵀ s 可能变成负的。

我在代码里加了一句ys <= 1e-12 * norm(y) * norm(s)就跳过更新。让 B 保持原样，有时候等于白跑一步，但总比 B 失去正定性后方向乱掉强。如果问题反复出现，可以考虑用阻尼 BFGS，也就是对 y 做修正，强行满足曲率条件。原理上相当于在 y 中加入一部分 B*s 的分量，让 yᵀ s > 0。不过对于大多数场景，直接跳过更新已经够用。

6.3 线搜索死循环

如果目标函数有噪声或者存在数值非光滑，Armijo 条件可能很难满足。线搜索会不停缩小 α，直到逼近浮点精度下限，然后一直在极小步长附近徘徊。

遇到这种情况，我有几个措施：

• 给最大回溯次数设一个上限，比如 30 次。
• 在函数内部对候选点做边界保护，比如限制坐标范围，防止变量溢出。
• 如果发现步长小到某个阈值以下仍不满足条件，直接返回当前步长并打印警告，让外层循环的停机准则决定是否终止。

在长期运行的大程序里，打印警告很关键。如果没有警告，你会在几小时之后才发现结果不对，定位起来非常痛苦。

6.4 停机阈值设得太严格

有同学会把梯度阈值设成 1e-10 甚至 1e-12，想着“越严格越精确”。但如果目标函数是数值计算出来的，本身就带噪声，梯度在极小点附近也会抖动。阈值设得太低，优化器会在极小点附近无限循环，永远无法满足停机条件。

我的经验是：如果目标函数来自仿真或数值积分，梯度阈值设在 1e-5 到 1e-6 已经非常理想了。如果用的是高精度解析函数，再往 1e-8 以下调才有意义。另外也可以同时用位移停机作为备选，上面主循环里已经有了这种写法。

6.5 MATLAB 矩阵维度的经典坑

MATLAB 里，'是共轭转置，.'是非共轭转置。对实数矩阵来说两者等价，但如果涉及复数，用混了就出问题。另外要特别小心函数返回的梯度是行向量还是列向量。

我见过不少人写g = gradient(fun, x)返回行向量，然后直接拿去做g' * p。行向量乘列向量没问题，但p = -B \ g就会因为 B 是 n×n、g 是 1×n 而报维度错误。解决办法很简单：在进入主函数时统一转成列向量：

x = x0(:);

这个方法也适用于梯度，在调用 gfun 后强制g = gfun(x); g = g(:);或者干脆在测试脚本里就让 gfun 返回列向量。

7. 扩展：L-BFGS、非光滑场景与约束问题的思路

如果你想把这个优化器用在更大规模的实际问题上，下面几个扩展方向很值得了解。

7.1 大规模场景：从 BFGS 到 L-BFGS

BFGS 需要存储一个完整的 n×n 矩阵 B。当 n = 1000 时，B 就有 100 万个元素，内存占用约 8 MB，还能接受；当 n = 100000 时，内存占用直接飙到 80 GB，完全不可行。

L-BFGS（Limited-Memory BFGS）的基本思想是：不显式存储 B，而是存储最近 m 步的 (s, y) 对。计算搜索方向时通过两循环递归反推，达到近似 BFGS 的效果。m 通常取 5 到 20，存储量从 n² 降到 m*n，内存问题就解决了。

MATLAB 里你自己实现 L-BFGS 也不难，核心是那个两循环递归。实际工程里，处理几万个参数的机器学习模型，L-BFGS 是标配之一。

7.2 非光滑目标函数的情况

BFGS 和 Armijo 条件都依赖梯度存在。如果你的目标函数里有绝对值、L1 范数之类不可导的地方，直接套用会出问题。

工程上常用两种解决思路：

• 把目标函数中的非光滑部分替换为光滑近似，比如用 Huber 损失代替 L1 损失。
• 用近端梯度法处理非光滑项，BFGS 只负责光滑部分的二次近似。

对于简单的不可导点，也可以在每次迭代后做一个投影操作，把变量投影回可行域。

7.3 有约束问题怎么借用这套思路

对于带约束的优化问题，最常见的手段是罚函数法，把约束乘以一个很大的惩罚系数加到目标函数上，然后当作无约束问题用 BFGS 求解。这种方法实现简单，适合约束不复杂、惩罚系数能选好的场景。

如果约束比较严格，可以换用增广拉格朗日方法。外层迭代更新拉格朗日乘子，内层用 BFGS 求解子问题，效果比简单罚函数稳定得多。

我自己在实际项目中用得更多的是：无约束 + 边界约束的组合。对边界约束，可以直接在每次 BFGS 更新后用投影操作把变量拉回边界内。对于内部的光滑目标函数，BFGS + Armijo 仍然非常有效。

说实话，BFGS 这套东西数学上已经非常成熟，真正拉开差距的是你在工程实现中的细节处理。比如对曲率条件的保护、对停机准则的合理设计、对线搜索参数的调优，以及后续扩展到大规模问题的能力。把这些细节处理好了，一个自己写的优化器在很多场景下并不比商业工具箱差，而且调试起来思路清晰得多。希望这篇记录能给正在折腾非线性优化的人一些实在的参考。