Compact representation of the full Broyden class of quasi-Newton updates

TL;DR

This paper develops a compact representation for matrices in the Broyden class of quasi-Newton updates, enabling efficient computation of inverses, solutions to linear systems, and eigenvalues, thus facilitating sensitivity analysis.

Contribution

It extends previous work by providing a compact representation for the full Broyden class, including rank-one and rank-two updates, with practical algorithms for inverse and linear system computations.

Findings

Accurately represents Broyden class matrices using the compact form.

Efficiently computes inverses and solves linear systems with these matrices.

Enables eigenvalue, condition number, and sensitivity analysis for Broyden matrices.

Abstract

In this paper, we present the compact representation for matrices belonging to the the Broyden class of quasi-Newton updates, where each update may be either rank-one or rank-two. This work extends previous results solely for the restricted Broyden class of rank-two updates. In this article, it is not assumed the same Broyden update is used each iteration; rather, different members of the Broyden class may be used each iteration. Numerical experiments suggest that a practical implementation of the compact representation is able to accurately represent matrices belonging to the Broyden class of updates. Furthermore, we demonstrate how to compute the compact representation for the inverse of these matrices, as well as a practical algorithm for solving linear systems with members of the Broyden class of updates. We demonstrate through numerical experiments that the proposed linear solver is able to efficiently solve linear systems with members of the Broyden class of matrices to high accuracy. As an immediate consequence of this work, it is now possible to efficiently compute the eigenvalues of any limited-memory member of the Broyden class of matrices, allowing for the computation of condition numbers and the ability perform sensitivity analysis.

Tables

Table 1. Table 1. The values of ϕ i subscript italic-ϕ 𝑖 \phi_{i} for 0 ≤ i ≤ 4 0 𝑖 4 0\leq i\leq 4 for each experiment. The choice of ϕ 1 = 1 subscript italic-ϕ 1 1 \phi_{1}=1 corresponds to the BFGS update while ϕ 3 = 0 subscript italic-ϕ 3 0 \phi_{3}=0 corresponds to the DFP update. Note that Experiment 1 does not use SR1 updates.

Experiment	${\varphi }_0$	${\varphi }_1$	${\varphi }_2$	${\varphi }_3$	${\varphi }_4$
1		1		0	${\varphi }_4>1$
2		1	${\varphi }_2^{SR11}$	0	${\varphi }_4>1$
3		1	${\varphi }_2^{SR11}$	${\varphi }_3^{SR11}$	${\varphi }_4>1$
4	${\varphi }_0^{SR11}$	1	${\varphi }_2^{SR11}$	0	${\varphi }_4>1$

Table 2. Table 2. Average relative error over ten different trials for each experiment with n = 100 , 1 , 000 , 𝑛 100 1 000 n=100,1,000, and 10 , 000 10 000 10,000 .

$n$	Exp. 1	Exp. 2	Exp. 3	Exp. 4
100	1.1315e-13	1.3383e-11	1.6749e-12	2.2855e-14
1,000	3.2039e-14	1.1225e-14	5.4247e-15	1.0155e-15
10,000	1.3426e-13	8.5453e-14	1.9969e-13	2.8354e-16

Table 3. Table 3. Average relative error over ten different trials for each experiment with n = 100 , 1 , 000 , 𝑛 100 1 000 n=100,1,\!000, and 10 , 000 10 000 10,\!000 .

$n$	Exp. 1	Exp. 2	Exp. 3	Exp. 4
100	4.0158e-13	1.342e-10	1.3065e-09	2.8160e-14
1,000	1.518e-14	7.6460e-14	6.1744e-14	1.8431e-13
10,000	2.4175e-12	1.6079e-12	4.3284e-12	1.8795e-14

Table 4. Table 4. Average computational times for solving B k + 1 ​ r = z subscript 𝐵 𝑘 1 𝑟 𝑧 B_{k+1}r=z using the inverse compact representation (ICR) of the inverse and the MATLAB “backslash” command with n = 100 , 1 , 000 𝑛 100 1 000 n=100,1,000 , and 10 , 000 10 000 10,000 .

	$n=100$		$n=1,000$		$n=10,000$
Exp.	ICR	MATLAB	ICR	MATLAB	ICR	MATLAB
1	9.9e-04	3.3e-04	1.1e-03	2.9e-02	3.7e-03	1.1e+01
2	7.2e-04	3.7e-04	1.1e-03	2.9e-02	3.5e-03	1.1e+01
3	6.7e-04	3.1e-04	1.1e-03	3.0e-02	3.2e-03	1.1e+01
4	5.8e-04	3.5e-04	9.8e-04	3.0e-02	2.8e-03	1.1e+01

Equations

$$s_k\stackrel{△}{=}{x}_{k+1}-x_k\text{and}{y}_k\stackrel{△}{=}\nabla f(x_{k+1})-\nabla f(x_k).$$

$$B_{k+1}=B_k-\frac{1}{s_k^T{B}_k{s}_k}{B}_k{s}_k{s}_k^T{B}_k+\frac{1}{y_k^T{s}_k}{y}_k{y}_k^T+{\varphi }_k(s_k^T{B}_k{s}_k)w_k{w}_k^T,$$

$$w_k=\frac{y_k}{y_k^T{s}_k}-\frac{B_k{s}_k}{s_k^T{B}_k{s}_k}.$$

$$S_k$$

$$Y_k$$

$$S_k^T{Y}_k=L_k+D_k+R_k,$$

$$B_{k+1}=B_0+{\Psi }_k{M}_k{\Psi }_k^T,$$

$${\Psi }_k\stackrel{△}{=}\left(\begin{array}{cc}{B}_0{S}_k& Y_k\end{array}\right)\text{and}{M}_k\stackrel{△}{=}-{\left(\begin{array}{cc}{S}_k^T{B}_0{S}_k& L_k\\ L_k^T& -D_k\end{array}\right)}^{-1},$$

$${\Psi }_k\stackrel{△}{=}\left(\begin{array}{cc}{B}_0{S}_k& Y_k\end{array}\right)\text{and}{M}_k={\left(\begin{array}{cc}-S_k^T{B}_0{S}_k+\varphi {\Lambda }_k& -L_k+\varphi {\Lambda }_k\\ -L_k^T+\varphi {\Lambda }_k& D_k+\varphi {\Lambda }_k\end{array}\right)}^{-1},$$

$${\lambda }_i\stackrel{△}{=}{\left(-\frac{1-\varphi }{s_i^T{B}_i{s}_i}-\frac{\varphi }{s_i^T{y}_i}\right)}^{-1}.$$

$${\Psi }_k=Y_k-B_0{S}_k\text{and}{M}_k=(D_k+L_k+L_k^T-S_k^T{B}_0{S}_k{)}^{-1}.$$

$$B_{k+1}=B_0+{\Psi }_k{M}_k{\Psi }_k^T=B_0+QR{M}_k{R}^T{Q}^T.$$

$$B_{k+1}=B_0+QVD{V}^T{Q}^T=\gamma I+QVD{V}^T{Q}^T=QV(\gamma I+\hat{D})V^T{Q}^T,$$

$$B_{k+1}=B_k+\left(B_k{s}_k{y}_k\right)O_k\left(\begin{array}{c}(B_k{s}_k{)}^T\\ y_k^T\end{array}\right),$$

$$O_k=\left(\begin{array}{cc}-\frac{(1-{\varphi }_k)}{s_k^T{B}_k{s}_k}& -\frac{{\varphi }_k}{y_k^T{s}_k}\\ -\frac{{\varphi }_k}{y_k^T{s}_k}& \left(1+{\varphi }_k\frac{s_k^T{B}_k{s}_k}{y_k^T{s}_k}\right)\frac{1}{y_k^T{s}_k}\end{array}\right).$$

$$\det (O_k)$$

$$-\frac{(1-{\varphi }_k)}{s_k^T{B}_k{s}_k}-\frac{{\varphi }_k}{y_k^T{s}_k}=0;$$

$$O_k^{-1}=\left(\begin{array}{cc}-s_k^T{B}_k{s}_k+\frac{{\varphi }_k}{{\alpha }_k+{\beta }_k}& \frac{{\varphi }_k}{{\alpha }_k+{\beta }_k}\\ \frac{{\varphi }_k}{{\alpha }_k+{\beta }_k}& y_k^T{s}_k+\frac{{\varphi }_k}{{\alpha }_k+{\beta }_k}\end{array}\right),$$

$${\Pi }_k=\left(\begin{array}{cccc}{I}_k& 0& 0& 0\\ 0& 0& I_k& 0\\ 0& 1& 0& 0\\ 0& 0& 0& 1\end{array}\right),$$

$${\Xi }_k=\left(\begin{array}{cc}{\Pi }_{k-1}^T{\Xi }_{k-1}& 0\\ 0& E_k\end{array}\right),\text{where }{E}_k=\{\begin{array}{ll}(-11{)}^T& \text{if ϕk=ϕkSR11}\\ I_2& \text{otherwise,}\end{array}$$

$$\Gamma_{k}=\underset{0\leq j\leq k}{\text{diag}}\big{(}\gamma_{j}\big{)},\qquad\text{where \ }\gamma_{j}=\begin{cases}\phi_{j}\left(-\frac{1-\phi_{j}}{s_{j}^{T}B_{j}s_{j}}-\frac{\phi_{j}}{s_{j}^{T}y_{j}}\right)^{-1}&\text{if}\ \phi_{j}\neq\phi_{j}^{{\small SR1}1}\\ \qquad\quad\ 0&\text{otherwise.}\end{cases}$$

$$B_{k+1}=B_0+{\hat{\Psi }}_k{\hat{M}}_k{\hat{\Psi }}_k^T,$$

$${\hat{M}}_k={\left({\Xi }_k^T{\Pi }_k\left(\begin{array}{cc}-S_k^T{B}_0{S}_k+{\Gamma }_k& -L_k+{\Gamma }_k\\ -L_k^T+{\Gamma }_k& D_k+{\Gamma }_k\end{array}\right){\Pi }_k^T{\Xi }_k\right)}^{-1},$$

$${\hat{\Psi }}_k={\Psi }_k{\Pi }_k^T{\Xi }_k.$$

$${\hat{M}}_0={\left({\Xi }_0^T{\Pi }_0\left(\begin{array}{cc}-s_0^T{B}_0{s}_0+{\gamma }_0& {\gamma }_0\\ {\gamma }_0& d_0+{\gamma }_0\end{array}\right){\Pi }_0^T{\Xi }_0\right)}^{-1}.$$

$$M_0\stackrel{△}{=}\left(\begin{array}{cc}-\frac{(1-{\varphi }_0)}{s_0^T{B}_0{s}_0}& -\frac{{\varphi }_0}{y_0^T{s}_0}\\ -\frac{{\varphi }_0}{y_0^T{s}_0}& \left(1+{\varphi }_0\frac{s_0^T{B}_0{s}_0}{y_0^T{s}_0}\right)\frac{1}{y_0^T{s}_0}\end{array}\right).$$

$$M_0=\left(\begin{array}{cc}{\alpha }_0& {\beta }_0\\ {\beta }_0& {\delta }_0\end{array}\right),$$

$${\alpha }_0=-\frac{(1-{\varphi }_0)}{s_0^T{B}_0{s}_0},{\beta }_0=-\frac{{\varphi }_0}{y_0^T{s}_0},\text{and}{\delta }_0=\left(1+{\varphi }_0\frac{s_0^T{B}_0{s}_0}{y_0^T{s}_0}\right)\frac{1}{y_0^T{s}_0}.$$

$$M_0=\left(\begin{array}{cc}-{\beta }_0& {\beta }_0\\ {\beta }_0& -{\beta }_0\end{array}\right)=(-{\beta }_0)\left(\begin{array}{cc}1& -1\\ -1& 1\end{array}\right)=-{\beta }_0{\Xi }_0{\Xi }_0^T.$$

$${\hat{M}}_0={\left({\Xi }_0^T\left(\begin{array}{cc}-s_0^T{B}_0{s}_0+{\gamma }_0& {\gamma }_0\\ {\gamma }_0& s_0^T{y}_0+{\gamma }_0\end{array}\right){\Xi }_0\right)}^{-1}=\frac{1}{s_0^T{y}_0-s_0^T{B}_0{s}_0}=-{\beta }_0,$$

$$M_0^{-1}=\left(\begin{array}{cc}-s_0^T{B}_0{s}_0+{\gamma }_0& {\gamma }_0\\ {\gamma }_0& s_0^T{y}_0+{\gamma }_0\end{array}\right)={\hat{M}}_0^{-1}.$$

Full text

Research supported in part by NSF grants CMMI-1334042 and CMMI-1333326.

1. Introduction

Quasi-Newton methods for minimizing a continuously differentiable function $f:\Re^{n}\rightarrow\Re$ generate a sequence of iterates $\{x_{k}\}$ such that $f$ is strictly decreasing at each iterate. Crucially, at each iteration a quasi-Newton matrix is used to approximate $\nabla^{2}f(x_{k})$ that is assumed to be either too computationally expensive to compute or unavailable. The approximation to the Hessian is updated each iteration using the most recently-computed iterate $x_{k+1}$ by defining a new quasi-Newton pair $(s_{k},y_{k})$ given by

[TABLE]

The quasi-Newton Broyden family of updates is given by

[TABLE]

where $\phi_{k}\in\Re$ and

[TABLE]

For $\phi_{k}\in[0,1]$, $B_{k+1}$ is said to be in the restricted or convex Broyden class of updates. Setting $\phi=0$ gives the Broyden-Fletcher-Goldfarb-Shanno (BFGS) update, arguably the most widely-used symmetric positive-definite update and a member of the restricted Broyden class. For $\phi\not\in[0,1]$, the sequence of quasi-Newton matrices generated by this update is not guaranteed to be positive definite. The most well-known update not in the restricted Broyden class is the symmetric rank-one (SR11) update, which is obtained by setting $\phi_{k}=(s_{k}^{T}y_{k})/(s_{k}^{T}y_{k}-s_{k}^{T}B_{k}s_{k})$.

Recently, there has been renewed interest in the entire Broyden class of updates, and in particular, in negative values of $\phi$. Research has shown that negative values of $\phi$ are desirable [7] and under some conditions, quasi-Newton methods based on negative values of $\phi$ exhibit superlinear convergence rates [7, 16]. There has also been empirical evidence that $\phi<0$ may lead to more efficient algorithms than BFGS [16, 14].

In this paper, we present the compact representation for the full Broyden class of quasi-Newton matrices, allowing $\phi$ to be negative and to change each iteration. We also demonstrate how to efficiently solve linear systems with any member of the Broyden class using the compact representation of its inverse. This paper can be viewed as an extension of the results found in [10, 11], which presented the compact representation for members of the restricted Broyden class and their inverses, as well as a practical method for solving linear systems involving only restricted Broyden class matrices (i.e., $\phi\in[0,1]$).

One important application of the compact representation is the ability to efficiently compute the eigenvalues and a partial eigenbasis when the number of stored pairs is small [10], which is the case in large-scale optimization with so-called limited-memory quasi-Newton updates. In this setting, only the most recently-computed $M$ quasi-Newton pairs $\{(s_{k},y_{k})\}$, $k=0,1,\ldots,M-1$, are stored and used to update $B_{k+1}$ using the recursive application of (1). Typically, in large-scale applications $M\leq 10$ regardless of $n$, i.e., $M\ll n$ (see, e.g.,  [8]). With the eigenvalues it is now possible to compute condition numbers, compute singular values, and perform sensitivity analysis.

This paper is organized in seven sections. In the second section, we review the compact formulation for the restricted Broyden class of updates ($\phi\in[0,1]$) as well as overview the efficient computation of their eigenvalues. The main result of the paper is in Section 3 where the compact representation is given for the entire Broyden class of updates that allows for $\phi$ to change each update. In this section, we also present a practical iterative method to compute the compact representation. In Section 4, we show how to perform linear solves with any member of the Broyden class using the compact representation of their inverse. Numerical experiments are reported in Section 5. Finally, Section 6 contains concluding remarks, and Section 7 includes acknowledgements for this work.

1.1. Notation and assumptions

Throughout this paper, we make use of the following matrices:

[TABLE]

Furthermore, we make use of the following decomposition of $S_{k}^{T}Y_{k}\in\Re^{(k+1)\times(k+1)}$:

[TABLE]

where $L_{k}$ is strictly lower triangular, $D_{k}$ is diagonal, and $R_{k}$ is strictly upper triangular. We assume that the matrix $B_{k}$ is nonsingular for each $k$. Finally, throughout the manuscript, $I_{j}$ denotes the $j\times j$ identity matrix.

2. Compact representation for the restricted Broyden class

Compact representations of matrices from the Broyden class of updates were first described by Byrd et al [8] as matrix decompositions of the form

[TABLE]

where $\Psi_{k}\in\Re^{n\times l}$, $M_{k}\in\Re^{l\times l}$, and $B_{0}$ is the initial matrix. The size of $l$ depends on the rank of the update; in the case of a rank-two update, $l=2(k+1)$, and in the case of a rank-one update, $l=k+1$. In the case of the BFGS update (i.e., $\phi=0$), $\Psi_{k}$ and $M_{k}$ are given in [8]:

[TABLE]

where $S_{k}$ and $Y_{k}$ are defined in (2). In [10], we presented the compact representation for any matrix in the restricted Broyden class (i.e, $\phi\in[0,1]$); in particular, for any matrix in the restricted Broyden class,

[TABLE]

where $L_{k}$ and $D_{k}$ are given in (4) and $\Lambda_{k}\in\Re^{(k+1)\times(k+1)}$ is the diagonal matrix $\Lambda_{k}=\text{diag}(\lambda_{i}),$ ($0\leq i\leq k$), given by

[TABLE]

To our knowledge, the only compact formulation known for a member of the Broyden class of updates outside the restricted class is for an SR11 matrix. As with the BFGS case, it is also given in [8]; in particular,

[TABLE]

Notice that $\Psi_{k}$ in the compact representation for SR11 matrices is half the size of that of $\Psi_{k}$ for the rank-two updates.

2.1. Applications of the compact representation

In this section, we briefly review how the eigenvalues of any quasi-Newton matrix that exhibits a compact representation can be efficiently computed. The first method to compute eigenvalues of limited-memory quasi-Newton matrices was proposed by Lu [15]. This method makes use of the singular value decomposition and an eigendecomposition of small matrices. An alternative approach, first described by Burdakov et al. [5], uses the QR factorization in lieu of the singular value decomposition. An overview of the method found in [5] follows below. For this section, we assume $k$ is small, as in the case of limited-memory quasi-Newton matrices; moreover, we assume $\Psi_{k}\in\Re^{n\times l}$ is full rank, where $l$ is either $l=2(k+1)$ or $l=k+1$. Finally, we assume $B_{0}=\gamma I$, where $\gamma\in\Re$.

Let $QR$ be the “thin” QR decomposition of $\Psi_{k}$, where $Q\in\Re^{n\times l}$ has orthonormal columns and $R\in\Re^{l\times l}$ is upper triangular (see, e.g., [12]). Then,

[TABLE]

The matrix $RM_{k}R^{T}$ is a real symmetric $l\times l$ matrix, whose spectral decomposition can be explicitly computed since $l$ is small. Letting $VDV^{T}$ be its spectral decomposition gives that

[TABLE]

where $\hat{D}$ is a diagonal matrix whose leading $l\times l$ block is $D$ while the rest of the matrix is zeros. Thus, the spectral decomposition of $B_{k+1}$ is given by (7). (Note that in practice, the matrices $Q$ and $V$ in (7) are not stored.) Note that the matrix $B_{k+1}$ has an eigenvalue of $\gamma$ with multiplicity $n-l$ and $l$ eigenvalues given by $\gamma+d_{i}$, where $1\leq i\leq l$ and $D_{i,i}\mathop{\,{\scriptstyle\stackrel{{\scriptstyle\triangle}}{{=}}}}\,d_{i}$. It also turns out that it is also possible to efficiently compute the eigenvectors associated with the nontrivial eigenvalues and only one eigenvector associated with the trivial eigenvalue $\gamma$. (For more details, see [5, 10].)

Generally speaking, computing the eigenvalues of $B_{k+1}$ directly is an $O(n^{3})$ process. In contrast, the above decomposition requires the QR factorization of $\Psi_{k}$ and the eigendecomposition of $RM_{k}R^{T}$, requiring $O(nl^{2})$ flops and $O(l^{3})$ flops, respectively. Since $l\ll n$, the proposed method’s runtimes should increase only linearly with $n$. (For some details regarding updating the (full) QR factorization after a new quasi-Newton pair is computed, see [10].) This efficient computation of eigenvalues and a partial eigenbasis appears in new methods for large-scale optimization [5, 4, 2, 1, 3].

The compact representation is also useful for solving linear systems with quasi-Newton matrices. In [6], Burke et al. use the compact formulation of a BFGS matrix to solve a linear system involving a diagonally-shifted BFGS matrix. In [11], the compact representation for the inverse of any member in the restricted Broyden class is given as well as a practical method to solve linear systems involving these matrices using this representation.

3. Compact Representation for any member of the Broyden class

The main result for this section is a theorem giving the compact representation for any member of the Broyden class. The representation allows $\phi$ to change each iteration and to be negative. In this section, we also present a practical algorithm for computing the compact representation.

We begin by observing that $B_{k+1}$ in (1) can be written as

[TABLE]

where

[TABLE]

We now state two lemmas about $O_{k}$; specifically, we provide the condition for $O_{k}$ when is singular as well as its inverse when it is nonsingular.

Lemma 1. The $2\times 2$ matrix $O_{k}$ is singular if and only if $\phi_{k}=(s_{k}^{T}y_{k})/(s_{k}^{T}y_{k}-s_{k}^{T}B_{k}s_{k})$.

Proof. The determinant of $O_{k}$ is given by

[TABLE]

Thus, $O_{k}$ is singular if and only if

[TABLE]

in other words, $\phi_{k}=(y_{k}^{T}s_{k})/(y_{k}^{T}s_{k}-s_{k}^{T}B_{k}s_{k})$. $\square$

Lemma 1 states that $O_{k}$ is singular if and only if the SR11 update is used. Special care will given to the SR11 case, since unlike other members of the Broyden class, this is a rank-one update. For the duration of this manuscript, we let $\phi_{k}^{{\small SR1}1}\mathop{\,{\scriptstyle\stackrel{{\scriptstyle\triangle}}{{=}}}}\,(y_{k}^{T}s_{k})/(y_{k}^{T}s_{k}-s_{k}^{T}B_{k}s_{k})$. For $\phi_{k}\neq\phi_{k}^{{\small SR1}1}$, $O_{k}$ is invertible and its inverse is given in Lemma 2. This result can be derived by using the formula for the inverse of a $2\times 2$ matrix.

Lemma 2. If $O_{k}$ is invertible, then

[TABLE]

where $\alpha_{k}=-(1-\phi_{k})/(s_{k}^{T}B_{k}s_{k})$ and $\beta_{k}=-\phi_{k}/(y_{k}^{T}s_{k}).$

We now state the main theorem of this manuscript that presents the compact representation for any member of the Broyden class, while allowing the parameter $\phi$ to vary at each iteration. After proving this theorem, we discuss several aspects of this compact representation as well as the key differences between the compact representation for the Broyden class of matrices (Theorem 1) and the compact representation of the restricted Broyden class reviewed in Section 2.

Theorem 1.

Let $\Psi_{k}=(B_{0}S_{k}\ \ Y_{k})\in\Re^{n\times 2(k+1)}$. Let $\Pi_{k}\in\Re^{2(k+1)\times 2(k+1)}$ be the permutation matrix

[TABLE]

with $\Pi_{0}\mathop{\,{\scriptstyle\stackrel{{\scriptstyle\triangle}}{{=}}}}\,I_{2}$. Additionally, let $\Xi_{k}$ be defined recursively as

[TABLE]

with $\Xi_{0}\mathop{\,{\scriptstyle\stackrel{{\scriptstyle\triangle}}{{=}}}}\,E_{0}$. Finally, let $\Gamma_{k}\in\Re^{(k+1)\times(k+1)}$ be a diagonal matrix such that

[TABLE]

If $B_{k+1}$ is a member of the Broyden class of updates, then

[TABLE]

where

[TABLE]

$L_{k}$ and $D_{k}$ are defined in (4), and

[TABLE]

Proof. This proof is by induction on $k$. For the base case ($k=0$), $D_{0}=y_{0}^{T}s_{0}$ with $L_{0}=R_{0}=0$, and $\Gamma_{0}$ is the scalar $\gamma_{0}$. Thus, $\widehat{M}_{0}$ defined in (15) reduces to

[TABLE]

By (8), $B_{1}$ is given by $B_{1}=B_{0}+\Psi_{0}M_{0}\Psi_{0}^{T}$ where $\Psi_{0}=(B_{0}s_{0}\ \ y_{0})$ and

[TABLE]

It remains to show that $M_{0}=\Pi_{0}^{T}\Xi_{0}\widehat{M}_{0}\Xi_{0}^{T}\Pi_{0}$. Since the initial permutation matrix is defined as $\Pi_{0}=I_{2}$, we only need to show $M_{0}=\Xi_{0}\widehat{M}_{0}\Xi_{0}^{T}.$ For simplicity, $M_{0}$ can be written as

[TABLE]

where

[TABLE]

From Lemma 1, $M_{0}$ is nonsingular if and only if $\phi_{0}\neq\phi_{0}^{{\small SR1}1}$. Thus, we consider the following two cases separately: (a) $\phi_{0}=\phi_{0}^{{\small SR1}1}$ and (b) $\phi_{0}\neq\phi_{0}^{{\small SR1}1}$.

Case (a): If $\phi_{0}=\phi_{0}^{{\small SR1}1}$, then $\Xi_{0}=E_{0}=(-1\ \ 1)^{T}$ by (12). By (10), $\alpha_{0}+\beta_{0}=0$, and thus, $M_{0}$ can be simplified as

[TABLE]

Finally, since $\phi_{0}=(y_{0}^{T}s_{0})/(s_{0}^{T}y_{0}-s_{0}^{T}B_{0}s_{0})$ and $\beta_{0}=-\phi_{0}/y_{0}^{T}s_{0}$, then

[TABLE]

and thus, $M_{0}=\Xi_{0}\widehat{M}_{0}\Xi_{0}^{T}$, as desired.

Case (b): If $\phi_{0}\neq\phi_{0}^{{\small SR1}1}$, then $M_{0}$ is nonsingular and $\Xi_{0}=I_{2}$. Thus, it remains to show $M_{0}=\widehat{M}_{0}$. By Lemma 1, $\alpha_{0}+\beta_{0}\neq 0$, making $\gamma_{0}=\phi_{0}/(\alpha_{0}+\beta_{0})$ well defined. By Lemma 2, the inverse of $M_{0}$ is given by

[TABLE]

Note that the last equality in (21) follows since $\Pi_{0}=I_{2}$.

For the induction step, assume

[TABLE]

where $\widehat{M}_{m-1}=\left(\Xi_{m-1}^{T}\Pi_{m-1}\Omega_{m-1}\Pi_{m-1}^{T}\Xi_{m-1}\right)^{-1}$ and

[TABLE]

From (8), we have

[TABLE]

where

[TABLE]

Multiplying (22) by $s_{m}$ on the right, we obtain

[TABLE]

Then, substituting this into (24) yields

[TABLE]

where $p_{m}\mathop{\,{\scriptstyle\stackrel{{\scriptstyle\triangle}}{{=}}}}\,\widehat{M}_{m-1}\widehat{\Psi}_{m-1}^{T}s_{m}$. Equivalently,

[TABLE]

where

[TABLE]

Note that $\mathcal{M}_{m}$ has the following decomposition:

[TABLE]

Thus, $\mathcal{M}_{m}$ is nonsingular if and only if $\alpha_{m}\delta_{m}-\beta_{m}^{2}\neq 0$; that is, $\mathcal{M}_{m}$ is nonsingular if and only if $\phi_{m}\neq\phi_{m}^{{\small SR1}1}$ (see Lemma 1). To complete the induction step, we will show that the last term in (27) is equal to $\widehat{\Psi}_{m}\widehat{M}_{m}\widehat{\Psi}_{m}^{T}$ by considering the following two cases separately: (i) $\phi_{m}=\phi_{m}^{{\small SR1}1}$ and (ii) $\phi_{m}\neq\phi_{m}^{{\small SR1}1}$.

Case (i): If $\phi_{m}=\phi_{m}^{{\small SR1}1}$, then by Lemma 1, $\alpha_{m}=-\beta_{m}\neq 0$. Then

[TABLE]

where

[TABLE]

We now show that $\widetilde{\mathcal{M}}_{m}=\widehat{M}_{m}$. By the inductive hypothesis, $\widehat{M}_{m-1}$ is nonsingular. Together with the fact that $\beta_{m}\neq 0$, it can be checked directly that

[TABLE]

The (2,2)-entry of $\widetilde{\mathcal{M}}_{m}^{-1}$ can be simplified by substituting in for $p_{m}$ and using the inductive step (22):

[TABLE]

Substituting this into (32) and using the inductive hypothesis gives:

[TABLE]

where

[TABLE]

Note that the middle matrix in (33) can be expressed as

[TABLE]

which is equivalent to

[TABLE]

Substituting this into (33), yields

[TABLE]

where $\Pi_{m}$ is defined in (11), replacing $k$ with $m$, and $\Omega_{m}$ is defined in (23), replacing $m-1$ with $m$. Thus, $\widetilde{\mathcal{M}}_{m}=\widehat{M}_{m}$, as defined in (15).

We finish this case of the proof by showing that the last term in (27) is equal to $\widehat{\Psi}_{m}\widehat{M}_{m}\widehat{\Psi}_{m}^{T}$. Substituting in (30) gives that the last term in (27) can be written as

[TABLE]

Using (35) this simplifies to

[TABLE]

or, in other words,

[TABLE]

which is exactly $\widehat{\Psi}_{m}\widehat{M}_{m}\widehat{\Psi}_{m}^{T}$. Thus, for $\phi_{m}=\phi_{m}^{{\small SR1}1}$, the inductive step is proven.

Case (ii): We consider the case that $\phi_{m}\neq\phi_{m}^{{\small SR1}1}$. We begin by showing that $\widehat{M}_{m}=\mathcal{M}_{m}$, given in (28). By Lemma 1, $\alpha_{m}+\beta_{m}\neq 0$. Second, $E_{m}=I_{2}$ (see (12)), and $\gamma_{m}=\phi_{m}/(\alpha_{m}+\beta_{m})$ is well-defined (see (13)). Then, the inverse of $\mathcal{M}_{m}$ can be computed using arguments similar to those found in [10]:

[TABLE]

where

[TABLE]

Simplifying the expressions in (37), yields

[TABLE]

We now simplify the entries of (36) using the same approach as in [10]. Since $p_{m}=\widehat{M}_{m-1}\widehat{\Psi}_{m-1}^{T}s_{m}$, then $\widehat{M}_{m-1}^{-1}p_{m}=\widehat{\Psi}_{m-1}^{T}s_{m}$, giving us an expression for the (1,2) and (2,1) entries. The (2,2) block entry is simplified by first multiplying (25) by $s_{m}^{T}$ on the left to obtain $s_{m}^{T}B_{m}s_{m}=s_{m}^{T}B_{0}s_{m}+p_{m}^{T}\widehat{M}_{m-1}^{-1}p_{m}$. Then,

[TABLE]

Thus, using (34), (36) can be written as

[TABLE]

proving that $\widehat{M}_{m}=\mathcal{M}_{m}$. Finally, using arguments similar to those in case (i), it can be shown that

[TABLE]

as desired. $\square$

There are two main differences in the compact representation for the full Broyden class (Theorem 1) and the restricted Broyden class (Section 2). First, in Theorem 1, $\Xi_{k}$ will always be the identity matrix for updates belonging to the restricted Broyden class. Second, in (14), the permutation matrices in the definitions of $\widehat{M}_{k}$ and $\widehat{\Psi}_{k}$, (equations (15) and (16), respectively) always cancel out in the restricted Broyden case. To emphasize that the permutation matrices do not cancel out for the general Broyden class updates, we use the notation $\widehat{M}_{k}$ and $\widehat{\Psi}_{k}$, in lieu of $M_{k}$ and $\Psi_{k}$ as in the restricted Broyden case.

Finally, we provide some insight regarding the permutation matrices (11). The permutation matrix $\Pi_{k}$ acts in the following manner:

[TABLE]

so that

[TABLE]

In other words, when applied on the right of $\Psi_{k}$, the product $\Pi_{k}^{T}\Xi_{k}$ permutes the columns of $\Psi_{k}$ and, using the matrices $\{E_{i}\}$, combines columns of $\Psi_{k}$ whenever the update is a rank-one update.

Unfortunately, computing $\widehat{M}_{k}$ is not straightforward. In particular, the diagonal matrix $\Gamma_{k}$ in Eq. (13) involves $s_{i}^{T}B_{i}s_{i}$ for each $i\in\{0,\ldots,k\}$, which requires $B_{i}$ for $0\leq i\leq k$. In the next section, we propose a recursive method for computing $\widehat{M}_{k}$ that does not require storing the matrices $B_{i}$ for $0\leq i\leq k$.

3.1. Computing $\widehat{M}_{k}$

In this section, we propose a recursive method for computing $\widehat{M}_{k}$ from $\widehat{M}_{k-1}$. This method is based on the method proposed in [11] for solving a linear system whose system matrix is generated using the restricted Broyden class of updates. In the proof of Theorem 1, we showed that

[TABLE]

which are given in (29) and (31). We now relate some of the entries in $\widehat{M}_{k}$ with other stored or computable quantities involving the pairs $\{s_{i},y_{i}\},i=0,\ldots,k$. The vector $p_{k}$ can be computed as

[TABLE]

Note that in (42), the vector $S_{k-1}^{T}B_{0}s_{k}$ is the first $k-1$ entries in the last column of $S_{k}^{T}B_{0}S_{k}$, and the vector $Y_{k-1}^{T}s_{k}$ is the first $k$ entries in the last column of $Y_{k}^{T}S_{k}$. Moreover, the entry $\alpha_{k}$, given by $\alpha_{k}=-(1-\phi_{k})/s_{k}^{T}B_{k}s_{k}$, can be computed from the following:

[TABLE]

In (43), the quantity $s_{k}^{T}B_{0}s_{k}$ is the $k$th diagonal entry in $S_{k}^{T}B_{0}S_{k}$, and $s_{k}^{T}\widehat{\Psi}_{k-1}p_{k}$ is the inner product of $p_{k}$ and $\widehat{\Psi}_{k-1}^{T}s_{k}$, the latter vector already having been computed in (42). Recall that the entry $\beta_{k}$ is given by $\beta_{k}=-\phi_{k}/y_{k}^{T}s_{k}$, where $y_{k}^{T}s_{k}$ is the $(k+1)$st diagonal entry in $S_{k}^{T}Y_{k}$. Finally, $\delta_{k}=(1+\phi_{k}s_{k}^{T}B_{k}s_{k}/y_{k}^{T}s_{k})/y_{k}^{T}s_{k}$, which uses the previously computed quantities $s_{k}^{T}B_{k}s_{k}$ and $y_{k}^{T}s_{k}$.

For the initialization of $\widehat{M}_{0}$, notice that $\widehat{M}_{0}$ in (17) can be written as

[TABLE]

where $\alpha_{0},\beta_{0},$ and $\delta_{0}$ are defined as in (20).

In Algorithm 1, we use the recursions described above to compute $\widehat{M}_{k}$ given in (41).

Note that the matrices $\Pi_{j-1}$ and $\Xi_{j-1}$ are not explicitly formed in Algorithm 1. Instead, (40) can be used to compute $\widehat{\Psi}_{j-1}$ in line 4 of Algorithm 1.

4. Solving linear systems

Given the compact representation of $B_{k+1}$, we can solve

[TABLE]

where $r,z\in\Re^{n}$, by computing the compact representation of the inverse of $B_{k+1}$. Intuitively speaking, computing the compact representation of the inverse is due to the fact that $H_{k+1}\mathop{\,{\scriptstyle\stackrel{{\scriptstyle\triangle}}{{=}}}}\,B_{k+1}^{-1}$ can also be written using a recursion relation [9]:

[TABLE]

where $H_{k}\mathop{\,{\scriptstyle\stackrel{{\scriptstyle\triangle}}{{=}}}}\,B_{k}^{-1}$, $v_{k}=s_{k}/(y_{k}^{T}s_{k})-(H_{k}y_{k})/(y_{k}^{T}H_{k}y_{k})$, and

[TABLE]

Note that when $\phi_{k}=\phi_{k}^{{\small SR1}1}$, then the corresponding $\Phi_{k}$ is given by

[TABLE]

In this section, we derive the compact representation of the inverse of a Broyden class member. This derivation is similar to the process of finding the inverse of a member of the restricted Broyden class presented in [11].

Applying the Sherman-Morrison-Woodbury formula (see, e.g., [13]) to the compact representation of $B_{k+1}$ given in (14), gives that

[TABLE]

For quasi-Newton matrices it is conventional to let $H_{i}$ denote the inverse of $B_{i}$ for each $i$; with this notation, the inverse of $B_{k+1}^{-1}$ is given by

[TABLE]

Using the definition of $\widehat{\Psi}_{k}$ in (16) gives that

[TABLE]

and thus,

[TABLE]

Substituting $H_{0}\widehat{\Psi}_{k}=H_{0}(B_{0}S_{k}\ Y_{k})\Pi_{k}^{T}\Xi_{k}=(S_{k}\ H_{0}Y_{k})\Pi_{k}^{T}\Xi_{k}$ into (48) gives the compact representation for the inverse of any member of the full Broyden class:

[TABLE]

where $\widetilde{\Psi}_{k}=(S_{k}\ H_{0}Y_{k})\Pi_{k}^{T}\Xi_{k}$ and

[TABLE]

Computing $\widetilde{M}_{k}$. Using an approach similar to how $\widehat{M}_{k}$ is computed, $\widetilde{M}_{k}$ can be computed as follows:

[TABLE]

where

[TABLE]

and $\tilde{p}_{k}=\widetilde{M}_{k-1}\widetilde{\Psi}_{k-1}^{T}y_{k}$. The initial matrix $\widetilde{M}_{0}$ is given by the following:

[TABLE]

where $\tilde{\alpha}_{0},\tilde{\beta}_{0},$ and $\tilde{\delta}_{0}$ are defined as in (52) with $k=0$. A practical iterative method to solve equations of the form (45) is given in Algorithm 2.

5. Numerical experiments

In this section we test the accuracy of Algorithm 1 to compute the compact representation by comparing it with the matrix obtained using the Broyden update formula (1). In addition, we demonstrate that solves with $B_{k+1}$ in (45) can be done efficiently using Algorithm 2 with respect to both accuracy and time. For these experiments, we used five (limited-memory) quasi-Newton pairs to compute $B_{k+1}$. To generate quasi-Newton pairs, we simulated a line-search method where the iterates are updated as follows:

[TABLE]

where $\alpha_{j}\in[0,1]$ was generated randomly. To initialize the process, we randomly generated initial points $x_{0}$ and $x_{1}$ so that $s_{0}=x_{1}-x_{0}$. The corresponding gradients, $g_{j}=\nabla f(x_{j})$ for $0\leq j\leq 5$, were also generated randomly in order to form $y_{j}=g_{j+1}-g_{j}$ for $0\leq j\leq 4$. The matrix $B_{0}$ was initially defined as $B_{0}=\gamma I$, where $\gamma>0$ was randomly generated. We considered four experiments where we vary the value of $\phi_{i}$ at each iteration $i$. In particular, we chose values of $\phi_{i}$ according to the scheme given in Table 1. We ran each experiment ten times with $n=10,100,$ and $1,000$ and report results.

5.1. Accuracy of the compact representation

To test the accuracy of the compact representation, we form each $B_{k+1}$ using (14) together with the proposed compact formulation given in Theorem 1. (In particular, we use Algorithm 1 to form $\widehat{M}_{k}$.) We denote the resulting matrix by $B_{k+1}^{\text{CR}}$. In Table 2, we report the average relative error of the compact representation in the Frobenius norm:

[TABLE]

where $B_{k+1}$ is computed using (1).

The small relative errors in Table 2 reflects the fact that the proposed compact representation for the full Broyden class of quasi-Newton matrices is correct; moreover, the relative errors suggest that Algorithm 1 provides a method to compute the compact representation to high accuracy.

5.2. Accuracy of the compact representation of the inverse

In these experiments, we test the accuracy of Algorithm 2 to solve linear systems of the form $B_{k+1}r=z$, where $r,z\in\Re^{n}$ and $B_{k+1}$ is a quasi-Newton matrix. The matrix $B_{k+1}$ is generated using five quasi-Newton pairs as described in the beginning of this section. Moreover, the righthand side $z$ is randomly generated for each experiment. In Table 3, we present the average residual error using the two-norm:

[TABLE]

where $r^{\text{ICR}}$ is the solution to $B_{k+1}r=z$ using the inverse compact representation computed by Algorithm 2. These results suggest that the compact representation of the inverse can be used to solve linear systems to high accuracy.

In addition, during the experiments, the computational time of the proposed method was recorded and compared to a similar solve using the MATLAB “backslash”. In particular, with the same quasi-Newton pairs, the backslash command was used to solve $B_{k+1}r=z$, where $B_{k+1}$ was formed using (1). The times required were averaged for each experiment and for each value of $n$. These results are given in Table 4 and do not include the time MATLAB required to form $B_{k+1}$. Note that the average computational times in Table 4 indicate that as $n$ increases using Algorithm 2 becomes significantly less computationally expensive than using the backslash command.

6. Conclusion

We derived the compact formulation for members of the full Broyden class of quasi-Newton updates. The compact representation allows for different $\phi_{k}$ at each iteration as well as different ranks of updates. With this compact formulation, we demonstrated how to solve linear systems defined by these limited-memory quasi-Newton matrices. Numerical results suggest that the compact representation can be computed to high accuracy and that we can solve (45) efficiently and accurately using the compact representation of the inverse of $B_{k+1}$. Future work includes integrating this linear solver inside large-scale optimization methods.

7. Acknowledgments

The authors would like to thank Lasith Adhikari and Johannes Brust for helpful discussions regarding this work. This research is supported by NSF grants CMMI-1334042 and CMMI-1333326.