Limited-Memory BFGS with Displacement Aggregation

TL;DR

This paper introduces a displacement aggregation strategy for limited-memory BFGS that achieves full-memory convergence properties with fewer stored pairs, improving optimization efficiency and performance.

Contribution

It proposes a novel displacement aggregation method that enables L-BFGS to attain superlinear convergence without storing all curvature pairs.

Findings

Displacement aggregation achieves full-memory BFGS convergence properties.

Numerical results show improved performance over standard L-BFGS.

The method guarantees local superlinear convergence under common assumptions.

Abstract

A displacement aggregation strategy is proposed for the curvature pairs stored in a limited-memory BFGS (a.k.a. L-BFGS) method such that the resulting (inverse) Hessian approximations are equal to those that would be derived from a full-memory BFGS method. This means that, if a sufficiently large number of pairs are stored, then an optimization algorithm employing the limited-memory method can achieve the same theoretical convergence properties as when full-memory (inverse) Hessian approximations are stored and employed, such as a local superlinear rate of convergence under assumptions that are common for attaining such guarantees. To the best of our knowledge, this is the first work in which a local superlinear convergence rate guarantee is offered by a quasi-Newton scheme that does not either store all curvature pairs throughout the entire run of the optimization algorithm or store an explicit (inverse) Hessian approximation. Numerical results are presented to show that displacement aggregation within an adaptive L-BFGS scheme can lead to better performance than standard L-BFGS.

Tables

Table 1. Table 1 : Numbers of iterations, function evaluations, and aggregations when algorithms are applied to solve problems from the CUTEst set with n ∈ [ 10 , 3000 ] 𝑛 10 3000 n\in[10,3000] .

		AggBFGS(5)			LBFGS(5)
Name	$n$	Iters.	Funcs.	Aggs.	Iters.	Funcs.
ARGLINA	200	3	12	1	3	12
ARGLINB	200	3	3	1	3	3
ARGLINC	200	3	3	1	3	3
ARGTRIGLS	200	1494	10220	0	1494	10220
ARWHEAD	1000	2	2	0	2	2
BA-L1LS	57	60	218	0	60	218
BA-L1SPLS	57	62	248	0	62	248
BDQRTIC	1000	160	734	28	399	2051
BOX	1000	16	76	6	---	---
BOXPOWER	1000	20	86	11	12	29
BROWNAL	200	3	4	0	3	4
BROYDN3DLS	1000	39	45	0	39	45
BROYDN7D	1000	1444	2996	0	1444	2996
BROYDNBDLS	1000	71	151	0	71	151
BRYBND	1000	71	151	0	71	151
CHAINWOO	1000	563	2791	11	506	2577
CHNROSNB	50	186	355	0	186	355
CHNRSNBM	50	533	1593	0	533	1593
COSINE	1000	112	467	0	112	467
CRAGGLVY	1000	247	1189	10	241	1153
CURLY10	1000	---	---	---	---	---
CURLY20	1000	---	---	---	---	---
CURLY30	1000	---	---	---	---	---
DIXMAANA	300	16	48	10	20	60
DIXMAANB	300	22	78	6	24	97
DIXMAANC	300	20	106	5	18	78
DIXMAAND	300	25	93	6	35	135
DIXMAANE	300	400	1865	5	378	1798
DIXMAANF	300	383	1906	6	332	1728
DIXMAANG	300	264	1428	8	319	1526
DIXMAANH	300	425	2148	34	526	2668
DIXMAANI	300	1763	9092	33	2054	10878
DIXMAANJ	300	307	1380	6	406	2006
DIXMAANK	300	333	1589	15	331	1618
DIXMAANL	300	302	1553	29	363	1860
DIXMAANM	300	2901	14233	22	5220	6061
DIXMAANN	300	701	3749	15	1623	8508
DIXMAANO	300	694	3423	48	720	3642
DIXMAANP	300	555	2662	52	677	3421
DIXON3DQ	1000	5247	5978	0	5247	5978
DMN15102LS	66	206	962	44	1835	13002
DMN15103LS	99	853	4535	111	1696	9250
DMN15332LS	66	1253	6574	194	978	5241
DMN15333LS	99	11244	39107	26	4339	16755
DMN37142LS	66	639	3314	105	710	3744
DMN37143LS	99	12	128	0	12	128
DQDRTIC	1000	40	191	21	23	75
DQRTIC	1000	311	1838	6	511	2973
EDENSCH	36	110	547	11	138	673
EG2	1000	5	5	3	5	5
EIGENALS	110	661	1898	0	661	1898
EIGENBLS	110	1239	3765	0	1239	3765
EIGENCLS	462	1898	5431	0	1898	5431
ENGVAL1	1000	50	74	1	50	74
ERRINROS	50	900	4597	92	3648	16633
ERRINRSM	50	662	3225	60	---	---
EXTROSNB	1000	301	1231	1	307	1350
FLETBV3M	1000	92	451	3	114	518
FLETCBV2	1000	1014	1022	0	1014	1022
FLETCBV3	1000	---	---	---	---	---
FLETCHBV	1000	---	---	---	---	---
FLETCHCR	1000	176	1625	0	176	1625
FMINSRF2	961	374	429	0	374	429
FMINSURF	961	245	309	0	245	309
FREUROTH	1000	23	64	0	23	64
GENHUMPS	1000	3368	9766	0	3368	9766
GENROSE	500	2299	13271	0	2299	13271
HILBERTA	10	31	145	9	55	240

Equations

$$\underset{x\in {\mathbb{R}}^n}{\min }f(x),$$

$$d_{k+1}\leftarrow \arg \underset{d\in {\mathbb{R}}^n}{\min }{m}_{k+1}(d),\text{where}{m}_{k+1}(d)=f_{k+1}+g_{k+1}^{T}d+\frac{1}{2}{d}^T{M}_{k+1}d$$

$$M_{k+1}{s}_k=g_{k+1}-g_k=:y_k.$$

$$\underset{W\in {\mathbb{S}}^n}{\min }\Vert W-W_k{\Vert }_{\mathcal{M}}s.t.W=W^T\text{and}W{y}_k=s_k,$$

$$S$$

$$\text{and}Y$$

$$\rho ={\left[\begin{array}{ccc}\frac{1}{s_1^T{y}_1}& \cdots & \frac{1}{s_m^T{y}_m}\end{array}\right]}^T\in {\mathbb{R}}_{>0}^m.$$

$$U_j$$

$$V_j$$

$$\text{and}\overline{W}$$

$$R\leftarrow \left[\begin{array}{ccc}{s}_1^T{y}_1& \cdots & s_1^T{y}_m\\ & \ddots & ⋮\\ & & s_m^T{y}_m\end{array}\right]\text{and}D\leftarrow \left[\begin{array}{ccc}{s}_1^T{y}_1& & \\ & \ddots & \\ & & s_m^T{y}_m\end{array}\right].$$

$$\overline{W}\leftarrow W+\left[\begin{array}{cc}S& WY\end{array}\right]\left[\begin{array}{cc}{R}^{-T}(D+Y^{T}WY)R^{-1}& -R^{-T}\\ -R^{-1}& 0\end{array}\right]\left[\begin{array}{c}{S}^T\\ Y^{T}W\end{array}\right].$$

$$d_k\leftarrow -W_k{g}_k\text{and}{x}_{k+1}\leftarrow x_k+{\alpha }_k{d}_k,$$

$$\tilde{S}$$

$$\text{and}\tilde{Y}$$

$$W_{1:j+1}=U_{j+1}^T{U}_j^T{W}_{1:j-1}{U}_j{U}_{j+1}+U_{j+1}^T{V}_j{U}_{j+1}+V_{j+1}.$$

$$U_j{U}_{j+1}$$

$$=$$

$$V_j{U}_{j+1}$$

$$=\left(\frac{1}{\tau s_{j+1}^T{y}_j}\right){\tau }^2{s}_{j+1}{s}_{j+1}^T(I-{\rho }_{j+1}{y}_{j+1}{s}_{j+1}^T)$$

$$=\left(\frac{1}{s_{j+1}^T{y}_j}\right)\tau \left(s_{j+1}{s}_{j+1}^T-{\rho }_{j+1}{s}_{j+1}(s_{j+1}^T{y}_{j+1})s_{j+1}^T\right)=0.$$

$$W_{1:j+1}=U_{j+1}^T{W}_{1:j-1}{U}_{j+1}+V_{j+1},$$

$$f(x_{(1)},x_{(2)})=100(x_{(2)}-x_{(1)}^2{)}^2+(1-x_{(1)}{)}^2.$$

$$(s_0,y_0)\in {\mathbb{R}}^n\times {\mathbb{R}}^n\text{with}{\rho }_0=1/s_0^T{y}_0>0$$

$$s_0=S_{1:m}\tau \text{for some}\tau \in {\mathbb{R}}^m.$$

$${\tilde{Y}}_{1:m}:=\left[\begin{array}{ccc}{\tilde{y}}_1& \cdots & {\tilde{y}}_m\end{array}\right]$$

$$\text{BFGS}(W,S_{0:m},Y_{0:m})=\text{BFGS}(W,S_{1:m},{\tilde{Y}}_{1:m}).$$

$${\tilde{Y}}_{1:m}=W^{-1}{S}_{1:m}\left[\begin{array}{cc}A& 0\end{array}\right]+y_0{\left[\begin{array}{c}b\\ 0\end{array}\right]}^T+Y_{1:m}$$

$${\tilde{R}}_{1:m}:=\left[\begin{array}{ccc}{s}_1^T{\tilde{y}}_1& \cdots & s_1^T{\tilde{y}}_m\\ & \ddots & ⋮\\ & & s_m^T{\tilde{y}}_m\end{array}\right]\text{and}{\tilde{D}}_{1:m}:=\left[\begin{array}{ccc}{s}_1^T{\tilde{y}}_1& & \\ & \ddots & \\ & & s_m^T{\tilde{y}}_m\end{array}\right],$$

$$R_{0:m}:=\left[\begin{array}{cc}{s}_0^T{y}_0& s_0^T{Y}_{1:m}\\ & R_{1:m}\end{array}\right]\text{and}{D}_{0:m}:=\left[\begin{array}{cc}{s}_0^T{y}_0& \\ & D_{1:m}\end{array}\right]$$

$$\text{BFGS}(W,S_{0:m},Y_{0:m})$$

$$\equiv$$

$$=$$

Full text

Limited-Memory BFGS with Displacement Aggregation

Albert S. Berahas E-mail: albertberahas@gmail.com Department of Industrial and Systems Engineering, Lehigh University

Frank E. Curtis E-mail: frank.e.curtis@gmail.com Department of Industrial and Systems Engineering, Lehigh University

Baoyu Zhou E-mail: baoyu.zhou@lehigh.edu Department of Industrial and Systems Engineering, Lehigh University

A displacement aggregation strategy is proposed for the curvature pairs stored in a limited-memory BFGS (a.k.a. L-BFGS) method such that the resulting (inverse) Hessian approximations are equal to those that would be derived from a full-memory BFGS method. This means that, if a sufficiently large number of pairs are stored, then an optimization algorithm employing the limited-memory method can achieve the same theoretical convergence properties as when full-memory (inverse) Hessian approximations are stored and employed, such as a local superlinear rate of convergence under assumptions that are common for attaining such guarantees. To the best of our knowledge, this is the first work in which a local superlinear convergence rate guarantee is offered by a quasi-Newton scheme that does not either store all curvature pairs throughout the entire run of the optimization algorithm or store an explicit (inverse) Hessian approximation. Numerical results are presented to show that displacement aggregation within an adaptive L-BFGS scheme can lead to better performance than standard L-BFGS.

1 Introduction

Quasi-Newton methods—in which one computes search directions using (inverse) Hessian approximations that are set using iterate and gradient displacement information from one iteration to the next—represent some of the most effective algorithms for minimizing nonlinear objective functions.111Quasi-Newton methods offer the ability to update Hessian and/or inverse Hessian approximations, which is why we state inverse parenthetically here. For ease of exposition throughout the remainder of the paper, we often drop mention of the inverse, although in many cases it is the approximation of the inverse, not the Hessian approximation, that is used in practice. The main advantages of such methods can be understood by contrasting their computational costs, storage requirements, and convergence behavior with those of steepest descent and Newton methods. Steepest descent methods require only first-order derivative (i.e., gradient) information, but only achieve a local linear rate of convergence. Newton’s method can achieve a faster (namely, quadratic) rate of local convergence, but at the added cost of forming and factoring second-order derivative (i.e., Hessian) matrices. Quasi-Newton methods lie between the aforementioned methods; they only require gradient information, yet by updating and employing Hessian approximations they are able to achieve a local superlinear rate of convergence. For many applications, the balance between cost/storage requirements and convergence rate offered by quasi-Newton methods makes them the most effective.

Davidon is credited as being the inventor of quasi-Newton methods [14]. For further information on their theoretical properties, typically when coupled with line searches to ensure convergence, see, e.g., [7, 9, 15, 16, 31, 32, 33, 34, 35]. Within the class of quasi-Newton methods, algorithms that employ Broyden–Fletcher–Goldfarb–Shanno (BFGS) approximations of Hessian matrices have enjoyed particular success; see [5, 18, 20, 38]. This is true when minimizing smooth objectives, as is the focus in the aforementioned references, but also when minimizing nonsmooth [4, 11, 12, 23, 26, 27, 39] or stochastic [1, 2, 6, 10, 22, 24, 28, 37, 40] functions.

For solving large-scale problems, i.e., minimizing objective functions involving thousands or millions of variables, limited-memory variants of quasi-Newton methods, such as the limited-memory variant of BFGS [30] (known as L-BFGS), have been successful for various applications. The main benefit of a limited-memory scheme is that one need not store nor compute matrix-vector products with the often-dense Hessian approximation; rather, one need only store a few pairs of vectors, known as curvature pairs, with dimension equal to the number of variables, and one can compute matrix-vector products with the corresponding Hessian approximation with relatively low computational cost. That said, one disadvantage of contemporary limited-memory methods is that they do not enjoy the local superlinear convergence rate properties achieved by full-memory schemes. Moreover, one does not know a priori what number of pairs should be maintained to attain good performance when solving a particular problem. Recent work [3] has attempted to develop an adaptive limited-memory BFGS method in which the number of stored curvature pairs is updated dynamically within the algorithm. However, while this approach has led to some improved empirical performance, it has not been coupled with any strengthened theoretical guarantees.

We are motivated by the following question: Is it possible to design an adaptive limited-memory BFGS scheme that stores and employs only a moderate number of curvature pairs, yet in certain situations can achieve the same theoretical convergence rate guarantees as a full-memory BFGS scheme (such as a local superlinear rate of convergence)? We have not yet answered this question to the fullest extent possible. That said, in this paper, we do answer the following question, which we believe is a significant first step: Is it possible to develop a limited-memory-type222By limited-memory-type BFGS algorithm, we mean one that stores and employs a finite set of curvature pairs rather than an explicit Hessian approximation. BFGS scheme that behaves exactly as a full-memory BFGS scheme—in the sense that the sequence of Hessian approximations from the full-memory scheme is represented through stored curvature pairs—while storing only a number of curvature pairs not exceeding the dimension of the problem? We answer this question in the affirmative by showing how one can use displacement aggregation such that curvature pairs computed by the algorithm may be modified adaptively to capture the information that would be stored in a full-memory BFGS approximation.

A straightforward application of the strategy proposed in this paper might not always offer practical advantages over previously proposed full-memory or limited-memory schemes in all settings. Indeed, to attain full theoretical benefits in some settings, one might have to store a number of curvature pairs up to the number of variables in the minimization problem, in which case the storage requirements and computational costs of our strategy might exceed those of a full-memory BFGS method. That said, our displacement aggregation strategy and corresponding theoretical results establish a solid foundation upon which one may design practically efficient approaches. We propose one such approach, leading to an adaptive L-BFGS scheme that we show to outperform a standard L-BFGS approach on a set of test problems. For this and other such adaptive schemes, we argue that displacement aggregation allows one to maintain more history in the same number of curvature pairs than in previously proposed limited-memory methods.

We conclude the paper with a discussion of how our approach can be extended to Davidon-Fletcher-Powell (DFP) quasi-Newton updating [14], the challenges of extending it to the entire Broyden class of updates [16], and further ways that displacement aggregation can be employed in practically efficient algorithms.

1.1 Contributions

We propose and analyze a displacement aggregation strategy for modifying the displacement pairs stored in an L-BFGS scheme such that the resulting method behaves equivalently to full-memory BFGS. In particular, we show that if a stored iterate displacement vector lies in the span of the other stored iterate displacement vectors, then the gradient displacement vectors can be modified in such a manner that the same Hessian approximation can be generated using the modified pairs with one pair (i.e, the one corresponding to the iterate displacement vector that is in the span of the others) being ignored. In this manner, one can iteratively throw out stored pairs and maintain at most a number of pairs equal to the dimension of the problem while generating the same sequence of Hessian approximations that would have been generated in a full-memory BFGS scheme. Employed within a minimization algorithm, this leads to an L-BFGS scheme that has the same convergence behavior as full-memory BFGS; e.g., it can achieve a local superlinear rate of convergence. To the best of our knowledge, this is the first work in which a local superlinear convergence rate guarantee is provided for a limited-memory-type quasi-Newton algorithm. We also show how our techniques can be employed within an adaptive L-BFGS algorithm—storing and employing only a small number of pairs relative to the dimension of the problem—that can outperform standard L-BFGS. We refer to our proposed “aggregated BFGS” scheme as Agg-BFGS.

Our proposal of displacement aggregation should be contrasted with the idea of merging information proposed in [25]. In this work, the authors prove—when minimizing convex quadratics only—that one can maintain an algorithm that converges finitely if previous pairs are replaced by linear combinations of pairs. This fact motivates a scheme proposed by the authors—when employing an L-BFGS-type scheme in general—of replacing two previous displacement pairs by a single “sum” displacement pair. In their experiments, the authors’ idea often leads to reduced linear algebra time, since fewer pairs are stored, but admittedly worse performance compared to standard L-BFGS in terms of function evaluations. By contrast, our approach does not use predetermined weights to merge two pairs into one. Our scheme aggregates information stored in any number of pairs using a strategy to ensure that no information is lost, even if the underlying function is not a convex quadratic. Our experiments show that our approach can often result in reduced iteration and function evaluations compared to standard L-BFGS.

1.2 Notation

Let $\mathbb{R}$ denote the set of real numbers (i.e., scalars), let $\mathbb{R}_{\geq 0}$ (resp., $\mathbb{R}_{>0}$) denote the set of nonnegative (resp., positive) real numbers, let $\mathbb{R}^{n}$ denote the set of $n$-dimensional real vectors, and let $\mathbb{R}^{m\times n}$ denote the set of $m$-by-$n$-dimensional real matrices. Let $\mathbb{S}^{n}$ denote the symmetric elements of $\mathbb{R}^{n\times n}$. Let $\mathbb{N}:=\{0,1,2,\dots\}$ denote the set of natural numbers. Let $\|\cdot\|:=\|\cdot\|_{2}$.

We motivate our proposed scheme in the context of solving

[TABLE]

where $f:\mathbb{R}^{n}\to\mathbb{R}$ is continuously differentiable. Corresponding to derivatives of $f$, we define the gradient function $g:\mathbb{R}^{n}\to\mathbb{R}^{n}$ and Hessian function $H:\mathbb{R}^{n}\to\mathbb{S}^{n}$. In terms of an algorithm for solving (1), we append a natural number as a subscript for a quantity to denote its value during an iteration of the algorithm; e.g., for each iteration number $k\in\mathbb{N}$, we denote $f_{k}:=f(x_{k})$. That said, when discussing BFGS updating for the Hessian approximations employed within an algorithm, we often simplify notation by referring to a generic set of curvature pairs that may or may not come from consecutive iterations within the optimization algorithm. In such settings, we clarify our notation at the start of each discussion.

1.3 Organization

In §2, we provide background on BFGS updating. In §3, we motivate, present, and analyze our displacement aggregation approach. The results of numerical experiments with Agg-BFGS are provided in §4. Concluding remarks are given in §5.

2 Background on BFGS

The main idea of BFGS updating can be described as follows. Let the $k$th iterate generated by an optimization algorithm be denoted as $x_{k}$. After an iterate displacement (or step) $s_{k}$ has been computed, one sets the subsequent iterate as $x_{k+1}\leftarrow x_{k}+s_{k}$. Then, in order to determine the subsequent step $s_{k+1}$, one uses the minimizer of the quadratic model $m_{k+1}:\mathbb{R}^{n}\to\mathbb{R}$ given by

[TABLE]

and $M_{k+1}\in\mathbb{R}^{n\times n}$ is a Hessian approximation. In a line search approach, for example, one computes $d_{k+1}$ by minimizing $m_{k+1}$, then computes $s_{k+1}\leftarrow\alpha_{k+1}d_{k+1}=-\alpha_{k+1}M_{k+1}^{-1}g_{k+1}$, where $\alpha_{k+1}$ is chosen by a line search for $f$ from $x_{k+1}$ along $d_{k+1}$. With $f_{k+1}$ and $g_{k+1}$ determined by evaluations at $x_{k+1}$, all that remains toward specifying the model $m_{k+1}$ is to choose $M_{k+1}$. In a quasi-Newton method [14], one chooses $M_{k+1}$ such that it is symmetric—i.e., like the exact Hessian $H_{k+1}$, it is an element of $\mathbb{S}^{n}$—and satisfies the secant equation

[TABLE]

Specifically in BFGS, one computes $W_{k+1}:=M_{k+1}^{-1}$ to solve

[TABLE]

where $\|\cdot\|_{\cal M}$ is a weighted Frobenius norm with weights defined by any matrix ${\cal M}$ satisfying the secant equation (or, for concreteness, one can imagine the weight matrix being the average Hessian between $x_{k}$ and $x_{k+1}$). If one chooses $M_{0}\succ 0$ and ensures that $s_{k}^{T}y_{k}>0$ for all $k\in\mathbb{N}$, as is often done by employing a (weak) Wolfe line search, then it follows that $M_{k}\succ 0$ for all $k\in\mathbb{N}$; see [31].

Henceforth, let us focus on the sequence of inverse Hessian approximations $\{W_{k}\}$. After all, this sequence, not $\{M_{k}\}$, is the one that is often computed in practice since, for all $k\in\mathbb{N}$, the minimizer of $m_{k}$ can be computed as $d_{k}\leftarrow-W_{k}g_{k}$. Moreover, all of our discussions about the inverse Hessian approximations $\{W_{k}\}$ have corresponding meanings in terms of $\{M_{k}\}$ since $W_{k}\equiv M_{k}^{-1}$ for all $k\in\mathbb{N}$.

2.1 Iterative and Compact Forms

As is well known, there are multiple ways to construct or merely compute a matrix-vector product with a BFGS inverse Hessian approximation [31]. For our purposes, it will be convenient to refer to two ways of constructing such approximations: an iterative form and a compact form [8]. For convenience, let us temporarily drop from our notation the dependence on the iteration number of the optimization algorithm and instead talk generically about constructing an inverse Hessian approximation $\mkern 1.5mu\overline{\mkern-1.5muW}\succ 0$. Regardless of whether one employs all displacements pairs since the start of the run of the optimization algorithm (leading to a full-memory approximation) or only a few recent pairs (leading to a limited-memory approximation), the approximation can be thought of as being constructed from some initial approximation $W\succ 0$ and a set of iterate and gradient displacements such that all iterate/gradient displacement inner products are positive, i.e.,

[TABLE]

where

[TABLE]

Algorithm 1 computes the BFGS inverse Hessian approximation from an initial matrix $W\succ 0$ and the displacement pairs in (3). In the context of an optimization algorithm using BFGS updating, this matrix would be set with a single update in each iteration rather than re-generated from scratch in every iteration using historical iterate/gradient displacements. However, we write the strategy in this iterative form for convenience of our analysis. The output $\mkern 1.5mu\overline{\mkern-1.5muW}$ represents the matrix to be employed in the step computation of the (outer) optimization algorithm. As a function of the inputs, we denote the output as $\mkern 1.5mu\overline{\mkern-1.5muW}=\textnormal{BFGS}(W,S,Y)$.

The updates performed in Algorithm 1 correspond to a set of projections and corresponding corrections. In particular, for each $j\in\{1,\dots,m\}$, the update projects out curvature information along the step/direction represented by $s_{j}$, then applies a subsequent correction based on the gradient displacement represented by $y_{j}$; see [13, Appendix B]. In this manner, one can understand the well-known fact that each update in a BFGS scheme involves a rank-two change of the matrix.

Rather than apply the updates iteratively, it has been shown that one can instead construct the BFGS approximation from the initial approximation by combining all low-rank changes directly. The scheme in Algorithm 2, which shows the compact form of the updates, generates the same output as Algorithm 1; see [8].

We note in passing that for computing a matrix-vector product with an BFGS (or L-BFGS) approximation without constructing the approximation matrix itself, it is well-known that one can use the so-called two-loop recursion [30].

2.2 Convergence and Local Rate Guarantees

An optimization algorithm using a BFGS scheme for updating Hessian approximations can be shown to converge and achieve a local superlinear rate of convergence under assumptions that are common for attaining such guarantees. This theory for BFGS relies heavily on so-called bounded deterioration of the updates or the self-correcting behavior of the updates, and, for attaining a fast local rate of convergence, on the ability of the updating scheme to satisfy the well-known Dennis-Moré condition for superlinear convergence. See [15, 16] for further information.

We show in the next section that our Agg-BFGS approach generates the same sequence of matrices as full-memory BFGS. Hence, an optimization algorithm that employs our updating scheme maintains all of the convergence and convergence rate properties of full-memory BFGS. For concreteness, we state one such result as the following theorem, which follows from Theorems 6.5 and 6.6 in [31]. We cite this result later in the paper to support our claim that our Agg-BFGS scheme is a limited-memory-type quasi-Newton approach that can be used to achieve local superlinear convergence for an optimization algorithm.

Theorem 2.1**.**

Suppose that $f$ is twice continuously differentiable and that one employs an algorithm that generates a sequence of iterates $\{x_{k}\}$ according to

[TABLE]

where $\{W_{k}\}$ is generated using the BFGS updating scheme and, for all $k\in\mathbb{N}$, the stepsize $\alpha_{k}\in\mathbb{R}_{>0}$ is computed from a line search (initialized with a unit stepsize) to satisfy the Armijo-Wolfe conditions; see equation (3.6) in [31]. In addition, suppose that the algorithm converges to a point $x_{*}\in\mathbb{R}^{n}$ at which the Hessian is positive definite and Lipschitz continuous. Then, $\{x_{k}\}$ converges to $x_{*}$ at a superlinear rate.

Such a result cannot be proved for L-BFGS [30]. Common theoretical results for L-BFGS merely show that if the pairs employed have $s_{k}^{T}y_{k}$ sufficiently positive and bounded above for all $k\in\mathbb{N}$, then the Hessian approximations are sufficiently positive definite and bounded and one can achieve a local linear rate of convergence, i.e., no better than the rate offered by a steepest descent method.

3 Displacement Aggregation

We begin this section by proving a simple, yet noteworthy result about a consequence that occurs when one makes consecutive BFGS updates with iterate displacements that are linearly dependent. We use this result and other empirical observations to motivate our proposed approach, which is stated in this section. We close this section by proving that our approach is well defined, and we discuss how to implement it in an efficient manner.

3.1 Motivation: Parallel Consecutive Iterate Displacements

The following theorem shows that if one finds in BFGS that an iterate displacement is a multiple of the previous one, then one can skip the update corresponding to the prior displacement and obtain the same inverse Hessian approximation.

Theorem 3.1**.**

Let $(S,Y)$ be defined as in (3), with $\rho$ as in (4), and suppose that $s_{j}=\tau s_{j+1}$ for some $j\in\{1,\dots,m-1\}$ and some nonzero $\tau\in\mathbb{R}$. Then, with

[TABLE]

Algorithm 1 $($and 2$)$ yields $\textnormal{BFGS}(W,S,Y)=\textnormal{BFGS}(W,\tilde{S},\tilde{Y})$ for any $W\succ 0$.

Proof.

Let $W\succ 0$ be chosen arbitrarily. For any $j\in\{1,\dots,m-1\}$, let $W_{1:j}:=\textnormal{BFGS}(W,S_{1:j},Y_{1:j})$ where $S_{1:j}:=\begin{bmatrix}s_{1}&\cdots&s_{j}\end{bmatrix}$ and $Y_{1:j}:=\begin{bmatrix}y_{1}&\cdots&y_{j}\end{bmatrix}$. By (5),

[TABLE]

Since $s_{j}=\tau s_{j+1}$, it follows that

[TABLE]

and that

[TABLE]

Combining these facts with (8), one finds that

[TABLE]

meaning that one obtains the same matrix by applying the updates corresponding to $(s_{j},y_{j})$ and $(s_{j+1},y_{j+1})$ as when one skips the update for $(s_{j},y_{j})$ and only applies the one for $(s_{j+1},y_{j+1})$. The result now follows by the fact that, if one starts with the same initial matrix $W_{1:j+1}$, then applying the updates defined by (5) with the same pairs $\{(s_{j+2},y_{j+2}),\dots,(s_{m},y_{m})\}$ yields the same result. ∎∎

Theorem 3.1 is a consequence of the over-writing process that signifies BFGS. That is, with the update associated with each curvature pair, the BFGS update over-writes curvature information along the direction corresponding to each iterate displacement. The theorem shows that if two consecutive iterate displacement vectors are linearly dependent, then the latter update completely over-writes the former—regardless of the gradient displacement vectors—and as a result the same matrix is derived if the former update is skipped.

What if a previous iterate displacement is not parallel with a subsequent displacement, but is in the span of multiple subsequent displacements? One might suspect that the information in the previous displacement might get over-written and can be ignored. It is not so simple. We illustrate this in two ways.

• •

Suppose that one accumulates curvature pairs $\{(s_{j},y_{j})\}_{j=0}^{k}$ and compares the corresponding BFGS approximations with those generated by an L-BFGS scheme with a history length of $n$. (Hereafter, we refer to the latter scheme as L-BFGS($n$).) Suppose also that for $k>n$ one finds that the latest $n$ iterate displacements span $\mathbb{R}^{n}$. If this fact meant that the updates corresponding to these pairs would over-write all curvature information from previous pairs, then one would find no difference between BFGS and L-BFGS($n$) approximations for $k>n$. However, one does not find this to be the case. In Figure 1(a), we plot the maximum absolute difference between corresponding matrix entries of the BFGS and L-BFGS($2$) inverse Hessian approximations divided by the largest absolute value of an element of the BFGS inverse Hessian approximation when an algorithm is employed to minimize the Rosenbrock function [36]

[TABLE]

The pairs are generated by running an optimization algorithm using BFGS updating with an Armijo-Wolfe line search; the different approximations were computed as side computations. It was verified that for $k>2$ the iterate displacements used to generate the L-BFGS(2) matrices were linearly independent. One finds that the differences between the approximations is large for $k>n$. For the purposes of comparison, we also plot the differences between the BFGS inverse Hessian approximations and those generated by our Agg-BFGS scheme described later; these are close to machine precision.

• •

Another manner in which one can see that historical information contained in a BFGS approximation is not always completely over-written once subsequent iterate displacements span $\mathbb{R}^{n}$ is to compare BFGS approximations with those generated by a BFGS scheme that starts $j$ iterations late. Let us refer to such a scheme as BFGS$(k-j)$. For example, for $j=1$, BFGS$(k-1)$ approximations employ all pairs since the beginning of the run of the algorithm except the first one. By observing the magnitude of the differences in the approximations as $k$ increases, one can see how long the information from the first pair lingers in the BFGS approximation. In Figure 1(b), we plot the maximum absolute difference between corresponding matrix entries of the BFGS and BFGS$(k-j)$ inverse Hessian approximations divided by the largest absolute value of an element of the BFGS inverse Hessian approximation for various values of $j$ using the same pairs generated by the algorithm from the previous bullet. We only plot once the corresponding matrices start to deviate, i.e., for the BFGS and BFGS$(k-j)$ approximations the matrices only start to differ at the $j+1$st iteration. One finds, e.g., that the information from the first pair lingers—in the sense that it influences the BFGS approximation—for iteration numbers beyond $n=2$.

On the other hand, one might expect that the information stored in a full-memory BFGS matrix could be contained in a set of at most $n$ curvature pairs, although not necessarily a subset of the pairs that are generated by the optimization algorithm. For example, this might be expected by recalling the fact that if a BFGS method with an exact line search is used to minimize a strongly convex quadratic function, then the Hessian of the quadratic can be recovered exactly if $n$ iterations are performed; see [16, 25, 31]. Letting $W$ represent a BFGS matrix—perhaps computed from more than $n$ pairs—this means that one could run an auxiliary BFGS method to minimize $x^{T}Wx$ to re-generate $W$ using at most $n$ pairs. By noting the equivalence between the iterates generated by this auxiliary BFGS scheme and the conjugate gradient (CG) method, one finds that the modified iterate displacements would lie in a certain Krylov subspace determined by the initial point and the matrix $W$ [16, 31].

Practically speaking, it would not be computationally efficient to run an entire auxiliary BFGS scheme (for minimizing $W$) in order to capture $W$ using a smaller number of curvature pairs. Moreover, we are interested in a scheme that can be used to reduce the number of curvature pairs that need to be stored even when an iterate displacement lies in the span of, say, only $m<n$ subsequent displacements. Our Agg-BFGS scheme is one efficient approach for achieving this goal.

3.2 Basics of Agg-BFGS

The basic building block of our proposed scheme can be described as follows. Suppose that, in addition to the iterate/gradient displacement information $(S,Y)\equiv([s_{1}\ \cdots\ s_{m}],[y_{1}\ \cdots\ y_{m}])=:(S_{1:m},Y_{1:m})$ as defined in (3) with inner products yielding (4), one also has a previous curvature pair

[TABLE]

such that

[TABLE]

(As in the proof of Theorem 3.1, we have introduced the subscript “$1:m$” for $S$ and $Y$ to indicate the dependence of these matrices on quantities indexed $1$ through $m$. We make use of similar notation throughout the rest of the paper.) Given the linear dependence of the iterate displacement $s_{0}$ on the iterate displacements in $S_{1:m}$, our goal is to determine aggregated gradient displacements

[TABLE]

such that, for a given initial matrix $W\succ 0$, one finds

[TABLE]

That is, our goal is to determine $\tilde{Y}_{1:m}$ such that the matrix $\textnormal{BFGS}(W,S_{0:m},Y_{0:m})$ is equivalently generated by ignoring $(s_{0},y_{0})$ and employing $(S_{1:m},\tilde{Y}_{1:m})$.

Remark 3.1**.**

For simplicity, we are presenting the basics of Agg-BFGS using indices for the iterate and gradient displacement vectors from 0 to $m$. However, these need not correspond exactly to the iterates of the outer optimization algorithm. In other words, our strategy and corresponding theoretical results apply equally for $(S,Y)\equiv([s_{k_{1}}\ \cdots\ s_{k_{m}}],[y_{k_{1}}\ \cdots\ y_{k_{m}}])$ where $\{k_{1},\dots,k_{m}\}$ are some $($not necessarily consecutive$)$ iteration numbers in the outer optimization algorithm. We remark on this generality further when presenting the implementation of our scheme within a complete optimization algorithm; see §3.4.

A basic view of our approach is stated as Algorithm 3, wherein we invoke the compact form of BFGS updates (recall (7)). However, while Algorithm 3 provides an understanding of the intent of our displacement aggregation scheme, it does not specify the manner in which $\tilde{Y}_{1:m}$ yielding (13) can be found. It also does not address the fact that multiple such matrices might exist. Toward a specific scheme, one can expand the compact forms of $\textnormal{BFGS}(W,S_{0:m},Y_{0:m})$ and $\textnormal{BFGS}(W,S_{1:m},\tilde{Y}_{1:m})$, specify that the aggregated displacements have the form

[TABLE]

for some $A\in\mathbb{R}^{m\times(m-1)}$ and $b\in\mathbb{R}^{m-1}$, and compare like terms in order to derive three key equations that ensure that (13) holds. The zero blocks in the variable matrices in (14) are motivated by Theorem 3.1, since in the case of $m=1$ one can set $\tilde{Y}_{1:1}=Y_{1:1}$; otherwise, $A$ and $b$ need to be computed to satisfy (13).

Following these steps, one obtains the more detailed view of our scheme that is stated in Algorithm 4, where the three key equations are stated as (16). Algorithm 4 does not specify a precise scheme for computing $(A,b)$ in order to satisfy (16). We specify such a scheme in §3.4 after first showing in the following subsection that real values for $(A,b)$ always exist to satisfy (16).

3.3 Existence of Real Solutions for Agg-BFGS

Our goal now is to prove the following theorem for our Agg-BFGS scheme.

Theorem 3.2**.**

Suppose one has $W\succ 0$ along with:

• •

$(S_{1:m},Y_{1:m})$* as defined in (3) with the columns of $S_{1:m}$ being linearly independent and the vector $\rho_{1:m}$ as defined in (4), and*

• •

$(s_{0},y_{0},\rho_{0})$* defined as in (10) such that (11) holds for some $\tau\in\mathbb{R}^{m}$.*

Then, there exist $A\in\mathbb{R}^{m\times(m-1)}$ and $b\in\mathbb{R}^{m-1}$ such that, with $\tilde{Y}_{1:m}\in\mathbb{R}^{n\times m}$ defined as in (14), the equations (16) hold. Consequently, for this $\tilde{Y}_{1:m}$, one finds that

[TABLE]

and $\textnormal{BFGS}(W,S_{0:m},Y_{0:m})=\textnormal{BFGS}(W,S_{1:m},\tilde{Y}_{1:m})$.

Proof.

First, observe that there are $(m+1)(m-1)$ unknowns in the formula (14) for $\tilde{Y}_{1:m}$; in particular, there are $m(m-1)$ unknowns in $A$ and $m-1$ unknowns in $b$, yielding $m(m-1)+(m-1)=(m+1)(m-1)$ unknowns in total. To see the number of equations effectively imposed by (16), first notice that (14) imposes $\tilde{y}_{m}=y_{m}$. Hence, by ignoring the last column of $R_{1:m}$ to define the submatrix

[TABLE]

and similarly defining $\tilde{P}$ with size $m\times(m-1)$ as a submatrix of $\tilde{R}_{1:m}$, the key equations in (16) can be simplified to have the following form:

[TABLE]

Observing the number of nonzero entries in (19a) (recall (18)) and recognizing the symmetry in (19c), one finds that the effective number of equations are:

[TABLE]

Hence, (19) (and so (16)) is a square system of linear and quadratic equations to be solved for the unknowns in the matrix $A$ and vector $b$.

Equation (19b) represents a formula for $b\in\mathbb{R}^{m-1}$. Henceforth, let us assume that $b$ is equal to the right-hand side of this equation, meaning that all that remains is to determine that a real solution for $A$ exists. Let us write

[TABLE]

Using this notation, one finds from (14) that the $j$th column of (19a) requires

[TABLE]

hence, (19a) reduces to the system of affine equations

[TABLE]

For each $j\in\{1,\dots,m-1\}$, let us write

[TABLE]

with $a_{j,1}$ having length $j$ and $a_{j,2}$ having length $m-j$. (Notice that $Q$ is positive definite under the conditions of the theorem, which requires that $S_{1:m}$ has full column rank.) Then, in order for (20) to be satisfied, one must have

[TABLE]

Moreover, with this value for $a_{j,1}$, it follows that (20), and hence (19a), is satisfied for any $a_{j,2}$. Going forward, our goal is to show the existence of $a_{j,2}\in\mathbb{R}^{m-j}$ for all $j\in\{1,\dots,m-1\}$ such that (19c) holds, completing the proof.

Observing (19c), one finds with (14) that it may be written as

[TABLE]

where

[TABLE]

One may rewrite equation (23) as

[TABLE]

Let us now rewrite the equations in (25) in a particular form that will be useful for the purposes of our proof going forward. Consider the matrix $QA+\Omega$ in (25). By the definitions of $a_{j,1}$, $a_{j,2}$, $Q$, and $\Omega$ in (21)–(24), as well as of $P$ from (18), the $j$th column of this matrix is given by

[TABLE]

where $0_{j}$ is a vector of zeros of length $j$. Letting $L\in\mathbb{R}^{m\times m}$ be any matrix such that $L^{T}L=Q^{-1}$ (whose existence follows since $Q^{-1}\succ 0$), defining $Z:=[z_{1}\ \cdots\ z_{m-1}]\in\mathbb{R}^{(m-1)\times(m-1)}$ as any matrix such that $Z^{T}Z=\omega\omega^{T}+\Omega^{T}Q^{-1}\Omega$ (whose existence follows since $\omega\omega^{T}+\Omega^{T}Q^{-1}\Omega\succeq 0$), and defining, for all $j\in\{1,\dots,m-1\}$,

[TABLE]

it follows that the $(i,j)\in\{1,\dots,m-1\}\times\{1,\dots,m-1\}$ element of (25) is

[TABLE]

Using the notation from the preceding paragraph, let us use an inductive argument to prove the existence of a real solution of (25). This induction will follow the indices $\{1,\dots,m-1\}$ in reverse order. As a base case, consider the index $m-1$, in which case one has the one-dimensional unknown $a_{m-1,2}$. One finds with

[TABLE]

that

[TABLE]

Hence, letting $a_{m-1,2}=a_{m-1,2}^{*}+\lambda_{m-1}$, where $\lambda_{m-1}$ is one-dimensional, one finds that the left-hand side of the $(i,j)=(m-1,m-1)$ equation in (27) is $\|\phi_{m-1}(a_{m-1,2})\|_{2}^{2}$. This is a strongly convex quadratic in the unknown $\lambda_{m-1}$. Since $\|\phi_{m-1}(a_{m-1,2}^{*})\|_{2}=0$ and $z_{m-1}^{T}z_{m-1}\geq 0$, there exists $\lambda_{m-1}^{*}\in\mathbb{R}$ such that $a_{m-1,2}=a_{m-1,2}^{*}+\lambda_{m-1}^{*}\in\mathbb{R}$ satisfies the $(i,j)=(m-1,m-1)$ equation in (27).

Now suppose that there exists real $\{a_{\ell+1,2},\dots,a_{m-1,2}\}$ such that (27) holds for all $(i,j)$ with $i\in\{\ell+1,\dots,m-1\}$ and $j\in\{i,\dots,m-1\}$. (By symmetry, these values also satisfy the $(j,i)$ elements of (27) for these same values of the indices $i$ and $j$.) To complete the inductive argument, we need to show that this implies the existence of $a_{\ell,2}\in\mathbb{R}^{m-\ell}$ satisfying (27) for all $(i,j)$ with $i\in\{j,\dots,m-1\}$ and $j=\ell$, i.e., solving the following system for $a_{\ell,2}$:

[TABLE]

Notice that (28a)–(28b) are affine equations in $a_{\ell,2}$, whereas (28c) is a quadratic equation in $a_{\ell,2}$. For all $t\in\{\ell+1,\dots,m-1\}$, let

[TABLE]

so that, by (26), one may write

[TABLE]

Our strategy is first to find $a_{\ell,2}^{*}\in\mathbb{R}^{m-\ell}$ satisfying (28a)–(28b) such that

[TABLE]

and (cf. (28c))

[TABLE]

Once this is done, we will argue the existence of a nonzero vector $\mkern 1.5mu\overline{\mkern-1.5mua}_{\ell,2}\in\mathbb{R}^{m-\ell}$ such that $a_{\ell,2}^{*}+\lambda_{\ell}\mkern 1.5mu\overline{\mkern-1.5mua}_{\ell,2}$ satisfies (28a)–(28b) for arbitrary one-dimensional $\lambda_{\ell}$. From here, it will follow by the fact that the left-hand side of (28c) is a strongly convex quadratic in the unknown $\lambda_{\ell}$ and the fact that (30) holds that we can claim that there exists $\lambda_{\ell}^{*}\in\mathbb{R}$ such that $a_{\ell,2}=a_{\ell,2}^{*}+\lambda_{\ell}^{*}\mkern 1.5mu\overline{\mkern-1.5mua}_{\ell,2}\in\mathbb{R}^{m-\ell}$ satisfies (28).

To achieve the goals of the previous paragraph, first let $c$ be the column rank of $[\psi_{\ell+1,\ell+1}\ \cdots\ \psi_{\ell+1,m-1}]$ so that there exists $\{t_{1},\dots,t_{c}\}\subseteq\{\ell+1,\dots,m-1\}$ with

[TABLE]

For completeness, let us first consider the extreme case when $c=0$. In this case,

[TABLE]

Hence, by our induction hypothesis, it follows from (27) and (32) that

[TABLE]

Consequently from (32) and (33), the affine equations (28a)–(28b) are satisfied by any $a_{\ell,2}^{*}\in\mathbb{R}^{m-\ell}$. In particular, one can choose

[TABLE]

and find by (26) that

[TABLE]

which shows that this choice satisfies (30). Now consider the case when $c>0$. For $a_{\ell,2}^{*}$ to satisfy (29), it follows with (31) that we must have

[TABLE]

where $\beta_{\ell}$ has length $c$. Choosing

[TABLE]

it follows with (35) that, for any $t\in\{t_{1},\dots,t_{c}\}$, one finds

[TABLE]

We shall now prove that, for any $t\in\{\ell+1,\dots,m-1\}\backslash\{t_{1},\dots,t_{c}\}$, one similarly finds that $\phi_{t}(a_{t,2})^{T}\phi_{\ell}(a_{\ell,2}^{*})=z_{t}^{T}z_{\ell}$. Toward this end, first notice that for any such $t$ it follows from (31) that $\psi_{\ell+1,t}=[\psi_{\ell+1,t_{1}}\ \cdots\ \psi_{\ell+1,t_{c}}]\gamma_{\ell,t}$ for some $\gamma_{\ell,t}\in\mathbb{R}^{c}$. Combining the relationship (31) along with the inductive hypothesis that, for any pair $(i,j)$ with $i\in\{\ell+1,\dots,m-1\}$ and $j\in\{i,\dots,m-1\}$, one has

[TABLE]

it follows with the positive definiteness of $Q^{-1}=L^{T}L$ that

[TABLE]

From (39), it follows that for any $t\in\{\ell+1,\dots,m-1\}\backslash\{t_{1},\dots,t_{c}\}$ there exists $\bar{\gamma}_{\ell,t}\in\mathbb{R}^{c}$ such that $z_{t}=[z_{t_{1}}\ \cdots\ z_{t_{c}}]\bar{\gamma}_{\ell,t}$. Combining the definitions of $\gamma_{\ell,t}$ and $\bar{\gamma}_{\ell,t}$ along with (38), it follows for any such $t$ that

[TABLE]

from which it follows that $\gamma_{\ell,t}=\bar{\gamma}_{\ell,t}$. Thus, with (37) and the definitions of $\gamma_{\ell,t}$ and $\bar{\gamma}_{\ell,t}$, it follows for any $t\in\{\ell+1,\dots,m-1\}\backslash\{t_{1},\dots,t_{c}\}$ that

[TABLE]

Combining this with (37), it follows that $a_{\ell,2}^{*}$ from (35) with $\beta_{\ell}$ from (36) satisfies (28a)–(28b), as desired. Let us show now that this $a_{\ell,2}^{*}$ also satisfies (30). Indeed, by (35), (36), and (38), it follows that

[TABLE]

where the last inequality comes from the fact that the eigenvalues of

[TABLE]

are all in $\{0,1\}$. (As an aside, one finds that the inequality above is strict if $z_{\ell}\notin\operatorname{span}\{z_{t_{1}},\dots,z_{t_{c}}\}$.) Hence, we have shown that $a_{\ell,2}^{*}$ from (35) satisfies (30).

As previously mentioned (in the text following (30)), our goal now is to show that there exists a nonzero $\mkern 1.5mu\overline{\mkern-1.5mua}_{\ell,2}\in\mathbb{R}^{m-\ell}$ such that $a_{\ell,2}^{*}+\lambda_{\ell}\mkern 1.5mu\overline{\mkern-1.5mua}_{\ell,2}$ satisfies (28a)–(28b) for arbitrary $\lambda_{\ell}$. From (28a)–(28b), such an $\mkern 1.5mu\overline{\mkern-1.5mua}_{\ell,2}\in\mathbb{R}^{m-\ell}$ must satisfy

[TABLE]

Since

[TABLE]

it follows that this matrix has a null space of dimension at least $\ell+1$, i.e., there exist at least $\ell+1$ linearly independent vectors in $\mathbb{R}^{m}$ belonging to the null space of this matrix. Let $N_{\ell+1}\in\mathbb{R}^{m\times(\ell+1)}$ be a matrix whose columns are $\ell+1$ linearly independent vectors in $\mathbb{R}^{m}$ lying in the null space of (41). Since this null space matrix has $\ell+1$ columns, there exists a nonzero vector $\zeta_{\ell+1}\in\mathbb{R}^{\ell+1}$ such that the first $\ell$ elements of $N_{\ell+1}\zeta_{\ell+1}$ are zero. Letting

[TABLE]

one finds that $\mkern 1.5mu\overline{\mkern-1.5mua}_{\ell,2}$ satisfies (40), as desired. Consequently, as stated in the text following (30), it follows by the fact that the left-hand side of (28c) is a strongly convex quadratic in the unknown $\lambda_{\ell}$ and the fact that (30) holds that we can claim that there exists $\lambda_{\ell}^{*}\in\mathbb{R}$ such that $a_{\ell,2}=a_{\ell,2}^{*}+\lambda_{\ell}^{*}\mkern 1.5mu\overline{\mkern-1.5mua}_{\ell,2}\in\mathbb{R}^{m-\ell}$ satisfies (28).

Combining all previous aspects of the proof, we have shown the existence of $A\in\mathbb{R}^{m\times(m-1)}$ and $b\in\mathbb{R}^{m-1}$ such that, with $\tilde{Y}_{1:m}\in\mathbb{R}^{n\times m}$ defined as in (14), the equations (16) hold. The remaining desired conclusions, namely, that (17) holds and that $\textnormal{BFGS}(W,S_{0:m},Y_{0:m})=\textnormal{BFGS}(W,S_{1:m},\tilde{Y}_{1:m})$, follow from the existence of $\tilde{Y}_{1:m}\in\mathbb{R}^{n\times m}$ (that we have proved), the fact that the equations in (17) are a subset of the equations in (19a), and the fact that (16) were derived explicitly to ensure that, with $\tilde{Y}_{1:m}$ satisfying (14), one would find that (13) holds. ∎∎

3.4 Implementing Agg-BFGS

We now discuss how one may implement our Agg-BFGS scheme to iteratively aggregate displacement information in the context of an optimization algorithm employing BFGS approximations. We also discuss the dominant computational costs of applying the scheme, and comment on certain numerical considerations that one should take into account in a software implementation. The procedures presented in this section are guided by our proof of Theorem 3.2.

As will become clear in our overall approach, in contrast to a traditional limited-memory scheme in which one always maintains the most recent curvature pairs to define a Hessian approximation, the pairs used in our approach might come from a subset of the previous iterations, with the gradient displacements potentially having been modified through our aggregation mechanism. For concreteness, suppose that during the course of the run of an optimization algorithm for solving (1), one has accumulated a set of curvature pairs, stored in the sets

[TABLE]

where $\{k_{i}\}_{i=0}^{m-1}\subset\mathbb{N}$ with $k_{i}<k_{i+1}$ for all $i\in\{0,\dots,m-2\}$, such that the vectors in the former set (i.e., the iterate displacements) are linearly independent. (Here, the elements of ${\cal Y}$ are not necessarily the gradient displacements computed in iterations $\{k_{0},\dots,k_{m-1}\}$, but, for simplicity of notation, we denote them as ${\cal Y}$ even though they might have been modified during a previous application of our aggregation scheme.) Then, suppose that a new curvature pair $(s_{k_{m}},y_{k_{m}})$ for $k_{m}\in\mathbb{N}$ with $k_{m-1}<k_{m}$ is available. Our goal in this section is to show how one may add and, if needed, aggregate the information in these pairs in order to form new sets

[TABLE]

such that

1. (i)

the sets $\tilde{\cal S}$ and $\tilde{\cal Y}$ have the same cardinality, which is either $m$ or $m+1$, 2. (ii)

the vectors in the set $\tilde{\cal S}$ are linearly independent, and 3. (iii)

the BFGS inverse Hessian approximation generated from the data in $({\cal S}\cup\{s_{k_{m}}\},{\cal Y}\cup\{y_{k_{m}}\})$ is the same as the approximation generated from $(\tilde{\cal S},\tilde{\cal Y})$.

As in the previous section, we henceforth simplify the subscript notation and refer to the “previously stored” displacement vectors as those in the sets $\{s_{0},\dots,s_{m-1}\}$ and $\{y_{0},\dots,y_{m-1}\}$, and refer to the “newly computed” curvature pair as $(s_{m},y_{m})$.

Once the newly computed pair $(s_{m},y_{m})$ is available, there are three cases.

1. Case 1.

The set of vectors $\{s_{0},s_{1},\dots,s_{m}\}$ is linearly independent. In this case, one simply adds the new curvature pair, which leads to the $(m+1)$-element sets

[TABLE]

Notice that if $m=n$, then this case is not possible. 2. Case 2.

The new iterate displacement vector $s_{m}$ is parallel to the most recently stored iterate displacement vector, i.e., $s_{m-1}=\tau s_{m}$ for some $\tau\in\mathbb{R}$. In this case, one should discard the most recently stored pair and replace it with the newly computed one, which leads to the $m$-element sets

[TABLE]

This choice is justified by Theorem 3.1. 3. Case 3.

For some $j\in\{0,\dots,m-2\}$, an iterate displacement vector $s_{j}$ lies in the span of the subsequent iterate displacements vectors, i.e., $s_{j}\in\operatorname{span}\{s_{j+1},\dots,s_{m}\}$. In this case, one should apply our aggregation scheme to determine the vectors $\{\tilde{y}_{j+1},\dots,\tilde{y}_{m}\}$, then remove the pair $(s_{j},y_{j})$, leading to the $m$-element sets

[TABLE]

This choice is justified by Theorem 3.2.

Computationally, the first step is to determine which of the three cases occurs. One way to do this efficiently is to maintain a Cholesky factorization of an inner product matrix corresponding to the previously stored iterate displacement vectors, then attempt to add to it a new row/column corresponding to the newly computed iterate displacement, checking whether the procedure breaks down. Specifically, before considering the newly computed pair $(s_{m},y_{m})$, suppose that one has a lower triangular matrix $\Theta\in\mathbb{R}^{m\times m}$ with positive diagonal elements such that

[TABLE]

which exists due to the fact that the vectors in $\{s_{0},\dots,s_{m-1}\}$ are linearly independent. A Cholesky factorization of an augmented inner product matrix would consist of a scalar $\mu\in\mathbb{R}_{>0}$, vector $\delta\in\mathbb{R}^{m}$, and lower triangular $\Delta\in\mathbb{R}^{m\times m}$ with

[TABLE]

As is well known, equating terms in (42) and (43) one must have $\mu=\|s_{m}\|$, $\delta^{T}=[s_{m}^{T}s_{m-1}\ \cdots\ s_{m}^{T}s_{0}]/\mu$, and $\Delta\Delta^{T}=\Theta\Theta^{T}-\delta\delta^{T}$, meaning that $\Delta$ can be obtained from $\Theta$ through a rank-one downdate; e.g., see [19]. If this downdate does not break down—meaning that all diagonal elements of $\Delta$ are computed to be positive—then one is in Case 1 and the newly updated Cholesky factorization has been made available when yet another curvature pair is considered (after a subsequent optimization algorithm iteration). Otherwise, if the downdate does break down, then it is due to a computed diagonal element being equal to zero. This means that, for some smallest $i\in\{1,\dots,m\}$, one has found a lower triangular matrix $\Xi\in\mathbb{R}^{i\times i}$ with positive diagonal elements and a vector $\xi\in\mathbb{R}^{i}$ such that

[TABLE]

Letting $\tau\in\mathbb{R}^{i}$ be the unique vector satisfying $\Xi^{T}\tau=\xi$, one finds that the vector $[\tau^{T},-1]^{T}$ lies in the null space of (44), from which it follows that

[TABLE]

where the first element of $\tau$ must be nonzero since $\{s_{m-1},\dots,s_{m-i}\}$ is a set of linearly independent vectors. If the breakdown occurs for $i=1$, then one is in Case 2. If the breakdown occurs for $i>1$, then one is in Case 3 with the vector $\tau$ that one needs to apply our aggregation scheme to remove the pair $(s_{m-i},y_{m-i})$.

Notice that if the breakdown occurs in the rank-one downdate as described in the previous paragraph, then one can continue with standard Cholesky factorization updating procedures in order to have the factorization of

[TABLE]

available in subsequent iterations. For brevity and since it is outside of our main scope, we do not discuss this in detail. Overall, the computational costs so far are ${\cal O}(mn)$ (for computing the inner products $\{s_{m}^{T}s_{m-1},\dots,s_{m}^{T}s_{0}\}$) plus ${\cal O}(m^{2})$ (for updating the Cholesky factorization and, in Case 2 or Case 3, computing $\tau$).

If Case 1 occurs, then no additional computation is necessary; one simply adds the new curvature pair as previously described. Similarly, if Case 2 occurs, then again no additional computation is necessary; one simply replaces the most recent stored pair with the newly computed one. Therefore, we may assume for the remainder of this section that Case 3 occurs, we have identified an index $j$ ($=m-i$ using the notation above) such that $s_{j}\in\operatorname{span}\{s_{j+1},\dots,s_{m}\}$, and we have $\tau\in\mathbb{R}^{m-j}$ such that $s_{j}=S_{j+1:m}\tau$. Our goal then is to apply our aggregation scheme to modify the gradient displacement vectors $\{y_{j+1},\dots,y_{m}\}$ to compute

[TABLE]

where $W_{0:j-1}$ represents the BFGS inverse Hessian approximation defined by some initial positive definite matrix $W\succ 0$ and the curvature pairs $\{s_{i},y_{i}\}_{i=0}^{j-1}$, and where $A\in\mathbb{R}^{(m-j)\times(m-j-1)}$ and $b\in\mathbb{R}^{m-j-1}$ are the unknowns to be determined.

For simplicity in the remainder of this section, let us suppose that $j=0$ so that the pair $(s_{0},y_{0})$ is to be removed as in the notation of §3.2 and §3.3. As one might expect, this is the value of $j$ for which the computational costs of computing $A\in\mathbb{R}^{m\times(m-1)}$ and $b\in\mathbb{R}^{m-1}$ are the highest. For all other values of $j$, there is a cost for computing $W_{0:j-1}^{-1}S_{j+1:m}$ as needed in (45). This matrix can be computed without forming the BFGS Hessian approximation $W_{0:j-1}^{-1}$; it can be constructed, say, by computing matrix-vector products with a compact representation of this approximation for a total cost of ${\cal O}(j(m-j)n)\leq{\cal O}(m^{2}n)$; see §7.2 in [31].

Let us now describe how one may implement our aggregation scheme such that, given $S_{1:m}$ with full column rank, $Y_{1:m}$, $\rho_{1:m}>0$, $\tau\in\mathbb{R}^{m}$ satisfying $s_{0}=S_{1:m}\tau$, $y_{0}$, and $\rho_{0}>0$, one may compute $A\in\mathbb{R}^{m\times(m-1)}$ and $b\in\mathbb{R}^{m-1}$ in order to obtain $\tilde{Y}_{1:m}$ as in (14). By Theorem 3.2, it follows that real values for $A$ and $b$ exist to satisfy (16). The computation of the vector $b$ is straightforward; it can be computed by the formula (19b). Assuming that the products in $S_{1:m}^{T}Y_{1:m-1}$ have already been computed (which, using previously computed inner products and recomputing any as needed if/when gradient displacements have been modified by aggregation, costs ${\cal O}((m-j)^{2}n)$), the cost of computing $b$ is ${\cal O}(m^{2})$.

For computing $A$, let us first establish some notation since the elements of this matrix will be computed with a specific order. As in the proof of Theorem 3.2, let $A=[a_{1}\ \cdots\ a_{m-1}]$ where $a_{\ell}\in\mathbb{R}^{m}$ for all $\ell\in\{1,\dots,m-1\}$, and, as in (21), let

[TABLE]

Here and going forward, our computations require products with $Q^{-1}$. Rather than form this matrix explicitly, one may maintain a Cholesky factorization of $Q$ and add/delete rows/columns—as described previously for the inner product matrix corresponding to the iterate displacements—as the iterate displacement set is updated throughout the run of the (outer) optimization algorithm. With such a factorization, products with $Q^{-1}$ are obtained by triangular solves in a standard fashion. If $W\succ 0$ is diagonal or defined by a limited memory approximation, then the cost of updating this factorization in each instance is ${\cal O}(mn)$ (for computing $W^{-1}s_{m}$) plus ${\cal O}(m^{2})$ for updating the factorization. Each backsolve costs ${\cal O}(m^{2})$.

An approach for computing $A$ is now summarized as Algorithm 5. As explained in the proof of Theorem 3.2, the values of $\{a_{\ell,1}\}_{\ell=1}^{m-1}$ are set to ensure that (19a) is satisfied. Assuming that $S_{1:m}^{T}y_{0}$ has been computed at cost ${\cal O}(mn)$, the cost of computing $\{a_{\ell,1}\}_{\ell=1}^{m-1}$ is ${\cal O}(m^{2})$. The values for $\{a_{\ell,2}\}_{\ell=1}^{m-1}$ are then computed in reverse order to solve the system of linear and quadratic equations comprising (19c). Assuming (as has been mentioned) that a factorization of $Q$ is available, and assuming that the elements of the right-hand side of (25) have been pre-computed (at cost ${\cal O}(m^{3})$), the cost of computing $a_{m-1,2}$ is ${\cal O}(1)$. As for computing the remaining vectors, the presented scheme follows the proof of Theorem 3.2. The most expensive operation in each iteration of this scheme is the QR factorization of the matrix in (41), which for each $\ell$ is of size $(m-(\ell+1))\times m$. Summing the cost of these for $\ell=m-2$ to $\ell=1$, the total cost is found to be ${\cal O}(m^{4})$.

This completes our description of a manner in which our Agg-BFGS scheme can be implemented. Observing the computational costs that have been cited, one finds that a conservative estimate of the total cost is ${\cal O}(m^{2}n)+{\cal O}(m^{4})$. Upon closer inspection, one finds that for $j\approx 0$ the cost is dominated by the cost of computing/updating $S_{1:m}^{T}\tilde{Y}_{1:m}$ after aggregation has been performed, whereas for $j\approx m$ the cost is dominated by the cost of computing $W_{0:j-1}^{-1}S_{j+1:m}$ before aggregation is performed. In either case, if $m\ll n$ (say with $3\leq m\leq 10$, as is typical in practice), then the overall cost is not too dissimilar from ${\cal O}(4mn)$ per iteration, which is the cost of other limited-memory schemes such as L-BFGS [8].

We end this section by stating the following result, for closure.

Theorem 3.3**.**

If one applies Agg-BFGS as described in this section, then one need only store at most $m\leq n$ curvature pairs such that the corresponding BFGS inverse Hessian approximations are equivalent to those in a full-memory BFGS scheme. Consequently, under the same conditions as in Theorem 2.1, the resulting optimization algorithm produces $\{x_{k}\}$ that converges to $x_{*}$ at a superlinear rate.

4 Numerical Demonstrations

Our goals in this section are to provide additional numerical demonstrations of applications of Agg-BFGS. (Recall that a preview demonstration has been provided in Figure 1(a).) Our first goal is to show empirically that limited-memory-type BFGS inverse Hessian approximations provided by Agg-BFGS accurately represent the approximations provided by full-memory BFGS. We demonstrate that while numerical errors might accumulate to some degree as Agg-BFGS is applied over a sequence of iterations, the inverse Hessian approximations provided by Agg-BFGS are not necessarily poor after multiple iterations. Our second goal is to provide the results of numerical experiments with an adaptive L-BFGS method that uses our Agg-BFGS scheme and show that over a diverse test set it can outperform a standard L-BFGS approach. These experiments are run in the practical regime when the number of pairs stored and employed is small relative to $n$.

4.1 Simulated Data

We implemented Agg-BFGS in MATLAB and ran two sets of experiments. First, for $(n,m)\in\{4,8,16,32,64,128\}^{2}$ with $m\leq n$, we generated data to show the error of applying Agg-BFGS to aggregate a single curvature pair. For each $(n,m)$, we generated 100 datasets using the following randomized procedure. First, MATLAB’s built-in sprandsym function was used to generate a random positive definite matrix with condition number approximately $10^{4}$. This matrix defined a quadratic function. Second, a random fixed point was determined using MATLAB’s built-in randn routine. From this point, a mock optimization procedure for minimizing the generated quadratic was run to generate $\{(s_{k},y_{k})\}_{k=1}^{m}$; in particular, for $k\in\{1,\dots,m\}$, starting with the fixed point, a descent direction was chosen as the negative gradient plus noise (specifically, the norm of the gradient divided by 10 times a random vector generated with randn) and a stepsize was chosen by an exact line search to compute the subsequent iterate and gradient displacement pair. Third, a vector $\tau\in\mathbb{R}^{m}$ was generated randomly using randn in order to define $s_{0}=[s_{1}\ \cdots\ s_{m}]\tau$. The corresponding gradient displacement $y_{0}$ was determined by stepping backward from the fixed point along $s_{0}$. In this manner, we obtained $\{(s_{k},y_{k})\}_{k=0}^{m}$ from a mock optimization procedure from some initial point in such a way that $s_{0}$ was guaranteed to lie the span of the subsequent iterate displacements.

For each dataset starting from $W_{0}=I$, we computed the BFGS inverse Hessian approximation from $\{(s_{k},y_{k})\}_{k=0}^{m}$ and constructed the inverse Hessian approximation corresponding to the set of pairs when Agg-BFGS was used to aggregate the information into $\{(s_{k},\tilde{y}_{k})\}_{k=1}^{m}$. Figure 2 shows box plots for relative errors between each pair of inverse Hessian approximations, where error is measured in terms of the maximum component-wise absolute difference between matrix entries divided by the largest absolute value of an element of the BFGS inverse Hessian approximation. The results show that while the errors are larger as $n$ increases and as $m$ is closer to $n$, they remain accurate relative to machine precision for all $(n,m)$.

As a second experiment, for $(n,m)\in\{(8,8),(32,32),(128,128)\}$, we aimed to investigate how errors might accumulate as Agg-BFGS is applied over a sequence of iterations. For these experiments, we generated data using a similar mock optimization procedure as in the previous experiment. Specifically, from some randomly generated starting point, we computed a random step as in the aforementioned procedure, continuing until $n+8$ iterations were performed. Figure 3 shows the errors for each iteration beyond $k=n$. As in Figure 1(a), these results show that the errors do not accumulate too poorly as $k$ increases. We conjecture that part of the reason for this is that errors that result from each application of Agg-BFGS can be over-written eventually, at least to some extent, similar to the manner in which curvature information is ultimately over-written in full-memory BFGS.

4.2 A practical adaptive L-BFGS method using Agg-BFGS

To demonstrate the use of Agg-BFGS in a minimization algorithm, we compare the results of three algorithms using quasi-Newton inverse Hessian approximations for computing the search directions, where for each algorithm the same weak Wolfe line search was used for computing stepsizes. The first algorithm employed L-BFGS approximations. The second algorithm also employed L-BFGS approximations, but with Agg-BFGS used to aggregate information when deemed appropriate in the manner described below. Since this second algorithm does not always aggregate the oldest pair, but rather uses a particular scheme for choosing which historical pair to use in the aggregation approach (as explained later on), we also considered a third algorithm. This algorithm uses the same scheme as the second to determine which pair to consider, but simply removes this pair rather than perform aggregation with it. Consequently, the comparison between the second and third algorithms shows the effect of aggregation itself, rather than only the difference between removing the oldest versus some other historical pair in an L-BFGS scheme.

We performed experiments with problems from the CUTEst collection [21]. We chose all problems from the collection for which the number of variables could be chosen in the interval $[10,3000]$. For most problems, $n\gg m$, meaning that the costs of performing aggregation was negligible compared to the costs of computing search directions, which were the same for all algorithms using the standard two-loop recursion for L-BFGS. (See Tables 4.2 and 4.2 for $n$ for all test problems used.) All algorithms were run until an iteration, call it $k\in\mathbb{N}$, was reached with $\|g_{k}\|_{\infty}\leq 10^{-6}\max\{1,\|g_{0}\|_{\infty}\}$, or until an iteration limit of $10^{5}$ was surpassed.

For the second algorithm, Agg-BFGS was employed in the following manner. Suppose in iteration $k\in\mathbb{N}$ that the algorithm had in storage the displacement pairs $\{(s_{k_{j}},\tilde{y}_{k_{j}})\}_{j=1}^{\mkern 1.5mu\overline{\mkern-1.5mum}}$ for some indices $\{k_{1},\dots,k_{\mkern 1.5mu\overline{\mkern-1.5mum}}\}\subseteq\{1,\dots,k-1\}$ with $\mkern 1.5mu\overline{\mkern-1.5mum}\leq m$ and the new pair $(s_{k},y_{k})$ has been computed. Using techniques described in §3.4, starting with $j=\mkern 1.5mu\overline{\mkern-1.5mum}$ and iterating down to $j=1$, we determined if $s_{k_{j}}$ lay approximately in $\operatorname{span}\{s_{k_{j+1}},\dots,s_{k_{\mkern 1.5mu\overline{\mkern-1.5mum}}},s_{k}\}$ (see next paragraph). If so, then we applied Agg-BFGS to aggregate the information in the pair with index $k_{j}$; otherwise, if no such iterate displacement was determined, then either the pair $(s_{k},y_{k})$ was simply added to the set of pairs in storage (if $\mkern 1.5mu\overline{\mkern-1.5mum}<m$) or the pair with index $k_{1}$ was dropped (if $\mkern 1.5mu\overline{\mkern-1.5mum}=m$, as in standard L-BFGS). In this manner, the number of pairs stored and employed remained less than or equal to $m$, as usual for L-BFGS.

To determine if $s_{k_{j}}$ lay approximately in $\operatorname{span}\{s_{k_{j+1}},\dots,s_{k_{\mkern 1.5mu\overline{\mkern-1.5mum}}},s_{k}\}$, we first computed the orthogonal projection of $s_{k_{j}}$ onto this subspace, call it $\hat{s}_{k_{j}}$. Computing this projection is inexpensive given a Cholesky factorization of the inner product matrix $[s_{k_{j+1}}\ \cdots\ s_{k_{\mkern 1.5mu\overline{\mkern-1.5mum}}}\ s_{k}]^{T}[s_{k_{j+1}}\ \cdots\ s_{k_{\mkern 1.5mu\overline{\mkern-1.5mum}}}\ s_{k}]$, which can be updated with only ${\cal O}(mn)$ cost in each outer iteration as described in §3.4. If the condition

[TABLE]

held for some $\texttt{tol}\in(0,1)$, then it was determined that $s_{k_{j}}$ lay approximately in the subspace, since in this case the norm of the projection $\hat{s}_{k_{j}}$ was sufficiently large compared to the orthogonal component $s_{k_{j}}-\hat{s}_{k_{j}}$. For $j=\{\mkern 1.5mu\overline{\mkern-1.5mum},\dots,2\}$, we used $\texttt{tol}=10^{-8}$, while for $j=1$ we loosened the tolerance to $\texttt{tol}=10^{-4}$ to promote aggregation of the oldest pair. If (46) held for some index $k_{j}$, then aggregation was performed with the pair $(\hat{s}_{k_{j}},\tilde{y}_{k_{j}})$. Notice that this ensures that the aggregation is performed with an iterate displacement vector (namely, the projected vector $\hat{s}_{k_{j}}$) that lies in the subspace defined by subsequent displacements, as required in our algorithms and theoretical results in §3. (We remark in passing that, using this scheme, it is possible that aggregation could be triggered multiple times in a single iteration. Our software allows for this, but it did not occur in our experiments.)

One might expect that, under these conditions, Agg-BFGS might not perform aggregation often. However, on the contrary, our results show that aggregation was performed quite often.333This provides evidence for the belief, held by some optimization researchers, that when solving certain large-scale problems one often observes that consecutive steps lie approximately in low-dimensional subspaces. Tables 4.2 and 4.2 show detailed results for the first two algorithms in our experiments for $m=5$. (We do not include the third algorithm since, as shown below, it was inferior to the first and second algorithms.) For conciseness, we refer to the algorithm that employed a standard LBFGS(5) strategy as LBFGS(5), and we refer to the strategy that employed Agg-BFGS as AggBFGS(5). For AggBFGS(5), the table reports the number of aggregations performed, which in many cases was significant compared to the number of iterations.