Limit-memory BFGS Method

Slide: Starting with the classical Newton's method, we will examine its limitations and the improvements introduced by quasi-Newton methods. Finally, we will delve into the L-BFGS algorithm, which is particularly suited for large-scale problems. Inspired by Dr. Zhepei Wang's Lecture "Numerical Optimization for Robotics".

Introduction

Unconstrained Optimization

Consider a smooth and twice-differentiable unconstrained optimization problem:

$$\underset{\mathbf{x}}{min}f(\mathbf{x})$$

Descent methods provide an iterative solution:

$${\mathbf{x}}^{k+1}={\mathbf{x}}^k+{\alpha }^k\cdot {\mathbf{d}}^k$$

where ${\mathbf{d}}^k$ is the direction, and ${\alpha }^k$ is the step size.

Newton's Method

By second-order Taylor expansion,

$$f(\mathbf{x})\approx f({\mathbf{x}}^k)+\mathrm{\nabla }f({\mathbf{x}}^k{)}^{\mathrm{\top }}(\mathbf{x}-{\mathbf{x}}^k)+\frac{1}{2}(\mathbf{x}-{\mathbf{x}}^k{)}^{\mathrm{\top }}{\mathrm{\nabla }}^{2}f({\mathbf{x}}^k)(\mathbf{x}-{\mathbf{x}}^k)$$

Minimizing quadratic approximation

$${\mathrm{\nabla }}^{2}f({\mathbf{x}}^k)(\mathbf{x}-{\mathbf{x}}^k)+\mathrm{\nabla }f({\mathbf{x}}^k)=0$$

For ${\mathrm{\nabla }}^{2}f({\mathbf{x}}^k)\succ 0$

$${\mathbf{x}}^{k+1}={\mathbf{x}}^k-[{\mathrm{\nabla }}^{2}f({\mathbf{x}}^k){]}^{-1}\mathrm{\nabla }f({\mathbf{x}}^k)$$

Courtesy: Ardian Umam

Damped Newton Method

For ${\mathrm{\nabla }}^{2}f({\mathbf{x}}^k)\nsucc 0$, the direction ${\mathbf{d}}^k$ cannot be directly solved from ${\mathrm{\nabla }}^{2}f({\mathbf{x}}^k){\mathbf{d}}^k=-\mathrm{\nabla }f({\mathbf{x}}^k)$. In such cases, a PD matrix ${\mathbf{M}}^k$ must be constructed to approximate the Hessian.

If the function is convex, ${\mathrm{\nabla }}^{2}f({\mathbf{x}}^k)$ may be singular. Adding a regularization term ensures positive definiteness:

$${\mathbf{M}}^k={\mathrm{\nabla }}^{2}f({\mathbf{x}}^k)+\lambda \mathbf{I}$$

$\lambda >0$ starts small and grows until Cholesky decomposition works.

If the function is nonconvex, ${\mathrm{\nabla }}^{2}f({\mathbf{x}}^k)$ may be indefinite. To handle this, the Bunch-Kaufman decomposition is applied to obtain $\mathbf{L}\mathbf{D}{\mathbf{L}}^{\mathrm{\top }}$ and $\tilde{\mathbf{D}}$.

$${\mathbf{M}}^k=\mathbf{L}\tilde{\mathbf{D}}{\mathbf{L}}^{\mathrm{\top }}$$

Finally, direction is solved from ${\mathbf{M}}^k{\mathbf{d}}^k=-\mathrm{\nabla }f({\mathbf{x}}^k)$.

Practical Newton Method

Moreover, we can select ${\alpha }^k$ by backtracking line search to ensure sufficient decrease in the objective function, satisfying the Armijo condition:

$$f({\mathbf{x}}^k+{\alpha }^k{\mathbf{d}}^k)\le f({\mathbf{x}}^k)+c_1\cdot {\alpha }^k\mathrm{\nabla }f({\mathbf{x}}^k{)}^{\mathrm{\top }}{\mathbf{d}}^k$$

where $c_1\in (0,1)$ is a small constant.

Courtesy: Cornell University

Quasi-Newton Methods

Newton's Method: Limitations

• High Cost: Computing the Hessian and its inverse requires $\mathcal{O}(n^3)$ operations, impractical for large problems.
• Indefinite Hessian: In nonconvex cases, the Hessian may lead to steps toward saddle points.
• Ill-Conditioning: Poorly conditioned Hessians amplify errors and hinder convergence.
• Inaccurate Model: Local quadratic approximations may fail for complex functions, causing inefficiency or divergence.

Quasi-Newton Approximation

Newton approximation:

$$f(\mathbf{x})\approx f({\mathbf{x}}^k)+(\mathbf{x}-{\mathbf{x}}^k{)}^{\mathrm{\top }}{\mathbf{g}}^k+\frac{1}{2}(\mathbf{x}-{\mathbf{x}}^k{)}^{\mathrm{\top }}{\mathbf{H}}^k(\mathbf{x}-{\mathbf{x}}^k)$$$${\mathbf{H}}^k{\mathbf{d}}^k=-{\mathbf{g}}^k$$

Quasi-Newton approximation:

$$f(\mathbf{x})\approx f({\mathbf{x}}^k)+(\mathbf{x}-{\mathbf{x}}^k{)}^{\mathrm{\top }}{\mathbf{g}}^k+\frac{1}{2}(\mathbf{x}-{\mathbf{x}}^k{)}^{\mathrm{\top }}{\mathbf{B}}^k(\mathbf{x}-{\mathbf{x}}^k)$$$${\mathbf{B}}^k{\mathbf{d}}^k=-{\mathbf{g}}^k$$

The matrix ${\mathbf{B}}^k$ should:

• Avoid full second-order derivatives.
• Have a closed-form solution for linear equations.
• Retain first-order curvature information.
• Preserve the descent direction.

BFGS Method

Descent Direction

Search direction $d^k$ should make acute angle with negative gradient.

$$({\mathbf{g}}^k{)}^{\mathrm{\top }}{\mathbf{d}}^k=-({\mathbf{g}}^k{)}^{\mathrm{\top }}({\mathbf{B}}^k{)}^{-1}{\mathbf{g}}^k<0$$

${\mathbf{B}}^k$ must be PD since for all non-negative ${\mathbf{g}}^k$, $({\mathbf{g}}^k{)}^{\mathrm{\top }}({\mathbf{B}}^k{)}^{-1}{\mathbf{g}}^k>0$.

Courtesy: Active Calculus

Curvature Information

At the point ${\mathbf{x}}^{k+1}$, the gradient is ${\mathbf{g}}^{k+1}$. We want ${\mathbf{B}}^{k+1}$ to satisfy:

$${\mathbf{B}}^{k+1}({\mathbf{x}}^{k+1}-{\mathbf{x}}^k)\approx {\mathbf{g}}^{k+1}-{\mathbf{g}}^k$$$${\mathbf{B}}^{k+1}{\mathbf{s}}^k={\mathbf{y}}^k$$

The Optimal ${\mathbf{B}}^{k+1}$?

Infinitely many ${\mathbf{B}}^{k+1}$ satisfy the secant condition, how to choose the best one. We define the following weighted least square problem

$$\underset{\mathbf{B}}{min}\Vert \mathbf{B}-{\mathbf{B}}^k{\Vert }_{\mathbf{W}}^2\text{subject to}\mathbf{B}={\mathbf{B}}^{\mathrm{\top }},\mathbf{B}{\mathbf{s}}^k={\mathbf{y}}^k$$

In BFGS, weight matrix is select to be:

$$\mathbf{W}={\int }_0^1{\mathrm{\nabla }}^{2}f[(1-\tau ){\mathbf{x}}^k+\tau {\mathbf{x}}^{k+1}]\mathrm{d}\tau$$

Closed-form Solution for ${\mathbf{B}}^{k+1}$

To derive the optimal ${\mathbf{B}}^{k+1}$, we construct the Lagrangian function:

$$\mathcal{L}(\mathbf{B},\mathbf{\Lambda })=\frac{1}{2}\Vert \mathbf{B}-{\mathbf{B}}^k{\Vert }_{\mathbf{W}}^2+\text{tr}\left[{\mathbf{\Lambda }}^{\mathrm{\top }}(\mathbf{B}{\mathbf{s}}^k-{\mathbf{y}}^k)\right]$$

Taking the derivative of the Lagrangian with respect to $\mathbf{B}$ and setting it to zero gives:

$$\frac{\partial \mathcal{L}}{\partial \mathbf{B}}=\mathbf{W}(\mathbf{B}-{\mathbf{B}}^k)\mathbf{W}+\mathbf{\Lambda }({\mathbf{s}}^k{)}^{\mathrm{\top }}=0$$

Rearranging the terms, we express $\mathbf{B}$ as:

$$\mathbf{B}={\mathbf{B}}^k-{\mathbf{W}}^{-1}\mathbf{\Lambda }({\mathbf{s}}^k{)}^{\mathrm{\top }}{\mathbf{W}}^{-1}$$

Substituting this expression into the secant condition $\mathbf{B}{\mathbf{s}}^k={\mathbf{y}}^k$, we obtain:

$$({\mathbf{B}}^k-{\mathbf{W}}^{-1}\mathbf{\Lambda }({\mathbf{s}}^k{)}^{\mathrm{\top }}){\mathbf{s}}^k={\mathbf{y}}^k$$

Solving for $\mathbf{\Lambda }$, we find:

$$\mathbf{\Lambda }=\mathbf{W}({\mathbf{y}}^k-{\mathbf{B}}^k{\mathbf{s}}^k){(({\mathbf{s}}^k{)}^{\mathrm{\top }}{\mathbf{W}}^{-1}{\mathbf{s}}^k)}^{-1}$$

Finally, substituting $\mathbf{\Lambda }$ back, the closed-form solution for ${\mathbf{B}}^{k+1}$ is:

$${\mathbf{B}}^{k+1}={\mathbf{B}}^k+\frac{{\mathbf{y}}^k({\mathbf{y}}^k{)}^{\mathrm{\top }}}{({\mathbf{s}}^k{)}^{\mathrm{\top }}{\mathbf{y}}^k}-\frac{{\mathbf{B}}^k{\mathbf{s}}^k({\mathbf{B}}^k{\mathbf{s}}^k{)}^{\mathrm{\top }}}{({\mathbf{s}}^k{)}^{\mathrm{\top }}{\mathbf{B}}^k{\mathbf{s}}^k}$$

Broyden-Fletcher-Goldfarb-Shanno (BFGS)

Given the initial value ${\mathbf{B}}^0=\mathbf{I}$, the updates are performed iteratively:

$${\mathbf{B}}^{k+1}={\mathbf{B}}^k+\frac{{\mathbf{y}}^k({\mathbf{y}}^k{)}^{\mathrm{\top }}}{({\mathbf{s}}^k{)}^{\mathrm{\top }}{\mathbf{y}}^k}-\frac{{\mathbf{B}}^k{\mathbf{s}}^k({\mathbf{B}}^k{\mathbf{s}}^k{)}^{\mathrm{\top }}}{({\mathbf{s}}^k{)}^{\mathrm{\top }}{\mathbf{B}}^k{\mathbf{s}}^k}$$

with ${\mathbf{s}}^k={\mathbf{x}}^{k+1}-{\mathbf{x}}^k$ and ${\mathbf{y}}^k={\mathbf{g}}^{k+1}-{\mathbf{g}}^k$.

For computational efficiency, we often work with the inverse of ${\mathbf{B}}^k$ directly. The iterative update for $({\mathbf{B}}^k{)}^{-1}$ is given by:

$${\mathbf{C}}^{k+1}=(I-\frac{{\mathbf{s}}^k({\mathbf{y}}^k{)}^{\mathrm{\top }}}{({\mathbf{s}}^k{)}^{\mathrm{\top }}{\mathbf{y}}^k}){\mathbf{C}}^k(I-\frac{{\mathbf{y}}^k({\mathbf{s}}^k{)}^{\mathrm{\top }}}{({\mathbf{s}}^k{)}^{\mathrm{\top }}{\mathbf{y}}^k})+\frac{{\mathbf{s}}^k({\mathbf{s}}^k{)}^{\mathrm{\top }}}{({\mathbf{s}}^k{)}^{\mathrm{\top }}{\mathbf{y}}^k}$$

Guaranteeing PD of ${\mathbf{B}}^{k+1}$

To ensure that ${\mathbf{B}}^{k+1}$ remains positive definite (PD), the following curvature condition must hold:

$$({\mathbf{y}}^k{)}^{\mathrm{\top }}{\mathbf{s}}^k>0$$

For any nonzero vector $\mathbf{z}$, using the Cauchy-Schwarz inequality:

$$\begin{array}{rl}{\mathbf{z}}^{\mathrm{\top }}{\mathbf{B}}^{k+1}\mathbf{z}& ={\mathbf{z}}^{\mathrm{\top }}{\mathbf{B}}^k\mathbf{z}+\frac{({\mathbf{z}}^{\mathrm{\top }}{\mathbf{y}}^k{)}^2}{({\mathbf{y}}^k{)}^{\mathrm{\top }}{\mathbf{s}}^k}-\frac{({\mathbf{z}}^{\mathrm{\top }}{\mathbf{B}}^k{\mathbf{s}}^k{)}^2}{({\mathbf{s}}^k{)}^{\mathrm{\top }}{\mathbf{B}}^k{\mathbf{s}}^k}\\ & \ge \frac{{\mathbf{z}}^{\mathrm{\top }}{\mathbf{B}}^k\mathbf{z}({\mathbf{s}}^k{)}^{\mathrm{\top }}{\mathbf{B}}^k{\mathbf{s}}^k-({\mathbf{z}}^{\mathrm{\top }}{\mathbf{B}}^k{\mathbf{s}}^k{)}^2}{({\mathbf{s}}^k{)}^{\mathrm{\top }}{\mathbf{B}}^k{\mathbf{s}}^k}\ge 0\end{array}$$

Equalities hold only when ${\mathbf{z}}^{\mathrm{\top }}{\mathbf{y}}^k=0$ and $\mathbf{z}\parallel {\mathbf{s}}^k$. Given that $({\mathbf{y}}^k{)}^{\mathrm{\top }}{\mathbf{s}}^k>0$, these conditions cannot hold simultaneously. Therefore, if ${\mathbf{B}}^k\succ 0$, it follows that ${\mathbf{B}}^{k+1}\succ 0$.

Guranteeing $({\mathbf{y}}^k{)}^{\mathrm{\top }}{\mathbf{s}}^k>0$

Armijo Condition (AC) cannot gurantee curvature, we need curvature condition (CC):

$$({\mathbf{d}}^k{)}^{\mathrm{\top }}\mathrm{\nabla }f({\mathbf{x}}^k+{\alpha }^k{\mathbf{d}}^k)\ge c_2\cdot ({\mathbf{d}}^k{)}^{\mathrm{\top }}\mathrm{\nabla }f({\mathbf{x}}^k)$$

Typically, $c_1={10}^{-4},c_2=0.9$.

Courtesy: Ján Kopačka

L-BFGS Method

Lewis-Overton Line Search

The Lewis-Overton line search is a sophisticated backtracking line search designed specifically for quasi-Newton methods:

1. Given search direction ${\mathbf{d}}^k$, current point ${\mathbf{x}}^k$ and gradient ${\mathbf{g}}^k$
2. Initialize trial step $\alpha :=1$, ${\alpha }_l:=0$, ${\alpha }_r:=+\mathrm{\infty }$
3. Repeat
   1. Update bounds:
      • If $\text{AC}(\alpha )$ fails, set ${\alpha }_r:=\alpha$
      • Else if $\text{CC}(\alpha )$ fails, set ${\alpha }_l:=\alpha$
      • Otherwise, return $\alpha$
   2. Update $\alpha$:
      • If ${\alpha }_r<+\mathrm{\infty }$, set $\alpha :=\text{CubicInterpolate}({\alpha }_l,{\alpha }_r)$ or $\alpha :=({\alpha }_l+{\alpha }_r)/2$
      • Otherwise, set $\alpha :=\text{CubicExtrapolate}({\alpha }_l,{\alpha }_r)$ or $\alpha :=2{\alpha }_l$
   3. Ensure $\alpha \in [{\alpha }_{min},{\alpha }_{max}]$

Cautious Update

Sometimes, when line search is inexact or the function is poorly conditioned, $({\mathbf{y}}^k{)}^{\mathrm{\top }}{\mathbf{s}}^k>0$ cannot gurantee. To ensure numerical stability and maintain the PD Hessian approximation, L-BFGS employs a cautious update strategy:

Skip update condition: If the curvature condition $({\mathbf{y}}^k{)}^{\mathrm{\top }}{\mathbf{s}}^k>ϵ|{\mathbf{s}}^k{|}^2$ is not satisfied, where $ϵ$ is a small positive constant (e.g., ${10}^{-6}$), skip the update for this iteration: ${\mathbf{B}}^{k+1}={\mathbf{B}}^k$.

Powell's Damping: If the curvature condition $({\mathbf{y}}^k{)}^{\mathrm{\top }}{\mathbf{s}}^k\ge \eta ({\mathbf{s}}^k{)}^{\mathrm{\top }}{\mathbf{B}}^k{\mathbf{s}}^k$ is not satisfied, where $\eta$ is a small positive constant (e.g., $0.2$ or $0.25$).

$${\tilde{\mathbf{y}}}^k=\theta {\mathbf{y}}^k+(1-\theta ){\mathbf{B}}^k{\mathbf{s}}^k,\theta =\frac{(1-\eta )\cdot ({\mathbf{s}}^k{)}^{\mathrm{\top }}{\mathbf{B}}^k{\mathbf{s}}^k}{({\mathbf{s}}^k{)}^{\mathrm{\top }}{\mathbf{B}}^k{\mathbf{s}}^k-({\mathbf{y}}^k{)}^{\mathrm{\top }}{\mathbf{s}}^k}$$

Cautious updates guaranteed to have its iterates converge to a critical point if the function has bounded sublevel sets and a Lipschitz continuous gradient.

Two-Loop Recursion

L-BFGS uses a two-loop recursion to compute the search direction without explicitly forming the Hessian approximation. The algorithm maintains a history of the most recent $m$ pairs ${({\mathbf{s}}^i,{\mathbf{y}}^i)}_{i=k-m}^{k-1}$, where typically $m$ is between 5 and 20.

1. Initialize an empty array $\mathcal{A}$ of length $m$, $\mathbf{d}={\mathbf{g}}^k$
2. For $i=k-1,k-2,\dots ,k-m$:
   1. ${\mathcal{A}}^{i+m-k}:=⟨{\mathbf{s}}^i,\mathbf{d}⟩/⟨{\mathbf{s}}^i,{\mathbf{y}}^i⟩$
   2. $\mathbf{d}:=\mathbf{d}-{\mathbf{L}}^{i+m-k}{\mathbf{y}}^i$
3. $\mathbf{d}:=\mathbf{d}\cdot ⟨{\mathbf{s}}^{k-1},{\mathbf{y}}^{k-1}⟩/⟨{\mathbf{y}}^{k-1},{\mathbf{y}}^{k-1}⟩$
4. For $i=k-m,k-m+1,\dots ,k-1$:
   1. $a:=⟨{\mathbf{y}}^i,\mathbf{d}⟩/⟨{\mathbf{s}}^i,{\mathbf{y}}^i⟩$
   2. $\mathbf{d}:=\mathbf{d}+{\mathbf{s}}^i({\mathcal{A}}^{i+m-k}-a)$
5. Return $d$

This approach reduces the storage requirement from $\mathcal{O}(n^2)$ to $\mathcal{O}(mn)$ and the computational cost per iteration from $\mathcal{O}(n^2)$ to $\mathcal{O}(mn)$.

Algorithm Summary

The complete L-BFGS algorithm with cautious update and Lewis-Overton line search:

1. Initialize ${\mathbf{x}}^0,{\mathbf{g}}^0:=\mathrm{\nabla }f({\mathbf{x}}^0)$, choose $m$
2. For $k=0,1,2,\dots$ until convergence:
   1. Compute search direction: ${\mathbf{d}}^k$ using L-BFGS two-loop recursion
   2. Find step size ${\alpha }^k$ using Lewis-Overton line search
   3. Update: ${\mathbf{x}}^{k+1}={\mathbf{x}}^k+{\alpha }^k{\mathbf{d}}^k$
   4. Compute ${\mathbf{s}}^k={\mathbf{x}}^{k+1}-{\mathbf{x}}^k$, ${\mathbf{y}}^k=\mathrm{\nabla }f({\mathbf{x}}^{k+1})-\mathrm{\nabla }f({\mathbf{x}}^k)$
   5. Apply cautious update to $(\{{\mathbf{s}}^k\},\{{\mathbf{y}}^k\})$

Open Source Implementation

• https://github.com/chokkan/liblbfgs
• https://github.com/hjmshi/PyTorch-LBFGS
• https://github.com/ZJU-FAST-Lab/LBFGS-Lite
• https://github.com/yixuan/LBFGSpp