Trust Region Algorithms and Timestep Selection.
Unconstrained optimization problems are closely related to systems of ordinary differential equations (ODEs) with gradient structure. In this work, we prove results that apply to both areas. We analyze the convergence properties of a trust region, or Levenberg--Marquardt, algorithm for optimization. The algorithm may also be regarded as a linearized implicit Euler method with adaptive timestep for gradient ODEs. From the optimization viewpoint, the algorithm is driven directly by the Levenberg--Marquardt parameter rather than the trust region radius. This approach is discussed, for example, in [R. Fletcher, Practical Methods of Optimization, 2nd ed., John Wiley, New York, 1987], but no convergence theory is developed. We give a rigorous error analysis for the algorithm, establishing global convergence and an unusual, extremely rapid, type of superlinear convergence. The precise form of superlinear convergence is exhibited---the ratio of successive displacements from the limit point is bounded above and below by geometrically decreasing sequences. We also show how an inexpensive change to the algorithm leads to quadratic convergence. From the ODE viewpoint, this work contributes to the theory of gradient stability by presenting an algorithm that reproduces the correct global dynamics and gives very rapid local convergence to a stable steady state.
Introduction
. This work involves ideas from two areas of numerical anal-
ysis: optimization and the numerical solution of ordinary differential equations
(odes). We begin by pointing out a connection between the underlying mathematical
problems.
Given a smooth function f R, an algorithm for unconstrained optimization
seeks to find a local minimizer ; that is, a point x ? such that f(x ?
for all x in some neighborhood of x ? . The following standard result gives necessary
conditions and sufficient conditions for x ? to be a local minimizer. Proofs may be
found, for example, in [5, 6, 7].
Theorem 1.1. The conditions rf(x ? positive semi-definite
are necessary for x ? to be a local minimizer, whilst the conditions rf(x ?
positive definite are sufficient.
On the other hand, given a smooth function F
may consider the ode system
Now suppose that F in (1.1) has the form F(x) j \Gammarf (x). By the Chain Rule, if
solves (1.1) then
d
dt
@f
dt
From (1.2) we see that along any solution of the ode the quantity f(x(t)) decreases
in Euclidean norm as t increases. Moreover, it strictly decreases unless
Hence, solving the ode up to a large value of t may be regarded as an attempt to
compute a local minimum of f . The conditions given in Theorem 1.1 may now be
interpreted as necessary conditions and sufficient conditions for x ? to be a linearly
stable fixed point of the ode.
If it is possible to write F(x) in the form \Gammarf (x) then the ode (1.1) is said
to have a gradient structure; see, for example, [19]. Several authors have noted the
Department of Mathematics, University of Strathclyde, Glasgow, G1 1XH, UK. Supported by
the Engineering and Physical Sciences Research Council of the UK under grant GR/K80228. This
manuscript appears as University of Strathclyde Mathematics Research Report 3 (1998).
connection between optimization and gradient odes. Schropp [18] examined fixed
timestep Runge-Kutta (rk) methods from a dynamical systems viewpoint, and
found conditions under which the numerical solution of the gradient ode converges
to a stationary point of f . Schropp also gave numerical evidence to suggest that
there are certain problem classes for which the ode formulation is preferable to
the optimization analogue. The book [11] shows that many problems expressible in
optimization terms can also be written as odes, often with gradient structure. Chu
has exploited this idea in order to obtain theoretical results and numerical methods
for particular problems; see [3] for a review. In the optimization literature, the
gradient ode connection has also been mentioned; see, for example, the discussion
on unconstrained optimization in [17]. Related work [1, 2] has looked at the use of
ode methods to solve systems of nonlinear algebraic equations.
The study of numerical methods applied to odes in gradient form has lead to
the concept of gradient stability [14, 20, 21]. The gradient structure arises in many
application areas, and provides a very useful framework for analysis of ode algo-
rithms. (In contrast with the classical linear and strictly contractive test problems,
gradient systems allow multiple equilibria.)
In [14, 20] positive results were proved about the ability of rk methods to preserve
the gradient structure, and hence to capture the correct long term dynamics,
for small, fixed timesteps. Most of these results require an extra assumption on
F that imposes either a one-sided Lipschitz condition or a form of dissipativity.
Adaptive rk methods, that is, methods that vary the timestep dynamically, were
analyzed in [21]. Here the authors considered a very special class of rk formula
pairs and showed that a traditional error control approach forces good behavior for
sufficiently small values of the error tolerance, independently of the initial data.
This would be regarded as a global convergence proof in the optimization litera-
ture. These results require a one-sided Lipschitz condition on F. A similar result
was proved in [12] for general ode methods that successfully control the local error-
per-unit-step. In this case the error tolerance must be chosen in a way that depends
on the initial data.
The work presented here has two main contributions.
ffl First, we note a close similarity between a trust region, or Levenberg-
Marquadt, algorithm for optimization and an adaptive, linearized, implicit
Euler method for gradient odes. We analyze the optimization algorithm
and establish a new result about its convergence properties. This also adds
to the theory of gradient stability for odes. Under a mild assumption on
f we show that the method is globally convergent and enjoys a very rapid
form of superlinear convergence. (The notion of the rate of convergence
to equilibrium is widely studied in optimization, but appears not to have
been considered in the gradient ode context. It is easily seen that any fixed
timestep rk formula that approaches equilibrium will do so at a generically
linear rate, in terms of the timestep number.)
ffl Second, we use the ideas from the gradient analysis to construct a timestepping
method for general odes that gives rapid superlinear local convergence
to a stable fixed point.
The presentation is organized as follows. In the next section we introduce New-
ton's Method and some simple numerical ode methods. Section 3 is concerned with
a specific trust region algorithm for unconstrained optimization. The algorithm,
which is essentially the same as one found in [6], is defined in x3.1. A non-rigorous
discussion of the convergence properties is given in x3.2, and the main convergence
theorems are proved in x3.3. The algorithm may also be regarded as a timestepping
process for a gradient ode algorithm and the analogous results are stated in x4. In
x5 we develop a timestepping scheme for general odes that gives superlinear local
convergence to stable fixed points.
2. Numerical Methods. Most numerical methods for finding a local minimizer
of f begin with an initial guess x 0 and generate a sequence fx k g. Similarly,
one-step methods for the ode (1.1) produce a sequence fx k g with x k - x(t k ).
The time-levels ft k g are determined dynamically by means of the timestep \Deltat k :=
The Steepest Descent method for optimization has the form
where ff k is a scalar that may arise, for example, from a line search. This is equivalent
to the explicit Euler Method applied to the corresponding gradient ode with
timestep \Deltat k j ff k . We note in passing that the poor performance of Steepest
Descent in the presence of steep-sided narrow valleys is analogous to the poor performance
of Euler's Method on stiff problems. Indeed, Figure 4j in [7] and Figure 1.2
in [10] illustrate essentially the same behavior, viewed from these two different perspectives

Newton's Method for optimization is based on the local quadratic model
Note that q k (ffi) is the quadratic approximation to f(x k + ffi) that arises from a
Taylor series expansion about x k . If r 2 f(x k ) is positive definite then q k (ffi) has
the unique minimizer
We thus arrive at Newton's
Method
The following result concerning the local quadratic convergence of Newton's
Method may be found, for example, in [5, 6, 7].
Theorem 2.1. Suppose that f 2 C 2 and that r 2 f satisfies a Lipschitz condition
in a neighborhood of a local minimizer x ? . If x 0 is sufficiently close to x ? and
positive definite, then Newton's method is well defined for all k and
converges at second order.
The Implicit Euler Method applied to (1.1) with F(x) j \Gammarf (x) using a
timestep of \Deltat k produces the equation
This is generally a nonlinear equation that must be solved for x k+1 . Applying
one interation of Newton's Method (that is, Newton's Method for solving nonlinear
equations) with initial guess x
This method is sometimes referred to as the Linearized Implicit Euler Method; see,
for example, [22]. Note that for large values of \Deltat k we have
and the ode method looks like Newton's Method (2.3). On the other hand, for
small \Deltat k we have
which corresponds to a small step in the direction of steepest descent (2.1). Hence,
at the extremes of large and small \Deltat k , the ode method behaves like well-known
optimization methods. However, we can show much more: for any value of \Deltat k ,
the method (2.5) can be identified with a trust region process in optimization.
This connection was pointed out by Goldfarb in the discussion on unconstrained
optimization in [17]. The relevant optimization theory is developed in the next
section.
3. A Trust Region Algorithm.
3.1. The Algorithm. We have seen that Newton's Method is based on the
idea of minimizing the local quadratic model q k (ffi) in (2.2) on each step. Since the
model is only valid locally, it makes sense to restrict the increment; that is, to seek
an increment ffi that minimizes q k (ffi) subject to some constraint kffik - h k . Here h k
is a parameter that reflects how much trust we are prepared to place in the model.
Throughout this work we use k \Delta k to denote the Euclidean vector norm and
the corresponding induced matrix norm. In this case a solution to the locally-
constrained quadratic model problem can be characterized. The following Lemma is
one half of [6, Theorem 5.2.1]; a weaker version was proved in [8]. For completeness,
we give a proof here.
Lemma 3.1. Given G 2 R m\Thetam and g 2 R m , if, for some - 0,
\Gammag
and G+ -I is positive semi-definite, then b ffi is a solution of
min
subject to kffik - k b ffik:
Furthermore, if G+ -I is positive definite, then b ffi is the unique solution of (3.2).
Proof. In the case where G + -I is positive semi-definite, it is straightforward
to show that b ffi minimizes
Hence, for all ffi we have b
solves the problem (3.2). When G + -I is positive definite, the inequality is
strict for ffi 6= b ffi, and hence the solution is unique.
Note that Lemma 3.1 does not show how to compute an increment b ffi given
a trust region constraint kffik - h k . Such an increment may be computed or approximated
using an iterative technique; see, for example, [6, pages 103-107] or [5,
pages 131-143]. However, as mentioned in [6], it is reasonable to regard - in (3.1)
as a parameter that drives the algorithm-having chosen a value for - and checked
that G+-I is positive definite, we may solve the linear system (3.1) and a posteriori
obtain a trust region radius h k := k b ffik. It easily shown that if G + -I is positive
definite then increasing - in (3.1) decreases k b ffik.
These remarks motivate Algorithm 3.2 below. We use - min (M) to denote the
smallest eigenvalue of a symmetric matrix M and let ffl ? 0 be a small constant.
Given x 0 and - 0 ? 0 a general step of the trust region algorithm proceeds as follows.
Algorithm 3.2.
Compute
Solve
Compute
Compute
Compute
using (3.3)
else
set r
If r k - 0
set x
else
set x
The algorithm involves the function
Note that r k records the ratio of the reduction in f from x k to x and the
reduction that is predicted by the local quadratic model. If r k is significantly less
than 1 then the model has been over-optimistic. This information is used in (3.3) to
update the trust region parameter -. In the case where the local quadratic model
has performed poorly, we double the - parameter, which corresponds to reducing
the trust region radius on the next step. If the performance is reasonable, we retain
the same value for -. In the case of good performance we halve the value of -,
thereby indirectly increasing the trust region radius.
We emphasize that Algorithm 3.2 is a trust region algorithm in the sense that
on each step ffi k solves the local restricted problem
min
subject to kffik - kffi k k:
Also, we remark that the algorithm is essentially the same as that described in [6,
pages 102-103]. The underlying idea of adding a multiple of the identity matrix to
ensure positive definiteness was first applied to the case where f has sum-of-squares
form, leading to the Levenberg-Marquadt algorithm. Goldfeld et al. [8] extended
the approach to a general objective function, and gave some theoretical justification.
Theorems 5.1.1 and 5.1.2 of [6] provide a general convergence theory for a wide
class of trust region methods. However, these results do not apply immediately to
Algorithm 3.2, since the algorithm does not directly control the radius h k := kffi k k,
but, rather, controls it indirectly via adaption of - k . In fact, we will see that the
behavior established in Theorem 5.1.2 of [6], local quadratic convergence, does not
hold for Algorithm 3.2. We are not aware of any existing convergence analysis that
applies directly to Algorithm 3.2, except for general results of the form encapsulated
in the Dennis-Mor'e Characterization Theorem for superlinear convergence [4, 5, 6]
and the "strongly consistent approximation to the Hessian" theory given in [16].
These references are discussed further in the remarks that follow Theorem 3.4.
3.2. Motivation for the Convergence Analysis. The proofs in x3.3 and
the appendix are rather technical, and hence, to help orient the reader, we give
below a heuristic discussion of the key points.
Theorem 3.3 establishes global convergence, and the proof uses arguments that
are standard in the optimization literature. Essentially, global convergence follows
from the fact that when the local quadratic model is inaccurate the algorithm
chooses a direction that is close to that of steepest descent. Perhaps of more interest
is the rate of local convergence. Suppose that x k
positive definite, and suppose that for k - b k we have r k ? 3=4, and
hence - It follows that, for some constant C 1 ,
Note also that G k and G \Gamma1
are bounded for large k.
given a large k, let ffi
Newt
k denote the correction that would arise from
Newton's method applied at x k , so that we have
Newt
Expanding (3.6), using (3.7),
Newt
Letting d k := x
Hence, in (3.8)
Newt
Using (3.5) we find that
Newt
for some constant C 2 .
Now, since x
Newt
k is the Newton step from x k , we have, from Theorem 2.1,
Newt
for some constant C 3 . The triangle inequality gives
Newt
Newt
and inserting (3.9) and (3.10) we arrive at the key inequality
for some constant C 4 . The first term on the right-hand side of (3.11) distinguishes
the algorithm from Newton's Method, and dominates the rate of convergence. To
proceed, it is convenient to consider a shifted sequence; let b e k := e k+N , for some
fixed N to be determined. Then from (3.11),
be k
Choosing N so that 2 N ? C 4 , we have
Now, neglecting the O(be 2
leads to
If, in addition to ignoring the O(be 2
in (3.13) we also assume that equality
holds, then we get equality in (3.14) and
be 0
but
be 0
We see from (3.15) that the error sequence is not quadratically convergent. However,
(3.16) corresponds to a very rapid form of superlinear convergence. Although this
analysis used several simplifying assumptions, the main conclusions can be made
rigorous, as we show in the next subsection. The type of superlinear convergence
that we establish is likely to be as good as quadratic convergence in practice. This
matter is discussed further after the proof of Theorem 3.4.
3.3. Convergence Analysis of the Trust Region Algorithm. The following
theorem shows that Algorithm 3.2 satisfies a global convergence result. The
structure of the proof is similar to that of [6, Theorem 5.1.1].
Theorem 3.3. Suppose that Algorithm 3.2 produces an infinite sequence such
that x k 2 B ae R m and g k 6= 0 for all k, where B is bounded and f 2 C 2 on B.
Then there is an accumulation point x 1 that satisfies the necessary conditions for
a local minimizer in Theorem 1.1.
Proof. Any sequence in B must have a convergent subsequence. Hence, we have
collects the indices in the convergent subsequence. It
is convenient to distinguish between two cases:
(i) sup
Case (i): From the form of V (r; -) in (3.3), there must be an infinite subsequence
whose indices form a set b
4 . Also, using the
boundedness of G k and g k , we have
and hence
Suppose that the gradient limit there exists a
descent direction s, normalized so that ksk = 1, such that
since ffi k solves the local restricted subproblem (3.4) we have q k (kffi k ks) -
ks T
Also, a Taylor expansion of f(x
We conclude from (3.18), (3.20) and (3.21) that r
S , which
contradicts r k ! 1. Hence,
Now suppose that G 1 := G(x 1 ) is not positive semi-definite; so there is a
direction v, with
Pick
S, and
S with k - b k. Then, since
solves the local restricted subproblem (3.4), we have
and hence
It follows from (3.18), (3.21) and (3.23) that r
S , which
contradicts
4 . Hence, G 1 is positive semi-definite.
Case (ii): From the form of V (r; -) in (3.3), there must be an infinite subsequence
whose indices form a set -
If
and hence
where Gmax := sup x2B
This gives
Hence, removing the earlier indices from -
S if necessary, we have, with h k := kffi k k,
min
we have \Deltaf
S. From r k - 1
4 it follows that \Deltaq k ! 0. Let
kffik - h and set -
x
Hence, is feasible on the subproblem that is solved by ffi k , and so
Letting
S , it follows from (3.25) that q k ( - ffi) - f
also minimizes q 1 (ffi) on kffik - h, and since the constraint is inactive, the necessary
conditions of Theorem 1.1 must be satisfied. Hence, g 1 6= 0 is contradicted.
Now, with case (ii), we have
as
S . Suppose G 1 is not positive semi-definite. Then the arguments
giving (3.22)-(3.23) may be applied, and we conclude that r
in -
S . It then follows from (3.3) that - k ! 0, and since - min must
have G 1 positive semi-definite. This gives the required contradiction.
Note that, as mentioned in [6], since the algorithm computes a non-increasing
sequence f k , the bounded region B required in this theorem will exist if any level
set is bounded.
In Theorem 3.3 we assume that g k 6= 0 for all k. If g b
the algorithm essentially terminates, giving x
However, in this case we cannot conclude that r 2 f(x k ) is positive semi-definite for
The next theorem quantifies the local convergence rate of Algorithm 3.2. The
first part of the proof is based on that of [6, Theorem 5.1.2].
Theorem 3.4. If the accumulation point x 1 of Theorem 3.3 also satisfies
the sufficient conditions for a local minimizer in Theorem 1.1, then for the main
sequence 1. Further, the displacement error e k :=
for some constant C, and if e k ? 0 for all k,
e
for constants e
C, but the ratio e k+1 =e 2
k is unbounded.
Proof. First, we show that case (i) of (3.17) in the proof of Theorem 3.3 can
be ruled out. Suppose that case (i) arises. Then r
k !1 in b
S.
positive definite, the matrix G k is also positive definite for large
k in b
S . In this case the Newton correction, ffi
Newt
Newt
\Gammag k , is
well defined and gives a global minimum of the local quadratic model q k . Define ff
by ffkffi Newt
and note that since ffi k solves the local restricted subproblem
(3.4), we have ff - 1. Then
Newt
Newt
Newt
Newt
Newt
Hence, using f
Newt
Newt
where - min ? 0 is a lower bound for the smallest eigenvalue of G k for large k in b
S .
It follows that
We may now conclude from (3.21) that r k ! 1 as
S . Hence, case (i)
cannot arise.
For case (ii), we have
as k !1 with k 2 -
S . Further, since
is a lower bound for the smallest eigenvalue of G k for large k in b
S .
It follows from (3.21) that as k !1 in -
S we must have
Having established that - k ! 0, we now know that the correction used in the
algorithm looks like the Newton correction ffi
Newt
k , which satisfies G k ffi
Newt
\Gammag k .
Let x Newt
Newt
k . Also, let d k := x
k !1 in -
S, and, by the triangle inequality,
The quadratic convergence property of Newton's Method given in Theorem 2.1
implies that for x k is sufficiently close to x 1
for some constant A 1 .
Expanding the other term in (3.31), we find
Newt
k ), we find that
Using (3.32) and (3.33) in (3.31) gives, for large k 2 -
e k+1 --
where A 2 is a constant.
Repeating the arguments that generated the inequalities (3.29) and (3.30), we
can show that there is a neighborhood N around x 1 such that if x
then r k - 3=4, so that - =2. Hence, from (3.34), there is some - k 2 -
S for
which x- k 2 N and the main sequnce lies in N for k -
k. So in the main sequence
we have x large k.
Hence, (3.34) may be extended to the bound
where A 3 and A 4 are constants. Lemma A.1 now gives (3.26).
To obtain a lower bound on e k+1 we use the triangle inequality in the form
From (3.32) and (3.33) we have
A 5
for constants A 5 ? 0 and A 6 . Lemma A.1 gives the required result.
We now list a number of remarks about Theorem 3.4.
1. The theorem shows that Algorithm 3.2 does not achieve a quadratic local
convergence rate. This is caused by the fact that - k does not approach
zero quickly enough. We have which is reflected in the first
term on the right-hand side of (3.35). A straightforward adaptation of the
proof shows that by increasing the rate at which - k ! 0, it is possible to
make the second term on the right-hand side of (3.35) significant, so that
quadratic convergence is recovered. For example, this occurs if we alter the
strategy for changing - k so that -
(and - otherwise). However, as explained in item 4
below, we would not expect this change to improve performance in practice.
Quadratic convergence is also discussed in item 5 below.
2. The power k 2 =3 appearing in (3.26) and (3.27) has been chosen partly
on the basis of simplicity-it is clear from the proofs of Lemma A.1 and
Theorem 3.4 that it can be replaced by ak 2 , for any a ! 1=2. (This will, of
course, cause the constant C to change.)
3. It is also clear from the proof that the result is independent of the precise
numerical values appearing in the algorithm. The values 1=4 and 3=4 in
(3.3) can be replaced by any ff and fi, respectively, with
the factor 2 in (3.3) can be replaced by any factor greater than unity. If
the factor 1=2 in (3.3) is replaced by 1=K, for K ? 1, then the statement of
the theorem remains true with powers of 2 replaced by powers of K. (The
changes mentioned here will, of course, alter the constants C, e
C and b
C.)
4. Theorem 3.4 shows that e k+1 =e k ! 0, and hence the convergence rate is
superlinear. However, the geometrically decreasing upper and lower bounds
on e k+1 =e k in (3.28) give us much more information. Asymptotically, whilst
Newton's Method gives twice as many bits of accuracy per step, the bound
(3.28) corresponds to k more bits of accuracy on the kth step. In both cases,
the asymptotic regime where e k is small enough to make the convergence
rate observable, but not so small that rounding errors are significant, is
likely to consist of only a small number of steps.
5. Several authors have found conditions that are sufficient, or necessary and
sufficient, for superlinear convergence of algorithms for optimization or
rootfinding. The most comprehensive result of this form is the Dennis-Mor'e
Characterization Theorem [4], [5, Theorem 8.2.4] and [6, Theorem 6.2.3].
Also, section 11.2 of [16] analyzes a class of rootfinding algorithms that employ
"consistent approximations to the Hessian", and this approach may
be used to establish superlinear convergence of Algorithm 3.2. However,
these references, which cover general classes of algorithms, do not derive
sharp upper and lower bounds on the rate of superlinear convergence of the
type given in Theorem 3.4. In the terminology of [16, x11.2], Algorithm 3.2
uses a strongly consistent approximation to the Hessian and superlinear
convergence is implied by - k ! 0. It also follows from [16, Result 11.2.7]
that quadratic convergence arises if we ensure that - k - Ckg k k and convergence
at R-order at least (1
5)=2 occurs if - k - Ckx
some constant C.
4. Timestepping on Gradient Systems. If we identify the trust region
parameter - k with the inverse of the timestep \Deltat k , then the Linearized Implicit
Euler Method (2.5) is identical to the updating formula in Algorithm 3.2. Hence,
Algorithm 3.2 can be regarded as an adaptive Linearized Implicit Euler Method for
gradient odes, and the convergence analysis of x3 applies. For completeness, we
re-write Algorithm 3.2 as a timestepping algorithm.
Given \Deltat 0 ? 0 and x 0 (= x init ), a general step of the algorithm for the gradient
system (1.1) with F(x) j \Gammarf (x) proceeds as follows.
Algorithm 4.1.
Compute
Solve
Compute
Compute
Compute
using (4.1)
else
set r
If r k - 0
set x
else
set x
The appropriate analogue of (3.3) is the function
\Deltat; 1- r - 3;
The following result is a re-statement of Theorems 3.3 and 3.4 in this context.
Theorem 4.2. Suppose that Algorithm 4.1 for (1.1) with F(x) j \Gammarf (x)
produces an infinite sequence such that x
is bounded and f 2 C 2 on B. Then there is an accumulation point x 1 that satisfies
the necessary conditions for a local minimizer in Theorem 1.1.
If the accumulation point x 1 also satisfies the sufficient conditions for a local
minimizer in Theorem 1.1, then for the main sequence
1. Further, the displacement error e k :=
for some constant C, and if e k ? 0 for all k,
e
for constants e
C, but the ratio e k+1 =e 2
k is unbounded.
In addition to the remarks at the end of x3, the following points should be
noted.
1. Algorithm 4.1 requires a check on the positive definiteness of the symmetric
. This is an unusual requirement for a timestepping
algorithm; however, we point out that an inexpensive and numerically stable
test can be performed in the course of a Cholesky factorization [13,
page 225]. If the test - min is omitted from Algorithm 4.1
then the local convergence rate is unaffected, but the global convergence
proof breaks down.
2. The rule for changing timestep is different in spirit to the usual local error
control philosophy for odes [9, 10]. This is to be expected, since the aim
of reaching equilibrium as quickly as possible is at odds with the aim of
following a particular trajectory accurately in time. The timestep control
policy in Algorithm 4.1 is based on a measurement of closeness to linearity
of the ode, and this idea is generalized in the next section. We also
note that local error control algorithms typically involve a user-supplied
tolerance parameter, with the understanding that a smaller choice of tolerance
produces a more accurate solution. Algorithm 4.1 on the other hand,
involves fixed parameters.
3. In [12] it is shown that, under certain assumptions, the use of local error
control on gradient odes forces the numerical solution close to equilibrium.
(Typically the solution remains within O(-) of an equilibrium point, where
- is the tolerance parameter.) This suggests that local error control may
form an alternative to the positive definiteness test as a means of ensuring
global convergence. Having driven the solution close to equilibrium with
local error control, the closeness to linearity test could be used to give
superlinear convergence.
5. Timestepping to a General, Stable Steady State. Motivated by x4,
we now develop an algorithm that gives rapid local convergence to stable equilibrium
for a general ode. We let F 0 denote the Jacobian of F, and define F k := F(x k ) and
k is not symmetric in general. The ratio
l k :=
indicates how close F is to behaving linearly in a region containing x k and x k+1 .
Given \Deltat 0 ? 0 and x 0 (= x init ), a general step of the algorithm for (1.1)
proceeds as follows.
Algorithm 5.1.
Compute
Solve
Compute using (5.1)
\Deltat
else
arbitrary
Since we are concerned only with local convergence properties, the action taken
when does not affect the analysis.
In the theorem below, B(fl; z) denotes the open ball of radius fl ? 0 about
z flg.
Theorem 5.2. Suppose that F(x ? in a neighborhood
of x ? and F 0 strictly negative real parts. Then given any
Algorithm 5.1 there is a fl ? 0 such that for any x 0 2
Further, the displacement error e k :=
for some constant C, and if e k ? 0 for all k,
e
for constants e
C, but the ratio e k+1 =e 2
k is unbounded.
Proof. There exists a b fl ? 0 such that F 0 (x) is nonsingular for x 2
Letting D 1 be an upper bound for kF
Now, from (5.1),
l k :=
and it follows from (5.5) that by reducing b fl if necessary we have jl
and F
k ), we have, for small e k ,
Now there exists a \Deltat ? ? 0 such that
\Deltat
By continuity, and by reducing b fl if necessary, we have
\Deltat
Hence, for large \Deltat we have a contraction in (5.6). We now show that for x 0
sufficiently close to x ? , \Deltat k increases beyond \Deltat ? while x k remains in B(bfl; x ? ).
Let
\Deltat-\Deltat
\Deltat
From (5.6), for x k 2 B(bfl; x ? ) we have
Let b k be such that 2 b k \Deltat 0 - \Deltat ? . From (5.8) we may choose sufficiently small
so that
Then, from (5.6) and (5.7), since \Deltat k - \Deltat ? for k - b k,
and hence e k ! 0 as k !1, and \Deltat \Deltat 0 for all k. We then have
F
\Deltat k
F
Since both F 0
k and F
are bounded, (5.6) gives
k and e k+1 - D 5
for constants completes the result.
It is straightforward to show that any fixed timestep rk or linear multistep
method can produce only a linear rate of convergence to equilibrium in general.
From Theorem 5.2, we see that Algorithm 5.1 provides a systematic means of increasing
the timestep in order to achieve a rapid form of superlinear convergence.
This has many applications, particularly in the area of computational fluid dynam-
ics, where it is common to solve a discretized steady partial differential equation
by introducing an artificial time derivative and driving the solution to equilibrium;
see, for example [22].
It is clear from the proof of Theorem 5.2 that for sufficiently large \Deltat 0 the
algorithm permits local convergence to an unstable fixed point. This can be regarded
as a consequence of the fact that the Implicit Euler Method is over-stable, in the
sense that the absolute stability region contains the infinite strip fz
in the right-half of the complex plane; see, for example, [15, page 229]. Another
explanation is that Newton's Method for optimizing f is identical to Newton's
Method for algebraic equations applied to see, for example, [5, page 100].
Hence, unless other measures are taken, there is no reason why stable fixed points
should be preferred. In Algorithm 4.1 for gradient odes we check that - min
which helps to force the numerical solution to a stable fixed
point. It is likely that traditional ode error control would also direct the solution
away from unstable fixed points, and hence the idea of combining optimization and
ode ideas forms an attractive area for future work.

Acknowledgements

This work has benefited from my conversations with a
number of optimizers and timesteppers; most notably Roger Fletcher and David
Griffiths.


Appendix

A. Convergence Rate Lemma.
Lemma A.1. Let
k for all k:
Then
for some constant C:
Further, if e k ? 0 for all k then
and if, in addition,
R
k for all k;
then
e
but the ratio e k+1 =e 2
k is unbounded.
Proof. Choose
We first prove a result under restricted circumstances, and then generalize to
the full result. We assume that
3:
Our induction hypothesis is
Note that, from (A.7), this holds for 3. If (A.8) is true for
using (A.1),
using (A.6) and (A.7),
Therefore, by induction, (A.8) is true for all k, if (A.7) holds.
Now, consider the shifted sequence b e k := e k+N , for some fixed N . We have
it is possible to choose N such that
3:
From (A.9) and (A.10), the result (A.8) holds for this shifted sequence, so
; for all k:
Translating this into a result for the original sequence, we find that,
Relabelling C as C=2 N 2 +N and letting b
Now
N . Hence,
Clearly, by increasing C, if necessary, the result will also hold for the finite sequence
N . Hence, (A.2) is proved. The inequality (A.3) follows after dividing
by e k in (A.1) and using (A.2).
From (A.2), for sufficiently large k we have 2 k T e k - R=2, so that
R=:
C:
Clearly, by reducing -
C, if necessary, this result must hold for all k. Now, reduce -
if necessary, so that
1. From (A.12),
letting e
e
e
Inequalities (A.12) and (A.13) give (A.5), as required.
Finally, using (A.2) and (A.5) we find that
e



--R

Fast local convergence with single and multistep methods for nonlinear equations.
The solution of nonlinear systems of equations by A-stable integration tech- niques
A list of matrix flows with applications.


Practical Methods of Optimization.
Practical Optimization.
Maximisation by quadratic hill-climbing
Solving Ordinary Differential Equations I
Solving Ordinary Differential Equations II
Optimization and Dynamical Systems.
Analysis of the dynamics of local error control via a piecewise continuous residual.
Accuracy and Stability of Numerical Algorithms.

Numerical Methods for Ordinary Differential Systems.
Iterative Solution of Nonlinear Equations in Several Variables.
Nonlinear Optimization
Using dynamical systems methods to solve minimization problems.
Nonlinear Dynamics and Chaos.
Model problems in numerical stability theory for initial value problems.
The essential stability of local error control for dynamical systems.
Global asymptotic behaviour of iterative implicit schemes.
--TR
steady state;global convergence;superlinear convergence;gradient system;levenberg-marquardt;unconstrained optimization;quadratic convergence