拟牛顿法

拟牛顿法的基本思想是用不包含二阶导数的矩阵近似牛顿法中的 Hessian 矩阵的逆矩阵，从而避免计算二阶导。拟牛顿法具有二次终止性，且对于一般情形具有 n 步二级收敛速率。缺点是所需存储量较大。是求解无约束最优化问题最有效的一类方法。

一、拟牛顿条件

  牛顿法的迭代公式为：

\[\boldsymbol{x}^{(k+1)}=\boldsymbol{x}^{(k)}+\lambda_k\boldsymbol{d}^{(k)} \]

其中， \(\lambda_k\) 是沿着牛顿方向搜索的步长， \(\boldsymbol{d}^{(k)}\) 是在点 \(\boldsymbol{x}^{(k)}\) 处的牛顿方向：

\[\boldsymbol{d}^{(k)}=-\nabla^2f(\boldsymbol{x}^{(k)})^{-1}\nabla f(\boldsymbol{x}^{(k)}) \]

  为构造 \(\nabla^2f(\boldsymbol{x}^{(k)})^{-1}\) 的近似矩阵 \(\bold{H}_k\) ，先分析 \(\nabla^2f(\boldsymbol{x}^{(k)})^{-1}\) 与一阶导数的关系。假设在第 \(k\) 次迭代后得到点 \(\boldsymbol{x}^{(k+1)}\) ，在该点处将目标函数进行二阶 Taylor 近似：

\[\begin{aligned} f(\boldsymbol{x}) \approx & f\left({\boldsymbol{x}}^{(k+1)}\right)+\nabla f\left({\boldsymbol{x}}^{(k+1)}\right)^{\mathrm{T}}\left({\boldsymbol{x}}-{\boldsymbol{x}}^{(k+1)}\right) \\ &+\frac{1}{2}\left({\boldsymbol{x}}-{\boldsymbol{x}}^{(k+1)}\right)^{\mathrm{T}} \nabla^{2} f\left({\boldsymbol{x}}^{(k+1)}\right)\left({\boldsymbol{x}}-{\boldsymbol{x}}^{(k+1)}\right) \end{aligned} \]

因此在 \(\boldsymbol{x}^{(k+1)}\) 附近有（移项然后按求导定义变换）：

\[\nabla f(\boldsymbol{x}) \approx \nabla f\left(\boldsymbol{x}^{(k+1)}\right)+\nabla^{2} f\left(\boldsymbol{x}^{(k+1)}\right)\left(\boldsymbol{x}-\boldsymbol{x}^{(k+1)}\right) \]

令 \(\boldsymbol{x}=\boldsymbol{x}^{(k)}\) ，则

\[\nabla f\left({\boldsymbol{x}}^{(k)}\right) \approx \nabla f\left({\boldsymbol{x}}^{(k+1)}\right)+\nabla^{2} f\left({\boldsymbol{x}}^{(k+1)}\right)\left({\boldsymbol{x}}^{(k)}-{\boldsymbol{x}}^{(k+1)}\right) \]

记：

\[\begin{aligned} &\boldsymbol{p}^{(k)}={\boldsymbol{x}}^{(k+1)}-{\boldsymbol{x}}^{(k)} \\ &\boldsymbol{q}^{(k)}=\nabla f\left({\boldsymbol{x}}^{(k+1)}\right)-\nabla f\left({\boldsymbol{x}}^{(k)}\right) \end{aligned} \]

则有：

\[\boldsymbol{q}^{(k)} \approx \nabla^{2} f\left(\boldsymbol{x}^{(k+1)}\right) \boldsymbol{p}^{(k)} \]

假设 Hessian 矩阵 $\nabla^{2} f\left(\boldsymbol{x}^{(k+1)}\right) $ 可逆，则有（Hessian矩阵的逆与一阶导的关系式）：

\[\boldsymbol{p}^{(k)} \approx \nabla^{2} f\left({\boldsymbol{x}}^{(k+1)}\right)^{-1} \boldsymbol{q}^{(k)} \]

因此当计算出 \(\boldsymbol{p}^{(k)}\) 和 \(\boldsymbol{q}^{(k)}\) 后，就可以根据上式估计出点 \(\boldsymbol{x}^{(k+1)}\) 处的 Hessian 矩阵的逆，因此可用不包含二阶导数的矩阵 \(\bold{H}_{k+1}\) 替代 \(\nabla^{2} f\left({\boldsymbol{x}}^{(k+1)}\right)^{-1}\) ，需要满足：

\[\boldsymbol{p}^{(k)}=\bold{H}_{k+1}\boldsymbol{q}^{(k)} \]

称之为拟牛顿条件。

二、秩 1 校正

原理

  当 \(\nabla^{2} f\left({\boldsymbol{x}}^{(k)}\right)^{-1}\) 是 \(n\) 阶对称正定矩阵时，满足拟牛顿条件的矩阵 \(\bold{H}_{k}\) 也应该是 \(n\) 阶对称正定矩阵。这种近似矩阵的一般构造策略是：初始的 \(\bold{H}_{k}\) （\(k=1\)）取任意一个 \(n\) 阶对称正定矩阵（一般取单位矩阵 \(I\) ），然后通过修正 \(\bold{H}_{k}\) 给出 \(\bold{H}_{k+1}\) ，令

\[\bold{H}_{k+1}=\bold{H}_{k}+\Delta\bold{H}_{k} \]

其中，\(\Delta\bold{H}_{k}\) 称为校正矩阵。

  确定 \(\Delta\bold{H}_{k}\) 的方法之一是令

\[\Delta\bold{H}_{k}=\alpha_{k} \boldsymbol{z}^{(k)}\left(\boldsymbol{z}^{(k)}\right)^{\mathrm{T}} \]

其中，\(\alpha_k\) 是一个常数，\(z^{(k)}\) 是 \(n\) 维列向量。这样得到的 \(\Delta\bold{H}_{k}\) 是一个秩为1的对称矩阵，\(\boldsymbol{z}^{(k)}\) 的选择应使得拟牛顿条件得到满足，因此

\[\boldsymbol{p}^{(k)}=\bold{H}_{k} \boldsymbol{q}^{(k)}+\alpha_{k} \boldsymbol{z}^{(k)} \boldsymbol{z}^{(k) \mathrm{T}} \boldsymbol{q}^{(k)} \]

等式两端同乘 \(\boldsymbol{q}^{(k)\mathrm{T}}\) ，整理可得到

\[\boldsymbol{q}^{(k) \mathrm{T}}\left(\boldsymbol{p}^{(k)}-\bold{H}_{k} \boldsymbol{q}^{(k)}\right)=\alpha_{k}\left(\boldsymbol{z}^{(k) \mathrm{T}} \boldsymbol{q}^{(k)}\right)^{2} \]

另外

\[\boldsymbol{z}^{(k)}=\frac{\boldsymbol{p}^{(k)}-\bold{H}_{k} \boldsymbol{q}^{(k)}}{\alpha_{k} \boldsymbol{z}^{(k) \mathrm{T}} \boldsymbol{q}^{(k)}} \]

从而可得到 秩1校正公式：

\[\bold{H}_{k+1}=\bold{H}_{k}+\frac{\left(\boldsymbol{p}^{(k)}-\bold{H}_{k} \boldsymbol{q}^{(k)}\right)\left(\boldsymbol{p}^{(k)}-\bold{H}_{k} \boldsymbol{q}^{(k)}\right)^{\mathrm{T}}}{\boldsymbol{q}^{(k) T}\left(\boldsymbol{p}^{(k)}-\bold{H}_{k} \boldsymbol{q}^{(k)}\right)} \]

用法

在利用该公式进行极小化目标函数 \(f(x)\) 时，第 \(k\) 次迭代中令搜索方向为

\[\boldsymbol{d}^{(k)}=-\bold{H}_{k} \nabla f\left(\boldsymbol{x}^{(k)}\right) \]

然后沿着该搜索方向搜索，计算步长 \(\lambda_k\) ，从而确定后继点 \(x^{(k+1)}\) 。进而计算后继点 \(x^{(k+1)}\) 处的梯度 \(\nabla f(\boldsymbol{x}^{(k+1)})\) 以及 \(\boldsymbol{p}^{(k)}\) 和 \(\boldsymbol{q}^{(k)}\) ，再利用 秩1校正公式 计算 \(\bold{H}_{k+1}\) ，接着继续计算搜索方向 \(\boldsymbol{d}^{(k+1)}\)。以此类推直至满足算法的终止条件（最大迭代次数或精度要求）。

优缺点

• 优点

  1. 不需要计算 Hessian 矩阵并求逆；
  2. 具有二次终止性，收敛性好。

• 缺点

  1. 需要满足一定的条件；

     仅当 \(\boldsymbol{q}^{(k) \mathrm{T}}\left(\boldsymbol{p}^{(k)}-\bold{H}_{k} \boldsymbol{q}^{(k)}\right)>0\) ，所求 \(\bold{H}_{k+1}\) 才具有正定性。

  2. 并且上一点是无法保证的，即使满足条件，也可能由于舍入误差的影响，导致 \(\Delta\bold{H}_k\) 无界，数值计算困难。

三、DFP算法（变尺度法）

校正矩阵

\[\Delta \bold{H}_{k}=\frac{\boldsymbol{p}^{(k)} \boldsymbol{p}^{(k) T}}{\boldsymbol{p}^{(k) T} \boldsymbol{q}^{(k)}}-\frac{\bold{H}_{k} \boldsymbol{q}^{(k)} \boldsymbol{q}^{(k) \mathrm{T}} \bold{H}_{k}}{\boldsymbol{q}^{(k) T} \bold{H}_{k} \boldsymbol{q}^{(k)}} \]

DFP公式

这样定义的校正矩阵能够使得 \(\bold{H}_{k+1}\) 满足拟牛顿条件，记为 DFP 公式，表示为

\[\bold{H}_{k+1}=\bold{H}_{k}+\frac{\boldsymbol{p}^{(k)} \boldsymbol{p}^{(k) T}}{\boldsymbol{p}^{(k) T} \boldsymbol{q}^{(k)}}-\frac{\bold{H}_{k} \boldsymbol{q}^{(k)} \boldsymbol{q}^{(k) \mathrm{T}} \bold{H}_{k}}{\boldsymbol{q}^{(k) T} \bold{H}_{k} \boldsymbol{q}^{(k)}} \]

其计算步骤与秩1校正公式近似：

1. 选取初始点 \(x^{(1)}\) ，给定终止精度要求 \(\varepsilon>0\) 。

2. 置 \(\bold{H}_{1}=\bold{I}_{n}\) （单位矩阵），计算在点 \(x^{(1)}\) 处的梯度

   \[\boldsymbol{g}_1=\nabla f(\boldsymbol{x}^{(1)}) \]

   置 \(k=1\) 。

3. 搜索方向 \(\boldsymbol d^{(k)}=-\bold H_{k} \boldsymbol g_{k}\) 。

4. 从 \(\boldsymbol x^{(k)}\) 出发，沿着 \(\boldsymbol d^{(k)}\) 方向搜索，求出步长 \(\lambda_k\) ，使得 \(\lambda_k=\arg\min _\limits{\lambda \geq 0} f\left(\boldsymbol x^{(k)}+\lambda \boldsymbol d^{(k)}\right)\) 。然后迭代搜索点 \(\boldsymbol x^{(k+1)}=\boldsymbol x^{(k)}+\lambda_k\boldsymbol d^{(k)}\) 。

5. 检验是否达到终止条件（\(\left\|\nabla f\left(\boldsymbol{x}^{(k+1)}\right)\right\| \leqslant \varepsilon\)），若达到则终止迭代；否则进入第6步。

6. 若 \(k=n\) ，则令 \(\boldsymbol x^{(1)}=\boldsymbol x^{(k+1)}\) ，返回步骤2；否则进入步骤7。

7. 令 \(\boldsymbol g_{(k+1)}=\nabla f(\boldsymbol{x}^{(k+1)}),\boldsymbol p^{(k)}=\boldsymbol x^{(k+1)}-\boldsymbol x^{(k)}, \boldsymbol q^{(k)}=\boldsymbol g_{k+1}-\boldsymbol g_{k}\) 。利用 DFP 公式计算 \(\bold{H}_{k+1}\) ，置 \(k=k+1\) ，返回步骤3。

优缺点

• 用 DFP 方法构造出来的 \(\bold{H}_{k}\) 均为对称正定矩阵。
• 具有二次终止性，且对于一般情形具有 n 步二级收敛速率。

四、BFGS公式及Broyden族

另一种形式的拟牛顿条件（DFP公式的对偶公式）。

参考《最优化理论与算法》P314页。

参考

• 《最优化理论与算法》陈宝林著。