BFGS with Update Skipping and Varying Memory.
We give conditions under which limited-memory quasi-Newton methods with exact line searches will terminate in n steps when minimizing n-dimensional quadratic functions. We show that although all Broyden family methods terminate in n steps in their full-memory versions, only BFGS does so with limited-memory.  Additionally, we show that full-memory Broyden family methods with exact line searches terminate in at most n steps when p matrix updates are skipped. We introduce new limited-memory BFGS variants and test them on nonquadratic minimization problems.
Introduction
. The quasi-Newton family of algorithms remains a standard
workhorse for minimization. Many of these methods share the properties of finite
termination on strictly convex quadratic functions, a linear or superlinear rate of
convergence on general convex functions, and no need to store or evaluate the second
derivative matrix. In general, an approximation to the second derivative matrix is
built by accumulating the results of earlier steps. Descriptions of many quasi-Newton
algorithms can be found in books by Luenberger [16], Dennis and Schnabel [7], and
Golub and Van Loan [11].
Although there are an infinite number of quasi-Newton methods, one method surpasses
the others in popularity: the BFGS algorithm of Broyden, Fletcher, Goldfarb,
and Shanno; see, e.g., Dennis and Schnabel [7]. This method exhibits more robust behavior
than its relatives. Many attempts have been made to explain this robustness,
but a complete understanding is yet to be obtained [23]. One result of the work in this
paper is a small step toward this understanding, since we investigate the question of
how much and which information can be dropped in BFGS and other quasi-Newton
methods without destroying the property of quadratic termination.
We answer this question in the context of exact line search methods, those that
find a minimizer on a one-dimensional subspace at every iteration. (In practice,
inexact line searches that satisfy side conditions such as those proposed by Wolfe, see
x4.3, are substituted for exact line searches.) We focus on modifications of well-known
quasi-Newton algorithms resulting from limiting the memory, either by discarding
the results of early steps (x2) or by skipping some updates to the second derivative
approximation (x3). We give conditions under which quasi-Newton methods will
terminate in n steps when minimizing quadratic functions of n variables. Although
all Broyden family methods (see x2) terminate in n steps in their full-memory versions,
we show that only BFGS has n-step termination under limited-memory. We also show
that the methods from the Broyden family terminate in n steps even if p updates
are skipped, but termination is lost if we both skip updates and limit the memory.
y Applied Mathematics Program, University of Maryland, College Park, MD 20742.
gibson@math.umd.edu. This work was supported in part by the National Physical Science Con-
sortium, the National Security Agency, and the University of Maryland.
z Department of Computer Science and Institute for Advanced Computer Studies, University of
Maryland, College Park, MD 20742. oleary@cs.umd.edu. This work was supported by the National
Science Foundation under grant NSF 95-03126.
x Department of Pure and Applied Mathematics, Washington State University, Pullman, WA
99164. nazareth@amath.washington.edu.
T. Gibson, D. P. O'Leary, L. Nazareth
In x4, we report the results of experiments with new limited-memory BFGS variants
on problems taken from the CUTE [3] test set, showing that some savings in
time can be achieved.
Notation. Matrices and vectors are denoted by boldface upper-case and lower-case
letters respectively. Scalars are denoted by Greek or Roman letters. The superscript
"T" denotes transpose. Subscripts denote iteration number. Products are always
taken from left to right. The notation spanfx 1 denotes the subspace
spanned by the vectors x Whenever we refer to an n-dimensional strictly
convex quadratic function, we assume it is of the form
where A is a positive definite n \Theta n matrix and b is an n-vector.
2. Limited-Memory Variations of Quasi-Newton Algorithms. In this
section we characterize full-memory and limited-memory methods that terminate in
n iterations on n-dimensional strictly convex quadratic minimization problems using
exact line searches. Most full-memory versions of the methods we will discuss are
known to terminate in n iterations. Limited-memory BFGS (L-BFGS) was shown by
Nocedal [22] to terminate in n steps. The preconditioned conjugate gradient method,
which can be cast as a limited-memory quasi-Newton method, is also known to terminate
in n iterations; see, e.g., Luenberger [16]. Little else is known about termination
of limited-memory methods.
Let f(x) denote the strictly convex quadratic function to be minimized, and let
g(x) denote the gradient of f . We define is the kth iterate. Let
denote the change in the current iterate and
denote the change in gradient.
Let x 0 be the starting point, and let H 0 be the initial inverse Hessian approximation.
For
1. Compute
2. Choose ff k ? 0 such that f(x k
3. Set s
4. Set x
5. Compute
7. Choose H k+1 .
Fig. 2.1. General Quasi-Newton Method
We present a general result that characterizes quasi-Newton methods, see Figure 2.1,
that terminate in n iterations. We restrict ourselves to methods with an update of
the form
mk
Here,
L-BFGS Variations 3
1. H 0 is an n \Theta n symmetric positive definite matrix that remains constant for
all k, and fl k is a nonzero scalar that can be thought of as an iterative rescaling of
2. P k is an n \Theta n matrix that is the product of projection matrices of the form
where is an n \Theta n matrix
that is the product of projection matrices of the same form where u is any n-vector
3. m k is a nonnegative integer, w ik n-vector, and z ik
vector in spanfs g.
We refer to this form as the general form. The general form fits many known
quasi-Newton methods, including the Broyden family and the limited-memory BFGS
method. We do not assume that these quasi-Newton methods satisfy the secant
condition,
nor that H k+1 is positive definite and symmetric. Symmetric positive definite updates
are desirable since this guarantees that the quasi-Newton method produces descent
directions. Note that if the update is not positive definite, we may produce a d k such
that d T
which case we choose ff k over all negative ff rather than all positive
ff.
Example. The method of steepest descent [16] fits the general form (2.1). For
each k we define
Note that neither w nor z vectors are specified since m
Example. The 1)st update for the conjugate gradient method with preconditioner
fits the general form (2.1) with
Example. The L-BFGS update, see Nocedal [22], with limited-memory constant
m can be written as
k\Gammam k+1;k
Y
L-BFGS fits the general form (2.1) if at iteration k we choose
4 T. Gibson, D. P. O'Leary, L. Nazareth
Observe that P k ; Q k and z ik all obey the constraints imposed on their construction.
BFGS is related to L-BFGS is the following way: if we were to use every (s; y)
pair in the formation of each update (i.e. we have unlimited memory), we would be
creating the same updates as BFGS. In practice, however, one would never do that
because it would take more memory than storing the BFGS matrix.
Example. We will define limited-memory DFP (L-DFP). Our definition is consistent
with the definition of limited-memory BFGS given in Nocedal [22]. Let m - 1
mg. In order to define the L-DFP update we need to create
a sequence of auxiliary matrices for
where
U DFP (H; s;
ss T
The matrix -
k+1 is the result of applying the DFP update m k times to the matrix
H 0 with the m k most recent (s; y) pairs. Thus, the 1)st L-DFP matrix is given
by
To simplify our description, note that -
k+1 can be rewritten as
k\Gammam k+i
k\Gammam k+i
k\Gammam k+i
k\Gammam k+i
y k\Gammam k+i
k\Gammam k+j
k\Gammam k+j
y k\Gammam k+j
Y
y k\Gammam k+l
k\Gammam k+l
Thus H k+1 can be written as
mk
ik
k\Gammam k+i
k\Gammam k+i
y k\Gammam k+i
where
mk
Y
y k\Gammam k+j
H (j \Gamma1)
k\Gammam k+j
H (j \Gamma1)
L-BFGS Variations 5
Equation (2.7) looks very much like the general form given in (2.1). L-DFP fits the
general form with the following choices:
k\Gammam k+i
y k\Gammam k+i ); and z
Except for the choice of P k , it is trivial to verify that the choices satisfy the general
form. To prove that P k satisfies the requirements, we need to show
Proposition 2.1. For limited-memory DFP, the following two conditions hold
for each value of k:
Proof. We will prove this via induction. Suppose We have
(Recall that spanfs 0 g is trivially equal to spanfH 0 g 0 g.) Furthermore,
So we can conclude,
and so the base case holds.
Assume that
We will use induction on i to show (2.10) for the 1)st case. For
Using the induction assumptions from the induction on k, we get that
6 T. Gibson, D. P. O'Leary, L. Nazareth
Assume that -
(induction assumption for i). Next,
For values of
maps any vector v into
and so -
is in
Using the induction assumptions for both i and k, we get
and we can continue the induction on i. If
so
Hence the induction on i is complete and this proves (2.10) in the (k 1)st case.
consider
mk
k\Gammam k+i
k\Gammam k+i
Using the structure of V jk and (2.10) we see that
Hence, (2.11) also holds in the (k 1)st case.
Example. The Broyden Class or Broyden Family is the class of quasi-Newton
methods whose matrices are linear combinations of the DFP and BFGS matrices:
see, e.g., Luenberger [16, Chap. 9]. The parameter OE is usually restricted to values
that are guaranteed to produce a positive definite update, although recent work with
SR1, a Broyden Class method, by Khalfan, Byrd and Schnabel [14] may change this
practice. No restriction on OE is necessary for the development of our theory. The
Broyden class update can be expressed as
Variations 7
We sketch the explanation of how the full-memory version fits the general form
given in (2.1). The limited-memory case is similar. We can rewrite the Broyden Class
update as follows:
Hence,
where
hi
y
s
y
It is left to the reader to show that H k y k is in spanfs thus the
Broyden Class updates fit the form in (2.1).
2.1. Termination of Limited-Memory Methods. In this section we show
that methods fitting the general form (2.1) produce conjugate search directions (The-
orem 2.2) and terminate in n iterations (Corollary 2.3) if and only if P k maps
spanfy into spanfy for each n. Furthermore, this
condition on P k is satisfied only if y k is used in its formation (Corollary 2.4).
Theorem 2.2. Suppose that we apply a quasi-Newton method (Figure 2.1) with
an update of the form (2.1) to minimize an n-dimensional strictly convex quadratic
function. Then for each k before termination (i.e. g k+1 6= 0),
As
if and only if
Proof. (() Assume that (2.15) holds. We will prove (2.12)-(2.14) by induction.
Since the line searches are exact, g 1 is orthogonal to s 0 . Using the fact that P 0 y
8 T. Gibson, D. P. O'Leary, L. Nazareth
from (2.15), and the fact that z i0 2 spanfs 0 g implies g T
see that s 1 is conjugate to s 0 since
z i0 w T
0:
Lastly, spanfs 0 g, and so the base case is established.
We will assume that claims (2.12)-(2.14) hold for
that they also hold for
The vector g - k+1 is orthogonal to s - k since the line search is exact. Using the
induction hypotheses that g -
k is orthogonal to fs
is conjugate to
g, we see that for
Hence, (2.12) holds for
To prove (2.13), we note that
As
so it is sufficient to prove that g T
We will use the
following facts:
- k+1 since the v in each of the projections used to form Q - k is in
k+1 is orthogonal to that span.
since each z i - k is in spanfs
is orthogonal to that span.
(iii) Since we are assuming that (2.15) holds true, for each
there
exists can be expressed as
P-
(iv) For is orthogonal to H 0 y i because g -
k+1 is orthogonal
to spanfs
from (2.14).
Thus,
Variations 9
0:
Thus, (2.13) holds for
k.
Lastly, using (i) and (ii) from above,
maps any vector v into spanfv; s by construction, there exist
Hence,
so
To show equality of the sets, we will show that H 0 g - k+1 is linearly independent of
g. (We already know that the basis fH 0 is linearly
independent since it spans the same space as the linearly independent set fs
and has the same number of elements.) Suppose that H 0 g - k+1 is not linearly indepen-
dent. Then there exist OE
k , not all zero, such that
Recall that g - k+1 is orthogonal to fs
g. By our induction assumption, this
implies that g -
k+1 is also orthogonal to fH 0 g. Thus for any j between 0
and - k,
positive definite and g j is nonzero, we conclude that OE j must be zero.
Since this is true for every j between zero and k, we have a contradiction. Thus, the
set fH 0 is linearly independent. Hence, (2.14) holds for
k.
Assume that (2.12)-(2.14) hold for all k such that g k+1 6= 0 but that (2.15)
does not hold; i.e., there exist j and k such that g k+1 6= 0, j is between 0 and k, and
(2.
T. Gibson, D. P. O'Leary, L. Nazareth
This will lead to a contradiction. By construction of P k , there exist -
that
By assumption (2.16), - k must be nonzero. From (2.13), it follows that g T
Using facts (i), (ii), and (iv) from before, (2.14) and (2.17), we get
mk
z ik w T
ik
mk
Thus since neither fl k nor - k is zero, we must have
but this is a contradiction since H 0 is positive definite and g k+1 was assumed to be
nonzero.
When a method produces conjugate search directions, we can say something about
termination.
Corollary 2.3. Suppose we have a method of the type described in Theorem 2.2
satisfying (2.15). Suppose further that H j Then the scheme
reproduces the iterates from the conjugate gradient method with preconditioner H 0 and
terminates in no more than n iterations.
Proof. Let k be such that are all nonzero and such that H i
we have a method of the type described in Theorem 2.2 satisfying
(2.15), conditions (2.12) - (2.14) hold. We claim that the (k 1)st subspace of
search directions, spanfs is equivalent to the 1)st Krylov subspace,
g.
From (2.14), we know that spanfs g. We will
show via induction that spanfH 0 g.
This base case is trivial, so assume that
for some
L-BFGS Variations 11
and from (2.14) and the induction hypothesis,
which implies that
Hence, the search directions span the Krylov subspace. Since the search directions
are conjugate (2.13) and span the Krylov subspace, the iterates are the same as those
produced by conjugate gradients with preconditioner H 0 .
Since we produce the same iterates as the conjugate gradient method and the
conjugate gradient method is well-known to terminate within n iterations, we can
conclude that this scheme terminates in at most n iterations.
Note that we require that H j g j be nonzero whenever g j is nonzero; this requirement
is necessary since not all the methods produce positive definite updates and it
is possible to construct an update that maps g j to zero. If this were to happen, we
would have a breakdown in the method.
The next corollary defines the role that the latest information (s k and y k ) plays
in the formation of the kth H-update.
Corollary 2.4. Suppose we have a method of the type described in Theorem
2.2 satisfying (2.15). Suppose further that at the kth iteration P k is composed
of p projections of the form in (2.2). Then at least one of the projections must have
is a single projection (p = 1), then v
must be of the form
Proof. Consider the case of p = 1. We have
where k+1g. We will assume that
for some scalars oe i and ae i . By (2.15), there exist -
Then
and so
(2.
12 T. Gibson, D. P. O'Leary, L. Nazareth
From (2.13), the set fs is conjugate and thus linearly independent. Since we
are working with a quadratic, y As i for all i; and since A is symmetric positive
definite, the set fy is also linearly independent. So the coefficient of the y k
on the left-hand side of (2.18) must match that on the right-hand side, thus
Hence,
and y k must make a nontrivial contribution to P k .
Next we will show that ae Assume that j is between 0
As j
As j
Now s j As j is nonzero because A is positive definite. If ae j is nonzero then the coefficient
of u is nonzero and so y k must make a nontrivial contribution to P k y j , implying
that g. This is a contradiction. Hence, ae
To show that ae k 6= 0, consider P k y k . Suppose that ae
As k+1
and so
This contradicts P k y k 2 spanfy must be nonzero.
Now we will discuss that p ? 1 case. Label the u-components of the p projections
as
for some scalars fl 1 through fl p . We know that
L-BFGS Variations 13
and that
Thus
and since u we can conclude that at least one u i
must have a nontrivial contribution from y k .
2.2. Examples of Methods that Reproduce the CG Iterates. Here are
some specific examples of methods that fit the general form, satisfy condition (2.15)
of Theorem 2.2, and thus terminate in at most n iterations.
Example. The conjugate gradient method with preconditioner H 0 , see (2.4),
satisfies condition (2.15) of Theorem 2.2 since
Example. Limited-memory BFGS, see (2.6), satisfies condition (2.15) of Theorem
2.2 since
ae 0 for
Example. DFP (with full memory), see (2.8), satisfies condition (2.15) of Theorem
2.2. Consider P k in the full memory case. We have
Y
For full-memory DFP, H i y 1. Using this fact, one can easily
verify that P k y Therefore, full-memory DFP satisfies condition
(2.15) of Theorem 2.2. The same reasoning does not apply to the limited-memory
case as we shall show in x2.3.
The next corollary gives some ideas for other methods that are related to L-BFGS
and terminate in at most n iterations on strictly convex quadratics.
Corollary 2.5. The L-BFGS (2.5) method will terminate in n iterations on
an n-dimensional strictly convex quadratic function even if any combination of the
following modifications is made to the update:
1. Vary the limited-memory constant, keeping m k - 1.
2. Form the projections used in V k from the most recent along with
any set of other pairs from f(s
3. Form the projections used in V k from the most recent along with
linear combinations of pairs from f(s 0 ; y
4. Iteratively rescale H 0 .
Proof. For each variant, we show that the method fits the general form in (2.1),
satisfies condition (2.15) of Theorem 2.2 and hence terminates by Corollary 2.3.
1. Let m ? 0 be any value which may change from iteration to iteration, and
define
Y
14 T. Gibson, D. P. O'Leary, L. Nazareth
Choose
These choices fit the general form. Furthermore,
so this variation satisfies condition (2.15) of Theorem 2.2.
2. This is a special case of the next variant.
3. At iteration k, let (- s (i)
y (i)
k ) denote the ith choice of any
linear combination from the span of the set
and let (- s (m)
Y
(- y (i)
Choose
These choices satisfy the general form (2.1). Furthermore,
ae
k for some i; and
Hence, this variation satisfies condition (2.15) of Theorem 2.2.
4. Let fl k in (2.1) be the scaling constant, and choose the other vectors and
matrices as in L-BFGS (2.6).
Combinations of variants are left to the reader.
Remark. Part 3 of the previous corollary shows that the "accumulated step"
method of Gill and Murray [10] terminates on quadratics.
Remark. Part 4 of the previous corollary shows that scaling does not affect
termination in L-BFGS. In fact, for any method that fits the general form, it is easy
to see that scaling will not affect termination on quadratics.
2.3. Examples of Methods that Do Not Reproduce the CG Iterates.
We will discuss several methods that fit the general form given in (2.1) but do not
satisfy the conditions of Theorem 2.2.
L-BFGS Variations 15
Example. Steepest descent, see (2.3), does not satisfy condition (2.15) of Theorem
2.2 and thus does not produce conjugate search directions. This fact is well-
known; see, e.g., Luenberger [16].
Example. Limited-memory DFP, see (2.8), with does not satisfy the condition
on P k (2.15) for all k, and so the method will not produce conjugate directions.
For example, suppose that we have a convex quadratic with
Using a limited-memory constant of exact arithmetic, it can be seen that
the iteration does not terminate within the first 20 iterations of limited-memory DFP
with I. The MAPLE notebook file used to compute this example is available
on the World Wide Web [9].
Remark. Using the above example, we can easily see that no limited-memory
Broyden class method except limited-memory BFGS terminates within the first n
iterations.
3. Update-Skipping Variations for Broyden Class Quasi-Newton Algo-
rithms. The previous section discussed limited-memory methods that behave like
conjugate gradients on n-dimensional strictly convex quadratic functions. In this sec-
tion, we are concerned with methods that skip some updates in order to reduce the
memory demands. We establish conditions under which finite termination is preserved
but delayed for the Broyden Class.
3.1. Termination when Updates are Skipped. It was shown by Powell [26]
that if we skip every other update and take direct prediction steps (i.e. steps of length
one) in a Broyden class method, then the procedure will terminate in no more than
2n+1 iterations on an n-dimensional strictly convex quadratic function. An alternate
proof of this result is given by Nazareth [21].
We will prove a related result. Suppose that we are doing exact line searches using
a Broyden Class quasi-Newton method on a strictly convex quadratic function and
decide to "skip" p updates to H (i.e. choose H occasions). Then, the
algorithm terminates in no more than n iterations. In contrast to Powell's result,
it does not matter which updates are skipped or if multiple updates are skipped in a
row.
Theorem 3.1. Suppose that a Broyden Class method using exact line searches
is applied to an n-dimensional strictly convex quadratic function and p updates are
skipped. Let
the update at iteration j is not skippedg:
Then for all
As
Furthermore, the method terminates in at most n iterations at the exact minimizer

Proof. We will use induction on k to show (3.1) and
T. Gibson, D. P. O'Leary, L. Nazareth
Then (3.2) follows easily since for all j 2 J(k),
As
0:
be the least value of k such that J(k) is nonempty; i.e., J(k 0 g.
Then g k0+1 is orthogonal to s k0 since line searches are exact, and H k0+1 y
since all members of the Broyden Family satisfy the secant condition. Hence, the base
case is true. Now assume that (3.1) and (3.3) hold for all values of
We will show that they also hold for
Case I. Suppose that - k 62 J( - k). Then H -
k and J( -
any j 2 J( - k),
As j
and
Case II. Suppose that - k 2 J( - k). Then H - k+1 satisfies the secant condition and
kg. Now g - k+1 is orthogonal to s k since the line searches are exact,
and it is orthogonal to the older s j by the argument in (3.4). The secant condition
guarantees that H - k+1 y
k , and for
we have
!/
As j
!/
As j
In either case, the induction result follows.
Suppose that we skip p updates. Then the set J(n cardinality n.
Without loss of generality, assume that the set fs i g i2J(n\Gamma1+p) has no zero elements.
From (3.2), the vectors are linearly independent. By (3.1),
and so gn+p must be zero. This implies that xn+p is the exact minimizer of f .
L-BFGS Variations 17
3.2. Loss of Termination for Update Skipping with Limited-Memory.
Unfortunately, updates that use both limited-memory and repeated update-skipping
do not produce n conjugate search directions for n-dimensional strictly convex qua-
dratics, and the termination property is lost. We will show a simple example, limited-memory
skipping every other update. Note that according to
Corollary 2.4, we would still be guaranteed termination if we used the most recent
information in each update.
Example. Suppose that we have a convex quadratic with
We apply limited-memory BFGS with limited-memory constant
and skip every-other update to H. Using exact arithmetic in MAPLE, we observe
that the process does not terminate even after 100 iterations [9].
4. Experimental Results. The results of x2 and x3 lead to a number of ideas
for new methods for unconstrained optimization. In this section, we motivate, de-
velop, and test these ideas. We describe the collection of test problems in x4.2. The
test environment is described in x4.3. Section 4.4.1 outlines the implementation of the
L-BFGS method (our base for all comparisons) and xx4.4.2-4.4.7 describe the varia-
tions. Pseudo-code for L-BFGS and its variations is given in Appendix B. Complete
numerical results, many graphs of the numerical results, and the original FORTRAN
code are available [9].
4.1. Motivation. So far we have only given results for convex quadratic func-
tions. While termination on quadratics is beautiful in theory, it does not necessarily
yield insight into how these methods will do in practice.
We will not present any new results relating to convergence of these algorithms
on general functions; however, many of these can be shown to converge using the
convergence analysis presented in x7 of [15]. In [15], Liu and Nocedal show that
a limited-memory BFGS method implemented with a line search that satisfies the
strong Wolfe conditions (see x4.3 for a definition) is R-linearly convergent on a convex
function that satisfies a few modest conditions.
4.2. Test Problems. For our test problems, we used the Constrained and Unconstrained
Testing Environment (CUTE) by Bongartz, Conn, Gould and Toint. The
package is documented in [3] and can be obtained via the world wide web [2] or via ftp
[1]. The package contains a large collection of test problems as well as the interfaces
necessary for using the problems. The test problems are stored as "SIF" files. We
chose a collection of 22 unconstrained problems. The problems ranged in size from
to 10,000 variables, but each took L-BFGS with limited-memory constant
at least 60 iterations to solve. Table 4.1 enumerates the problems, giving the SIF file
name, the dimension (n), and a description for each problem. The CUTE package
also provides a starting point
4.3. Test Environment. We used FORTRAN77 code on an SGI Indigo 2 to
run the algorithms, with FORTRAN BLAS routines from NETLIB. We used the
compiler's default optimization level.

Figure

2.1 outlines the general quasi-Newton implementation that we followed.
For the line search, we use the routines cvsrch and cstep written by Jorge J. Mor'e
T. Gibson, D. P. O'Leary, L. Nazareth
No. SIF Name n Description & Reference
Extended Rosenbrock function (nonseparable
version) [30, Problem 10].
problem [17, Problem 20].
3 TOINTGOR 50 Toint's operations research problem [29].
4 TOINTPSP 50 Toint's PSP operations research problem [29].
5 CHNROSNB 50 Chained Rosenbrock function [29].
6 ERRINROS 50 Nonlinear problem similar to CHNROSNB [28].
7 FLETCHBV 100 Fletcher's boundary value problem [8, Problem
1].
8 FLETCHCR 100 Fletcher's chained Rosenbrock function [8, Problem
2].
9 PENALTY2 100 Second penalty problem [17, Problem 24].
Problem 5].
11 BDQRTIC 1000 Quartic with a banded Hessian with band-
diagonal variant of the Broyden tridiagonal
system with a band away from diagonal [29].
First penalty problem [17, Problem 23].
14 POWER 1000 Power problem by Oren [25].
MSQRTALS 1024 The dense matrix square root problem by Nocedal
and Liu (case 0) seen as a nonlinear equation
problem [4, Problem 204].
MSQRTBLS 1025 The dense matrix square root problem by Nocedal
and Liu (case 1) seen as a nonlinear equation
problem [4, Problem 201].
17 CRAGGLVY 5000 Extended Cragg & Levy problem [30, Problem
test problem [5, Problem
57].
19 POWELLSG 10000 Extended Powell singular function [17, Problem
13].
Another function with nontrivial groups and repetitious
elements [12].
tridiagonal matrix square root
problem [4, Problem 151].
22 TRIDIA 10000 Shanno's TRIDIA quadratic tridiagonal problem
[30, Problem 8].

Table
Test problem collection. Each problems was chosen from the CUTE package.
and David Thuente from a 1983 version of MINPACK. This line search routine finds an
ff that meets the strong Wolfe conditions,
see, e.g., Nocedal [23]. We used 0:9. Except for the first
iteration, we always attempt a step length of 1.0 first and only use an alternate value
if 1.0 does not satisfy the Wolfe conditions. In the first iteration, we initially try a
step length equal to kg . The remaining line search parameters are detailed in


Appendix

A.
We generate the matrix H k by either the limited-memory update or one of the
variations described in x4.4, storing the matrix implicitly in order to save both memory
and computation time.
We terminate the iterations if any of the following conditions are met at iteration
L-BFGS Variations 19
k:
1. The inequality
is satisfied,
2. the line search fails, or
3. the number of iterations exceeds 3000.
We say that the iterates have converged if the first condition is satisfied. Otherwise,
the method has failed.
4.4. L-BFGS and its variations. We tried a number of variations to the standard
L-BFGS algorithm. L-BFGS and these variations are described in this subsection
and summarized in Table 4.2.
4.4.1. L-BFGS: Algorithm 0. The limited-memory BFGS update is given in
(2.5) and described fully by Byrd, Nocedal and Schnabel [22]. Our implementation
and the following description come essentially from [22].
Let H 0 be symmetric and positive definite and assume that the m k pairs
each satisfy s T
We will let
and m is some positive integer. We will assume that
I and that H 0 is iteratively rescaled by a constant fl k as is commonly done
in practice. Then, the matrix H k obtained by k applications of the limited-memory
BFGS update can be expressed as
\GammaU
where U k and D k are the m k \Theta m k matrices given by
ae
and
We will describe how to compute d k in the case that k ? 0. Let x k
be the current iterate. Let m Given s , the matrices
and the vectors S T
1. Update the n \Theta m k\Gamma1 matrices S k\Gamma1 and Y k\Gamma1 to get the n \Theta m k matrices
using s k\Gamma1 and y
2. Compute the m k -vectors S T
3. Compute the m k -vectors S T
by using the fact that
We already know components of S k g k\Gamma1 from S k\Gamma1 g k\Gamma1 , and likewise for
. We need only compute s T
and do the subtractions.
20 T. Gibson, D. P. O'Leary, L. Nazareth
No. Reference Brief Description
x4.4.1 L-BFGS with no options.
Allow m to vary iteratively basing the choice of m of kgk
and not allowing m to decrease.
Allow m to vary iteratively basing the choice of m of kgk
and allowing m to decrease.
Allow m to vary iteratively basing the choice of m of kg=xk
and not allowing m to decrease.
Allow m to vary iteratively basing the choice of m of kg=xk
and allowing m to decrease.
5 x4.4.3 Dispose of old information if the step length is greater than
one.
6 x4.4.4, Variation 1 Back-up if the current iteration is odd.
7 x4.4.4, Variation 2 Back-up if the current iteration is even.
8 x4.4.4, Variation 3 Back-up if a step length of 1.0 was used in the last iteration.
9 x4.4.4, Variation 4 Back-up if kg k k ? kg
Back-up if a step length of 1.0 was used in the last iteration
and we did not back-up on the last iteration.
and we did not back-up on the
last iteration.
neither of the two vectors to be merged is itself
the result of a merge and the 2nd and 3rd most recent steps
taken were of length 1.0.
13 x4.4.5, Variation 2 Merge if we did not do a merge the last iteration and there
are at least two old s vectors to merge.
14 x4.4.6, Variation 1 Skip update on odd iterations.
update on even iterations.
Alg. 5 & Alg. 8 Dispose of old information and back-up on the next iteration
if the step length is greater than one.
Alg. 13 & Alg. 1 Merge if we did not do a merge the last iteration and there
are at least two old s vectors to merge, and allow m to vary
iteratively basing the choice of m of kgk and not allowing
m to decrease.
19 Alg. 13 & Alg. 3 Merge if we did not do a merge the last iteration and there
are at least two old s vectors to merge, and allow m to
vary iteratively basing the choice of m of kg=xk and not
allowing m to decrease.
Alg. 13 & Alg. 2 Merge if we did not do a merge the last iteration and there
are at least two old s vectors to merge, and allow m to vary
iteratively basing the choice of m of kgk and allowing m to
decrease.
Alg. 13 & Alg. 2 Merge if we did not do a merge the last iteration and there
are at least two old s vectors to merge, and allow m to vary
iteratively basing the choice of m of kg=xk and allowing m
to decrease.

Table
Description of Numerical Algorithms
4. Compute
k . Rather than recomputing U
k , we update the matrix from
the previous iteration by deleting the leftmost column and topmost row if m
and appending a new column on the right and a new row on the bottom. Let ae
1=s T
be the (m lower right submatrix of U
and let (S T
be the upper
L-BFGS Variations 21
Note that s T
so is already computed.
5. Assemble Y T
We have already computed all the components.
6. Update D k using D k\Gamma1 and s T
7. Compute
Note that both y T
8. Compute two intermediate values
9. Compute
The storage costs for this are very low. In order to reconstruct H k , we need to
store
diagonal matrix) and a few m-vectors. This requires
only 2mn Assuming m !! n, this is much less storage than
the n 2 storage required for typical implementation of BFGS.
Step Operation Count
9

Table
Operations Count for Computation of H k g k . Steps with no operations are not shown.
The computation of Hg takes at most O(mn) operations assuming n ?? m. (See

Table

4.3.) This is much less than the O(n 2 normally needed to compute Hg
when the whole matrix H is stored.
We are using L-BFGS as our basis for comparison. For information on the performance
of L-BFGS see Liu and Nocedal [15] and Nash and Nocedal [19].
4.4.2. Varying m iteratively: Algorithms 1-4. In typical implementations
of L-BFGS, m is fixed throughout the iterations: once m updates have accumulated,
m updates are always used. We considered the possibility of varying m iteratively,
preserving finite termination on convex quadratics. Using an argument similar to that
presented in [15], we can also prove that this algorithm has a linear rate of convergence
on a convex function that satisfies a few modest conditions.
We tried four different variations on this theme. All were based on the following
linear formula that scales m in relation to the size of kgk. Let m k be the number of
iterates saved at the kth iteration, with Here, think of m as the maximum
allowable value of m k . Let the convergence test be given by kg k k=kx k k ! ffl. Then
the formula for m k at iteration k is
ae
log
oe
22 T. Gibson, D. P. O'Leary, L. Nazareth
Alg. No.

Table
The number of failures of the algorithms on the 22 test problems. An algorithm is said to have
"failed" on a particular problem if a line search fails or the maximum allowable number of iterations
(3000 in our case) is exceeded.
We have two choices for ffi k , and a choice of whether or not we will allow m k to decrease
as well as increase. The four variations are
1.
2.
3.
4.
We used four values of m: 5,10,15 and 50, for each algorithm. The results are
summarized in Tables 4.4 - 4.8. More extensive results can be obtained [9].

Table

4.4 shows that these algorithms had roughly the same number of failures
as L-BFGS.

Table

4.5 compares each algorithm to L-BFGS in terms of function evaluations.
For each algorithm and each value of m, the number of times that the algorithm
used as few or fewer function evaluations than L-BFGS is listed relative to the total
number of admissible problems. Problems are admissible if at least one of the two
methods solved it. We observe that in all but three cases, the algorithm used as few
or fewer function evaluations than L-BFGS for over half the test problems.

Table

4.6 compares each algorithm to L-BFGS in terms of time. The entries are
similar to those in Table 4.5. Observe that Algorithms 1-4 did very well in terms of
time, doing as well or better than L-BFGS in nearly every case.
For each problem in each algorithm, we computed the ratio of the number of
function evaluations for the algorithm to the number of function evaluations for L-
BFGS.

Table

4.7 lists the means of these ratios. A mean below 1.0 implies that
the algorithm does better than L-BFGS on average. The average is better for the
algorithms in 6 out of 16 cases for the first four algorithms. Observe, however, that
all the means are close to one.
L-BFGS Variations 23
Alg. No. m= 5
5 19/22 20/22 20/22 21/22
7 8/22 12/22 10/22 10/22
8 12/22 14/22 12/22 15/22
9 6/22 13/22 12/22 16/22
13 3/22 4/22 4/22 4/22
14 2/21 2/22 2/22 2/21
19 2/22 3/22 4/22 4/22

Table
Function Evaluations Comparison. The first number in each entry is the number of times the
algorithm did as well as or better than normal L-BFGS in terms of function evaluations. The second
number is the total number of problems solved by at least one of the two methods (the algorithm
and/or L-BFGS).
We experience savings in terms of time for the first four algorithms. These algorithms
will tend save fewer vectors than L-BFGS since m k is typically less than m;
and so less work is done computing H k g k in these algorithms. Table 4.8 gives the
mean of the ratios of time to solve for each value of m in each algorithm. Note that
most of the ratios are far below one in this case.
These variations did particularly well on problem 7. See [9] for more information.
4.4.3. Disposing of old information: Algorithm 5. We may decide that we
are storing too much old information and that we should stop using it. For example,
we may choose to throw away everything except for the most recent information
whenever we take a big step, since the old information may not be relevant to the
new neighborhood. We use the following test: If the last step length was bigger than
1, dispose of the old information.
The algorithm performed nearly the same as L-BFGS. There was substantial
deviation on only one or two problems for each value of m, and this seemed evenly
divided in terms of better and worse. From Table 4.4, we see that this algorithm
successfully converged on every problem. Table 4.5 shows that it almost always did
as well or better than L-BFGS in terms of function evaluations. However, Table 4.7
shows that the differences were minor. In terms of time, we observe that the algorithm
generally was faster than L-BFGS (Table 4.6), but again, considering the mean ratios
of time (Table 4.8), the differences were minor. The method also does particularly
well on problem 7 [9].
4.4.4. Backing Up in the Update to H: Algorithms 6-11. As discussed
in x2.2, if we always use the most recent s and y in the update, we preserve quadratic
termination regardless of which older values of s and y we use.
T. Gibson, D. P. O'Leary, L. Nazareth
Alg. No.
5 15/22 13/22 14/22 15/22
7 11/22 11/22 10/22 7/22
9 9/22 10/22 7/22 8/22
13 5/22 10/22 13/22 17/22
14 2/21 2/22 2/22 3/21
19 11/22 11/22 17/22 19/22

Table
Time Comparison. The first number in each entry is the number of times the algorithm did as
well as or better than normal L-BFGS in terms of time. The second number is the total number of
problems solved by at least one of the two methods (the algorithm and/or L-BFGS).
Using this idea, we created some algorithms. Under certain conditions, we discard
the next most recent values of s and y in the H although we still use the most recent
s and y vectors and any other vectors that have been saved from previous iterations.
We call this "backing up" because it as if we back-up over the next most recent values
of s and y. These algorithms used the following four tests to trigger backing up:
1. The current iteration is odd.
2. The current iteration is even.
3. A step length of 1.0 was used in the last iteration.
4. kg k k ? kg
In two additional algorithms, we varied situations 3 and 4 by not allowing a back-up
if a back-up was performed on the previous iteration.
The backing up strategy seemed robust in terms of failures. In 4 out of the 6
variations we did for this algorithm, there were no failures at all. See Table 4.4 for
more information.
It is interesting to observe that backing up on odd iterations (Algorithm
backing up on even iterations (Algorithm 7) caused very different results. Backing
up on odd iterations seemed to have almost no effect on the number of function
evaluations (Table 4.7) and little effect on the time (Table 4.8). However, backing up
on even iterations causes much different behavior from L-BFGS. It does worse than
L-BFGS on most problems, but better on a few.
Algorithms were two variations of the same idea: backing up if the
previous step length was one. This wipes out the data from the previous iteration
after it has been used in one update. Both show improvement over L-BFGS in terms
of function evaluations; in fact, these two algorithms have the best function evaluation
ratio for the case (Table 4.7). Unfortunately, these algorithms did not
compete with L-BFGS in terms of time (Table 4.8). There is little difference between
L-BFGS Variations 25
Alg. No. m= 5
6 1.000 1.000 1.000 1.000
9 1.035 1.371 1.005 0.947
14 7.521 7.917 8.288 8.502
19 1.212 1.959 1.242 1.387

Table
Mean function evaluations ratios for each algorithm compared to L-BFGS. Problems for which
either method failed are not used in this mean.
Algorithms probably because there were rarely many steps of length one
is a row.
Algorithms 9 and 11 are also two variations of the same idea: back-up on iteration
the norm of g k is bigger than the norm of g k+1 . There is a larger difference
between the results of 9 and 11 than there was between 8 and 10. In terms of function
evaluation ratios (Table 4.7), Algorithm 11 did better, indicating that it may not be
wise to back-up twice in a row. Both of these did poorly in terms of time as compared
with L-BFGS (Table 4.8).
4.4.5. Merging s and y information in the update: Algorithms 12 and
13. Yet another idea is to "merge" s data so that it takes up less storage and computation
time. By merging, we mean forming some linear combination of various s
vectors. The y vectors would be merged correspondingly. Corollary 2.5 shows that
as long as the most recent s and y are used without merge, old s vectors may be
replaced by any linear combination of the old s vectors in L-BFGS.
We used this idea in the following way: if certain criteria were met, we replaced
the second and third newest s vectors in the collection by their sum, and did similarly
for the y vectors. We used various tests to determine when we would do a merge:
1. Neither of the two vectors to be merged is itself the result of a merge and
the second and third most recent steps taken were of length 1.0.
2. We did not do a merge the last iteration and there are at least two old s
vectors to merge.
The first variation (Algorithm 12) performs almost identically to L-BFGS, especially
in terms of time (Table 4.5). Occasionally it did worse in terms of time

Table

4.6). These observations are also reflected in the other results in Table 4.7 and

Table

4.8. It is likely that very few vectors were merged.
The second variation (Algorithm 13) makes gains in terms of time, especially for
26 T. Gibson, D. P. O'Leary, L. Nazareth
Alg. No.
6 1.007 0.983 0.977 0.995
9 1.032 1.220 1.043 1.173
14 4.585 3.703 3.228 2.417

Table
Mean time ratios for each algorithm compared to L-BFGS. Problems for which either method
failed are not used in this mean.
the larger values of m (Table 4.6 and Table 4.8). Unfortunately, this reflects only a
saving in the amount of linear algebra required. The number of function evaluations
generally is larger for this algorithm than L-BFGS (Table 4.5 and Table 4.7).
4.4.6. Skipping Updates to H: Algorithms 14-16. If every other update
to H is skipped and a step length of one is always chosen, BFGS will terminate in
2n iterations on a strictly convex quadratic function. The same holds true when
doing an exact line search. (See x3.) Unfortunately, neither property holds in the
limited-memory case. We will, however, try some algorithms motivated by this idea.
The idea is that, every so often, we do not use the current s and y to update H,
and instead just use the old H. There are three variations on this theme.
1. Skip update on odd iterations.
2. Skip update on even iterations.
3. Skip update if kg k+1k ? kg k k.
As with the algorithms that did back-ups, the results of the skipping on odd or
even iterations were quite different. Skipping on odd updates (Algorithm 14) did
extremely well for every value of m on only two problems: 1 and 12. Otherwise, it did
very badly. Skipping on even updates (Algorithm 15) performed somewhat better. It
did extremely well on problem 7 but not on problems 1 and 12. It also did better than
L-BFGS in terms of time on more occasions than Algorithm 14 (Table 4.6). Neither
did well in terms of function evaluations, but the mean ratios for function evaluations

Table

4.7) and time (Table 4.8) were usually far greater than one.
Skipping the update if the norm of g increased (Algorithm 16) did not do well at
all. It only did better in terms of function evaluations for one problem for each value
of m (

Table

4.5) and rarely did better in terms of time (Table 4.6). It ratios were
very bad for function evaluations (Table 4.7) and time (Table 4.8)
L-BFGS Variations 27
4.4.7. Combined Methods: Algorithms 17-21. We did some experimentation
with combinations of methods described in the previous sections.
In Algorithm 17, we combined Algorithms 5 and 8: we dispose of old information
and back-up on the next iterations if the step length is greater than one. Essentially
we are assuming that we have stepped out of the region being modeled by the quasi-Newton
matrix if we take a long step and we should thus rid the quasi-Newton matrix
of that information. This algorithm did well in terms of function evaluations, having
mean ratios of less than one for three values of m (Table 4.7), but it did not do as
well in terms of time.
In Algorithms 19-21, we combined merging and varying m. These algorithms did
well in terms of time for larger m (Table 4.8) but not in terms of function evaluations

Table

4.7).
5. Conclusions. There is a spectrum of quasi-Newton methods, ranging from
those that require the storage of an n \Theta n approximate Hessian (e.g. the Broyden fam-
ily) to those that require only the storage of a few vectors (e.g. conjugate gradients).
Limited-memory quasi-Newton methods fall in between these extremes in terms of
performance and storage. There are other methods that fall into the middle ground;
for example, conjugate gradient methods such as those proposed by Shanno [27] and
Nazareth [20], the truncated-Newton method [24, 6] and the partitioned quasi-Newton
method [13].
We have characterized which limited-memory quasi-Newton methods fitting a general
form (2.1) have the property of producing conjugate search directions on convex
quadratics. We have shown that limited-memory BFGS is the only Broyden family
member that has a limited-memory analog with this property. We also considered
update-skipping, something that may seem attractive in a parallel environment. We
show that update skipping on quadratic problems is acceptable for full-memory Broyden
class members in that it only delays termination, but that we lose the property
of finite termination if we both limit memory and skip updates.
We have also introduced some simple-to-implement modifications of the standard
limited-memoryBFGS algorithm that seem to behave well on some practical problems.


Appendix

A. Line Search Parameters. Table A.1 give the line search parameters
used for our code. Note that in the first iteration, the initial steplength is
rather than 1.0.
Variable Value Description
STP 1.0 Step length to try first.
\Gamma4 Value of ! 1 in Wolfe conditions.
GTOL 0.9 Value of ! 2 in Wolfe conditions.
Relative width of interval of uncertainty.
Maximum number of function evaluations.

Table

A.1
Line Search Parameters


Appendix

B. Pseudo-Code.
B.1. L-BFGS: Algorithm 0. The pseudo-code for the computation of d
\GammaH k g k at iteration k for L-BFGS is given in Figure B.2. The update of H is also
handled implicitly in this computation.
28 T. Gibson, D. P. O'Leary, L. Nazareth
if (sze ==
else
% This is needed for Step 3 before we overwrite Stg and Ytg
Fig. B.1. MATLAB pseudo-code for the computation of d = Hg in L-BFGS. sze is the number
of s vectors available for the update this iteration and oldsze is the number of s vectors that were
available the previous iteration. For L-BFGS, sze is chosen as the minimum of oldsze
(the limited-memory constant).
B.2. Varying m iteratively: Algorithms 1-4. Suppose that m k denotes the
number of (s; y) pairs to be used in the kth update. Then simply chose sze as the
minimum of oldsze computing d k .
B.3. Disposing of old information: Algorithm 5. If the disposal criterion
is met at iteration k, set oldsze to zero and sze to one before computing d k .
B.4. Backing Up in the Update to H: Algorithms 6-11. If we are to
back-up at iterations k, set oldsze to the one less than the previous value of sze and
set sze as the minimum of oldsze m, as usual.
B.5. Merging s and y information in the update: Algorithms 12 and 13.
Merging is the most complicated variation to handle. Before we determine the newest
sze and before we compute d k , we execute the pseudo-code given in Figure B.1. We
then set oldsze to one less than the previous value of sze and set sze as the minimum
of oldsze m, as usual. We are assuming we are at iteration k, but that the
Variations 29
newest values of s and y have not yet been added to S and Y.
Execute before choosing new value for sze and before computing d
Fig. B.2. MATLAB pseudo-code for the merge variation. This fixes the values of the components
that are used in the computation of d k .
B.6. Skipping Updates to H: Algorithms 14-16. To skip the update at
iteration k, set sze to oldsze. Compute Stg and Ytg before Step 0 and then skip to
Step 8 and continue.



--R

ftp://thales.


Test functions for unconstrained minimization
Performance of a multifrontal scheme for partially separable optimization

Numerical Methods for Unconstrained Optimization and Nonlinear Equations
An optimal positive definite update for sparse Hessian matrices
http://www.

Matrix Computations
Private communication to authors of
Partitioned variable metric updates for large structured optimization problems
A theoretical and experimental study of the symmetric rank one update
On the limited memory BFGS method for large scale optimization
Linear and Nonlinear Programming


A numerical study of the limited memory BFGS method and the truncated-Newton method for large scale optimization
A relationship between BFGS and conjugate gradient algorithms and its implications for new algorithms

Updating quasi-Newton matrices with limited storage

A discrete Newton algorithm for minimizing a function of many variables

Quadratic termination properties of minimization algorithms I.
Conjugate gradient methods with inexact line searches
An error in specifying problem CHNROSNB.


--TR

--CTR
Sun Linping, Updating the Self-Scaling Symmetric Rank One Algorithm with Limited Memory for Large-Scale Unconstrained Optimization, Computational Optimization and Applications, v.27 n.1, p.23-29, January 2004
Adi Ditkowski , Gadi Fibich , Nir Gavish, Efficient Solution of A, Using A-1, Journal of Scientific Computing, v.32 n.1, p.29-44, July      2007
M. Al-Baali, Extra-updates criterion for the limited memory BFGS Algorithm for large scale nonlinear optimization, Journal of Complexity, v.18 n.2, p.557-572, June 2002
update skipping;limited-memory;broyden family;BFGS;minimization;quasi-Newton