10 - Quasi Newton Method

13 Jul 2014 Comments

Quasi Newton Method

前一篇文章中提到了 Classical Newton 的诸多问题，针对这些问题有了 Quasi Newton Method，其基本思想是在每步迭代计算 descent direction 时不用公式 $\b{d}^k = -(H^k)^{-1} \b{g}^k$，而是用 $\b{d}^k = -B^k \b{g}^k$，也就是找到一个矩阵 $B^k$ 去近似 $(H^k)^{-1}$，相当于在 $\b{x}^k$ 处选择下面的近似函数

利用这一近似我们就将求解 linear system ($H^k \b{d}^k = -\b{g}^k$) 的操作转变为 matrix vector multiplication ($-B^k \b{g}^k$) 操作，计算量从 $O(N^3)$ 降到了 $O(N^2)$

那如何得到一个好的对 $(H^k)^{-1}$ 的近似矩阵呢？Quasi Newton 对 $B^k$ 做了如下约束

综合上述三个约束，生成 $B^k$ 的问题可以表示成

$B^k$ 是一个 symmetric matrix，因此共包含 $\frac{n(n+1)}{2}$ 个变量，第二个对 leading principal minors 的约束对应 n 个不等式，最后一个约束对应 n 个等式，由于变量的个数多于等式和不等式的个数，所以上面的问题不止有一个解。根据对上述问题不同的解法，也就有了不同的 Quasi Newton Method，常见的包括如下 3 种

三种算法的核心区别就是如何得到 $B^k$，计算 descent direction 用的都是 $\b{d}^k = -B^k \b{g}^k$，其余包括 line search 什么的也都一样，下面分别介绍这三种算法

下面的介绍中我计算的是 $B^{k+1}$ 而不是 $B^k$，其本质没有任何区别，就是为了公式看上去能干净一些，如果用的是 $B^k$，则等号左边通常是一堆的 $k-1$ 上标，看上去有点乱

Rank One Correction

这种算法使用下面的公式求得 $B^{k+1}$

可以看到它在 $B^{k}$ 的基础上加了一个 rank 为 1 的 matrix，所以叫 rank one correction。这其中 $B^{k}$ 是已知的，我们需要确定的是 $a$ 和 $\b{u}$。

根据 secant equation

为了解最后一个 equation，Rank one correction 这么做，令 $a \b{u}^T \gamma^{k} = 1$，则有

这样计算 $B^{k+1}$ 的公式就是

下面分析一下 rank one correction

DFP Algorithm (Rank Two Correction)

DFP 是一个 rank two 的算法，最早由 Davidon 在 1959 年提出，后来 Fletcher 和 Powell 在 1963 年先后做了修改，所以算法取名为 DFP。DFP 计算 $B^{k+1}$ 的公式是

从公式可以看到，DFP 加了两个不同的 rank 为 1 的 matrix，比 rank one correction 多了一个。这里我们需要确定的变量有 4 个，分别是 $a, b, \b{u}, \b{v}$。另外，下面我就直接假设 $B^k$ 是 symmetric matrix。

根据 secant equation

显然最后一个 equation 是有很多解的，DFP 是这么解的，令

这样进一步可以求得 $a = \frac{1}{ {\delta^k}^T \gamma^k}, b = -\frac{1}{ {\gamma^k}^T B^k \gamma^k}$

因此计算 $B^{k+1}$ 的公式就是

在使用 exact line search 的情况下，如果 $B^k$ 是 symmetric positive definite matrix，则 $B^{k+1}$ 也是

BFGS Algorithm

BFGS 分别由 Broyden, Fletcher, Goldfarb, Shanno 四人于 1970 年独立提出，它也用到了 rank two correction，但不是用在 $B^{k+1}$ 上，而是用在 $(B^{k+1})^{-1}$ 上，令 $G^{k+1} = (B^{k+1})^{-1}$，则有

由于 $B^{k+1} \gamma^{k} = \delta^{k}$，因此 $G^{k+1} \delta^{k} = \gamma^{k}$，运用在 DFP 一节中给出的解题步骤，可得

可以看出，它和 DFP 中 $B^{k+1}$ 的迭代公式形式是完全一样的，不同的是 $\delta^k$ 和 $\gamma^k$ 的角色做了互换

有了 $G^{k+1}$，$B^{k+1} = (G^{k+1})^{-1}$，这里有一个 inverse operation，由于 $G^{k+1}$ 有相对简单的形式，所以对它做 inverse 直接有 closed form solution，下面给出相关的 Sherman-Morrison-Woodbury formula

Let $A$ be a nonsingular matrix. Let $\b{u}$ and $\b{v}$ be column vectors such that $1 + \b{v}^T A^{-1} \b{v} \neq 0$. Then $A + \b{u}\b{v}^T$ is nonsingular, and $$ (A + \b{u}\b{v}^T)^{-1} = A^{-1} - \frac{(A^{-1} \b{u})(\b{v}^T A^{-1})}{1 + \b{v}^T A^{-1} \b{v}}$$

连续两次在 $G^{k+1}$ 上应用这个公式可得

为了公式看上去简洁，对符号做了省略，其中 $B$ 就是 $B^k$，$\delta$ 就是 $\delta^k$，$\gamma$ 就是 $\gamma^k$

这个公式可以写成更简洁的形式

根据与 DFP 一节中给出的证明相同的证明，可以得出 $G^{k+1}$ 在 exact line search 的情况下一定是 positive definite matrix，因此 $B^{k+1}$ 也一定是 positive definite matrix。Powell 在 Some global convergence properties of a variable metric algorithm for minimization without exact line searches 这篇文章中进一步证明了对于 convex function，BFGS + Wolfe line search 可以达到 global convergence

limited-memory BFGS (lBFGS)

BFGS 虽然是个高效的算法，但其每步迭代要存储一个矩阵 $B^k$，对于参数规模较大的函数，这个空间上的消耗显然是不可接受的，也因此有了所谓的 limited-memory BFGS，其与 BFGS 的区别主要是以下两点

根据上面给出的公式

这是个递归公式，展开 $m$ 步可得 (方便起见，定义 $\rho^k = \frac{1}{ {\gamma^{k}}^T \delta^{k}}, \rho^k \in \mathbb{R}$)

最后一个式子中除了 $B^{k-m}$ 以外，其余的所有变量都能以 $\gamma$ 和 $\delta$ 表示。为了让 $B^k$ 能完全由前 $m$ 步的 $\gamma$ 和 $\delta$ 计算出来，lBFGS 在每一步迭代都选择一个矩阵去替换 $B^{k-m}$，这个矩阵每步都可以是不同的，但通常都是形式相对简单的矩阵，比如 diagonal matrix。令第 k 步选择的矩阵为 $B_0^k$，这样每步迭代中 lBFGS 计算

下面我们看看怎么实现这个公式使得它可以(几乎)完全由 vector operation 来完成

根据前面定义的 $\eta, \xi, \zeta$ 可以得出如下计算 $B^k \b{g}^k$ 的伪代码

$\eta_0 = \b{g}^k$
for $i$ from $1$ to $m$
  $\xi_i = \rho^{k-i} {\delta^{k-i}}^T \eta_{i-1}$
  $\eta_i = \eta_{i-1} - \xi_{i} \gamma^{k-i}$

$\zeta_{m+1} = B_0^k\eta_m$
for $i$ from $m$ to $1$
  $\zeta_{i} = \zeta_{i+1} + \delta^{k-i} (\xi_i - \rho^{k-i}({\gamma^{k-i}}^T\zeta_{i+1}))$

output $\zeta_1$

这是一个 $O(mn)$ 的算法。lBFGS 每轮迭代都运行上述代码得到 descent direction，并且如果 $k > m$ 还需要 update 保存 $\gamma, \delta$ 的队列，去掉最后一个，加入最新的一个，然后继续下一轮迭代，直至算法收敛

关于 $B_0^k$ 选择，一种被证明比较有效的方法是令 $B_0^k = \frac{ {\delta^{k-1}}^T \gamma^{k-1}}{ {\gamma^{k-1}}^T \gamma^{k-1}} I$

总结

一些关于 Quasi-Newton method 的理论分析这篇文章中并没有给出，下面给出一些有用的结论