Quasi-Newton and secant equation

Apply Taylor theorem to $\mathrm{\nabla }f$ at $x^{\ast }$ give us

$$\mathrm{\nabla }f(x^{\ast })=\mathrm{\nabla }{f}_k+{\mathrm{\nabla }}^2{f}_k\cdot (x^{\ast }-x_k)+o(\Vert x^{\ast }-x_k\Vert )$$

Ignoring high-order term, we get approximate

$$\begin{array}{rcl}\mathrm{\nabla }{f}_k+{\mathrm{\nabla }}^2{f}_k\cdot (x^{\ast }-x_k)& \approx & \mathrm{\nabla }f(x^{\ast })=0\\ x^{\ast }& \approx & x_k-({\mathrm{\nabla }}^2{f}_k{)}^{-1}\mathrm{\nabla }{f}_k\end{array}$$

Newton method uses the right hand side form to update $x_k$,

$$x_{k+1}=x_k-({\mathrm{\nabla }}^2{f}_k{)}^{-1}\mathrm{\nabla }{f}_k$$

However, the calculation of Hessian matrix’s inverse is usually expensive, and sometimes we cannot guarantee ${\mathrm{\nabla }}^2{f}_k$ is positive definite.

The idea of quasi-Newton method is to replace Hessian matrix with another positive definite symmetric matrix $B_k$, and refine it in every iteration with a much smaller computation cost.

From Newton method formula, we have

$${\mathrm{\nabla }}^2{f}_k\cdot (x_{k+1}-x_k)\approx \mathrm{\nabla }{f}_{k+1}-\mathrm{\nabla }{f}_k$$

Let $$\begin{array}{rcl}{s}_k& =& x_{k+1}-x_k\\ \text{(1)}& y_k& =& \mathrm{\nabla }{f}_{k+1}-\mathrm{\nabla }{f}_k\end{array}$$

The next matrix $B_{k+1}$ we are looking for should satisfy the above approximate,

$$\begin{array}{}\text{(2)}& B_{k+1}{s}_k=y_k\end{array}$$

This equation is called secant equation.

There is another approach to get the secant equation. Let

$$m_k(p)=f_k+\mathrm{\nabla }{f}_k^{T}p+\frac{1}{2}{p}^T{B}_{k}p$$

be the second-order approximate to $f$ at $x_k$. It is a function of $p$ with gradient

$$m_k^{\prime }(p)=\mathrm{\nabla }{f}_k+B_{k}p$$

When we iterate from $k$ to $k+1$, it’s easy to check that $m_{k+1}(0)=f_{k+1}$ and $m_{k+1}^{\prime }(0)=\mathrm{\nabla }{f}_{k+1}$. This means $m_{k+1}(p)$ equals $f$ with both function value and gradient at $x_{k+1}$. If we want this $m_{k+1}(p)$ to approximate $f$ even better, we may ask its gradient to be equal with $f$ at last step $x_k$ as well. At this point, $x_k=x_{k+1}-{\alpha }_k{p}_k$, So $p=-{\alpha }_k{p}_k=-s_k$, and we have

$$\begin{array}{rcl}{m}_{k+1}^{\prime }(-s_k)=\mathrm{\nabla }{f}_{k+1}-B_{k+1}{s}_k& =& \mathrm{\nabla }{f}_k\\ B_{k+1}{s}_k& =& \mathrm{\nabla }{f}_{k+1}-\mathrm{\nabla }{f}_k\\ B_{k+1}{s}_k& =& y_k\end{array}$$

Rank-two update

So the next question becomes how to iterate $B_k$. The Davidon-Fletcher-Powell (DFP) formula (Fletcher and Powell 1963) gives a rank-two matrix update approach.

$$\begin{array}{}\text{(3)}& B_{k+1}=B_k-\frac{B_k{s}_k{s}_k^T{B}_k}{s_k^T{B}_k{s}_k}+\frac{y_k{y}_k^T}{y_k^T{s}_k}\end{array}$$

In Nocedal and Wright (2006), the authors provides another interpretation. They view $B_{k+1}$ as the solution of the following optimization problem.

$$\begin{array}{r}\underset{B}{min}\Vert B-B_k{\Vert }_W\\ s.t.B=B^T,B{s}_k=y_k\end{array}$$

where $\Vert A{\Vert }_W$ norm is defined as

$$\Vert A{\Vert }_W=\Vert W^{1/2}A{W}^{1/2}{\Vert }_F$$

where $\Vert \cdot {\Vert }_F$ is Frobenius norm and $W$ is any matrix satisfying $W{y}_k=s_k$. This is a more intuitive way, however, I haven’t finish the detailed proof of this conclusion. I would revise this part later.

BFGS

The DFP formula maintains properties of $B_k$ like symmetric and positive definite. But we still have to calculate its inverse to get $p_k=-B_k^{-1}\mathrm{\nabla }{f}_k$. To avoid the inverse matrix calculation, we need to approximate matrix $H_k=B_k^{-1}$ directly.

Many materials say just applying Sherman–Morrison formula and we can get a iterate equation for $B_k^{-1}$. However, this matrix deduction is not trivial. After a little calculation, I think applying Woodbury matrix identity would be a more straightforward way.

Woodbury matrix identity says

$$\begin{array}{}\text{(4)}& (A+UCV{)}^{-1}=A^{-1}-A^{-1}U(C^{-1}+V{A}^{-1}U{)}^{-1}V{A}^{-1}\end{array}$$

Let $\rho =\frac{1}{y_k^T{s}_k}$ and drop subscript $k$ to rewrite DFP as

$$\begin{array}{rcl}{B}_{k+1}& =& B-\frac{Bs{s}^{T}B}{s^{T}Bs}+\frac{y{y}^T}{y^{T}s}\\ & =& B+(yBs)\left(\begin{array}{cc}\rho & 0\\ 0& -\frac{1}{s^{T}Bs}\end{array}\right)\left(\begin{array}{c}{y}^T\\ s^{T}B\end{array}\right)\\ \text{(5)}& & =& B+UCV\end{array}$$

So

$$\begin{array}{rcl}{C}^{-1}+V{A}^{-1}U& =& \left(\begin{array}{cc}\frac{1}{\rho }& 0\\ 0& -s^{T}Bs\end{array}\right)+\left(\begin{array}{c}{y}^T\\ s^{T}B\end{array}\right)H(y,Bs)\\ & =& \left(\begin{array}{cc}1/\rho +y^{T}Hy& 1/\rho \\ 1/\rho & 0\end{array}\right)\end{array}$$

We can easily check that

$$\left(\begin{array}{cc}a+b& a\\ a& 0\end{array}\right)\left(\begin{array}{cc}0& \frac{1}{a}\\ \frac{1}{a}& -\frac{a+b}{a^2}\end{array}\right)=I$$

Using this property, we can get

$$\begin{array}{rcl}(C^{-1}+V{A}^{-1}U{)}^{-1}& =& \left(\begin{array}{cc}0& \rho \\ \rho & -\rho -\rho y^{T}Hy\end{array}\right)\end{array}$$

Substitute this back to (4) with (5) (still ignore subscript $k$ for simplicity)

$$\begin{array}{rcl}{H}_{k+1}& =& H-H(yBs)\left(\begin{array}{cc}0& \rho \\ \rho & -\rho -\rho y^{T}Hy\end{array}\right)\left(\begin{array}{c}{y}^T\\ s^{T}B\end{array}\right)H\\ & =& H-\rho s{y}^{T}H-\rho Hy{s}^T+\rho s{y}^{T}Hy{s}^T+\rho s{s}^T\\ & =& (I-\rho s{y}^T)H(I-\rho y{s}^T)+\rho s{s}^T\end{array}$$

Add back subscript $k$ and we finally get the BFGS formula

$$\begin{array}{}\text{(6)}& H_{k+1}=(I-{\rho }_k{s}_k{y}_k^T)H_k(I-{\rho }_k{y}_k{s}_k^T)+{\rho }_k{s}_k{s}_k^T,{\rho }_k=\frac{1}{y_k^T{s}_k}\end{array}$$

Initial $H_0$

A simple choice of $H_0$ is setting it to identical matrix $I$. A similar strategy is to introduce a scalar $\beta$, then let $H_0=\beta I$.

There are many other approaches. For instance, we can also use other optimization method to “warm up”—do a few iterations in order to get a better approximate of ${\mathrm{\nabla }}^{2}f$, then change back to BFGS.

Convergence

I’m going to discuss the convergence theory of Newton and quasi-Newton method in another article later. Therefore, here I only mention the convergence result theorem.

Theorem Suppose that $f$ is twice continuously differentiable and that the iterates generated by the BFGS algorithm converge to a minimizer $x^{\ast }$ at which the Hessian matrix $G$ satisfies Lipschitz continuous in a neighborhood with some positive constant $L$,

$$\Vert G(x)-G(x^{\ast })\Vert \le L\Vert x-x^{\ast }\Vert$$

Suppose the sequence also satisfies

$$\sum _{k=1}^{\mathrm{\infty }}\Vert x_k-x^{\ast }\Vert <\mathrm{\infty }$$

Then $x_k$ converges to $x^{\ast }$ at a superlinear rate.

For L-BFGS and stochastic L-BFGS, some convergence discussion can be found in Mokhtari and Ribeiro (2015).

Wolfe condition

We introduce a helper function

$$\varphi (\alpha )=f(x_k+\alpha p_k),\alpha >0$$

The minimizer of $\varphi (\alpha )$ is what we need. However, solving this univariate minimum problem accurately could be too expensive. An inexact solution is acceptable as long as $\varphi (\alpha )=f(x_k+\alpha p_k)<f_k$.

However, a simple function value reduction may not be enough sometimes. The picture below show an example of this situation. We need a kind of sufficient decrease to avoid this.

Figure 1: Insufficient reduction¹

The following Wolfe condition is the formalization of this sufficient decrease.

$$\begin{array}{rcl}f(x_k+{\alpha }_k{p}_k)& \le & f(x_k)+c_1{\alpha }_k\mathrm{\nabla }{f}_k^T{p}_k\\ \text{(7)}& \mathrm{\nabla }f(x_k+{\alpha }_k{p}_k{)}^T{p}_k& \ge & c_2\mathrm{\nabla }{f}_k^T{p}_k\end{array}$$

where $0<c_1<c_2<1$.

The right hand side of the first inequality is a line $l(\alpha )$ drawn from $x_k$ with the same slope of $\varphi (\alpha )$. Thus the intuition behind the first inequality is that the function value at step ${\alpha }_k$ should below this line $l(\alpha )$. This condition is usually called Armijo condition, or sufficient decrease condition.

The intuition of the second inequality uses more information of $f$’s curvature. Since $p_k$ is the descent direction, we have $\mathrm{\nabla }{f}_k^T{p}_k<0$. Suppose the slope of $\varphi (\alpha )$ at step ${\alpha }_k$ is smaller (steeper) than the slope at the start point, i.e. $\mathrm{\nabla }{f}_k^T{p}_k$, then we can go further safely to reach a even low objective value. Therefore, a sufficient descent step should be at a point with a more gentle slope. This second condition is usually referred to as curvature condition.

Figure 2: The curvature condition²

There is a chance that step ${\alpha }_k$ goes too further that the slope becomes positive. To prevent this case, we can use strong Wolfe condition

$$\begin{array}{rcl}f(x_k+{\alpha }_k{p}_k)& \le & f(x_k)+c_1{\alpha }_k\mathrm{\nabla }{f}_k^T{p}_k\\ \text{(8)}& \Vert \mathrm{\nabla }f(x_k+{\alpha }_k{p}_k{)}^T{p}_k\Vert & \le & c_2\Vert \mathrm{\nabla }{f}_k^T{p}_k\Vert \end{array}$$

where $0<c_1<c_2<1$.

Existance and convergence

The next question is does these ${\alpha }_k$ exist; and if does, how to find them.

The following lemma and theorem (Nocedal and Wright 2006) guarantee the existence of ${\alpha }_k$ that satisfies Wolfe condition, and show us the superlinearly convergence under certain conditions.

Lemma Suppose that $f$ is continuously differentiable. Let $p_k$ be a descent direction at $x_k$, and assume that $f$ is bounded below along $\{x_k+\alpha p_k|\alpha >0\}$. For $0<c_1<c_2<1$, there exist intervals of step lengths that satisfying the Wolfe condition (7) and strong Wolfe condition (8).

Theorem Suppose that $f$ is twice continuously differentiable. Consider the iteration $x_{k+1}=x_k+{\alpha }_k{p}_k$ , where $p_k$ is a descent direction and ${\alpha }_k$ satisfies the Wolfe conditions (7) with $c_1\le 1/2$. If the sequence $\{x_k\}$ converges to a point $x^{\ast }$ such that $\mathrm{\nabla }f(x^{\ast })$ and ${\mathrm{\nabla }}^{2}f(x^{\ast })$ is positive definite, and if the search direction satisfies

$$\underset{k\to \mathrm{\infty }}{lim}\frac{\Vert \mathrm{\nabla }{f}_k+{\mathrm{\nabla }}^2{f}_k{p}_k\Vert }{\Vert p_k\Vert }=0$$

then

1. the step length ${\alpha }_k=1$ is admissible for all $k$ greater than a certain index $k_0$;
2. if ${\alpha }_k=1$ for all $k>k_0$, $\{x_k\}$ converges to $x^{\ast }$ superlinearly.

Linear search algorithm

We use linear search algorithm to locate a valid $\alpha$. The idea is to generate a sequence of monotonically increasing $\{{\alpha }_i\}$ in $(0,{\alpha }_{max})$. If ${\alpha }_i$ satisfies Wolfe condition, return that step; otherwise narrow the search interval.

We show the algorithm in Julia pseudo code below.

using Flux

function linear_search(ϕ, α_0=0, α_max=1, c1=0.0001, c2=0.9)
    α = Dict(0 => α_0)
    α[1] = choose(α_0, α_max)
    ϕ′(x) = gradient(ϕ, x)[1]
    i = 1
    while true
        y = ϕ(α[i])
        if y > ϕ(0) + c1 * α[i] * ϕ′(0) or (ϕ(α[i]) >= ϕ(α[i - 1]) and i > 1)
            return zoom(α[i - 1], α[i])
        end
        dy = ϕ′(α[i])
        if abs(dy) <= -c2 * ϕ′(0)
            return α[i]
        end
        if dy >= 0
            return zoom(α[i], α[i - 1])
        end
        α[i + 1] = choose(α[i], α_max)
        i = i + 1
    end
end

function zoom(ϕ, α_lo, α_hi, c2=0.9)
    while true
        # using quadratic, cubic, or bisection to find a trial step length α
        α = choose(α_lo, α_hi)
        y = ϕ(α)
        if y > ϕ(0) + c1 * α * ϕ′(0) or y >= ϕ(α_lo)
            α_hi = α
        else
            dy = ϕ′(α)
            if abs(dy) <= -c2 * ϕ′(0)
                return α
            end
            if dy * (α_hi - α_lo) >= 0
                α_hi = α_lo
            end
            α_lo = α
        end
    end
end

Notice that in function zoom(), the order of ${\alpha }_i$ and ${\alpha }_{i-1}$ could swap. ${\alpha }_{lo}$ always gives the smallest function value.

L-BFGS: vanilla iteration

Recall (1) and BFGS (6), and for simplicity, let

$$V_k=I-{\rho }_k{y}_k{s}_k^T$$

The BFGS formula becomes

$$\begin{array}{}\text{(9)}& H_{k+1}=V_k^T{H}_k{V}_k+{\rho }_k{s}_k{s}_k^T\end{array}$$

We can use a vanilla iterative way to calculate $H_k\mathrm{\nabla }{f}_k$ directly (let $q=\mathrm{\nabla }{f}_k$)

1. $V_k\mathrm{\nabla }{f}_k=q-{\rho }_k{y}_k{s}_k^{T}q$: calculate $\rho s^{T}q$ first in $n$ multiplications to get a scalar $\alpha$, then calculate $q-\alpha y$ in another $n$ multiplications.
2. Suppose $H_k^0$ is a diagonal matrix, so we can calculate $H_k^0{V}_{k}q$ in $n$ multiplications.
3. Multiply $V_k^T$ is analogous, which needs $2n$ multiplications as well.
4. Finally, $\rho s{s}^T$ and add the results together needs another $2n$ multiplications.

So this vanilla iteration requires $7n$ multiplications, leading to a total cost of $7nm$ multiplications.

Can we do better?

L-BFGS: two-loop recursion

L-BFGS has a two-loop recursion algorithm, which is quite brilliant.

function L_BFGS(H0, ∇f_k, s, y, ρ, k, m)
    q = ∇f_k
    for i = k - 1 : -1 : k - m
        α[i] = ρ[i] * transpose(s[i]) * q
        q = q - α[i] * y[i]
    end
    w = H0 * q
    for i = k - m : k - 1
        β = ρ[i] * transpose(y[i]) * w
        w = w + s[i] * (α[i] - β)
    end
    return w
end

The return value $w=H_k\mathrm{\nabla }{f}_k$ is what we need, and there are only $5nm$ multiplications in it.

The data manipulation in the two-loop recursion is not that easy to understand at first glimpse. To make it more clear, we expand (9) $m$ steps to get

Some key observations to help you understand the two-loop recursion:

• In the first loop, ${\alpha }_i={\rho }_i{s}_i^T{V}_{i+1}\cdots V_{k-1}\mathrm{\nabla }{f}_k$
• After the first loop and w = H0 * q, $q=H_k^0{V}_{k-m}\cdots V_{k-1}$
• In the second loop, after iteration $i$, the vector $w$ will be