企业官网建设流程全解析

上一篇写完，有读者留言说：“公式都对得上，代码也能跑，怎么换个目标函数就开始乱跳，有时候loss下到一半直接卡死？”我一看就明白，这是从“跑通示例”迈向“处理真实优化问题”的必经之路。BFGS算法的核心虽然只有那行矩阵更新，但在工程里，它默认了目标函数足够光滑、步长足够讲究、梯度足够准确。这三个前提，前几篇的实现里其实只解决了一小半。

这篇是BFGS系列的第四篇，我不再推公式，专门讲怎么在Python里装备一把能上生产线的BFGS：Wolfe条件线性搜索、数值梯度陷阱、矩阵更新保护，以及一个完整的Logistic回归拟合案例。适合已经写过一版BFGS、又觉得不顺手的人。我默认你已经理解了BFGS的基本思路：用梯度差值去逼近Hessian，再用这个近似矩阵引导搜索方向。如果你前几篇是跳着看的，也没关系，下面每一步代码都能直接落地。

1. 换个函数就翻车？先看“病态”到底是怎么吃掉收敛速度的

1.1 一个条件数1000的二次函数就够让人头疼

很多新手喜欢拿Rosenbrock函数当测试基准，跑通了就很开心。但真实目标函数往往比Rosenbrock更不讲道理：有的方向平坦得像平地，有的方向陡峭得像悬崖。BFGS算法之所以强大，就是因为它试图用曲率信息去适应这种“各向异性”，可一旦实现偷工减料，病态函数的杀伤力反而比梯度下降更明显。

我拿一个最简单的例子说明：

import numpy as np def f(x): return x[0] ** 2 + 1000.0 * x[1] ** 2 def grad(x): return np.array([2.0 * x[0], 2000.0 * x[1]])

极小点在原点，但在y方向上，曲率是x方向的1000倍。用固定步长做最朴素的迭代，步长取0.0001时，x方向的下降慢得让人抓狂，y方向又可能因为步长偏大而发生振荡；步长取0.001时，y方向直接发散。任何固定步长都没法同时照顾这两个方向，这就是“病态”的本质——特征值的跨度太大。

换到BFGS上，理想情况下，B矩阵会慢慢学习到Hessian的特征值比例，把坐标系拉成“圆球”，从而一步走出很远。但这个能力完全依赖于线搜索“喂”给它的步长是否合理。如果步长随便给，B矩阵学到的曲率就是一段被污染的记录，方向越走越偏，最终表现就是“换个函数就翻车”。

1.2 第一排查动作：不碰参数，先打印日志

遇到BFGS不收敛，我的习惯不是马上调代码，而是先看日志。至少打印三样东西：loss、梯度范数、当前步长。如果条件允许，再打印方向与梯度的夹角。

for it in range(max_iter): g = grad(x) p = -np.linalg.solve(B, g) # 或者用 H 矩阵直接乘 alpha = 0.01 x_new = x + alpha * p loss_new = f(x_new) if it % 5 == 0: print(f"iter {it:3d}, loss {loss_new:.6e}, ||g|| {np.linalg.norm(g):.3e}, alpha {alpha:.4f}")

你大概率会看到一种规律：loss前几步在下降，但下降速度越来越慢，甚至后面开始上升。这时候如果第一反应是“调小学习率”，就走错方向了。BFGS没有学习率这个旋钮，它需要的是“每一步都重新计算合适步长”的线性搜索。固定步长不是BFGS的一部分，是我为了演示故意加进去的“短路操作”。

实际上，我踩过几次坑之后发现，大多数所谓“手写BFGS不收敛”的问题，根因不是更新公式写错，而是步长策略太粗糙。BFGS从设计之初就假设每一步都有一个足够好的步长，这个“好”由Wolfe条件来定义。所以下一步就必须把线性搜索这个零件装上去。

2. 给BFGS配上Wolfe条件的线性搜索，不再裸奔

2.1 BFGS不是“有方向就能走”，需要满足下降要求和适量步长

先回忆一下方向是怎么来的。假设我们维护的是Hessian逆的近似矩阵H，方向就是：

p = -H @ g

如果H保持对称正定，p一定是一个下降方向。但“下降方向”只说明在当前点附近往下走，没说走多远。步长太小，迭代速度慢；步长太大，可能越过山谷，甚至发散。更关键的是，BFGS更新矩阵时需要用相邻两步的数据：

s = x_new - x y = grad_new - grad

这个数据要有效，必须满足曲率条件 y^T s > 0。如果步长选得不好，y^T s 可能小于等于0，那更新出来的H就会失去正定性，算法直接崩坏。

所以正确做法是让线搜索同时满足Armijo条件和曲率条件：

• Armijo条件：保证这一步确实有足够下降。
• 曲率条件：保证步长没有小到让曲率信息失效。

Armijo条件写出来是：

f(x + alpha * p) <= f(x) + c1 * alpha * grad^T p

c1通常取1e-4。这个条件很宽松，只要求下降率不太离谱。然后是曲率条件：

grad(x + alpha * p)^T p >= c2 * grad^T p

c2通常取0.9，这是弱Wolfe形式。更严格一点可以用绝对值形式，但弱Wolfe对BFGS已经够用。两个条件合起来，步长既不会大到让loss反弹，也不会小到让曲率信息丢失。

2.2 从0写一个“先大步长，再指数缩小”的搜索

对于拟牛顿算法，最自然的初始步长是1。为什么？因为H被看成Hessian逆的近似，理论上一步就想走到“牛顿步”的位置，所以alpha=1是第一步试探值。如果Armijo不满足，就把alpha乘以0.5往回缩。

def line_search(f, gradf, x, p, c1=1e-4, c2=0.9, rho=0.5, max_iter=50): alpha = 1.0 fx = f(x) g0 = gradf(x) slope0 = g0 @ p # 保护：如果p不是下降方向，强制回退到最速下降方向 if slope0 >= 0: p = -g0 slope0 = g0 @ p for _ in range(max_iter): x_new = x + alpha * p f_new = f(x_new) if f_new <= fx + c1 * alpha * slope0: g_new = gradf(x_new) if g_new @ p >= c2 * slope0: return alpha, x_new, f_new, g_new alpha *= rho return alpha, x + alpha * p, f(x + alpha * p), gradf(x + alpha * p)

这段代码看起来简单，但有一个容易被忽略的点：如果p不是下降方向，直接就返回最速下降方向。这个保护在实际工程里很重要，因为某些数值误差会让H矩阵在迭代中途失去正定性，导致p变成上山方向。没有这个保护，程序会马上炸掉。

2.3 曲率条件不过，宁可不更新矩阵

线搜索返回之后，更新的门槛是 y^T s 是否为正。注意，就算线搜索尽力了，真实代码里仍然可能遇到 y^T s 很小甚至为负的情况。原因可能是梯度计算有噪声，也可能是目标函数本身不光滑。这时候不要硬着头皮做更新，直接把旧的H保留下来，继续下一步。

sy = y @ s if sy > 1e-8 * np.linalg.norm(y) * np.linalg.norm(s): rho = 1.0 / sy I = np.eye(len(s)) A = I - rho * np.outer(s, y) B = I - rho * np.outer(y, s) H = A @ H @ B + rho * np.outer(s, s) else: # 跳过更新，保留旧H pass

这个阈值看起来是随便给的，其实有讲究。单纯判断 sy > 0 会被浮点噪声骗过去，但用 y 和 s 的范数做缩放，就与量纲无关了。这个经验是高维优化里非常实用的小技巧，常规教科书里不会写。

到这里，BFGS的核心零件就齐了：方向更新、线搜索、矩阵更新保护。但还有一个隐形大坑没聊——梯度到底算得准不准。

3. 数值梯度和BFGS是天生的冤家，怎么止损

3.1 三种求导方式的误差对比

很多人的第一版代码，目标函数太复杂，手动推梯度推到怀疑人生，于是直接用有限差分。这个做法能跑，但极容易把BFGS带沟里。原因在于BFGS的核心输入是梯度的差值 y = grad_new - grad，它相当于Hessian的一个“测量样本”。测量误差直接混进样本里，后面的矩阵更新越走越歪。

三类梯度求法：

前向差分的公式是：

grad_i ≈ (f(x + h*e_i) - f(x)) / h

中心差分则同时取左右两点：

grad_i ≈ (f(x + h*e_i) - f(x - h*e_i)) / (2*h)

中心差分截断误差更小，但计算量翻倍。更麻烦的是，h选大了截断误差大，选小了舍入误差大，两股力在打架。如果你用一个精度不合适的h，梯度误差会以肉眼可见的速度破坏BFGS。

3.2 用中心差分做梯度一致性校验

我自己的习惯是：能用解析梯度就用解析梯度；如果实在推不出来，先写一版中心差分做验证，再用自动微分工具替代。在动手训练之前，一定要做一次梯度校验。

def center_diff(f, x, h=1e-6): g = np.zeros_like(x) for i in range(x.size): e = np.zeros_like(x) e[i] = 1.0 g[i] = (f(x + h * e) - f(x - h * e)) / (2.0 * h) return g # 假设已经手推了 analytic_grad x0 = np.array([1.2, -0.8, 0.5, 2.0]) g_num = center_diff(f, x0) g_ana = analytic_grad(x0) print("最大误差:", np.max(np.abs(g_num - g_ana)))

如果最大误差在1e-5到1e-7之间，说明解析梯度大概率写对了。如果误差在1e-2以上，要么h太大，要么解析梯度本身有bug。注意，用数值梯度直接参与BFGS训练也不是不行，但必须接受一个事实：高维下这个梯度噪声会被放大，最终大概率比真梯度多花几倍迭代时间，甚至收敛到不理想的位置。

3.3 为什么BFGS对梯度误差比梯度下降更敏感

梯度下降对梯度误差有些容忍度，因为只要大方向对，多走几步也能到。但BFGS是个“闭环反馈系统”：梯度误差先污染方向p，p影响步长和x_new，x_new又产生新的s和y，y直接参与更新H，H再影响下一个方向。误差会一路累积，最后H的谱结构失真，表现为收敛慢、振荡甚至发散。

尤其要注意那些带条件分支的目标函数，比如：

if x[0] > 0: return x[0] ** 2 else: return x[0] ** 3

这种函数在决策边界处不可导，中心差分可能每次踩到的点落在不同的分支上，梯度数值时好时坏。BFGS遇到这种情况，常常会在边界附近来回横跳，看起来就是“loss卡住了”。先确认目标函数的光滑性，再来考虑是不是BFGS的问题。

4. 实战：用自研BFGS拟合Logistic回归，还真的收敛了

4.1 问题定义

前段时间我在做一个“模拟项目X”，需求是用Logistic回归拟合一组分类数据。特征维度30，样本量1000，训练一个二分类模型。本来用现成工具包一分钟就能搞定，但我当时就想验证一下手写的BFGS到底能不能扛住真实的数值任务，于是把整个优化循环都自己写了一遍。

数据大概是这个样子：每个样本有30个特征，标签是+1或者-1。注意，标签不能取0和1，原因下面会讲。目标函数是带L2正则的Logistic负对数似然：

def logistic_loss(beta, X, y, lam=1e-4): z = X @ beta loss = np.mean(np.logaddexp(0, -y * z)) + 0.5 * lam * np.dot(beta, beta) grad = -np.mean(y * X / (1.0 + np.exp(y * z)), axis=0) + lam * beta return loss, grad

这段代码里有几个细节值得说：

• 用logaddexp而不是直接写np.log(1 + np.exp(...))，是为了避免exp溢出。
• 梯度公式里的 y 必须是 ±1，不能是0和1。因为我的推导里用到了 sigmoid(-y z) = 1 / (1 + exp(y z)) 这个对称关系，换成0/1标签后公式就错了。
• 正则项系数取1e-4，对X的尺度很敏感，后面我做了特征标准化之后才稳定。

4.2 组装BFGS更新循环

BFGS有两种常见实现方式：维护B矩阵（Hessian近似）或者维护H矩阵（Hessian逆近似）。用H矩阵的好处是计算方向只需要一次矩阵乘向量，不需要解线性方程组。更新公式如下：

H_new = (I - rho * s y^T) H (I - rho * y s^T) + rho * s s^T rho = 1 / (y^T s)

写成代码就是：

def bfgs_update(H, s, y): sy = y @ s if sy <= 1e-12: return H rho = 1.0 / sy I = np.eye(len(s)) A = I - rho * np.outer(s, y) B = I - rho * np.outer(y, s) return A @ H @ B + rho * np.outer(s, s)

然后整个训练循环可以写得很短：

def train_bfgs(X, y, beta0, max_iter=100): beta = beta0.copy() d = len(beta) H = np.eye(d) for it in range(max_iter): loss, g = logistic_loss(beta, X, y) p = -H @ g alpha, beta_new, loss_new, g_new = line_search( lambda b: logistic_loss(b, X, y)[0], lambda b: logistic_loss(b, X, y)[1], beta, p ) if loss_new > loss - 1e-12: break s = beta_new - beta y = g_new - g H = bfgs_update(H, s, y) beta = beta_new if it % 10 == 0: print(f"iter {it:3d}, loss {loss_new:.6f}, ||g|| {np.linalg.norm(g_new):.4e}")

这里线搜索函数就是第二节写的那个。每次计算loss和grad时都调用了一次logistic_loss，实际工程里可以顺便把g_new缓存下来，避免重复算，不过我为了可读性没有做优化。

4.3 收敛结果和一个关键翻车点

训练日志大致长这样：

loss从0.693一路掉到0.07附近，梯度范数收敛到1e-4以下，参数估计也基本稳定了。我当时挺高兴，但比较早的一版代码其实翻过一次车：loss卡在0.693完全不动，梯度范数也几乎不下降。我查了很久，最后发现是标签用了0和1，导致梯度公式里的sigmoid对称性破坏，优化过程进入了一个“假平衡”。

这是个特别典型的坑，因为很多框架里二分类标签默认为0/1，但Logistic回归的经典推导用的是±1。以后你写这类代码，第一件事就要确认目标函数和标签约定一致，不然后面怎么调都白费。

4.4 特征缩放对BFGS的影响

这个实验中我还做了一次对比：X没有标准化的时候，loss下降明显慢很多，经常要走100步以上；对每列特征做零均值单位方差之后，50步就能达到同样的精度。原因仍然是条件数。

如果某个特征取值在0.01量级，另一个特征在10000量级，H矩阵初始状态下很难兼顾两者的尺度，导致前几次迭代步长波动剧烈，严重拖慢收敛。标准化之后，每个坐标的尺度接近，BFGS的H矩阵只需要拟合“形状”而不是同时拟合“单位”，省事得多。

所以，判断一个优化问题是不是BFGS“不行”，先检查特征缩放。特征都标准化了还不行，再回头查梯度。

5. B/H矩阵更新的工程保护：这层护甲必须穿上

5.1 sy > 0 只是最低门槛，量纲归一化更稳

这个前面提过，但我还是要单独说一次。更新保护的条件不能只写 sy > 0，因为当变量尺度很小或者梯度数值很大时，浮点误差会让一个本应该为负的 sy 变成微小的正数，更新后H矩阵直接污染。我建议统一的写法：

sy = y @ s if sy > 1e-8 * np.linalg.norm(y) * np.linalg.norm(s): # 正常更新 else: # 跳过更新

这个阈值选1e-8是比较通用的，如果问题特别病态，可以放大到1e-6。哪怕跳过了几次更新，也只是暂时退化成梯度下降，不会导致算法崩溃。

5.2 周期性检查H矩阵的对称性和特征值

理论上的H矩阵是严格对称正定的，但浮点运算会慢慢累积误差。迭代几百上千步之后，H可能变成一个“轻度不对称”的矩阵。这种说法其实不太严谨，因为每次更新都做了对称操作，但数值误差累计确实会让对称部分被微小破坏。

我习惯每隔一个周期做个检查：

def check_repair(H): err = np.linalg.norm(H - H.T, ord='fro') / np.linalg.norm(H, ord='fro') if err > 1e-8: H = (H + H.T) / 2.0 min_eig = np.linalg.eigvalsh(H)[0] if min_eig < 0: # 出现负特征值，说明前面的线搜索或梯度校验有漏洞 print("warning: H is not positive definite") H = H + (1e-8 - min_eig) * np.eye(H.shape[0]) return H

注意，负特征值本身就是一个诊断信号。如果你看到这个warning，不要只是“修一下矩阵”就继续跑，必须回头检查线搜索有没有真的满足Wolfe条件、梯度有没有噪声、目标函数是否光滑。矩阵修好了，病根还在，后面会反复犯。

5.3 初始H矩阵的尺度：别傻乎乎用单位矩阵

BFGS的初始H矩阵一般取单位矩阵，但这对真实问题来说太“天真”了。更实用的做法是先做一步最速下降试探，然后根据第一组s和y估算一个初始尺度。

我常用的计算公式是：

H = np.eye(d) # 第一次迭代后，用一维信息近似设置初始Hessian逆的尺度 if first_update_done: H0 = (y @ y) / (y @ s) * np.eye(d)

这个思想来自Barzilai-Borwein步长：在最优步长设为1的前提下，反推H的尺度应该约等于 y^T y / y^T s。用这个尺度初始化H，头几步的搜索质量会明显好过拍脑袋取单位矩阵。你不要指望第一组数据能给出精确的Hessian，但至少它能给出“整体大小”的量级。

5.4 高维度怎么办：L-BFGS是一句话的事，做起来是另一篇的事

当你把上面的代码都跑通之后，一定会遇到新问题：维度上去了，H矩阵存不下来。BFGS矩阵的存储和矩阵乘法都是O(d^2)代价。d=1000时，8字节的双精度数组就要8MB左右，这还勉强能忍；d=10000时就是800MB，直接让人崩溃。

应对方案就是L-BFGS：不显式维护H矩阵，只保存最近m组s和y，每次要算方向的时候用two-loop recursion还原出“H @ g”的效果。存储从O(d^2)降为O(md)，而且还可以处理非平滑目标的二阶信息。这块内容展开又是一篇的体量，我大概率会在系列第五篇里写。

不过从我的实操体会来说，如果你把这一篇讲的线性搜索、梯度和矩阵保护都做好了，绝大多数中低维优化问题已经不需要再折腾了。真正的瓶颈往往不在算法本身，而在数据尺度、步长策略和梯度质量这三件看似“不核心”的事情上。

真要说经验，我的习惯是：接到一个非凸、高维、带约束的优化需求时，先别急着秀BFGS。先做三件事：检查目标函数有没有解析梯度，用中心差分做一次校验；给线搜索加装Wolfe条件而不是固定步长；最后给H矩阵加更新保护。这三条都做到了，手写的BFGS基本能对标库实现。剩下那些收敛不稳定的情况，多半不是BFGS的问题，而是特征没缩放。