Bfgs and limited memory quasi newton methods

Fix an integer $n\ge 1$. For successive points $x,x^{+}\in {\mathbb{R}}^n$ and vectors $g,g^{+}\in {\mathbb{R}}^n$, define the step vector $s:=x^{+}-x$ and the gradient-difference vector $y:=g^{+}-g$. Their curvature is $s^{T}y$, using the Euclidean inner product of The Euclidean inner product $⟨x,y⟩={\sum }_{k<n}{x}_k{y}_k$ on ${\mathbb{R}}^n$.

A real matrix $B\in {\mathbb{R}}^{n\times n}$ is symmetric positive definite (SPD) when $B=B^T$ and $z^{T}Bz>0$ for every nonzero $z\in {\mathbb{R}}^n$. For an SPD matrix $B$, write $H=B^{-1}$ for its inverse; $I$ denotes the identity matrix of Rectangular matrix multiplication and the identity matrix $I_n$, including zero-sized shapes. If $s\mathrm{\ne }0$ and $s^{T}y>0$, set $\rho :=(s^{T}y{)}^{-1}$. The BFGS Hessian update and its inverse form are $B^{+}:=B-\frac{Bs{s}^{T}B}{s^{T}Bs}+\frac{y{y}^T}{s^{T}y},H^{+}:=(I-\rho s{y}^T)H(I-\rho y{s}^T)+\rho s{s}^T\mathrm{.}$ The denominators are positive under the stated hypotheses.

For Powell damping, first set $c:=s^{T}Bs$. If $s^{T}y\ge 0.2c$, set $\tilde{y}:=y$. If $s^{T}y<0.2c$, set $\theta :=\frac{0.8c}{c-s^{T}y},\tilde{y}:=\theta y+(1-\theta )Bs\mathrm{.}$ The damped BFGS update is the same formula with $y$ replaced by $\tilde{y}$ and $\rho =(s^T\tilde{y}{)}^{-1}$.

The secant equations are $B^{+}s=y$ and $H^{+}y=s$ for an undamped pair; for a damped pair they use $\tilde{y}$. An inverse pair $(B,H)$ is consistent when $BH=HB=I$.

The two update formulas are algebraically equivalent when $B$ and $H$ are consistent. A positive curvature pair is required for both updates. Powell damping changes the gradient-difference vector used in the secant equation; a convergence theorem for the undamped update does not automatically cover the damped branch.

Let $n\ge 1$, let $U\subseteq {\mathbb{R}}^n$ be open, and let $f:U\to \mathbb{R}$ have continuous first and second partial derivatives. Put $g(x):=\mathrm{\nabla }f(x)$ and $G(x):={\mathrm{\nabla }}^{2}f(x)$, using the conventions of The Jacobian matrix of partial derivatives and the gradient in the scalar-valued case and The Hessian matrix and critical points of a scalar field. Let $K\subseteq U$ be convex. Suppose that for every $x\in K$, $G(x)$ is symmetric and, for constants $0<\mu \le \mathrm{\Lambda }<\mathrm{\infty }$, $\mu \mathrm{\parallel }w{\mathrm{\parallel }}_2^2\le w^{T}G(x)w(w\in {\mathbb{R}}^n),\mathrm{\parallel }G(x){\mathrm{\parallel }}_2\le \mathrm{\Lambda },$ where $\mathrm{\parallel }A{\mathrm{\parallel }}_2:={\sup }_{\mathrm{\parallel }w{\mathrm{\parallel }}_2=1}\mathrm{\parallel }Aw{\mathrm{\parallel }}_2$. Then for all $u,v\in K$, with $d:=u-v$, $\mathrm{\parallel }g(u)-g(v){\mathrm{\parallel }}_2\le \mathrm{\Lambda }\mathrm{\parallel }d{\mathrm{\parallel }}_2,(g(u)-g(v){)}^{T}d\ge \mu \mathrm{\parallel }d{\mathrm{\parallel }}_2^2,$ and $\frac{\mu }{2}\mathrm{\parallel }d{\mathrm{\parallel }}_2^2\le f(u)-f(v)-g(v{)}^{T}d\le \frac{\mathrm{\Lambda }}{2}\mathrm{\parallel }d{\mathrm{\parallel }}_2^2\mathrm{.}$ Separately, let $W\subseteq U$ be convex, let $x_{\ast }\in W$, and suppose $G$ is $L_G$-Lipschitz on $W$ in operator norm. Then for all $u,v\in W$, with $d:=u-v$, $\mathrm{\parallel }g(u)-g(v)-G(x_{\ast })d{\mathrm{\parallel }}_2\le L_G\max \{\mathrm{\parallel }u-x_{\ast }{\mathrm{\parallel }}_2,\mathrm{\parallel }v-x_{\ast }{\mathrm{\parallel }}_2\}\mathrm{\parallel }d{\mathrm{\parallel }}_2\mathrm{.}$

Let $B\in {\mathbb{R}}^{n\times n}$ be SPD, let $H=B^{-1}$, and let $s,y\in {\mathbb{R}}^n$ satisfy $s\mathrm{\ne }0$ and $d:=s^{T}y>0$. Define $B^{+},H^{+}$ by Bfgs curvature pairs and Hessian/inverse updates. Then $B^{+}$ and $H^{+}$ are SPD, $B^{+}s=y$, $H^{+}y=s$, and $B^{+}{H}^{+}=I$. In particular, the two updates remain consistent inverses and satisfy their secant equations.

Let $f:D\to \mathbb{R}$ be $C^1$ on an open $D\subseteq {\mathbb{R}}^n$. At $x\in D$, let $p\mathrm{\ne }0$ be a strict descent direction, so $g^{T}p<0$ for $g:=\mathrm{\nabla }f(x)$. Suppose $\alpha >0$, $x+\alpha p\in D$, and $\alpha$ satisfies the Wolfe curvature condition of Wolfe and strong Wolfe line-search conditions. With $s:=\alpha p$ and $y:=\mathrm{\nabla }f(x+\alpha p)-\mathrm{\nabla }f(x)$, one has $s^{T}y\ge (1-c_2)\alpha \mathrm{\mid }{g}^{T}p\mathrm{\mid }>0.$

Let $B$ be SPD, let $s\mathrm{\ne }0$, and put $c:=s^{T}Bs>0$ and $q:=s^{T}y$ for an arbitrary $y$. Define $\tilde{y}$ by Powell damping in Bfgs curvature pairs and Hessian/inverse updates. Then $s^T\tilde{y}\ge 0.2c>0$. In the damped branch, the coefficient $\theta$ is strictly between $0$ and $1$ and $s^T\tilde{y}=0.2c$.

Input a real symmetric positive-definite matrix $B\in {\mathbb{R}}^{n\times n}$, $n\ge 1$, a vector $g\in {\mathbb{R}}^n$, and representable finite entries in the operationwise arithmetic model of Operationwise real and complex floating point arithmetic with rounded square roots and fused multiply-add. Output: either a computed direction $p$ or a failure status.

1. Call Failure-aware rounded Cholesky factorization on $B$. If it returns a failure status, return factorization-failure with that pivot and operation index.

2. If it returns a lower-triangular factor $\widehat{L}$ with positive normal diagonal, implement the recursions of Forward and backward substitution for triangular linear systems as follows, with indices $1,\dots ,n$. For the forward solve, process rows $i=1,\dots ,n$ in increasing order. Initialize $s:=0$; for $j=1,\dots ,i-1$ in increasing order, first compute $t:=\mathrm{fl}({\widehat{l}}_{ij}{z}_j)$, then $s:=\mathrm{fl}(s+t)$. After this inner loop, compute $r:=\mathrm{fl}((-g_i)-s)$, then $z_i:=\mathrm{fl}(r\mathrm{/}{\widehat{l}}_{ii})$. For the backward solve, process rows $i=n,\dots ,1$ in decreasing order. Initialize $s:=0$; for $j=i+1,\dots ,n$ in increasing order, first compute $t:=\mathrm{fl}({\widehat{l}}_{ji}{p}_j)$, then $s:=\mathrm{fl}(s+t)$. After this inner loop, compute $r:=\mathrm{fl}(z_i-s)$, then $p_i:=\mathrm{fl}(r\mathrm{/}{\widehat{l}}_{ii})$. Empty inner loops leave $s=0$, and sign changes are exact. Each displayed $\mathrm{fl}$ is a separate scalar operation in the given arithmetic model; do not reassociate or fuse operations. On the first range-failure, stop immediately and return triangular-solve-failure with the solve (forward or backward), row $i$, inner index $j$ when applicable, and operation (multiply, add, subtract, or divide).

3. Otherwise return $p$. The exact-real reference problem is $Bp=-g$ and has a unique solution by Cholesky solves Hermitian positive-definite systems and has about half the factorisation cost of LU. The returned vector is the computed solve from $\widehat{L}$; this definition asserts no backward-error bound for it.

The procedure performs a finite factorization loop and two finite triangular loops. Its use in the full BFGS algorithm includes a subsequent descent check; if the computed direction is not descending, that algorithm resets to a positive scalar multiple of the negative gradient.

Inputs are a $C^1$ objective $f:D\to \mathbb{R}$ on an open set $D\subseteq {\mathbb{R}}^n$ with $n\ge 1$, an initial point $x_0\in D$, a positive integer memory length $m$, a fallback scale ${\gamma }_0>0$, line-search constants $0<c_1<c_2<1\mathrm{/}2$, and finite expansion, zoom, and iteration budgets $E,Z,K\in \mathbb{N}$ (zero is allowed). All evaluations of $f$ and $\mathrm{\nabla }f$, and all scalar and vector arithmetic, are exact over $\mathbb{R}$. Initialize $k:=0$, the current point $x:=x_0$, the current gradient $g:=\mathrm{\nabla }f(x_0)$, and an empty history list. If $g=0$, return stationary. Otherwise, while $k<K$, store at most $m$ pairs $(s_i,y_i)$ in chronological order, retaining only pairs with $s_i^T{y}_i>0$. Let ${\rho }_i=(s_i^T{y}_i{)}^{-1}$ for each stored pair.

At each nonstationary iterate $x$ with current gradient $g$, all inner products are the Euclidean ones of The Euclidean inner product $⟨x,y⟩={\sum }_{k<n}{x}_k{y}_k$ on ${\mathbb{R}}^n$. Choose $\gamma ={\gamma }_0$ if the list is empty; otherwise choose $\gamma =\frac{s_{\mathrm{l}\mathrm{a}\mathrm{s}\mathrm{t}}^T{y}_{\mathrm{l}\mathrm{a}\mathrm{s}\mathrm{t}}}{y_{\mathrm{l}\mathrm{a}\mathrm{s}\mathrm{t}}^T{y}_{\mathrm{l}\mathrm{a}\mathrm{s}\mathrm{t}}}\mathrm{.}$ The denominator is positive because the retained pair has positive curvature, so $y_{\mathrm{l}\mathrm{a}\mathrm{s}\mathrm{t}}\mathrm{\ne }0$. Compute $r=H_{k}g$ by the following two loops, then set $p=-r$:

1. Set $q:=g$. Traverse the stored pairs from newest to oldest, setting ${\alpha }_i:={\rho }_i{s}_i^{T}q$ and then $q:=q-{\alpha }_i{y}_i$.
2. Set $r:=\gamma q$. Traverse the pairs from oldest to newest, setting ${\beta }_i:={\rho }_i{y}_i^{T}r$ and then $r:=r+s_i({\alpha }_i-{\beta }_i)$.

If $g=0$, return stationary. If $g^{T}p\ge 0$, clear the list, set $p:=-{\gamma }_{0}g$, and return direction-failure only if this fallback is not a strict descent direction.

Call Finite-budget bracketing and zoom for a strong Wolfe step with $(f,x,p,c_1,c_2,E,Z)$. If it returns a failure status, return that status with the current point. If it returns $\alpha >0$, set $x^{+}:=x+\alpha p$, compute $g^{+}:=\mathrm{\nabla }f(x^{+})$, and set $s:=x^{+}-x$, $y:=g^{+}-g$. If $s^{T}y>0$, append $(s,y)$ and $\rho =(s^{T}y{)}^{-1}$, dropping the oldest pair if the list would exceed $m$. If $s^{T}y\le 0$, skip that pair and preserve the existing valid list.

After either storing or skipping the new pair, commit the accepted step by setting $x:=x^{+}$ and $g:=g^{+}$, then increment $k:=k+1$. If the updated gradient is zero, return stationary at $x$; otherwise continue while $k<K$. When the iteration budget is exhausted, return the current point and iteration-budget-exhausted. A line-search failure returns its status with the current point, and the fallback branch returns direction-failure only if the fallback direction is not a strict descent direction. Finite line-search budgets may fail even when a strong Wolfe step exists.

Output the last point and one of stationary, the line-search status, direction-failure, or iteration-budget-exhausted.

The initial inverse-Hessian scale is $\gamma I$, an SPD matrix. Positive curvature is checked before storing a pair. The pair filter protects the recursion from a zero or negative denominator; it does not turn a failed finite search into an accepted step.

Let $B$ be SPD, $s\mathrm{\ne }0$, and $y$ satisfy $d:=s^{T}y>0$. Let $B^{+}$ be the Hessian update of Bfgs curvature pairs and Hessian/inverse updates. Define $m_s:=\frac{s^{T}y}{\mathrm{\parallel }s{\mathrm{\parallel }}_2^2},M_s:=\frac{\mathrm{\parallel }y{\mathrm{\parallel }}_2^2}{s^{T}y},q:=\frac{s^{T}Bs}{\mathrm{\parallel }s{\mathrm{\parallel }}_2^2},\eta :=\frac{s^{T}Bs}{\mathrm{\parallel }s{\mathrm{\parallel }}_2\mathrm{\parallel }Bs{\mathrm{\parallel }}_2}\mathrm{.}$ Then $\eta >0$, $\eta \le 1$, and $\mathrm{tr}(B^{+})=\mathrm{tr}(B)-\frac{\mathrm{\parallel }Bs{\mathrm{\parallel }}_2^2}{s^{T}Bs}+\frac{\mathrm{\parallel }y{\mathrm{\parallel }}_2^2}{s^{T}y},\det (B^{+})=\det (B)\frac{s^{T}y}{s^{T}Bs}\mathrm{.}$ The potential $\psi (B):=\mathrm{tr}(B)-\log \det (B)$ is well-defined and satisfies $\psi (B)\ge n$. With $h(t):=1-t+\log t$ for $t>0$, one has $h(t)\le 0$, with equality exactly at $t=1$; $h$ is strictly increasing on $(0,1)$ and strictly decreasing on $(1,\mathrm{\infty })$, and $\psi (B^{+})\le \psi (B)+(M_s-\log m_s-1)+h\text{\hspace{0.17em}⁣}\left(\frac{q}{{\eta }^2}\right)+\log {\eta }^2\mathrm{.}$

Inputs are a $C^1$ objective $f:D\to \mathbb{R}$ on an open $D\subseteq {\mathbb{R}}^n$ with $n\ge 1$, $x_0\in D$, consistent SPD matrices $B_0,H_0$ with $B_0{H}_0=H_0{B}_0=I$, a reset scale $\gamma >0$, constants $0<c_1<c_2<1\mathrm{/}2$, finite line-search and iteration budgets $E,Z,K\in \mathbb{N}$ (zero is allowed), a direction mode exact or rounded-Cholesky, and an update mode undamped or Powell. For rounded-Cholesky mode, also input a finite binary arithmetic format for the called direction procedure.

Initialize $k:=0$, $x_k:=x_0$, $B_k:=B_0$, $H_k:=H_0$, and $g_0:=\mathrm{\nabla }f(x_0)$. If $g_0=0$, return stationary at $x_0$. While $k<K$, if $g_k=0$, return stationary at $x_k$. Use the Euclidean inner product of The Euclidean inner product $⟨x,y⟩={\sum }_{k<n}{x}_k{y}_k$ on ${\mathbb{R}}^n$ in descent and curvature tests. In exact direction mode set $p_k:=-H_k{g}_k$. In rounded-Cholesky mode, first check that every entry of $B_k$ and $g_k$ is finite and representable in the selected format. If any entry is not, return direction-input-failure with $x_k$. Otherwise call Bfgs Hessian direction by failure-aware Cholesky solve with $(B_k,g_k)$ in that format; if the factorization or triangular solve fails, return its failure status with $x_k$. Otherwise set $p_k$ to its returned direction. If the computed direction is not descending, reset $B_k:={\gamma }^{-1}I$, $H_k:=\gamma I$, and set $p_k:=-\gamma g_k$.

Call Finite-budget bracketing and zoom for a strong Wolfe step with $(f,x_k,p_k,c_1,c_2,E,Z)$. If it returns a failure status, return that status and $x_k$. Otherwise take its positive strong Wolfe step ${\alpha }_k$, set $x_{k+1}:=x_k+{\alpha }_k{p}_k$, and compute $g_{k+1}:=\mathrm{\nabla }f(x_{k+1})$. Define $s_k:={\alpha }_k{p}_k,y_k:=g_{k+1}-g_k\mathrm{.}$ In undamped mode set ${\bar{y}}_k:=y_k$. In Powell mode put $c_k:=s_k^T{B}_k{s}_k$ and use the damping rule of Bfgs curvature pairs and Hessian/inverse updates to form ${\bar{y}}_k$. If $s_k^T{\bar{y}}_k\le 0$, skip the matrix update and return curvature-rejected with $x_{k+1}$. Otherwise update both matrices by $B_{k+1}=B_k-\frac{B_k{s}_k{s}_k^T{B}_k}{s_k^T{B}_k{s}_k}+\frac{{\bar{y}}_k{\bar{y}}_k^T}{s_k^T{\bar{y}}_k},H_{k+1}=(I-{\rho }_k{s}_k{\bar{y}}_k^T)H_k(I-{\rho }_k{\bar{y}}_k{s}_k^T)+{\rho }_k{s}_k{s}_k^T,{\rho }_k=(s_k^T{\bar{y}}_k{)}^{-1}\mathrm{.}$ If any required scalar or matrix operation is undefined, discard partial update results and return update-failure with $x_{k+1}$. After a successful update, commit $(x_{k+1},g_{k+1},B_{k+1},H_{k+1})$ as the state for the next iteration. If $g_{k+1}=0$, return stationary at $x_{k+1}$; if $k+1=K$, return iteration-budget-exhausted at $x_{k+1}$. Otherwise set $k:=k+1$ and continue from that committed state. A line-search failure returns its status with $x_k$. If the loop exits because $k=K$ without another return, return iteration-budget-exhausted with $x_k$; this covers $K=0$. The output is the last point and one of stationary, a line-search status, direction-input-failure, factorization/solve failure, curvature-rejected, update-failure, or iteration-budget-exhausted.