107

INEXACT NEWTON METHODS

and then F (x+ ) will be evaluated to test for termination. We also count

the evaluation of F (x0 ) as an additional cost in the evaluation of x1 , which

therefore has a minimum cost of three function evaluations. For each example

we compare a constant value of ·n with the choice given in (6.20) and used as

the default in nsolgm.

6.4.1. Chandrasekhar H-equation. We solve the H-equation on a 100-

point mesh with c = .9 using two schemes for selection of the parameter ·.

The initial iterate is the function identically one. We set the parameters in

(6.20) to γ = .9 and ·max = .25.

We used „r = „a = 10’6 in this example.

In Fig. 6.1 we plot the progress of the iteration using · = .1 with the solid

line and using the sequence given by (6.20) with the dashed line. We set the

parameters in (6.20) to γ = .9 and ·max = .25. This choice of ·max was the

one that did best overall in our experiments on this problem.

We see that the constant · iteration terminates after 12 function evalua-

tions, 4 outer iterations, and roughly 275 thousand ¬‚oating-point operations.

This is a slightly higher overall cost than the other approach in which {·n } is

given by (6.20), which terminated after 10 function evaluations, 3 outer itera-

tions, and 230 thousand ¬‚oating-point operations. The chord method with the

di¬erent (but very similar for this problem) l∞ termination condition for the

same problem, reported in § 5.6 incurred a cost of 3.3 million ¬‚oating-point

operations, slower by a factor of over 1000. This cost is entirely a result of the

computation and factorization of a single Jacobian.

0

10

-1

10

-2

10

Relative Nonlinear Residual

-3

10

-4

10

-5

10

-6

10

-7

10

-8

10

0 2 4 6 8 10 12

Function Evaluations

Fig. 6.1. Newton-GMRES for the H-equation, c = .9.

In Fig. 6.2 the results change for the more ill-conditioned problem with

c = .9999. Here the iteration with ·n = .1 performed slightly better and

Buy this book from SIAM at http://www.ec-securehost.com/SIAM/FR16.html.

Copyright ©1995 by the Society for Industrial and Applied Mathematics. This electronic version is for personal use and may not be duplicated or distributed.

108 ITERATIVE METHODS FOR LINEAR AND NONLINEAR EQUATIONS

0

10

-1

10

Relative Nonlinear Residual

-2

10

-3

10

-4

10

-5

10

-6

10

-7

10

0 5 10 15 20 25

Function Evaluations

Fig. 6.2. Newton-GMRES for the H-equation, c = .9999

0

10

-1

10

Relative Nonlinear Residual

-2

10

-3

10

-4

10

0 2 4 6 8 10 12 14 16 18

Function Evaluations

Fig. 6.3. Newton-GMRES for the PDE, C = 20.

terminated after 22 function evaluations and 7 outer iterations, and roughly

513 thousand ¬‚oating-point operations. The iteration with the decreasing {·n }

given by (6.20) required 23 function evaluations, 7 outer iterations, and 537

thousand ¬‚oating-point operations.

6.4.2. Convection-di¬usion equation. We consider the partial di¬eren-

tial equation

’ ∇2 u + Cu(ux + uy ) = f

(6.21)

Buy this book from SIAM at http://www.ec-securehost.com/SIAM/FR16.html.

Copyright ©1995 by the Society for Industrial and Applied Mathematics. This electronic version is for personal use and may not be duplicated or distributed.

109

INEXACT NEWTON METHODS

with homogeneous Dirichlet boundary conditions on the unit square (0, 1) —

(0, 1). f has been constructed so that the exact solution was the discretization

of

10xy(1 ’ x)(1 ’ y) exp(x4.5 ).

We set C = 20 and u0 = 0. As in § 3.7 we discretized (6.21) on a 31 — 31 grid

using centered di¬erences. To be consistent with the second-order accuracy of

the di¬erence scheme we used „r = „a = h2 , where h = 1/32.

We compare the same possibilities for {·n } as in the previous example and

report the results in Fig. 6.3. In (6.20) we set γ = .9 and ·max = .5. The

computation with constant · terminated after 19 function evaluations and 4

outer iterates at a cost of roughly 3 million ¬‚oating-point operations. The

iteration with {·n } given by (6.20) required 16 function evaluations, 4 outer

iterations, and 2.6 million ¬‚oating-point operations. In the computation we

preconditioned (6.21) with the fast Poisson solver fish2d.

In this case, the choice of {·n } given by (6.20) reduced the number of inner

iterations in the early stages of the iteration and avoided oversolving. In cases

where the initial iterate is very good and only a single Newton iterate might be

needed, however, a large choice of ·max or constant · may be the more e¬cient

choice. Examples of such cases include (1) implicit time integration of nonlinear

parabolic partial di¬erential equations or ordinary di¬erential equations in

which the initial iterate is either the converged solution from the previous

time step or the output of a predictor and (2) solution of problems in which

the initial iterate is an interpolant of the solution from a coarser mesh.

Preconditioning is crucial for good performance. When preconditioning

was not done, the iteration for C = 20 required more than 80 function

evaluations to converge. In Exercise 6.5.9 the reader is asked to explore this.

Buy this book from SIAM at http://www.ec-securehost.com/SIAM/FR16.html.

Copyright ©1995 by the Society for Industrial and Applied Mathematics. This electronic version is for personal use and may not be duplicated or distributed.

110 ITERATIVE METHODS FOR LINEAR AND NONLINEAR EQUATIONS

6.5. Exercises on inexact Newton methods

6.5.1. Verify (6.8).

6.5.2. Prove Proposition 6.1.1 by showing that if the standard assumptions

hold, 0 < < 1, and x is su¬ciently near x— then

F (x) /(1 + ) ¤ F (x— )e ¤ (1 + ) F (x) .

6.5.3. Verify (6.14).

6.5.4. Prove Theorems 6.2.2 and 6.2.1.

6.5.5. Prove Theorems 6.1.5 and 6.1.6.

6.5.6. Can anything like Proposition 6.2.1 be proved for a ¬nite-di¬erence Bi-

CGSTAB or TFQMR linear solver? If not, what are some possible

implications?

6.5.7. Give an example in which F (xn ) ’ 0 q-linearly but xn does not converge

to x— q-linearly.

6.5.8. In the DASPK code for integration of di¬erential algebraic equations,

[23], preconditioned Newton-GMRES is implemented in such a way

that the data needed for the preconditioner is only computed when

needed, which may be less often that with each outer iteration. Discuss

this strategy and when it might be useful. Is there a relation to the

Shamanskii method?

6.5.9. Duplicate the results in § 6.4. For the H-equation experiment with

various values of c, such as c = .5, .99999, .999999, and for the convection-

di¬usion equation various values of C such as C = 1, 10, 40 and how

preconditioning a¬ects the iteration. Experiment with the sequence of

forcing terms {·n }. How do the choices ·n = 1/(n + 1), ·n = 10’4

[32], ·n = 21’n [24], ·n = min( F (xn ) 2 , (n + 2)’1 ) [56], ·n = .05, and

·n = .75 a¬ect the results? Present the results of your experiments as

graphical iteration histories.

6.5.10. For the H-equation example in § 6.4, vary the parameter c and see how

the performance of Newton-GMRES with the various choices of {·n } is

a¬ected. What happens when c = 1? See [119] for an explanation of

this.

6.5.11. Are the converged solutions for the preconditioned and unpreconditioned

convection-di¬usion example in § 6.4 equally accurate?

6.5.12. For the convection-di¬usion example in § 6.4, how is the performance

a¬ected if GMRES(m) is used instead of GMRES?

Buy this book from SIAM at http://www.ec-securehost.com/SIAM/FR16.html.

Copyright ©1995 by the Society for Industrial and Applied Mathematics. This electronic version is for personal use and may not be duplicated or distributed.

111

INEXACT NEWTON METHODS

6.5.13. Apply Newton-CGNR, Newton-Bi-CGSTAB, or Newton-TFQMR to the

two examples in § 6.4. How will you deal with the need for a transpose in

the case of CGNR? How does the performance (in terms of time, ¬‚oating-

point operations, function evaluations) compare to Newton-GMRES?

Experiment with the forcing terms. The MATLAB codes fdkrylov,

which allows you to choose from a variety of forward di¬erence Krylov

methods, and nsola, which we describe more fully in Chapter 8, might

be of use in this exercise.

6.5.14. If you have done Exercise 5.7.25, apply nsolgm to the same two-point

boundary value problem and compare the performance.

Buy this book from SIAM at http://www.ec-securehost.com/SIAM/FR16.html.

Copyright ©1995 by the Society for Industrial and Applied Mathematics. This electronic version is for personal use and may not be duplicated or distributed.

112 ITERATIVE METHODS FOR LINEAR AND NONLINEAR EQUATIONS

Buy this book from SIAM at http://www.ec-securehost.com/SIAM/FR16.html.

Copyright ©1995 by the Society for Industrial and Applied Mathematics. This electronic version is for personal use and may not be duplicated or distributed.

Chapter 7

Broyden™s method

Quasi-Newton methods maintain approximations of the solution x— and the

Jacobian at the solution F (x— ) as the iteration progresses. If xc and Bc are

the current approximate solution and Jacobian, then

’1

(7.1) x+ = xc ’ Bc F (xc ).

After the computation of x+ , Bc is updated to form B+ . The construction of

B+ determines the quasi-Newton method.

The advantages of such methods, as we shall see in § 7.3, is that solution of

equations with the quasi-Newton approximation to the Jacobian is often much

cheaper than using F (xn ) as the coe¬cient matrix. In fact, if the Jacobian is

dense, the cost is little more than that of the chord method. The method we

present in this chapter, Broyden™s method, is locally superlinearly convergent,

and hence is a very powerful alternative to Newton™s method or the chord

method.

In comparison with Newton-iterative methods, quasi-Newton methods

require only one function evaluation for each nonlinear iterate; there is no cost

in function evaluations associated with an inner iteration. Hence, if a good

preconditioner (initial approximation to F (x— )) can be found, these methods

could have an advantage in terms of function evaluation cost over Newton-

iterative methods.

Broyden™s method [26] computes B+ by

(y ’ Bc s)sT F (x+ )sT

(7.2) B+ = Bc + = Bc + .

sT s sT s

In (7.2) y = F (x+ ) ’ F (xc ) and s = x+ ’ xc .

Broyden™s method is an example of a secant update. This means that the

updated approximation to F (x— ) satis¬es the secant equation

(7.3) B+ s = y.

In one dimension, (7.3) uniquely speci¬es the classical secant method. For

equations in several variables (7.3) alone is a system of N equations in N 2

113

Buy this book from SIAM at http://www.ec-securehost.com/SIAM/FR16.html.

Copyright ©1995 by the Society for Industrial and Applied Mathematics. This electronic version is for personal use and may not be duplicated or distributed.

114 ITERATIVE METHODS FOR LINEAR AND NONLINEAR EQUATIONS

unknowns. The role of the secant equation is very deep. See [62] and [64] for

discussion of this topic. In this book we focus on Broyden™s method and our

methods of analysis are not as general as those in [62] and [64]. In this chapter

we will assume that the function values are accurate and refer the reader to [66],

[65], and [114] for discussions of quasi-Newton methods with inaccurate data.

Locally superlinearly convergent secant methods can be designed to maintain

the sparsity pattern, symmetry, or positive de¬niteness, of the approximate

Jacobian, and we refer the reader to [28], [63], [62], and [64] for discussion and

analysis of several such methods. More subtle structural properties can also

be preserved by secant methods. Some examples can be found in [101], [96],

[95], [113], and [116].

Broyden™s method is also applicable to linear equations [29], [82], [84],

[104], [143], [130], Ax = b, where B ≈ A is updated. One can also apply the

method to linear least squares problems [84], [104]. The analysis is somewhat

clearer in the linear case and the plan of this chapter is to present results for

the linear case for both the theory and the examples. One interesting result,

which we will not discuss further, is that for linear problems Broyden™s method

converges in 2N iterations [29], [82], [84], [143].

We will express our results in terms of the error in the Jacobian

E = B ’ F (x— )

(7.4)

and the step s = x+ ’xc . While we could use E = B ’F (x), (7.4) is consistent

with our use of e = x ’ x— . When indexing of the steps is necessary, we will

write

sn = xn+1 ’ xn .

In this chapter we show that if the data x0 and B0 are su¬ciently good,

the Broyden iteration will converge q-superlinearly to the root.

7.1. The Dennis“Mor´ condition

e

In this section we consider general iterative methods of the form

’1

(7.5) xn+1 = xn ’ Bn F (xn )

where Bn = F (x— )+En ≈ F (x— ) is generated by some method (not necessarily

Broyden™s method).

Veri¬cation of q-superlinear convergence cannot be done by the direct

approaches used in Chapters 5 and 6. We must relate the superlinear

convergence condition that ·n ’ 0 in Theorem 6.1.2 to a more subtle, but

easier to verify condition in order to analyze Broyden™s method. This technique

is based on veri¬cation of the Dennis“Mor´ condition [61], [60] on the sequence

e

of steps {sn } and errors in the Jacobian {En }

En sn

(7.6) lim = 0.

sn

n’∞