Solving the Trust-Region Subproblem using the Lanczos Method.
The approximate minimization of a quadratic function within an ellipsoidal trust region is an important subproblem for many nonlinear programming methods. When the number of variables is large, the most widely used strategy is to trace the path of conjugate gradient iterates either to convergence or until it reaches the trust-region boundary. In this paper, we investigate ways of continuing the process once the boundary has been encountered. The key is to observe that the trust-region problem within the currently generated Krylov subspace has a very special structure which enables it to be solved very efficiently. We compare the new strategy with existing methods. The resulting software package is available as HSL_VF05 within the Harwell Subroutine Library.
Introduction
Trust-region methods for unconstrained minimization are blessed with both strong theoretical
convergence properties and a good reputation in practice. The main computational step in these
methods is to find an approximate minimizer of some model of the true objective function within
a "trust" region for which a suitable norm of the correction lies within a given bound. This
restriction is known as the trust-region constraint, and the bound on the norm is its radius. The
radius is adjusted so that successive model problems mimic the true objective within the trust
region.
The most widely-used models are quadratic approximations to the objective function, as these
are simple to manipulate and may lead to rapid convergence of the underlying method. From
a theoretical point of view, the norm which defines the trust region is irrelevant so long as it
"uniformly" related to the ' 2 norm. From a practical perspective, this choice certainly affects the
subproblem, and thus the methods one can consider when solving it. The most popular practical
choices are the ' 2 - and ' 1 -norms, and weighted variants thereof. In our opinion, it is important
that the choice of norm reflects the underlying geometry of the problem; simply picking the
may not be adequate when the problem is large, and the eigenvalues of the Hessian of
the model widely spread. We believe that weighting the norm is essential for many large-scale
problems.
In this paper, we consider the solution of the quadratic-model trust-region subproblem in a
weighted ' 2 -norm. We are interested in solving large problems, and thus cannot rely solely on
factorizations of the matrices involved. We thus concentrate on iterative methods. If the model of
the Hessian is known to be positive definite and the trust-region radius sufficiently large that the
trust-region constraint is inactive at the unconstrained minimizer of the model, the obvious way
to solve the problem is to use the preconditioned conjugate-gradient method. Note that the role
of the preconditioner here is the same as the role of the norm used for the trust-region, namely to
change the underlying geometry so that the Hessian in the rescaled space is better conditioned.
Thus, it will come as no surprise that the two should be intimately connected. Formally, we shall
require that the weighting in the ' 2 -norm and the preconditioning are performed by the same
matrix.
When the radius is smaller than a critical value, the unconstrained minimizer of the model
will no longer lie within the trust-region, and thus the required solution will lie on the trust-region
boundary. The simplest strategy in this case is to consider the piecewise linear path connecting
the conjugate-gradient iterates, and to stop at the point where this path leaves the trust region.
Such a strategy was first proposed independently by Steihaug (1983) and Toint (1981), and
we shall refer to the terminating point as the Steihaug-Toint point. Remarkably, it is easy to
establish the global convergence of a trust-region method based on such a simple strategy. The
key is that global convergence may be proved provided that the accepted estimate of the solution
has a model value no larger than at the Cauchy point (see Powell, 1975). The Cauchy point
is simply the minimizer of the model within the trust-region along the preconditioned steepest-descent
direction. As the first segment on the piecewise-linear conjugate-gradient path gives
N. I. M. Gould, S. Lucidi, M. Roma and Ph. L. Toint
precisely this point, and as the model value is monotonically decreasing along the entire path,
the Steihaug-Toint strategy ensures convergence.
If the model Hessian is indefinite, the solution must also lie on the trust-region boundary.
This case may also be simply handled using preconditioned conjugate gradients. Once again
the piecewise linear path is followed until either it leaves the trust-region, or a segment with
negative curvature is found (a vector p is a direction of negative curvature if the inner product
is the model Hessian). In the latter case, the path is continued downhill
along this direction of negative curvature as far as the constraint boundary. This variant was
proposed by Steihaug (1983), while Toint (1981) suggests simply returning to the Cauchy point.
As before, global convergence is ensured as either of these terminating points, as the objective
function values there are no larger than at the Cauchy point. For consistency with the previous
paragraph, we shall continue to refer to the terminating point in Steihaug's algorithm as the
Steihaug-Toint point, although strictly Toint's point in this case may be different.
The Steihaug-Toint method is basically unconcerned with the trust region until it blunders
into its boundary and stops. This is rather unfortunate, particularly as considerable experience
has shown that this frequently happens during the first few, and often the first, iteration(s) when
negative curvature is present. The resulting step is then barely, if at all, better than the Cauchy
direction, and this may lead to a slow but globally convergent algorithm in theory and a barely
convergent method in practice. In this paper, we consider an alternative which aims to avoid
this drawback by trying harder to solve the subproblem when the boundary is encountered, while
maintaining the efficiencies of the conjugate gradient method so long as the iterates lie interior.
The mechanism we use is the Lanczos method.
The paper is organized as follows. In Section 2 we formally define the problem and any
notation that we will use. The basis of our new method is given in Section 3, while in Section 4,
we will review basic properties of the preconditioned conjugate-gradient and Lanczos methods.
Our new method is given in detail in Section 5. Some numerical experiments demonstrating the
effectiveness of the approach are given in Section 6, and a number of conclusions and perspectives
are drawn in the final section.
Solving the trust-region subproblem using the Lanczos method 3
2 The trust-region subproblem and its solution
Let M be a symmetric positive-definite easily-invertible approximation to the symmetric matrix
H. Furthermore, define the M-norm of a vector as
where h\Delta; \Deltai is the usual Euclidean inner product. In this paper, we consider the M-norm trust-region
problem
minimize
subject to ksk M - \Delta; (2.1)
for some vector g and radius \Delta ? 0.
A global solution to the problem is characterized by the following result.
Theorem 2.1 (Gay, 1981, Sorensen, 1982) Any global minimizer s M of q(s) subject to
satisfies the equation
positive semi-definite, - M - 0 and - M (ks
is positive definite, s M is unique.
This result is the basis of a series of related methods for solving the problem which are appropriate
when forming factorizations of H(-) j H +-M for a number of different values of - is realistic.
For then, either the solution lies interior, and hence - g, or the solution lies
on the boundary and - M satisfies the nonlinear equation
denotes the pseudo-inverse of H. Equation (2.3) is straightforward to solve using a
safeguarded Newton iteration, except in the so-called hard case for which g lies in the null-space of
H(- M ). In this case, an additional vector in the range-space of H(- M ) may be required if a solution
on the trust-region boundary is sought. Goldfeldt, Quandt and Trotter (1966), Hebden (1973)
and Gay (1981) all proposed algorithms of this form. The most sophisticated algorithm to date,
by Mor'e and Sorensen (1983), is available as subroutine GQTPAR in the MINPACK-2 package, and
guarantees that a nearly optimal solution will be obtained after a finite number of factorizations.
While such algorithms are appropriate for large problems with special Hessian structure -
such as for band matrices - the demands of a factorization at each iteration limits their applicability
for general large problems. It is for this reason that methods which do not require
factorizations are of interest.
Throughout the paper, we shall denote the k by k identity matrix by I k , and its j-th column
by e j . A set of vectors fq i g are said to be M-orthonormal if hq the Kronecker delta,
4 N. I. M. Gould, S. Lucidi, M. Roma and Ph. L. Toint
and the matrix Q formed from these vectors is an M-orthonormal matrix. The set
of vectors fp i g are H-conjugate (or H-orthogonal) if hp
3 A new algorithm for large-scale trust-region subproblems
To set the scene for this paper, we recall that the Cauchy point may be defined as the solution
to the problem
minimize
subject to ksk M - \Delta; (3.1)
that is as the minimizer of q within the trust region where s is restricted to the 1-dimensional
subspace span \Phi
. The dogleg methods (see, Powell, 1970, Dennis and Mei, 1979) aim to
solve the same problem over a one-dimensional arc, while Byrd, Schnabel and Schultz (1985) do
the same over a two-dimensional subspace. In each of these cases the solution is easy to find as
the search space is small. The difficulty with the general problem (2.1) is that the search space
R n is large. This leads immediately to the possibility of solving a compromise problem
minimize
subject to ksk M - \Delta; (3.2)
where S is a specially chosen subspace of R n .
Now consider the Steihaug-Toint algorithm at an iteration k before the trust-region boundary
is encountered. In this case, the point s k+1 is the solution to (3.2) with the set
span
the Krylov space generated by the starting vector M \Gamma1 g and matrix M \Gamma1 H. That is, the
Steihaug-Toint algorithm gradually widens the search space using the very efficient preconditioned
conjugate gradient method. However, as soon as the Steihaug-Toint algorithm moves
across the trust-region boundary, the terminating point s k+1 no longer necessarily solves the
problem in (3.3), indeed it is very unlikely to do so when k ? 0. (As the iterates generated by
the method increase in M-norm, once an iterate leaves the trust region, the solution to (3.3),
and thus (2.1), must lie on the boundary. See, Steihaug, 1983, Theorem 2.1, for details). Can
we do better? Yes, by recalling that the preconditioned conjugate gradient and Lanczos methods
generate different bases for the same Krylov space.
Solving the trust-region subproblem using the Lanczos method 5
4 The preconditioned conjugate-gradient and Lanczos methods
The preconditioned conjugate-gradient and Lanczos methods may be viewed as efficient techniques
for constructing different bases for the same Krylov space, K k . The conjugate gradient
method aims for an H-conjugate basis, while the Lanczos method obtains an M-orthonormal
basis.
Algorithm 4.1: The preconditioned conjugate gradient method
perform the iteration,
Algorithm 4.2: Preconditioned Lanczos method
perform the iteration,
The conjugate gradient method generates the basis
from Algorithm 4.1, while the Lanczos method generates the basis
6 N. I. M. Gould, S. Lucidi, M. Roma and Ph. L. Toint
with Algorithm 4.2. The Lanczos iteration is often written in the more compact form
k+1 and (4.14)
is the matrix (q and the matrix
@
A
is tridiagonal. It then follows directly that
The two methods are intimately related. In particular, so long as the conjugate-gradient
iteration does not break down, the Lanczos vectors may be recovered from the conjugate-gradient
iterates as
while the Lanczos tridiagonal may be expressed as
ff
ff k\Gamma2
ff
A
The conjugate gradient iteration may breakdown if hp which can only occur if H is
not positive definite, and will stop if hg On the other hand, the Lanczos iteration can
only fail if K j is an invariant subspace for M \Gamma1 H.
If q(s) is convex in the manifold K j+1 , the minimizer s j+1 of q in this manifold satisfies
so long as the initial value s chosen. Thus this estimate is easy to recur from the conjugate-gradient
iteration. The minimizers in successive manifolds may also be easily obtained using the
Lanczos process, although the conjugate-gradient iteration is slightly less expensive, and thus to
be preferred.
Solving the trust-region subproblem using the Lanczos method 7
The vector g j+1 in the conjugate gradient method gives the gradient of q(s) at s j+1 . It is
quite common to stop the method as soon as this gradient is sufficiently small, and the method
naturally records the M \Gamma1 -norm of the gradient, kg This norm is also available
in the Lanczos method as
solves the tridiagonal linear system T k The last component, he k+1
of h k is available as a further by-product.
5 The truncated Lanczos approach
Rather than use the preconditioned conjugate gradient basis fp for S, we shall use
the equivalent Lanczos M-orthonormal basis fq g. The Lanczos basis has previously
been used by Nash (1984) - to convexify the quadratic model - and Lucidi and Roma (1997)
- to compute good directions of negative curvature - within linesearch based method for unconstrained
minimization. We shall consider vectors of the form
and seek
solves the problem
minimize
s 2S
subject to ksk M - \Delta: (5.2)
It then follows directly from (4.15), (4.17) and (4.18) that h k solves the problem
minimize
subject to khk 2 - \Delta: (5.3)
There are a number of crucial observations to be made here. Firstly, it is important to note that
the resulting trust-region problem involves the two-norm rather than the M-norm. Secondly,
as T k is tridiagonal, it is feasible to use the Mor'e-Sorensen algorithm to compute the model
minimizer even when n is large. Thirdly, having found h k , the matrix Q k is needed to recover
thus the Lanczos vectors will either need to be saved on backing store or regenerated.
As we shall see, we only need Q once we are satisfied that continuing the Lanczos process will
give little extra benefit. Fourthly, one would hope that as a sequence of such problems may be
solved, and as T k only changes by the addition of an extra diagonal and superdiagonal entry,
solution data from one subproblem may be useful for starting the next. We consider this issue in
Section 5.2.
The basic trust-region solution classification theorem, Theorem 2.1, shows that
8 N. I. M. Gould, S. Lucidi, M. Roma and Ph. L. Toint
I k+1 is positive semi-definite, What does this tell
us about s k ? Firstly, using (4.17), (4.18) and (5.4) we have
and additionally that
Comparing these with the trust-region classification theorem, we see that s k is the Galerkin
approximation to s M from the space spanned by Q k .
We may then ask how good the approximation is. In particular, what is the error (H
g? The simplest way of measuring this error would be to calculate h k and - k by
solving (5.3), then to recover s k as Q k h k and finally to substitute s k and - k into (H +-M)s+ g.
However this is inconvenient as it requires that we have easy access to Q k . Fortunately there is
a far better way.
Theorem 5.1
and
Proof. We have that
iw k+1 from (4.14)
iw k+1 from (5.4)
This then directly gives (5.6), and (5.7) follows from the M \Gamma1 -orthonormality of w k+1 . 2
Therefore we can indirectly measure the error (in the M \Gamma1 -norm) knowing simply fl k+1 and the
last component of h k , and we do not need s k or Q k at all. Observant readers will notice the strong
similarity between this error estimate and the estimate (4.22) for the gradient of the model in
the Lanczos method, but this is not at all surprising as the two methods are aiming for the same
point if the trust-region radius is large enough. An interpretation of (5.7) is also identical to that
of (4.22). The error will be small when either of fl k+1 or the last component of h k is small.
We now consider the problem (5.3) in more detail. We say that a symmetric tridiagonal matrix
is degenerate if one or more of its off-diagonal entries is zero; otherwise it is non-degenerate. We
then have the following preliminary result.
Solving the trust-region subproblem using the Lanczos method 9
Lemma 5.2 (See also, Parlett, 1980, Theorem 7.9.5) Suppose that the tridiagonal matrix
T is non-degenerate, and that v is an eigenvector of T . Then the first component of v is
nonzero.
Proof. By definition
for some eigenvector '. Suppose that the first component of v is zero. Considering the first
component of (5.8), we have that the second component of v is zero as T is tridiagonal and
non-degenerate. Repeating this argument for the i-th component of (5.8), we deduce that the
the 1-st component of v is zero for all i, and hence that contradicts the
assumption that v is an eigenvector, and so the first component of v cannot be zero. 2
This immediately yields the following useful result.
Theorem 5.3 Suppose that T k is non-degenerate. Then the hard case cannot occur for the
subproblem (5.3).
Proof. Suppose the hard case occurs. Then, by definition, fl 0 e 1 is orthogonal to v k , the
eigenvector corresponding to the leftmost eigenvalue, \Gamma' k , of T k . Thus, the first component of
v k is zero, which, following Lemma 5.2, contradicts the assumption that v k is an eigenvector.
Thus the hard case cannot occur. 2
This result is important as it suggests that the full power of the Mor'e and Sorensen (1983)
algorithm is not needed to solve (5.3). We shall return to this in Section 5.2. We also have an
immediate corollary.
Corollary 5.4 Suppose that T n\Gamma1 is non-degenerate. Then the hard case cannot occur for
the original problem (2.1).
Proof. When the columns of Q forms a basis for R n , Thus the problems (2.1)
and (5.2) are identical, and (5.2) and (5.3) are related through a nonsingular transformation.
The result then follows directly from Theorem 5.3 in the case
Thus, if the hard case occurs for (2:1), the Lanczos tridiagonal must degenerate at some stage.
N. I. M. Gould, S. Lucidi, M. Roma and Ph. L. Toint
Theorem 5.5 Suppose that T k is non-degenerate, that h k and - k satisfy (5.4) and that
I k+1 is positive semi-definite. Then I k+1 is positive definite.
Proof. Suppose that T k I k+1 is singular. Then there is a nonzero eigenvector v k for
which Hence, combining this with (5.4) reveals that
and hence that the first component of v k is zero. But this contradicts Lemma 5.2. Hence
I k+1 is both positive semi-definite and nonsingular, and thus positive definite. 2
This result implies that (5.4) has a unique solution. We now consider this solution.
Theorem 5.6 Suppose that he k+1
Proof. Suppose that T k is not degenerate. As the k 1-st component of h k is zero, then the
non-degeneracy of T k and the k 1-st equation of (5.4), we deduce that the k-th component
of h k is zero. Repeating this argument for the 1-st equation of (5.4), we deduce that the
i-th component of h k is zero for 1 - i - k, and hence that h contradicts the
first equation of (5.4), and thus T k must be degenerate. 2
Thus we see that of the two possibilities suggested by Theorem 5.1 for obtaining an s k for which
will be the possibility fl that occurs before he k+1
Theorem 5.7 Suppose that the hard case does not occur for (2.1), and that fl
Proof. If the Krylov space K k is an invariant subspace M \Gamma1 H, and by construction
the first basis element of this space is M \Gamma1 g. As the hard case does not occur for (2.1), this
space must also contain at least one eigenvector corresponding to the leftmost eigenvalue, \Gamma',
of M \Gamma1 H. Thus one of the eigenvalues of T k must be \Gamma', and - k - ' as T k +- k I k+1 is positive
semi-definite. But this implies that H positive semi-definite, which combines with
(5.1), (5.5) and Theorem 5.1 with to show that s k satisfies the optimality conditions
shown in Theorem 2.1. 2
Thus we see that in the easy case, the required solution will be obtained from the first non-degenerate
block of the Lanczos tridiagonal. It remains for us to consider the hard case. In
Solving the trust-region subproblem using the Lanczos method 11
view of Corollary 5.4, this case can only occur when T k is degenerate. Suppose therefore that T k
degenerates into ' blocks of the form
@
A
where each of the T k i
defines an invariant subspace for M \Gamma1 H and where the last block T k '
is
the first to yield the leftmost eigenvalue, \Gamma', of M \Gamma1 H. Then there are two cases to consider.
Theorem 5.8 Suppose that the hard case occurs for (2.1), that T k is as described by (5.9),
and the last block T k '
is the first to yield the leftmost eigenvalue, \Gamma', of M \Gamma1 H. Then,
1. if ' - k1 , a solution to (2.1) is given by s k1 solves the positive-definite
system
2. if ' ? - k1 , a solution to (2.1) is given by s
@
A
h is the solution of the nonsingular tridiagonal system
u is an eigenvector of T k '
corresponding to \Gamma', and ff is chosen so that
Proof. In case 1, H+- positive semi-definite as - k 1 - ', and the remaining optimality
conditions are satisfied as k1 +1 is positive
definite follows from Theorem 5.5. In case 2, H+'M is positive semi-definite. Furthermore, as
is easy to show that khk 2 ! kh k1 k 2 - \Delta, and hence that there is a root ff for which
ks \Delta. Finally, as each Q k i
defines an invariant subspace, HQ k i
Writing
u, we therefore have
N. I. M. Gould, S. Lucidi, M. Roma and Ph. L. Toint
and
all of the optimality conditions for (5.2). 2
Notice that to obtain s k as described in this theorem, we only require the Lanczos vectors corresponding
to blocks one and, perhaps, ' of T k .
We do not claim that to solve the problem as outlined in Theorem 5.8 is realistic, as it relies
on our being sure that we have located the left-most eigenvalue of M \Gamma1 H. With Lanczos-type
methods, one cannot normally guarantee that all eigenvalues, including the leftmost, will be
found unless one ensures that all invariant subspaces have been investigated, and this may prove
to be very expensive for large problems. In particular, the Lanczos algorithm, Algorithm 4.2,
terminates each time an invariant subspace has been determined, and must be restarted using
a vector q which is M-orthonormal to the previous Lanczos directions. Such a vector may be
obtained from the Gram-Schmidt process by re-orthonormalizing a suitable vector - a vector
with some component M-orthogonal to the existing invariant subspaces, perhaps a random vector
- with respect to the previous Lanczos directions, which means that these directions will have to
be regenerated or reread from backing store. Thus, while this form of the solution is of theoretical
interest, it is unlikely to be of practical interest if a cheap approximation to the solution is all
that is required.
5.1 The algorithm
We may now outline our algorithm, Algorithm 5.1, the generalized Lanczos trust region (GLTR)
method. We stress that, as our goal is merely to improve upon the value delivered by the
Steihaug-Toint method, we do not use the full power of Theorem 5.8, and are content just
to investigate the first invariant subspace produced by the Lanczos algorithm. In almost all
cases, this subspace contains the global solution to the problem, and the complications and costs
required to implement a method based on Theorem 5.8 are, we believe, prohibitive in our context.
Solving the trust-region subproblem using the Lanczos method 13
Algorithm 5.1: The generalized Lanczos trust region method
. Set the flag INTERIOR as
true. For convergence, perform the iteration,
using (4.20)
If INTERIOR is true, but ff k - 0 or ks k reset INTERIOR to false.
If INTERIOR is true
else
solve the tridiagonal trust-region subproblem (5.3) to obtain h k
If INTERIOR is true
test for convergence using the residual kg
else
test for convergence using the value fl k+1 jhe
If INTERIOR is false, recover s by rerunning the recurrences or obtaining Q k from
backing store.
When recovering s by rerunning the recurrences, economies can be made by saving the
during the first pass, and reusing them during the second. A potentially bigger saving
may be made if one is prepared to accept a slightly inferior value of the objective function. The
idea is simply to save the value of q at each iteration. On convergence, one looks back through
this list to find an iteration, ' say, for which a required percentage of the best value was obtained,
recompute h ' and then accept s as the required estimate of the solution. If the required
percentage occurs at an iteration before the boundary is encountered, both the final point before
the boundary and the Steihaug-Toint point are suitable and available without the need for the
second pass.
We note that we have used the conjugate-gradient method (Algorithm 4.1) to generate the
Lanczos vectors. If the inner-product hp k ; Hp k i proves to be tiny, it is easy to continue using the
Lanczos method (Algorithm 4.2) itself; the vectors
14 N. I. M. Gould, S. Lucidi, M. Roma and Ph. L. Toint
required to continue the Lanczos recurrence (4.11) are directly calculable from conjugate-gradient
method.
At each stage of both the Steihaug-Toint algorithm and our GLTR method (Algorithm 5.1),
we need to calculate ks k . This issue is not discussed by Steihaug as it is implicitly
assumed that M is available. However, it may be the case that all that is actually available is a
procedure which returns M \Gamma1 v for a given input v, and thus M is unavailable. Fortunately this
is not a significant drawback as it is possible to calculate ks k +ffp k k M from available information.
To see this, observe that
ks
ks
and thus that we can find ks k+1 k 2
M from ks k k 2
M so long as we already know hs k ; Mp k i and
M . But it is straightforward to show that these quantities may be calculated from the pair
of recurrences
and (5.12)
where, of course, hg k ; v k i has already been calculated as part of the preconditioned conjugate-gradient
method.
5.2 Solving the nondegenerate tridiagonal trust-region subproblem
In view of Theorem 5.3, the nondegenerate tridiagonal trust-region subproblem (5.3) is, in theory,
easier to solve than the general problem. This is so both because the Hessian is tridiagonal (and
thus very inexpensive to factorize), and because the hard case cannot occur. We should be
cautious here, because the so-called "almost" hard case - which occurs when g only has a tiny
component in the range-space of H(- M ) - may still happen, and the trust-region problem in
this case is naturally ill conditioned and thus likely to be difficult to solve.
The Mor'e and Sorensen (1983) algorithm is based on being able to form factorizations of the
model Hessian (which is certainly the case here as is tridiagonal), but does not try
to calculate the leftmost eigenvalue of the pencil H + -M . In the tridiagonal case, computing
the extreme eigenvalues is straightforward, particularly if a sequence of related problems are to
be solved. Thus, rather than using the Mor'e and Sorensen algorithm, we prefer the following
method.
We restrict ourselves to the case where the solution lies on the trust-region boundary - we
will only switch to this approach when the conjugate gradient iteration leaves the trust region.
The basic iteration is identical to that proposed by Mor'e and Sorensen (1983), namely to apply
Newton's method to
where
Solving the trust-region subproblem using the Lanczos method 15
Recalling that we denote the leftmost eigenvalue of T k by \Gamma' k , the main difference between
our approach and Mor'e and Sorensen's is that we always start from some point in the interval
this interval is characterized by both being positive definite and
as then the resulting Newton iteration is globally linearly, and asymptotically
quadratically, convergent without any further safeguards. The Newton iteration is performed
using Algorithm 5.2.
Algorithm 5.2: Newton's method to solve
1. Factorize are unit bidiagonal and diagonal
matrices, respectively.
2. Solve BDB T
3. Solve
4. Replace - by -
The Newton correction in Step 4 of this algorithm is given by
while the exact form given is obtained by using the identity
where w is as computed in Step 3. It is slightly more efficient to pick B to be unit upper-
bidiagonal rather than unit lower-bidiagonal, as then the Step 2 simplifies to B T
because of the structure of the right-hand side.
To obtain a suitable starting value, two possibilities are considered. Firstly, we attempt
to use the solution value - k\Gamma1 from the previous subproblem. Recall that T k is merely T
with an appended row and column. As we already have a factorization of T
trivial to obtain that of T and thus to determine if the latter is positive definite.
If turns out to be positive definite, h k (- k\Gamma1 ) is computed from (5.15) and if
is used to start the Newton iteration.
Secondly, if - k\Gamma1 is unsuitable, we monitor T k to see if it is indefinite. This is trivial, as for
instance, the matrix is positive definite so long as all of the ff i , generated by the
conjugate-gradient method are positive. If T k is positive definite, we start the Newton iteration
with the value which by assumption gives kh k (0)k 2 - \Delta as the unconstrained solution lies
outside the trust region. Otherwise, we determine the leftmost eigenvalue, \Gamma' k , of T k , and start
N. I. M. Gould, S. Lucidi, M. Roma and Ph. L. Toint
with positive number chosen to make T k numerically
"just" positive definite. By this we mean, that its BDB T factorization should exist, but that ffl
should be as small as possible. We have found that a value (1
is the unit
roundoff, is almost always suitable, but have added the precaution of multiplying this value by
increasing powers of 2 so long as the factorization fails.
If we need to compute the leftmost eigenvalue of T k , we use an iteration based upon the
last-pivot function proposed by Parlett and Reid (1981). The last-pivot function, ffi k ('), is simply
the value of the last diagonal entry of the BDB T factor D k (-) of T k \Gamma 'I k+1 . This value will be
zero, and the other diagonal entries positive, when
of uncertainty [' l ; ' u ] is placed around the required root. The initial interval is given by the
Gersgorin bounds on the leftmost eigenvalue. When it is known, the leftmost eigenvalue, \Gamma'
of T k\Gamma1 may be used to improve the lower bound, because of the Cauchy interlacing property
of the eigenvalues of T k\Gamma1 and T k (see, for instance, Parlett, 1980, Theorem 10.1.2). Given an
initial estimate of ' k , an improvement may be sought by applying Newton's method to ffi k ('); the
derivative of ffi k is easy to obtain by recurrence. However, as Parlett and Reid point out,
and thus has a pole at Hence it is better to choose the new point by fitting the model
to the function and derivative value at the current ', and then to pick the new iterate as the
larger root of ffi M
('). If the new iterate lies outside the interval of uncertainty, it is replaced by
the midpoint of the interval. The interval is then contracted by computing ffi k at the new iterate,
and replacing the appropriate endpoint by the iterate. The iteration is stopped if the length of
the interval or the value of
If ' k\Gamma1 is known, the initial iterate chosen as ' positive ffl -
successive iterates generated from (5.16), the iterates convergence globally, and asymptotically
superlinearly, from the left. If the Newton iteration is used, the required root is frequently
obscured, and the scheme resorts to interval bisection. Thus the Parlett and Reid scheme is to
be preferred.
Other means of locating the required eigenvalue, based on using the determinant det(T
were tried, but proved to be less reliable because of the huge numerical
range (and thus potential overflow) of the determinant.
Solving the trust-region subproblem using the Lanczos method 17
6 Numerical experiments
The algorithm sketched in Sections 5.1 and 5.2 has been implemented as a Fortran 90 module,
HSL VF05, within the Harwell Subroutine Library (1998).
As our main interest is in using the methods described in this paper within a trust-region
algorithm, we are particularly concerned with two issues. Firstly, can we obtain significantly
better values of the model by finding better approximations to its solution than the Steihaug-
Toint method? And secondly, do better approximations to the minimizer of the model necessarily
translate into fewer iterations of the trust-region method? In this section, we address these
outstanding questions.
Throughout, we will consider the basic problem of minimizing an objective f(x) of n real
variables x. We shall use the following standard trust-region method.
Algorithm 6.1: Standard Trust-Region Algorithm
An initial point x 0 and an initial trust-region radius \Delta 0 are given, as are constants ffl g ,
are required to satisfy the conditions
1. Stop if kr x f(x k )k 2 - ffl g .
2. Define a second-order Taylor series model q k and a positive-definite preconditioner
. Compute a step s k to "sufficiently reduce the model" q k within the trust-region
3. Compute the ratio
4. Set
Increment k by one and go to Step 1.
We choose the specific values ffl set an
upper limit of n iterations. The step s k in step 2 is computed using either Algorithm 5.1 or the
Steihaug-Toint algorithm. Convergence in both algorithms for the subproblem occurs as soon as
N. I. M. Gould, S. Lucidi, M. Roma and Ph. L. Toint
or if more than n iterations have been performed. In addition, of course, the Steihaug-Toint
algorithm terminates as soon as the boundary is crossed.
All our tests were performed on an IBM RISC System/6000 3BT workstation with 64 Mega-bytes
of RAM; the codes are all double precision Fortran 90, compiled under xlf90 with -O
optimization, and IBM library BLAS are used. The test examples we consider are the larger
examples from the CUTE test set (see Bongartz, Conn, Gould and Toint, 1995) for which negative
curvature is frequently encountered. Tests were terminated if more than thirty CPU minutes
elapsed.
6.1 Can we get much better model values than Steihaug-Toint?
We first consider problems of the form (2.1). Our test examples are generated by running Algorithm
6.1 on the CUTE set for 10 iterations, and taking the trust-region subproblem at iteration
as our example. The idea here is to simulate the kind of subproblems which occur in practice,
not those which result at the starting point for the algorithm as such points frequently have
special (favourable) properties.
Our aim is to see whether there is any significant advantage in continuing the minimization of
the trust-region subproblem once the boundary of the trust region has been encountered. We ran
HSL VF05 to convergence, stopping when kg more than
iterations had been performed.
In all of the experiments reported here, the best value found was in fact the optimum value
- a factorization of H was used to confirm that the matrix was positive semi-definite,
while the algorithm ensured that the remaining optimality conditions hold - although, of course,
there is no guarantee that this will always be the case. We measured the iteration (ST) and the
percentage (ratio) of the optimal value obtained at the point at which the Steihaug-Toint method
left the trust region, as well as the number of iterations taken to achieve 10%, 90% and 99% of
the optimal reduction (10%, 90%, 99% respectively).
The results of these experiments are summarized in Table 6.1. In this table we give the
name of each example used, along with its dimension n, and the statistics "ratio"(expressed in
the form x(y) as a shorthand for x \Theta 10 y ), "ST", "10%", "90%" and "99%" as just described.
Some of the problems had interior solutions, in which case the "ratio" and "ST" statistics are
absent (as indicated by a dash). We considered both the unpreconditioned method
and a variety of standard preconditioners - a band preconditioner with semi-bandwidth of 5,
and modified incomplete and sparse Cholesky factorizations, with the modifications as proposed
by Schnabel and Eskow (1991) - used by the LANCELOT package (see, Conn, Gould and Toint,
1992, Chapter 3). The Cholesky factorization methods both failed for the problem MSQRTALS for
which the Hessian matrix required too much storage.
We make a number of observations.
1. On some problems, the Steihaug-Toint point gives a model value which is a good approximation
to the optimal value.
Solving the trust-region subproblem using the Lanczos method 19
no preconditioner 5-band
example
BRYBND 1000 3(-5) 23 24 28
CHAINWOO 1000 4(-5) 15
COSINE 1000
GENROSE 1000
MSQRTALS 1024 1(-1) 12 11 23
Incomplete Cholesky Modified Cholesky
example
COSINE 1000
GENROSE 1000
MSQRTALS 1024 factorization failure factorization failure

Table

6.1: A comparison of the number of iterations required to achieve a given percentage of
the optimal model value for a variety of preconditioners. See the text for a key to the data.
N. I. M. Gould, S. Lucidi, M. Roma and Ph. L. Toint
2. On other problems, a few extra iterations beyond the Steihaug-Toint point pay handsome
dividends.
3. Getting to within 90% or even 99% of the best value very rarely requires many more
iterations than to find the Steihaug-Toint point.
In conclusion, based on these numbers, we suggest that a good strategy would be to perform
a few (say 5) iterations beyond the Steihaug-Toint point, and only accept the improved point
if its model value is significantly better (as this will cost a second pass to compute the Lanczos
vectors). We shall consider this further in the next section.
6.2 Do better values than Steihaug-Toint imply a better trust-region method?
We now consider how the methods we have described for approximately solving the trust-region
subproblem perform within a trust-region algorithm. Of particular interest is the question as to
whether solving the subproblem more accurately reduces the number of trust-region iterations,
or more particularly the cost of solving the problem - the number of iterations is of concern if
the evaluation of the objective function and its derivatives is the dominant cost as then there is
a direct correlation between the number of iterations and the overall cost of solving the problem.
In

Tables

6.2 and 6.3, we compare the Steihaug-Toint scheme with the GLTR algorithm
(Algorithm 5.1) run to high accuracy. We exclude the problem HYDC20LS for our reported results
as no method succeeded in solving the problem in fewer than our limit of n iterations, and
the problems BROYDN7D and SPMSRTLS as a number of different local minima were found. In
these tables, in addition to the name and dimension of each example, we give the number of
objective function ("#f") and derivative ("#g") values computed, the total number of matrix-vector
products ("#prod") required to solve the subproblems, and the total CPU time required
in seconds. We compare the same preconditioners M as we used in the previous section. We
indicate those cases where one or other method performs at least 10% better than its competitor
by highlighting the relevant figure in bold.
We observe the following.
1. The use of different M leads to radically different behaviour. Different preconditioners
appear to be particularly suited to different problems. Surprisingly, perhaps, the unpreconditioned
algorithm often performs the best overall.
2. In the unpreconditioned case, the model-optimum variant frequently requires significantly
fewer function evaluations than the Steihaug-Toint method. However, the extra algebraic
costs per iteration often outweigh the reduction in the numbers of iterations. The advantage
in function calls for the other preconditioners is less pronounced.
Ideally, one would like to retain the advantage in numbers of function calls, while reducing
the cost per iteration. As we noted in Section 6.1, one normally gets a good approximation to
the optimal model value after a modest number of iterations. Moreover, while the Steihaug-Toint
point often gives a significantly sub-optimal value, a few extra iterations usually suffices to give
Solving the trust-region subproblem using the Lanczos method 21
no preconditioner Steihaug-Toint model optimum
example
iterations 865 577 34419 145.02
DQRTIC 1000 43 43 83 0.3 43 43 91 0.3
FREUROTH 1000 17 17 34 0.4 17 17 34 0.4
GENROSE 1000 859 777 6092 28.8 773 642 24466 82.2
MSQRTALS 1024 44 34 7795 486.0
NONCVXUN 1000 492 466 177942 1017.9 ? 1800 seconds
SENSORS 100 20 19 37 6.4 20 19 140 8.8
SINQUAD 5000 182 114 363 24.3 161 106 382 24.6
5-band Steihaug-Toint model optimum
example
CHAINWOO 1000 146 99 145 4.8 191 123 196 6.3
COSINE 1000 21 15 20 0.4 21 15
CRAGGLVY 1000 22 22 21 1.1 22 22 21 1.1
DQRTIC 1000 54 54 53 0.9 54 54 53 1.0
EIGENALS 930 56 43 171 75.2 53 42 222 75.8
FREUROTH 1000 20
iterations ? n iterations
MANCINO 100 91 72 90 87.2 52 43 90 52.2
MSQRTALS 1024 88 62 9793 700.2 73 52 19416 1292.2
22 N. I. M. Gould, S. Lucidi, M. Roma and Ph. L. Toint
Incomplete Cholesky Steihaug-Toint model optimum
example
CHAINWOO 1000 174 115 173 8.1 183 121 309 10.3
COSINE 1000 22 17 26 0.8 22 19 49 1.2
CRAGGLVY 1000 22 22 21 1.5 22 22 21 1.5
DQRTIC 1000 54 54 53 0.9 54 54 53 1.1
iterations ? n iterations
GENROSE 1000 948 629 951 35.5 496 322 847 23.5
28
MSQRTALS 1024 factorization failure factorization failure
28 150 41.2
iterations ? n iterations
iterations ? n iterations
SINQUAD 5000 77 52 89 542.6 78 50 121 526.7
Modified Cholesky Steihaug-Toint model optimum
example
BRYBND 1000 15 15 14 2.2 59 37 61 7.7
CHAINWOO 1000 178 119 177 7.6 183 121 309 10.3
COSINE 1000
CRAGGLVY 1000 23 23 33 1.4 22 22 21 1.6
DQRTIC 1000 54 54 53 1.2 54 54 53 1.1
iterations ? n iterations
GENROSE 1000 462 332 463 16.5 496 322 847 23.4
MANCINO 100 31 28 28
MSQRTALS 1024 factorization failure factorization failure
Solving the trust-region subproblem using the Lanczos method 23
a large percentage of the optimum. Thus, we next investigate both of these issues in the context
of an overall trust-region method.
In

Tables

6.4 and 6.5, we compare the number of function evaluations (#f), and the CPU time
taken to solve the problem for the Steihaug-Toint ("ST") method with a number of variations
on our basic GLTR method (Algorithm 5.1). The basic requirement is that we compute a model
value which is at least 90% of the best value found during the first pass of the GLTR method. If
this value is obtained by an iterate before that which gives the Steihaug-Toint point, the Steihaug-
Toint point is accepted. Otherwise, a second pass is performed to recover the first point at which
90% of the best value was observed. The other ingredient is the choice of the stopping rule for
the first pass. One possibility is to stop this pass as soon as the test (6.4) is satisfied. We denote
this strategy by "90%best". The other possibility is to stop when either (6.4) is satisfied or at
most a fixed number of iterations beyond the Steihaug-Toint point have occurred. We refer to
this as "90%(ST+k)", where k gives the number of additional iterations allowed. We investigate
the cases Once again, we compare the same preconditioners M as we used in
the previous section. We highlight in bold those entries which are at least 10% better than the
competition.
The conclusions are broadly as before. Each method has its successes and failures, and there
is no clear overall best method or preconditioner, although the unpreconditioned version performs
surprisingly well. Restricting the number of iteration allowed after the Steihaug-Toint point has
been found appears to curb the worst behaviour of the unrestricted method.
N. I. M. Gould, S. Lucidi, M. Roma and Ph. L. Toint
no preconditioner ST 90%(ST+1) 90%(ST+5) 90%(ST+10) 90%best
example
CRAGGLVY 1000 19 1.0 19 0.9 19 0.9 19 0.9 19 1.0
DQRTIC 1000 43 0.3 43 0.3 43 0.3 43 0.3 43 0.3
GENROSE 1000 859 28.8 748 38.9 721 48.1 738 57.3 728 60.0
MSQRTALS 1024 44 486.0
NONCVXUN 1000 492 1017.9 368 861.3 ? 1800 secs. ? 1800 secs. 433 1198.6
SENSORS 100 20 6.4 23 7.3 21 8.1 21 8.0 21 8.1
SINQUAD 5000 182 24.3 152 20.8 152 21.7 152 21.4 152 21.5
5-band ST 90%(ST+1) 90%(ST+5) 90%(ST+10) 90%best
example
CHAINWOO 1000 146 4.8 159 5.1 159 5.1 159 5.2 159 5.1
COSINE 1000 21 0.4 21 0.5 21 0.4 21 0.4 21 0.5
CRAGGLVY 1000 22 1.1 22 1.0 22 1.1 22 1.1 22 1.1
DQRTIC 1000 54 0.9 54 0.9 54 1.0 54 1.0 54 1.0
MANCINO 100 91 87.2 52 51.8 52 51.8 52 52.0 52 51.8
MSQRTALS 1024 88 700.2 97 756.7 73 704.9 74 844.7 79 981.5
Solving the trust-region subproblem using the Lanczos method 25
Incomplete Cholesky ST 90%(ST+1) 90%(ST+5) 90%(ST+10) 90%best
example
BRYBND 1000 55 3.9 56 4.2 56 4.3 56 4.3 56 5.0
CHAINWOO 1000 174 8.1 199 9.7 199 10.1 199 10.2 199 10.1
CRAGGLVY 1000 22 1.5 22 1.6 22 1.6 22 1.5 22 1.6
DQRTIC 1000 54 0.9 54 1.0 54 1.1 54 1.1 54 1.1
EIGENALS 930 76 94.6 77 97.2 74 97.2 74 97.3 74 96.8
GENROSE 1000 948 35.5 500 22.4 499 23.0 499 23.0 499 23.0
MSQRTALS 1024 fact. failure fact. failure fact. failure fact. failure fact. failure
SINQUAD 5000 77 542.6 68 484.2 68 484.1 68 485.4 68 489.0
Modified Cholesky ST 90%(ST+1) 90%(ST+5) 90%(ST+10) 90%best
example
BRYBND 1000 15 2.2 15 2.2 15 2.3 15 2.2 15 2.2
CHAINWOO 1000 178 7.6 176 7.9 176 7.9 176 7.8 176 8.0
COSINE 1000 41 1.1 41 1.3 41 1.3 41 1.3 41 1.3
DQRTIC 1000 54 1.2 54 1.2 54 1.3 54 1.3 54 1.3
GENROSE 1000 462 16.5 434 18.8 434 19.3 434 19.1 434 19.1
MANCINO 100 31 129.3 64 232.3 77 275.9 77 275.5 77 275.6
MSQRTALS 1024 fact. failure fact. failure fact. failure fact. failure fact. failure
26 N. I. M. Gould, S. Lucidi, M. Roma and Ph. L. Toint
7 Perspectives and conclusions
We have considered a number of methods which aim to find a better approximation to the solution
of the trust-region subproblem than that delivered by the Steihaug-Toint scheme. These methods
are based on solving the subproblem within a subspace defined by the Krylov space generated
by the conjugate-gradient and Lanczos methods. The Krylov subproblem has a number of useful
properties which lead to its efficient solution. The resulting algorithm is available as a Fortran
90 module, HSL VF05, within the Harwell Subroutine Library (1998).
We must admit to being slightly disappointed that the new method did not perform uniformly
better than the Steihaug-Toint scheme, and were genuinely surprised that a more accurate approximation
does not appear to significantly reduce the number of function evaluations within
a standard trust-region method. While this may limit the use of the methods developed here,
it also calls into question a number of other recent eigensolution-based proposals for solving the
trust-region subproblem (see Rendl, Vanderbei and Wolkowicz, 1995, Rendl and Wolkowicz, 1997,
Sorensen, 1997, Santos and Sorensen, 1995). While these authors demonstrate that their methods
provide an effective means of solving the subproblem, they make no effort to evaluate whether
this is actually useful within a trust-region method. The results given in this paper suggest that
this may not in fact be the case. This also leads to the interesting question as to whether it is
possible to obtain useful low-accuracy solutions with these methods.
We should not pretend that the formulae given in this paper are exact or even accurate
in floating-point arithmetic. Indeed, it is well-known that the floating-point matrices Q k from
the Lanczos method quickly loose M-orthonormality (see, for instance, Parlett, 1980, Section
13.3). Despite this, the method as given appears to be capable of producing usable approximate
solutions to the trust-region subproblem. We are currently investigating why this should be so.
One further possibility, which we have not considered so far, is to find an estimate - using
the first pass of Algorithm 5.1, and then to compute the required s by minimizing the unconstrained
model
using the preconditioned conjugate gradient method.
The advantage of doing this is that any instability in the first pass does not necessarily reappear
in this auxiliary calculation. The disadvantages are that it may require more work than simply
using (5.1), and that - must be computed sufficiently large to ensure that H + -M is positive
semi-definite.

Acknowledgement

We would like to thank John Reid for his helpful advice on computing eigenvalues of tridiagonal
matrices, and Jorge Mor'e for his useful comments on the Mor'e and Sorensen (1983) method.
We are grateful to the British Council-MURST for a travel grant (ROM/889/95/53) which made
some of this research possible.
Solving the trust-region subproblem using the Lanczos method 27



--R

CUTE: Constrained and unconstrained testing environment.
A family of trust-region-based algorithms for unconstrained minimization with strong global convergence properties
LANCELOT: a Fortran package for large-scale nonlinear optimization (Release
Two new unconstrained optimization algorithms which use function and gradient values.
Computing optimal locally constrained steps.
Maximization by quadratic hill-climbing

A catalogue of subroutines (release 13).
An algorithm for minimization using exact second derivatives.
Numerical experience with new truncated Newton methods in large scale unconstrained optimization.
Computing a trust region step.

The Symmetric Eigenvalue Problem.
Tracking the progress of the Lanczos algorithm for large symmetric eigenproblems.
Convergence properties of a class of minimization algorithms.
A semidefinite framework for trust region subproblems with applications to large scale minimization.

A new matrix-free algorithm for the large-scale trust-region subproblem
A new modified Cholesky factorization.
Newton's method with a model trust modification.
Minimization of a large-scale quadratic function subject to a spherical constraint
SIAM Journal on Optimization
The conjugate gradient method and trust regions in large scale optimization.
Towards an efficient sparsity exploiting Newton method for minimization.

--TR

--CTR
Nicholas I. M. Gould , Philippe L. Toint, FILTRANE, a Fortran 95 filter-trust-region package for solving nonlinear least-squares and nonlinear feasibility problems, ACM Transactions on Mathematical Software (TOMS), v.33 n.1, p.3-es, March 2007
Giovanni Fasano , Massimo Roma, Iterative computation of negative curvature directions in large scale optimization, Computational Optimization and Applications, v.38 n.1, p.81-104, September 2007
Nicholas I. M. Gould , Dominique Orban , Philippe L. Toint, GALAHAD, a library of thread-safe Fortran 90 packages for large-scale nonlinear optimization, ACM Transactions on Mathematical Software (TOMS), v.29 n.4, p.353-372, December
Nicholas I. M. Gould , Dominique Orban , Philippe L. Toint, CUTEr and SifDec: A constrained and unconstrained testing environment, revisited, ACM Transactions on Mathematical Software (TOMS), v.29 n.4, p.373-394, December
Nicholas I. M. Gould , Philippe L. Toint, An iterative working-set method for large-scale nonconvex quadratic programming, Applied Numerical Mathematics, v.43 n.1-2, p.109-128, October 2002
trust-region subproblem;conjugate gradients;preconditioning;lanczos method