首页课程笔记自动化系凸优化

# 概述

拟牛顿方法（Quasi-Newton Methods）是一类重要的优化算法，它们通过近似 Hessian 矩阵来避免牛顿方法中计算和存储完整 Hessian 矩阵的高昂成本，同时保持较快的收敛速度。

# 拟牛顿方法的动机

  1. 牛顿方法的优势与劣势

    考虑光滑凸优化问题：

    min⁡xf(x)\min_x f(x)

    其中 ff 是凸函数，二次可微，且 dom(f)=Rn\text{dom}(f) = \mathbb{R}^n。

    梯度下降更新：

    x+=x−t∇f(x)x^+ = x - t \nabla f(x)

    牛顿方法更新：

    x+=x−t(∇2f(x))−1∇f(x)x^+ = x - t(\nabla^2 f(x))^{-1} \nabla f(x)

    比较：

    • 牛顿方法具有（局部）二次收敛，而梯度下降只有线性收敛
    • 但牛顿迭代的计算成本要高得多

  2. 牛顿方法的主要计算步骤

    牛顿迭代中的两个主要步骤：

    • 计算 Hessian 矩阵 ∇2f(x)\nabla^2 f(x)
    • 求解线性方程组 ∇2f(x)s=−∇f(x)\nabla^2 f(x) s = -\nabla f(x)

    这两个步骤都可能很昂贵（特别是对于大规模问题）。

  3. 拟牛顿方法的基本思想

    拟牛顿方法重复以下形式的更新：

    x+=x+tsx^+ = x + t s

    其中方向 ss 由线性方程组定义：

    Bs=−∇f(x)B s = -\nabla f(x)

    其中 BB 是 ∇2f(x)\nabla^2 f(x) 的某个近似。我们希望：

    • BB 易于计算
    • Bs=gB s = g 易于求解

  4. 拟牛顿算法模板

    设 x(0)∈Rnx^{(0)} \in \mathbb{R}^n，B(0)≻0B^{(0)} \succ 0。对于 k=1,2,3,…k = 1, 2, 3, \ldots，重复：

    Step 1：求解 B(k−1)s(k−1)=−∇f(x(k−1))B^{(k-1)} s^{(k-1)} = -\nabla f(x^{(k-1)})

    Step 2：更新 x(k)=x(k−1)+tks(k−1)x^{(k)} = x^{(k-1)} + t_k s^{(k-1)}

    Step 3：从 B(k−1)B^{(k-1)} 计算 B(k)B^{(k)}

    不同的拟牛顿方法在 Step 3 的实现上不同。通常我们可以从 (B(k−1))−1(B^{(k-1)})^{-1} 计算 (B(k))−1(B^{(k)})^{-1}。

    基本思想：由于 B(k−1)B^{(k-1)} 已经包含关于 Hessian 的信息，使用合适的矩阵更新来形成 B(k)B^{(k)}。

  5. 割线方程（Secant Equation）

    对 B(k)B^{(k)} 的合理要求（由割线法启发）：

    ∇f(x(k))=∇f(x(k−1))+B(k)s(k−1)\nabla f(x^{(k)}) = \nabla f(x^{(k-1)}) + B^{(k)} s^{(k-1)}

    等价地，可以写成：

    ∇f(x+)=∇f(x)+B+s\nabla f(x^+) = \nabla f(x) + B^+ s

    令 y=∇f(x+)−∇f(x)y = \nabla f(x^+) - \nabla f(x)，这变为：

    B+s=yB^+ s = y

    这称为割线方程（secant equation）。

  6. 对近似矩阵的要求

    除了割线方程，我们还希望：

    • B+B^+ 是对称的
    • B+B^+ "接近" BB
    • B≻0⇒B+≻0B \succ 0 \Rightarrow B^+ \succ 0（保持正定性）

# SR1, BFGS, DFP, Broyden 类

  1. SR1 更新（对称秩一更新，Symmetric Rank-One）

    更新形式

    尝试以下形式的更新：

    B+=B+auuTB^+ = B + a u u^T

    割线方程 B+s=yB^+ s = y 给出：

    (auTs)u=y−Bs(a u^T s) u = y - B s

    这只有当 uu 是 y−Bsy - B s 的倍数时才成立。令 u=y−Bsu = y - B s，求解得到 a=1/(y−Bs)Tsa = 1/(y - B s)^T s，从而得到：

    B+=B+(y−Bs)(y−Bs)T(y−Bs)TsB^+ = B + \frac{(y - B s)(y - B s)^T}{(y - B s)^T s}

    这称为对称秩一（SR1）更新。

    逆矩阵更新

    为了求解 B+s+=−∇f(x+)B^+ s^+ = -\nabla f(x^+)，除了传播 BB 到 B+B^+，我们也传播逆矩阵，即 C=B−1C = B^{-1} 到 C+=(B+)−1C^+ = (B^+)^{-1}。

    使用 Sherman-Morrison 公式：

    (A+uvT)−1=A−1−A−1uvTA−11+vTA−1u(A + u v^T)^{-1} = A^{-1} - \frac{A^{-1} u v^T A^{-1}}{1 + v^T A^{-1} u}

    对于 SR1 更新，逆矩阵也可以轻松更新：

    C+=C+(s−Cy)(s−Cy)T(s−Cy)TyC^+ = C + \frac{(s - C y)(s - C y)^T}{(s - C y)^T y}

    特点：

    • SR1 简单且计算成本低
    • 关键缺点：它不保持正定性

  2. BFGS 更新（Broyden-Fletcher-Goldfarb-Shanno）

    秩二更新

    现在尝试秩二更新：

    B+=B+auuT+bvvTB^+ = B + a u u^T + b v v^T

    割线方程 y=B+sy = B^+ s 给出：

    y−Bs=(auTs)u+(bvTs)vy - B s = (a u^T s) u + (b v^T s) v

    令 u=yu = y，v=Bsv = B s，求解 a,ba, b 得到：

    B+=B−BssTBsTBs+yyTyTsB^+ = B - \frac{B s s^T B}{s^T B s} + \frac{y y^T}{y^T s}

    这称为Broyden-Fletcher-Goldfarb-Shanno（BFGS）更新。

    逆矩阵更新

    使用 Woodbury 公式（Sherman-Morrison 的推广）：

    (A+UCV)−1=A−1−A−1U(C−1+VA−1U)−1VA−1(A + U C V)^{-1} = A^{-1} - A^{-1} U (C^{-1} + V A^{-1} U)^{-1} V A^{-1}

    应用到我们的情况，得到逆矩阵 CC 的秩二更新：

    C+=(I−syTyTs)C(I−ysTyTs)+ssTyTsC^+ = \left(I - \frac{s y^T}{y^T s}\right) C \left(I - \frac{y s^T}{y^T s}\right) + \frac{s s^T}{y^T s}

    特点：

    • BFGS 保持正定性（如果 yTs>0y^T s > 0）
    • 在实践中表现优异
    • 具有超线性收敛性

  3. DFP 更新（Davidon-Fletcher-Powell）

    DFP 更新是 BFGS 的对偶形式，更新逆 Hessian 近似 CC：

    C+=C−CyyTCyTCy+ssTyTsC^+ = C - \frac{C y y^T C}{y^T C y} + \frac{s s^T}{y^T s}

    对应的 Hessian 近似更新为：

    B+=B+(y−Bs)(y−Bs)T(y−Bs)Ts−BssTBsTBsB^+ = B + \frac{(y - B s)(y - B s)^T}{(y - B s)^T s} - \frac{B s s^T B}{s^T B s}

    特点：

    • DFP 也保持正定性
    • 与 BFGS 类似，但更新公式不同
    • 在实践中，BFGS 通常比 DFP 表现更好

  4. Broyden 类

    Broyden 类是一族拟牛顿更新，包含 BFGS 和 DFP 作为特殊情况：

    B+=B−BssTBsTBs+yyTyTs+ϕ(sTBs)(BssTBs−yyTs)(BssTBs−yyTs)TB^+ = B - \frac{B s s^T B}{s^T B s} + \frac{y y^T}{y^T s} + \phi (s^T B s) \left(\frac{B s}{s^T B s} - \frac{y}{y^T s}\right)\left(\frac{B s}{s^T B s} - \frac{y}{y^T s}\right)^T

    其中 ϕ\phi 是参数：

    • ϕ=0\phi = 0：BFGS 更新
    • ϕ=1\phi = 1：DFP 更新
    • 其他值：Broyden 类的其他成员

  5. 收敛性分析

    超线性收敛

    对于拟牛顿方法，在适当条件下可以证明超线性收敛。关键观察是：

    当序列 {x0,x1,…,xk,…}\{x_0, x_1, \ldots, x_k, \ldots\} 线性收敛到最优点 x∗x^* 时，当 kk 很大时，∥xk+1−xk∥\|x_{k+1} - x_k\| 已经非常小，而 xk,xk+1x_k, x_{k+1} 都很靠近最优点 x∗x^*。如果 Hessian 本身是充分连续的，那么我们基本上就有：

    ∇2f(x∗)(∇f(xk+1)−∇f(xk))≈xk+1−xk\nabla^2 f(x^*)(\nabla f(x_{k+1}) - \nabla f(x_k)) \approx x_{k+1} - x_k

    也就是说，每次更新 HkH_k 时，收集到的是准确的关于 ∇2f(x∗)\nabla^2 f(x^*) 的信息。因此可以证明：

    lim⁡k→∞∥(∇2f(x∗)−Hk−1)sk∥∥sk∥=0\lim_{k \to \infty} \frac{\|(\nabla^2 f(x^*) - H_k^{-1}) s_k\|}{\|s_k\|} = 0

    这一点再加上其他条件还可以推出超线性收敛。

  6. 性能比较

    以线性规划障碍问题为例（n=100,m=500n = 100, m = 500）：

    min⁡xcTx−∑i=1mlog⁡(bi−aiTx)\min_x \quad c^T x - \sum_{i=1}^m \log(b_i - a_i^T x)

    计算复杂度：

    • 牛顿更新：O(n3)O(n^3)
    • 拟牛顿更新：O(n2)O(n^2)

    迭代次数：

    • 拟牛顿收敛所需的迭代次数少于牛顿方法的 100 倍
    • 虽然每次迭代稍慢，但总体计算时间通常更短

  7. 割线法的几何解释

    拟牛顿方法本质上是割线法的推广。

    一元情况：

    在导函数上某一点做切线，这条切线的斜率就是二次导函数：

    f′′(xk)=−f′(xk)xk+1−xkf''(x_k) = -\frac{f'(x_k)}{x_{k+1} - x_k}

    化简得到牛顿法的一元形式：xk+1=xk−f′(xk)f′′(xk)x_{k+1} = x_k - \frac{f'(x_k)}{f''(x_k)}。

    考虑导函数之前两个点所形成的割线，得到：

    Bk(xk−xk−1)=f′(xk)−f′(xk−1)B_k(x_k - x_{k-1}) = f'(x_k) - f'(x_{k-1})

    这里的 BkB_k 就是这条割线的斜率。

    多元推广：

    设 sk=xk+1−xks_k = x_{k+1} - x_k，yk=∇f(xk+1)−∇f(xk)y_k = \nabla f(x_{k+1}) - \nabla f(x_k)，式子变成：

    Bksk−1=yk−1B_k s_{k-1} = y_{k-1}

    这就是拟牛顿法的本质——割线方程。

  8. 有限内存 BFGS（L-BFGS）

    动机

    对于大规模问题，拟牛顿更新可能变得过于昂贵。基本思想是：不显式计算和存储 CC，而是通过维护所有对 (y,s)(y, s) 来计算 CC 的隐式版本。

    BFGS 更新回顾

    回忆 BFGS 通过以下方式更新 CC：

    C+=(I−syTyTs)C(I−ysTyTs)+ssTyTsC^+ = \left(I - \frac{s y^T}{y^T s}\right) C \left(I - \frac{y s^T}{y^T s}\right) + \frac{s s^T}{y^T s}

    观察这导致：

    C+g=p+(α−β)sC^+ g = p + (\alpha - \beta) s

    其中：

    α=sTgyTs,q=g−αy,p=Cq,β=yTpyTs\alpha = \frac{s^T g}{y^T s}, \quad q = g - \alpha y, \quad p = C q, \quad \beta = \frac{y^T p}{y^T s}

    L-BFGS 算法

    我们看到 C+gC^+ g 可以通过两个长度为 kk 的循环计算（如果 C+C^+ 是 kk 次迭代后逆 Hessian 的近似）：

    1. 令 q=−∇f(x(k))q = -\nabla f(x^{(k)})
    2. 对于 i=k−1,…,0i = k-1, \ldots, 0：
      (a) 计算 αi=(s(i))Tq/((y(i))Ts(i))\alpha_i = (s^{(i)})^T q / ((y^{(i)})^T s^{(i)})
      (b) 更新 q=q−αyiq = q - \alpha y^i
    3. 令 p=C(0)qp = C^{(0)} q
    4. 对于 i=0,…,k−1i = 0, \ldots, k-1：
      (a) 计算 β=(y(i))Tp/((y(i))Ts(i))\beta = (y^{(i)})^T p / ((y^{(i)})^T s^{(i)})
      (b) 更新 p=p+(αi−β)s(i)p = p + (\alpha_i - \beta) s^{(i)}
    5. 返回 pp

    有限内存版本

    有限内存 BFGS（L-BFGS）简单地将每个循环的长度限制为 mm：

    1. 令 q=−∇f(x(k))q = -\nabla f(x^{(k)})
    2. 对于 i=k−1,…,k−mi = k-1, \ldots, k-m：
      (a) 计算 αi=(s(i))Tq/((y(i))Ts(i))\alpha_i = (s^{(i)})^T q / ((y^{(i)})^T s^{(i)})
      (b) 更新 q=q−αyiq = q - \alpha y^i
    3. 令 p=Cˉ(k−m)qp = \bar{C}^{(k-m)} q
    4. 对于 i=k−m,…,k−1i = k-m, \ldots, k-1：
      (a) 计算 β=(y(i))Tp/((y(i))Ts(i))\beta = (y^{(i)})^T p / ((y^{(i)})^T s^{(i)})
      (b) 更新 p=p+(αi−β)s(i)p = p + (\alpha_i - \beta) s^{(i)}
    5. 返回 pp

    在 Step 3 中，Cˉ(k−m)\bar{C}^{(k-m)} 是我们对 C(k−m)C^{(k-m)} 的猜测（不存储）。一个流行的选择是 Cˉ(k−m)=I\bar{C}^{(k-m)} = I，也存在更复杂的选择。

    特点：

    • 存储成本：O(mn)O(mn) 而不是 O(n2)O(n^2)
    • 计算成本：O(mn)O(mn) 而不是 O(n2)O(n^2)
    • 适用于大规模问题（nn 很大）
    • 通常 mm 取 5-20 就足够