拟牛顿法

拟牛顿法（quasi-Newton method）用相邻点的梯度变化估计曲率，再据此选择优化方向。它针对可微目标，不要求每次显式计算完整的二阶导数矩阵。这里以无约束最小化 $f:{ℝ}^n\to ℝ$ 的 BFGS 更新为主：矩阵如何学习曲率，比算法名称本身更重要。

从梯度变化到割线条件

设一次移动为 $s=x^{+}-x$，相应梯度差为 $y=\mathrm{\nabla }f(x^{+})-\mathrm{\nabla }f(x)$。若目标恰为二次函数 $f(x)=\frac{1}{2}{x}^{𝖳}Qx-b^{𝖳}x$，其中 $Q$ 对称，则 $y=Qs$。一般二阶连续可微函数则满足 $$y={\int }_0^1{\mathrm{\nabla }}^{2}f(x+ts)s\mathrm{d}t.$$ 这是沿线段对梯度逐分量使用微积分基本定理得到的。一次移动只探测了矩阵在方向 $s$ 上的作用，并没有测出整个矩阵。

若用 $B$ 近似 Hessian，新的估计希望满足 $B^{+}s=y$；若存储逆矩阵估计 $H$，相应条件为 $H^{+}y=s$。这称为割线条件。对称性、正定性和保留旧信息进一步限定了更新方式。BFGS 是其中一种选择，不能由割线条件单独唯一确定。

BFGS 的两个可验证性质

设 $H$ 实对称正定，且 $s^{𝖳}y>0$。定义 $\rho =1/(s^{𝖳}y)$，逆 BFGS 公式为 $$H^{+}=(I-\rho s{y}^{𝖳})H(I-\rho y{s}^{𝖳})+\rho s{s}^{𝖳}.$$ 先验证割线条件。右侧乘以 $y$，由于 $(I-\rho y{s}^{𝖳})y=0$，第一项消失；第二项为 $\rho s(s^{𝖳}y)=s$。

再取任意非零向量 $v$，令 $w=(I-\rho y{s}^{𝖳})v$，便有 $$v^{𝖳}{H}^{+}v=w^{𝖳}Hw+\rho (s^{𝖳}v{)}^2.$$ 两项均非负。若和为零，第二项迫使 $s^{𝖳}v=0$，从而 $w=v\ne 0$，与第一项为零矛盾。所以 $H^{+}$ 仍正定。于是当梯度 $g\ne 0$ 时，方向 $p=-Hg$ 满足 $g^{𝖳}p=-g^{𝖳}Hg<0$，确为局部下降方向。这证明了更新的代数性质，尚未证明任意步长下的迭代都会收敛。

旧估计为单位阵；BFGS 更新使 H⁺y 精确等于移动 s。

例如对 $Q=\mathrm{diag}(1,4)$，一次观测 $s=(1,1{)}^{𝖳}$ 给出 $y=(1,4{)}^{𝖳}$。从 $H=I$ 出发，$\rho =1/5$，直接相乘得到 $$H^{+}=\frac{1}{25}\left(\begin{array}{cc}37& -3\\ -3& 7\end{array}\right),H^{+}(\genfrac{}{}{0ex}{}{1}{4})=(\genfrac{}{}{0ex}{}{1}{1}).$$ 图中比较 $y$ 与更新前后施加逆曲率所得的向量；长度和方向的变化对应上述乘法。

步长、曲率与一次完整求解

迭代先取 $p_k=-H_k{g}_k$，经线搜索选择 ${\alpha }_k>0$，再置 $x_{k+1}=x_k+{\alpha }_k{p}_k$。弱 Wolfe 曲率条件为 $g_{k+1}^{𝖳}{p}_k\ge c_2{g}_k^{𝖳}{p}_k$，其中 $0<c_2<1$。结合下降性可得 $$s_k^{𝖳}{y}_k={\alpha }_k(g_{k+1}-g_k{)}^{𝖳}{p}_k\ge {\alpha }_k(c_2-1)g_k^{𝖳}{p}_k>0.$$ 这说明线搜索不仅控制目标值，也为正定更新提供曲率条件。只有充分下降条件的回溯并不自动保证此式。

继续用 $f(x)=\frac{1}{2}(x_1^2+4x_2^2)$，取 $x_0=(1,1)$、$H_0=I$。梯度为 $(1,4)$，第一方向为 $(-1,-4)$。沿这个方向的精确最小步长为 $${\alpha }_0=-\frac{g_0^{𝖳}{p}_0}{p_0^{𝖳}Q{p}_0}=\frac{17}{65},x_1=(\frac{48}{65},-\frac{3}{65}).$$ 由这次实际移动形成 $s_0,y_0$ 再更新，第二次精确线搜索到达原点。下面的代码逐步输出这条路径。二维严格凸二次模型和精确算术下的结果，不能照搬为任意非线性函数都只需两步。

圆点为 (1,1)、(48/65,−3/65) 与原点，辅助曲线是目标等高线。

曲率失效与有解变式

变式：在前述单次观测中改成 $s=(1,0)$、$y=(2,0)$，仍从 $H=I$ 开始。此时 $\rho =1/2$，两侧投影因子均为 $\mathrm{diag}(0,1)$，所以 $H^{+}=\mathrm{diag}(1/2,1)$：只学习到第一坐标的逆曲率。

若 $f(x)=-x^2/2$，取 $s=1$ 便有 $y=-1$。一维割线条件要求 $H^{+}=-1$，不可能同时正定。这是非凸负曲率的具体冲突。实际算法可跳过不合格曲率对、阻尼更新或改用信赖域；这些措施改变了算法，不能在未经说明时仍声称拥有原公式的全部结论。梯度很小也可能对应极大点或鞍点。

可运行的二维实现

代码仅用 Python 3 标准库，输入为两个实数坐标；返回迭代点和目标值，不改写输入。它专用于上述正定二次函数，精确线搜索公式依赖固定矩阵 $Q$。一般稠密 BFGS 每步存储和更新矩阵约需 $O(n^2)$；函数、梯度和线搜索评价另计。高维时可用有限内存BFGS算法保存少数曲率对。

Python 3

def bfgs_quadratic(start, steps=2):
    x = [float(v) for v in start]
    if len(x) != 2 or steps < 0:
        raise ValueError("two coordinates and nonnegative steps required")
    H = [[1.0, 0.0], [0.0, 1.0]]
    out = [(tuple(x), (x[0]**2 + 4*x[1]**2)/2)]
    for _ in range(steps):
        g = [x[0], 4*x[1]]
        if sum(v*v for v in g) < 1e-24:
            break
        p = [-sum(H[i][j]*g[j] for j in range(2)) for i in range(2)]
        a = -sum(g[i]*p[i] for i in range(2))/(p[0]**2+4*p[1]**2)
        z = [x[i]+a*p[i] for i in range(2)]
        s = [z[i]-x[i] for i in range(2)]
        y = [s[0], 4*s[1]]
        sy = sum(s[i]*y[i] for i in range(2))
        if sy <= 0:
            raise ArithmeticError("nonpositive curvature")
        rho = 1/sy
        V = [[float(i == j)-rho*s[i]*y[j] for j in range(2)] for i in range(2)]
        H = [[sum(V[i][a]*H[a][b]*V[j][b] for a in range(2) for b in range(2))
              + rho*s[i]*s[j] for j in range(2)] for i in range(2)]
        x = z
        out.append((tuple(x), (x[0]**2+4*x[1]**2)/2))
    return out

for x, f in bfgs_quadratic((1, 1)):
    print(tuple(round(v, 6) for v in x), round(f, 6))
assert bfgs_quadratic((0, 0)) == [((0.0, 0.0), 0.0)]

本例输出点依次约为 (1, 1)、(0.738462, -0.046154)、(0, 0)，目标值为 2.5、0.276923、数值零。最大迭代次数是计算预算；退出时应另看梯度范数，而非把预算耗尽等同于收敛。