Accelerated Inexact Newton Schemes for Large Systems of Nonlinear Equations.
Classical iteration methods for linear systems, such as Jacobi iteration, can be accelerated considerably by Krylov subspace methods like GMRES @. In this paper, we describe how inexact Newton methods for nonlinear problems can be accelerated in a similar way and how this leads to a general framework that includes many well-known techniques for solving linear and nonlinear systems, as well as new ones.  Inexact Newton methods are frequently used in practice to avoid the expensive exact solution of the large linear system arising in the (possibly also inexact) linearization step of Newton's process.  Our framework includes acceleration techniques for the "linear steps" as well as for the "nonlinear steps" in Newton's process.  The described class of methods, the accelerated inexact Newton (AIN) methods, contains methods like GMRES and GMRESR for linear systems, Arnoldi and JacDav{} for linear eigenproblems, and many variants of Newton's method, like damped Newton, for general nonlinear problems.  As numerical experiments suggest, the AIN{} approach may be useful for the construction of efficient schemes for solving nonlinear problems.
Introduction
. Our goal in this paper is twofold. A number of iterative solvers
for linear systems of equations, such as FOM [23], GMRES [26], GCR [31], Flexible
GMRES [25], GMRESR [29] and GCRO [7], are in structure very similar to iterative
methods for linear eigenproblems, like shift and invert Arnoldi [1, 24], Davidson [6, 24],
and Jacobi-Davidson [28]. We will show that all these algorithms can be viewed as
instances of an Accelerated Inexact Newton (AIN) scheme (cf. Alg. 3), when applied
to either linear equations or to linear eigenproblems. This observation may help us
in the design and analysis of algorithms by "transporting" algorithmic approaches
from one application area to another. Moreover, our aim is to identify efficient AIN
schemes for nonlinear problems as well, and we will show how we can learn from the
algorithms for linear problems.
To be more specific, we will be interested in the numerical approximation of the
solution u of the nonlinear equation
(1)
where F is some smooth (nonlinear) map from a domain in R n (or C n ) that contains
the solution u, into R n (or C n ), where n is typically large.
Some special types of systems of equations will play an important motivating role
in this paper.
The first type is the linear system of equations
(2)
Mathematical Institute, Utrecht University, P. O. Box 80.010, NL-3508 Utrecht, The Nether-
lands. E-mail: fokkema@math.ruu.nl, sleijpen@math.ruu.nl, vorst@math.ruu.nl.
D. R. FOKKEMA, G. L. G. SLEIJPEN, AND H. A. VAN DER VORST
where A is a nonsingular matrix and b, x are vectors of appropriate size; A and b are
given, x is unknown. The dimension n of the problem is typically large and A is often
sparse. With is equivalent to (1). This type will serve
as the main source of inspiration for our ideas.
The second type concerns the generalized linear eigenproblem
With we have that, for normalized v, with F (u) := Av \Gamma -Bv, equation (3)
is equivalent to (1). This type is an example of a mildly nonlinear system and will
serve as an illustration for the similarity between various algorithms, seen as instances
of AIN (see Section 6).
However, the AIN schemes, that we will discuss, will be applicable to more general
nonlinear problems, like, for instance, equations that arise from discretizing nonlinear
partial differential equations of the form
where\Omega is a domain in R 2 or R 3 ,
and u satisfies suitable boundary conditions. An example of (4) is, for instance
where\Omega is some domain in R 2 and
@\Omega (see also Section 8).
Guided by the known approaches for the linear system (cf. [25, 29, 7]) and the
eigenproblem (cf. [28, 27]) we will define accelerated Inexact Newton schemes for more
general nonlinear systems. This leads to a combination of Krylov subspace methods
for Inexact Newton (cf. [16, 4] and also [8]) with acceleration techniques (as in [2]) and
offers us an overwhelming choice of techniques for further improving the efficiency of
Newton type methods. As a side-effect this leads to a surprisingly simple framework
for the identification of many well-known methods for linear-, eigen-, and nonlinear
problems.
Our numerical experiments for nonlinear problems, like problem (5), serve as an
illustration for the usefulness of our approach.
The rest of this paper is organized as follows. In Section 2 we briefly review the
ideas behind the Inexact Newton method. In Section 3 we introduce the Accelerated
Inexact Newton methods (AIN). We will examine how iterative methods for linear
problems are accelerated and we will distinguish between a Galerkin approach and a
Minimal Residual approach. These concepts are then extended to the nonlinear case.
In Section 4 we make some comments on the implementation of AIN schemes. In
Section 5 we show how many well-known iterative methods for linear problems fit in
the AIN framework. In Section 6 and Section 7 we consider instances of AIN for the
mildly nonlinear generalized eigenproblem and for more general nonlinear problems.
In Section 8 we present our numerical results and some concluding remarks are in
Section 9.
ACCELERATED INEXACT NEWTON SCHEMES 3
1. choose an initial approximation u 0 .
2. Repeat until u k is accurate
(b) Compute the residual r
(c) Compute an approximation J k for the Jacobian F 0 (u k ).
(d) Solve the correction equation (approximately). Compute
an (approximate) solution p k of the correction equation
Update. Compute the new approximation:
Algorithm 1: Inexact Newton
2. Inexact Newton methods. Newton type methods are very popular for solving
systems of nonlinear equations as, for instance, represented by (4). If u k is the
approximate solution at iteration number k, Newton's method requires for the next
approximate solution of (1), the evaluation of the Jacobian J k := F 0 and the
solution of the correction equation
Unfortunately, it may be very expensive, or even practically impossible to determine
the Jacobian and/or to solve the correction equation exactly, especially for larger
systems.
In such situations one aims for an approximate solution of the correction equa-
tion, possibly with an approximation for the Jacobian (see, e.g., [8]). Alg. 1 is an
algorithmical representation of the resulting Inexact Newton scheme.
If, for instance, Krylov subspace methods are used for the approximate solution of
the correction equation (6), then only directional derivatives are required (cf. [8, 4]);
there is no need to evaluate the Jacobian explicitly. If v is a given vector then the
vector F 0 (u)v can be approximated using the fact that
The combination of a Krylov subspace method with directional derivatives combines
the steps 2c and 2d in Alg. 1.
For an initial guess u 0 sufficiently close to a solution, Newton's method has asymptotically
at least quadratic convergence behavior. However, this quadratic convergence
is usually lost if one uses inexact variants and often the convergence is not much faster
than linear. In the next section we make suggestions how this (linear) speed of convergence
may be improved.
Note: It is our aim to restore, as much as possible, the asymptotic convergence
behavior of exact Newton and we do not address the question of global convergence.
4 D. R. FOKKEMA, G. L. G. SLEIJPEN, AND H. A. VAN DER VORST
1. choose an initial approximation x 0 .
2. Repeat until x k is accurate
(c)
(d) x
Algorithm 2: Jacobi Iteration
3. Accelerating Inexact Newton methods. Newton's method is a one step
method, that is, in each step, Newton's method updates the approximate solution
with information from the previous step only. However, in the computational process,
subspaces have been built gradually, that contain useful information concerning the
problem. This information may be exploited to improve the current approximate
solution, and this is what we propose to do. More precisely, we will consider alternative
update strategies for step 2e of the Inexact Newton algorithm Alg. 1.
3.1. Acceleration in the linear case. The linear system can be written
as F (x) := and \GammaA is the Jacobian of F . When the approximate solution p k
is computed as preconditioning matrix (approximating
A), then the Inexact Newton algorithm, Alg. 1, reduces to a standard Richardson-
type iteration process for the splitting R. For instance, the choice
leads to Jacobi iteration (see Alg. 2).
One may improve the convergence behavior of standard iterations schemes by
ffl using more sophisticated preconditioners M , and/or
ffl applying acceleration techniques in the update step.
Different preconditioners and different acceleration techniques lead to different algo-
rithms, some of which are well-known.
Examples of iterations schemes that use more sophisticated preconditioners are,
for instance, Gauss-Seidel iteration, where with L is the strict lower
triangular part of A and is a
relaxation parameter.
Examples of iterations schemes that use acceleration techniques are algorithms
that take their updates to the approximate solution as a linear combination of previous
directions p j . Preferable updates ~
are those for which b \Gamma Ax k+1 ,
where x
k , is minimal in some sense: e.g., kb \Gamma Ax k+1 k 2 is minimal, as in
GMRES [26] and GCR [31], or b \Gamma Ax k+1 is orthogonal to the p j for j - k, as in FOM
or GENCG [23], or b \Gamma Ax k+1 is "quasi-minimal", as in Bi-CG [17], and QMR [11].
Of course the distinction between preconditioning and acceleration is not a clear
one. Acceleration techniques with a limited number of steps can be seen as a kind
of dynamic preconditioning as opposed to the static preconditioning with fixed M .
In this view one is again free to choose an acceleration technique. Examples of such
iteration schemes are Flexible GMRES [25], GMRESR [29] and GCRO [7].
All these accelerated iteration schemes for linear problems construct approximations
the solution of a smaller or an
ACCELERATED INEXACT NEWTON SCHEMES 5
easier projected problem. For example, GMRES computes y k such that
equivalently (AV k )   (b \Gamma
such that V
as the
solution of a larger tri-diagonal problem obtained with oblique projections.
For stability (and efficiency) reasons one usually constructs another basis for the
span of V k with certain orthogonality properties, depending on the selected approach.
3.2. Acceleration in the nonlinear case. We are interested in methods for
finding a zero of a general nonlinear mapping F . For the linear case, the methods
mentioned above are, apart from the computation of the residual, essentially a mix of
two components: (1) the computation of a new search direction (which involves the
residual), and (2) the update of the approximation (which involves the current search
direction and possibly previous search directions, and the solution y k of a projected
problem). The first component may be interpreted as preconditioning, while the
second component is the acceleration.
Looking more carefully on how y k in the linear case is computed, we can distinguish
between two approaches based on two different conditions. With G k (y) := F
the FOM and the other "oblique" approaches lead to methods that compute y such
that (for appropriate W k
the GMRES
approach leads to methods that compute y such that kG k (y)k 2 is minimal (a Minimal
Residual condition).
From these observations for the linear case we now can formulate iteration schemes
for the nonlinear case.
The Inexact Newton iteration can be accelerated in a similar way as the standard
linear iteration. This acceleration can be accomplished by updating the solution by a
correction ~
k in the subspace spanned by all correction directions p j (j - k).
To be more precise, the update ~
k for the approximate solution is given by ~
y, where of the search space V k spanned by
Furthermore, with G k we propose to determine y by
ffl a Galerkin condition on G k (y): y is a solution of
where W k is some matrix of the same dimensions as V k ,
ffl or a Minimal Residual (MR) condition on G k (y): y is a solution of
min y
ffl or a mix of both, a Restricted Minimal Residual (RMR) condition on G k (y): y
is a solution of
min y
Equation (7) generalizes the FOM approach, while equation (8) generalizes the GMRES
approach.
Solving (7) means that the component of the residual r k+1 in the subspace W k
(spanned by W k ) vanishes. For W k one may choose, for instance, W (as in
6 D. R. FOKKEMA, G. L. G. SLEIJPEN, AND H. A. VAN DER VORST
FOM), or W is the component of J k p k orthogonal to W
(as in GMRES: w
linear equations the Minimal Residual
and the Galerkin approach coincide for the last choice.
As is known from the linear case, a complication of the Galerkin approach is
that equation (7) may have no solution, which means that this approach may lead to
breakdown of the method. In order to circumvent this shortcoming to some extend
we have formulated the Restricted Minimal Residual approach (9). Compared to (7)
this formulation is also attractive for another reason: one can apply standard Gauss-Newton
[9] schemes for solving general nonlinear least squares problems to it. One
might argue that a drawback of a Gauss-Newton scheme is that it may converge slowly
(or not at all). However, for least squares problems with zero residual solutions,
the asymptotic speed of convergence of a Gauss-Newton method is that of Newton's
method. This means, that if the Galerkin problem (7) has a solution, a Gauss-Newton
scheme applied to (9) will find it fast and efficient (see also Section 8).
Note that equations (7)-(9) represent nonlinear problems in only k variables,
which may be much easier to solve than the original problem. If these smaller non-linear
problems can be formulated cheaply, then the costs for an update step may be
considered as being relatively small.
Note also that since equations (7)-(9) are nonlinear, they may have more than one
solution. This fact may be exploited to steer the computational process to a specific
preferable solution of the original problem.
Accelerated Inexact Newton. For the Galerkin approach, step 2e in the Inexact
Newton's algorithm, Alg. 1, is replaced by four steps in which
ffl the search subspace V k\Gamma1 is expanded by an approximate "Newton correction"
and a suitable basis is constructed for this subspace,
ffl a shadow space W k is selected on which we project the original problem,
ffl the projected problem (7) is solved,
ffl and the solution is updated.
This is represented, by the steps 3e-3h in Alg. 3. The Minimal Residual approach
and the Restricted Minimal Residual approach can be represented in a similar way.
4. Computational considerations. In this section we make some comments on
implementation details that mainly focus on limiting computational work and memory
space.
4.1. Restart. For small k, problems (7)-(9) are of small dimension and may
often be solved at relatively low computational costs (e.g., by some variant of Newton's
method).
For larger k they may become a serious problem in itself. In such a situation,
one may wish to restrict the subspaces V and W to subspaces of smaller dimension
(see Alg. 3, step 3i). Such an approach limits the computational costs per iteration,
but it may also have a negative effect on the speed of convergence.
For example, the simplest choice, restricting the search subspace to a 1-dim. sub-space
leads to Damped Inexact Newton methods, where, for instance, the damping
parameter ff is the solution of min ff kG k (ff)k 2 , where G k
ACCELERATED INEXACT NEWTON SCHEMES 7
1. choose an initial approximation u 0 .
2.
3. Repeat until u k is accurate
(b) Compute the residual r
(c) Compute an approximation J k for the Jacobian F 0 (u k ).
(d) Solve the correction equation (approximately). Compute an
(approximate) solution p k for the correction equation
Expand the search space. Select a v k in the span(V that is
linearly independent of V k\Gamma1 and update
(f) Expand the shadow space. Select a w k that is linearly
independent of W k\Gamma1 and update W
(g) Solve the projected problem. Compute nontrivial solutions y of
the projected system
Update. Select a y k (from the set of solutions y) and update the
approximation:
(i) Restart. large, select an '
select ' \Theta ' 0 matrices R V and RW and compute
suitable combinations of the columns of
Algorithm 3: Accelerated Inexact Newton
Of course, a complete restart is also feasible, say after each mth step (cf. step 3i
of Alg. 3):
The disadvantage of a complete restart is that we have to rebuild subspace information
again. Usually it leads to a slower speed of convergence.
It seems like an open door to suggest that parts of subspaces better be retained at
restart, but in practical situations it is very difficult to predict what those parts should
be. A meaningful choice would depend on spectral properties of the Jacobian as well
as on the current approximation. When solving linear equations with GMRESR [29],
good results have been reported in [30] when selecting a number of the first and the
last columns (cf. step 3i of Alg. 3); e.g.,
8 D. R. FOKKEMA, G. L. G. SLEIJPEN, AND H. A. VAN DER VORST
In [7], a variant of GMRESR, called GCRO, is proposed, which implements another
choice. For subspaces of dimension l + m, the first l columns are retained, together
with a combination of the last m columns. This combination is taken such that
the approximate solution, induced by a minimal residual condition, is the same for both
the subspaces of dimension l +m and l + 1. To be more specific, if u
where y k solves min y replaced by [V k l ; V km y km (denoting
4.2. Update. In the update step (step 3h of Alg. 3), a solution y k of the projected
problem has to be selected from the set of solutions y. Selection may be necessary
since many nonlinear problems have more than one solution. Sometimes, this may
be the reason for poor convergence of (Inexact) Newton: the sequence of approximate
solutions "wavers" between different exact solutions. For larger search subspaces, the
search subspace may contain good approximations for more than one solution. This
may be exploited to steer the sequence of approximate solution to the wanted solution
may help to avoid wavering convergence behavior.
The selection of y k should be based on additional properties of the solution u k+1 .
For instance, we may look for the solution largest in norm, or as in the case of eigenvalue
problems, for a solution of which one component is close to some specific value
(for instance, if one is interested in eigenvalues close to, say, 0, the Ritz vector with
Ritz value closest to 0 will be chosen).
4.3. The projected problem. Even though problems of small dimension can
be solved with relatively low computational costs, step 3g in Alg. 3 is not necessarily
inexpensive. The projected problem is embedded in the large subspace and it may
require quite some computational effort to represent the problem in a small subspace
(to which y belongs) of dimension ' := dim(span(V k In the case of linear
equations (or linear eigenvalue problems) the computation of an ' \Theta '-matrix as W
products. For this type of problems, and for many others as well, one
may save on the computational costs by re-using information from previous iterations.
4.4. Expanding the search subspace. The AIN algorithm breaks down if
the search subspace is not expanded. This happens when p k belongs to the span of
(or, in finite precision arithmetic, if the angle between p k and this subspace is
very small). Similar as for GMRES, one may then replace p k by J k v ' , where v ' is the
last column vector of the matrix V k\Gamma1 .
With approximate solution of the correction equation, a breakdown will also occur
if the new residual r k is equal to the previous residual r k\Gamma1 . We will have such a
situation if y instead of modifying the expansion process in iteration
number k, one may also take measures in iteration number in order to avoid
In [29] a few steps by LSQR are suggested when the linear solver is a Krylov
subspace method:
may already cure the stagnation.
5. How linear solvers fit in the AIN framework. In this section we will
show how some well-known iterative methods for the solution of linear systems fit
in the AIN framework. The methods that follow from specific choices in AIN are
equivalent to well-known methods only in the sense that, at least in exact arithmetic,
ACCELERATED INEXACT NEWTON SCHEMES 9
they produce the same basis vectors for the search spaces, the same approximate
solutions, and the same Newton corrections (in the same sense as in which GMRES
and ORTHODIR are equivalent).
With the linear equation (2) is equivalent to the one in (1)
and J \GammaA. In this section, M denotes a preconditioning matrix for A (i.e., for a
vector v, M \Gamma1 v is easy to compute and approximates A \Gamma1 v).
5.1. GCR. With the choice,
Alg. 3 (without restart) is equivalent to preconditioned
GCR [31].
5.2. FOM and GMRES. The choice
algorithms that are related to FOM and GMRES [26]. With the
additional choice w Alg. 3 is just FOM, while the choice
gives an algorithm that is equivalent to GMRES.
5.3. GMRESR. Taking w k such that w k is perpendicular to
GCR and taking p k as an approximate solution of the equation
Alg. 3 is equivalent to the GMRESR algorithms [29]. One
might compute p k by a few steps of GMRES, for instance.
6. AIN schemes for mildly nonlinear problems. In this section we will discuss
numerical methods for the iterative solution of the generalized eigenproblem (3).
We will show that they also fit in the general framework of the AIN Alg. 3 methods.
As already mentioned these AIN methods consist of two parts. In one part
an approximate solution of the correction equation (cf. step 3d in Alg. 3) is used
to extend the search space. In the other part a solution of the projected problem
(cf. step 3g in Alg. 3) is used to construct an update for the approximate solution.
We will start with the derivation of a more suitable form for the (Newton) correction
equation for the generalized eigenproblem. After that, we will make some
comments on how to solve the projected problem.
The correction equation. In order to avoid some of the complications that go
with complex differentiation, we will mainly focus on the numerical computation of
eigenvectors with a fixed component in some given direction (rather then on the computation
of eigenvectors with a fixed norm).
First, let ~
u be a fixed vector with a nontrivial component in the direction of the
desired eigenvector x. We want to compute approximations u k for x with a normalized
component in the ~
u-direction: (u; ~
We will select # k such that the
residual r k := is orthogonal to w, where w is another fixed nontrivial
vector, i.e., the approximate eigenvalue # k is given by # k := w   Au k =w   Bu k .
Consider the map F given by
F (u) :=
and u belongs to the hyper-plane fy 2 C 1g. The Jacobian J
then given by
D. R. FOKKEMA, G. L. G. SLEIJPEN, AND H. A. VAN DER VORST
and the correction equation reads as
u; such that
equation is equivalent to
~ u   0
\Gammar k#
that is, p is the solution of (10) if and only if p is the solution of (11).
The projected problem. For the Generalized eigenvalue problem we are in the
fortunate position that all the solutions of problems of moderate size can be computed
by standard methods such as for instance the QZ [19] method. However, before we
can apply these methods we have to reformulate the projected problem because of the
exceptional position of u k in W   F
The key to this reformulation is the observation that in the methods we consider
the affine subspace u k +span(V k ) is equal to V k because V k contains u k by itself. Now,
as an alternative to step 3g in Alg. 3, we may also compute all the solutions y of
This problem can now be solved by for instance the QZ method, and after selecting
y k a new approximation u k+1 is given by u
6.1. Arnoldi's method. We consider the simplified case where I , i.e., the
standard eigenproblem. If we do only one step of a Krylov subspace method (Krylov
dimension 1) for the solution of the correction equation (10), then we obtain for the
correction
Hence, . Note that this may be a poor (very) approximation,
because, in general, r k 6? ~
u. The approximate eigenvector u k belongs to the search
subspace span(V expanding the search subspace by the component of p k orthogonal
to span(V k\Gamma1 ) is equivalent to expanding this space with the orthogonal component
of Au k , which would be the "expansion" vector in Arnoldi's method. Hence,
the search subspace is precisely the Krylov subspace generated by A and u 0 . Ap-
parently, Arnoldi's method is an AIN method (with a ``very inexact Newton step'')
without restart.
The choice W corresponds to the standard one in Arnoldi and produces
#'s that are called Ritz values, while the choice W leads to Harmonic Ritz
values [22].
6.2. Davidson's method. As in Arnoldi's method, Davidson's method [6] also
carries out only one step of a Krylov subspace method for the solution of the correction
equation. However, in contrast to Arnoldi's method, Davidson also incorporates a
preconditioner.
ACCELERATED INEXACT NEWTON SCHEMES 11
He suggests to solve (10) approximately by p k with
where M is the inverse of the diagonal of A \Gamma # k B. Other choices have been suggested
as well (cf. e.g., [5, 21]). Because of the preconditioner, even if I , the search
space is not simply the Krylov subspace generated by A and u 0 . This may lead to an
advantage of Davidson's method over Arnoldi's method.
For none of the choices of the preconditioner, proper care has been taken of the
projections (see (10)): the preconditioner should approximate the inverse of the projected
matrix (see (10)) as a map from ~
rather than of A \Gamma # k B.
However, if M is the diagonal of A \Gamma #B, and we choose ~ u and w equal to the
same, arbitrary standard basis vector (as Davidson does [6]) then
. Note that p ? because M is diagonal
and r k ? w. Therefore, for this particular choice of w (and ~
u), the diagonal M may
be expected to be a good preconditioner for the correction equation (including the
projections) in the cases where M is a good preconditioner for A \Gamma # k B. Observe that
this argument does not hold for non-diagonal preconditioners M .
6.3. Jacobi-Davidson. Davidson methods with a non-diagonal preconditioner
do not take care properly of the projections in the correction equation (10). This
observation was made in [28], and a new algorithm was proposed for eigenproblems by
including the projections in the Davidson scheme. In addition, these modified schemes
allow for more general approximate solutions p k than . For instance, the
use of ' steps of a preconditioned Krylov subspace method for the correction equation
is suggested, leading to Arnoldi type of methods in which the variable polynomial pre-conditioning
is determined efficiently and the projections have been included correctly.
The new methods have been called Jacobi-Davidson methods (Jacobi took proper care
of the projections, but did not build a search subspace as Davidson did (see [28] for
details and further references)).
The analysis and results in [3, 27] show that these Jacobi-Davidson methods can
also be effective for solving generalized eigenproblems, even without any matrix inversion

The Jacobi-Davidson methods allow for a variety of choices that may improve efficiency
of the steps and speed of convergence and are good examples of AIN methods
in which the projected problem (7) is used to steer the computation.
For an extensive discussion, we refer to [27].
7. AIN schemes for general nonlinear problems. In this section we summarize
some iterative methods for the solution of nonlinear problems that have been
proposed by different authors, and we show how these methods fit in the AIN frame-work

Brown and Saad [4] describe a family of methods for solving nonlinear problems.
They refer to these methods as nonlinear Krylov subspace projection methods. Their
12 D. R. FOKKEMA, G. L. G. SLEIJPEN, AND H. A. VAN DER VORST
modifications to Newton's method are intended to enhance robustness and are heavily
influenced by ideas presented in [9]. One of their methods is a variant of Damped
Inexact Newton, in which they approximate the solution of the correction equation
by a few steps of Arnoldi or GMRES and determine the damping parameter ff by
a "linesearch backtracking technique". So this is just another AIN scheme, with a
special 1-dimensional subspace acceleration. They also propose a model trust region
approach, where they take their update to the approximation from the Krylov subspace
e
Vm generated by m steps of (preconditioned) Arnoldi or GMRES as
Vm ~
~ y k is the point on the dogleg curve for which k~y k k the trust region size: ~ y k is an
approximation for min y
Vm ~
This could be considered as a block version
of the previous method.
In [2] Axelsson and Chronopoulos propose two nonlinear versions of a (truncated)
Generalized Conjugate Gradient type of method. Both methods fit in the AIN frame-
work. The first method, NGCG, is a Minimal Residual AIN method with
and V k orthonormal; in other words the correction equation is not solved. The second
method, NNGCG, differs from NGCG in that p k is now computed as an approximate
solution (by some method) of the correction equation (6), where the accuracy
is such that non-increasing sequence
[8]). So the method NNGCG is a Minimal Residual AIN
method. It can be viewed as generalization of GMRESR [29]. Under certain conditions
on the map F they prove global convergence.
In [15], Kaporin and Axelsson propose a class of nonlinear equation solvers, GNKS,
in which the ideas presented in [4] and [2] are combined. There, the direction vectors
are obtained as linear combinations of the columns of e
Vm and V k . To be more
precise,
This problem
is then solved by a special Gauss-Newton iteration scheme, which avoids excessive
computational work, by taking into account the acute angle between r k and J k p k , and
the rate of convergence. The method generalizes GCRO [7].
8. Numerical Experiments. In this section we test several AIN schemes and
present results of numerical experiments on three different nonlinear problems. For
tests and test results with methods for linear- and eigenproblems we refer to their
references. The purpose of this presentation is to show that acceleration may be
useful also in the nonlinear case. By useful, we mean that additional computational
cost is compensated for by faster convergence.
Different AIN schemes distinguish themselves by the way they (approximately)
solve the correction equation and the projected problem (cf. Section 3.2 and 7). Out
of the overwhelming variety of choices we have selected a few possible combinations,
some of which lead to AIN schemes that are equivalent to already proposed methods
and some of which lead to new methods. We compare the following (existing) Minimal
Residual AIN schemes:
ffl linesearch, the backtracking linesearch technique [4, pp. 458];
ffl dogleg, the model trust region approach as proposed in [4, pp. 462];
ffl nngcg, a variant of the method proposed in [2], solving (8) by the Levenberg-Marquardt
algorithm [20];
ACCELERATED INEXACT NEWTON SCHEMES 13
ffl gnks, the method proposed in [15];
and the (new) Restricted Minimal Residual AIN schemes:
ffl rmr a, choosing W
ffl rmr b, choosing W is the component of J k p k orthogonal
to W k\Gamma1 .
For these last two schemes, the minimization problem solved by the Gauss-Newton
variant described in [15]. The necessary subspaces for the direction p k or the
projected problem were obtained by 10 steps of GMRES, or (in the third example)
also by at most 50 iterations of the generalized CGS variant CGS2 [10].
In all cases the exact Jacobian was used. Furthermore, we used orthonormal
matrices V k and W k , obtained from a modified Gram-Schmidt process and restricted
to the last 10 columns in an attempt to save computational work. The computations
were done on a Sun Sparc 20 in double precision and the iterations were stopped when
method failed, either when the convergence was too slow, i.e., when
or when to number of nonlinear iterations (per step)
exceeded 200.
Since the computational cost of the methods is approximately proportional to the
costs of the number of function evaluations and matrix multiplications, the following
counters are given in the tables:
ffl ni, the number of nonlinear iterations;
ffl fe, the number of function evaluations;
ffl mv, the number of multiplications by the Jacobian;
ffl pre, the number of applications of the preconditioner;
ffl total, the sum of fe, mv and pre.
8.1. A 1D Burgers' Equation. As a first test problem we consider the following
1D Burgers' Equation [14]
@t
@x
We discretized the spatial variable x with finite differences with 64 grid points
and for the time derivative we used
\Deltat
with denotes the solution at time t n\Deltat. For this
test the solution u n was computed for and as an initial guess to u n+1
we took u n . No preconditioning was used.
In table Tab. 1 we show the results for problem (12) with
14 D. R. FOKKEMA, G. L. G. SLEIJPEN, AND H. A. VAN DER VORST
Method ni fe mv total
linesearch 594 1818 4853 6671
dogleg 644 3559 6140 9699
nngcg 229 4426 11769 16195
gnks 187 852 7067 7919
rmr a 230 926 4471 5397

Table

1: Results for Burgers' Equation
A plot of the solutions is given in Fig. 1. The table shows the cumulative
value of the counters for each method after completing the computation of u
If we look at the number of nonlinear iterations (ni), we see that acceleration
indeed reduces this number. However, in the case of gnks this does not result in less
work, because the number of matrix multiplications (mv) increases too much. Here
both the Galerkin approaches rmr a and rmr b are less expensive than all the other
methods. rmr a being the winner.
x

Figure

1: Solution of Burgers' Equation
8.2. The Bratu problem. As a second test problem we consider the Bratu
problem [12, 4]. We seek a solution (u; -) of the nonlinear boundary value problem:
For\Omega we took the unit square and we discretized with finite differences on a 31 \Theta 31
regular grid. It is known, cf. [12], that there exist a critical value -   such that for
problem (13) has two solutions and for -   problem (13) has no
solutions. In order to locate this critical value we use the arc length continuation
method as described in [12, section 2.3 and 2.4]. Problem (13) is replaced by a problem
ACCELERATED INEXACT NEWTON SCHEMES 15
Method ni fe mv pre total
linesearch 391 1013 3732 3421 8166
dogleg 381 2664 3010 3010 8684
nngcg 361 1297 4243 3091 8631
gnks 358 1056 6896 2780 10732
rmr a 389 539 4005 3399 7943

Table

2: Results for the Bratu Problem, solved by the arc length continuation
method.
Method ni fe mv pre total
linesearch 29 85 336 308 729
dogleg
gnks 38 119 1806 370 2295
rmr a 6 13 77 55 145

Table

3: Single solve of the Bratu Problem, u
of the form
where ', a scalar valued function, is chosen such that s is some arc length on the
solution branch and u s is the solution of (13) for -(s). We preconditioned GMRES
by ILU(0) [18] of the discretized Laplace operator \Delta.
The first table Tab. 2 shows the results after a full continuation run: starting
from the smallest solution (u; -) with solution branch is followed along the
(discretized) arc with s
we see that acceleration may be useful, in spite of the fact that there is little room for
it, because on the average approximately only 4.5 Newton iterations where necessary
to compute the solution per continuation step. In this example rmr b performs better
than rmr a.

Table

Tab. 3 shows the results for the case where we solve (13) for fixed
(near the critical value). In this case Galerkin acceleration is even more useful and
the differences are more pronounced.
The sup norm of the solution for the different values of - are plotted in Fig. 2.
The two solutions at - 4 along the diagonal of the unit square are shown in Fig. 3.
8.3. The Driven Cavity Problem. In this Section we present test results for
the classical driven cavity problem from incompressible fluid flow. We follow closely
D. R. FOKKEMA, G. L. G. SLEIJPEN, AND H. A. VAN DER VORST

Figure

2: Sup norms of the solution u
along the arc.

Figure

3: Solutions u at - 4 along
the diagonal of the domain.
the presentations in [12, 4]. In stream function-vorticity formulation the equations are
@/
@n
where\Omega is the unit square and the viscosity - is the reciprocal of the Reynolds number
Re. In terms of / alone this can be written as
subject to the same boundary conditions. This equation was discretized with finite
differences on a 25 \Theta 25 grid, see Fig. 4. The grid lines are distributed as the roots
of the Chebychev polynomial of degree 25. As preconditioner we used the Modified
[13] decomposition of the biharmonic operator . Starting from the solution
for computed several solutions, using the the arc length continuation
method (cf. the previous example, and [12]) with step sizes \Deltas = 100 for 0 - Re -
Tab. 4 shows the results of this test when using 10 steps of GMRES and CGS2 [10]
for the correction equation. In the case of CGS2 we approximately solved the correction
equation to a relative residual norm precision of 2 \Gammak , where k is the current
Newton step [8], with a maximum of 50 steps. Clearly, the methods using (the basis
produced steps of GMRES perform very poorly for this example. Only gnks
is able to complete the full continuation run, but requires a large number of Newton
steps. If we look at the results for the AIN schemes for which CGS2 is used, we see
that, except for the linesearch method, these methods perform much better. The
Restricted Minimal Residual methods are again the most efficient ones.
ACCELERATED INEXACT NEWTON SCHEMES 17
GMRES
Method ni fe mv pre total
linesearch fails at Re = 400 after a total of 545.
dogleg fails at Re = 100 after a total of 113.
nngcg fails at Re = 2200 after a total of 19875
gnks 641 2315 30206 6210 38731
rmr a fails at Re = 2000 after a total of 13078
rmr b fails at Re = 800 after a total of 7728
Method ni fe mv pre total
linesearch fails at Re = 1300 after a total of 4342.
rmr a 137 297 6266 5969 12532

Table

4: Results for the Driven Cavity problem, solved by the arc length continuation
method for Re

Figure

4: Grid for the Driven Cavity
problem, (25 \Theta 25).
Figure

5: Stream lines of the Driven
Cavity problem, Re = 100.
This test also reveals a possible practical drawback of methods like dogleg and
gnks. These methods exploit an affine subspace to find a suitable update for the
approximation. This may fail, when the problem is hard or when the preconditioner
is not good enough. In that case the dimension of the affine subspace must be large,
which may be, because of storage requirements and computational overhead, not fea-
sible. For the schemes that use approximate solutions of the correction equation,
delivered by some arbitrary, iterative method, e.g., CGS2, one can easily adapt the
precision, which leaves more freedom.
D. R. FOKKEMA, G. L. G. SLEIJPEN, AND H. A. VAN DER VORST
Figure

lines of the Driven
Cavity problem, Re = 400.
Figure

7: Stream lines of the Driven
Cavity problem, Re = 1600.
Figure

8: Stream lines of the Driven
Cavity problem, Re = 2000.
Figure

9: Stream lines of the Driven
Cavity problem, Re = 3000.
Plots of the stream lines for the values
0:0; 0:0025; 0:001; 0:0005; 0:0001; 0:00005
(cf. [12]) are given in Fig. 5-9. The plots show virtually the same solutions as in [12].
9. Conclusions. We have shown how the classical Newton iteration scheme for
nonlinear problems can be accelerated in a similar way as standard Richardson-type
iteration schemes for linear equations. This leads to the AIN framework in which
ACCELERATED INEXACT NEWTON SCHEMES 19
many well-known iterative methods for linear-, eigen-, and general nonlinear problems
fit. From this framework an overwhelming number of possible iterations schemes can
be formulated. We have selected a few and shown by numerical experiments that
especially the Restricted Minimal Residual methods can be very useful for further
reducing computational costs.



--R

The principle of minimized iterations in the solution of the matrix eigenvalue problem
On nonlinear generalized conjugate gradient methods
te Riele
Hybrid Krylov methods for nonlinear systems of equations
The Davidson method
The iterative calculation of a few of the lowest eigenvalues and corresponding eigenvectors of large real symmetric matrices
Nested Krylov methods and preserving the orthogonality


Generalized conjugate gradient squared
QMR: A quasi minimal residual method for non-Hermitian linear systems

A class of first order factorizations methods
The partial differential equation u t
On a class of nonlinear equation solvers based on the residual norm reduction over a sequence of affine subspaces
Acceleration techniques for decoupling algorithms in semiconductor simulation
Solution of systems of linear equations by minimized iteration
An iterative solution method for linear systems of which the coefficient matrix is a symmetric M-matrix
An algorithm for generalized matrix eigenvalue problems

Generalizations of Davidson's method for computing eigenvalues of large non-symmetric matrices
Approximate solutions and eigenvalue bounds from Krylov subspaces
Krylov subspace method for solving large unsymmetric linear systems


GMRES: A generalized minimum residual algorithm for solving nonsymmetric linear systems

A Jacobi-Davidson iteration method for linear eigenvalue problems
GMRESR: A family of nested GMRES methods
Further experiences with GMRESR
Generalized conjugate-gradient acceleration of nonsymmetrizable iterative methods
--TR

--CTR
Keith Miller, Nonlinear Krylov and moving nodes in the method of lines, Journal of Computational and Applied Mathematics, v.183 n.2, p.275-287, 15 November 2005
P. R. Graves-Morris, BiCGStab, VPAStab and an adaptation to mildly nonlinear systems, Journal of Computational and Applied Mathematics, v.201 n.1, p.284-299, April, 2007
Heng-Bin An , Ze-Yao Mo , Xing-Ping Liu, A choice of forcing terms in inexact Newton method, Journal of Computational and Applied Mathematics, v.200 n.1, p.47-60, March, 2007
Heng-Bin An , Zhong-Zhi Bai, A globally convergent Newton-GMRES method for large sparse systems of nonlinear equations, Applied Numerical Mathematics, v.57 n.3, p.235-252, March, 2007
Monica Bianchini , Stefano Fanelli , Marco Gori, Optimal Algorithms for Well-Conditioned Nonlinear Systems of Equations, IEEE Transactions on Computers, v.50 n.7, p.689-698, July 2001
nonlinear problems;iterative methods;newton's method;inexact Newton