Local q-linear convergence. Theorem 7.2.3 states that Broyden™s method

is locally q-linearly convergent. The proof is more complex than similar results

in Chapters 5 and 6. We must explore the relation between the size of E0

needed to enforce q-linear convergence and the q-factor.

Theorem 7.2.3. Let the standard assumptions hold. Let r ∈ (0, 1) be

given. Then there are δ and δB such that if x0 ∈ B(δ) and E0 2 < δB the

Broyden sequence for the data (F, x0 , B0 ) exists and xn ’ x— q-linearly with

q-factor at most r.

Proof. Let δ and δ1 be small enough so that the conclusions of Theo-

rem 5.4.3 hold. Then, reducing δ and δ1 further if needed, the q-factor is at

most

KA (δ + δ1 ) ¤ r,

where KA is from the statement of Theorem 5.4.3.

We will prove the result by choosing δB and reducing δ if necessary so that

E0 2 ¤ δB will imply that En 2 ¤ δ1 for all n, which will prove the result.

To do this reduce δ if needed so that

γ(1 + r)δ

(7.33) δ2 = < δ1

2(1 ’ r)

and set

(7.34) δB = δ1 ’ δ2 .

We show that En 2 ¤ δ1 by induction. Since E0 2 < δB < δ1 we may

begin the induction. Now assume that Ek 2 < δ1 for all 0 ¤ k ¤ n. By

Theorem 7.2.2,

En+1 ¤ En + γ( en+1 + en 2 )/2.

2 2 2

Since En < δ1 , en+1 ¤ r en and therefore

2 2 2

+ γ(1 + r)rn δ/2.

En+1 ¤ En + γ(1 + r) en 2 /2 ¤ En

2 2 2

We can estimate En 2 , En’1 2 , . . . in the same way to obtain

n

(1 + r)γδ

rj

En+1 ¤ E0 2 +

2

2 j=0

(1 + r)γδ

¤ δB + ¤ δ1 ,

2(1 ’ r)

which completes the proof.

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.

123

BROYDEN™S METHOD

Veri¬cation of the Dennis“Mor´ condition. The ¬nal task is to verify

e

the Dennis“Mor´ condition.

e

Theorem 7.2.4. Let the standard assumptions hold. Then there are δ and

δB such that if x0 ∈ B(δ) and E0 2 < δB the Broyden sequence for the data

(F, x0 , B0 ) exists and xn ’ x— q-superlinearly.

Proof. Let δ and δB be such that the conclusions of Theorem 7.2.3 hold.

By Theorem 7.1.1 we need only show that Dennis“Mor´ condition (7.6) holds

e

to complete the proof.

As in the proof of Theorem 7.2.1 let

sn sT

n

Pn = .

sn 22

Set

1

(F (xn + tsn ) ’ F (x— )) dt.

∆n =

0

Let φ ∈ RN be arbitrary. Note that

En+1 φ = (I ’ Pn )En φ + Pn ∆T φ.

T T

n

We wish to apply Lemma 7.2.1 with

T

= Pn ∆T φ.

ψn = En φ, ·n = sn / sn 2 , and n n

The hypothesis in the lemma that

<∞

n2

n

holds since Theorem 7.2.3 and the standard assumptions imply

¤ γ(1 + r)rn δ/2.

∆n 2

Hence Lemma 7.2.1 implies that

(En φ)T sn

T En sn

T

= φT

·n ψn = ’ 0.

sn 2 sn 2

This, as in the proof of Theorem 7.2.1, implies the Dennis“Mor´ condition and

e

completes the proof.

7.3. Implementation of Broyden™s method

The implementation of Broyden™s method that we advocate here is related

to one that has been used in computational ¬‚uid mechanics [73]. Some

such methods are called limited memory formulations in the optimization

literature [138], [142], [31]. In these methods, after the storage available for

the iteration is exhausted, the oldest of the stored vectors is replaced by the

most recent. Another approach, restarting, clears the storage and starts over.

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.

124 ITERATIVE METHODS FOR LINEAR AND NONLINEAR EQUATIONS

This latter approach is similar to the restarted GMRES iteration. Our basic

implementation is a nonlinear version of the implementation for linear problems

in [67]. We use a restarted approach, though a limited memory approach could

be used as well.

We begin by showing how the initial approximation to F (x— ), B0 , may be

regarded as a preconditioner and incorporated into the de¬ntion of F just as

preconditioners were in the implementation of Newton-GMRES in § 6.4.

Lemma 7.3.1. Assume that the Broyden sequence ({xn }, {Bn }) for the data

(F, x0 , B0 ) exists. Then the Broyden sequence ({yn }, {Cn }) exists for the data

’1

(B0 F, x0 , I) and

(7.35) xn = yn , Bn = B0 Cn

for all n.

Proof. We proceed by induction. Equation (7.35) holds by de¬nition when

n = 0. If (7.35) holds for a given n, then

’1 ’1

yn+1 = yn ’ Cn B0 F (yn )

= xn ’ (B0 Cn )’1 F (xn ) = xn ’ Bn F (xn )

’1

(7.36)

= xn+1 .

Now sn = xn+1 ’ xn = yn+1 ’ yn by (7.36). Hence

F (xn+1 )sT

n

Bn+1 = Bn + 2

sn 2

’1

B0 F (yn+1 )sT

n

= B0 Cn + B0 = B0 Cn+1 .

2

sn 2

This completes the proof.

The consequence of this proof for implementation is that we can assume

that B0 = I. If a better choice for B0 exists, we can incorporate that into the

’1

function evaluation. Our plan will be to store Bn in a compact and easy to

apply form. The basis for this is the following simple formula [68], [176], [177],

[11].

Proposition 7.3.1. Let B be a nonsingular N — N matrix and let

u, v ∈ RN . Then B + uv T is invertible if and only if 1 + v T B ’1 u = 0. In

this case,

(B ’1 u)v T

T ’1

B ’1 .

(7.37) (B + uv ) = I ’ T B ’1 u

1+v

The expression (7.37) is called the Sherman“Morrison formula. We leave

the proof to Exercise 7.5.2.

In the context of a sequence of Broyden updates {Bn } we have for n ≥ 0,

T

Bn+1 = Bn + un vn ,

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.

125

BROYDEN™S METHOD

where

un = F (xn+1 )/ sn and vn = sn / sn 2 .

2

Setting

’1 T ’1

wn = (Bn un )/(1 + vn Bn un )

we see that, if B0 = I,

’1 = (I ’ w T T T

Bn n’1 vn’1 )(I ’ wn’2 vn’2 ) . . . (I ’ w0 v0 )

(7.38)

n’1 T

= (I ’ wj vj ).

j=0

Since the empty matrix product is the identity, (7.38) is valid for n ≥ 0.

’1

Hence the action of Bn on F (xn ) (i.e., the computation of the Broyden

n’1

step) can be computed from the 2n vectors {wj , vj }j=0 at a cost of O(N n)

¬‚oating-point operations. Moreover, the Broyden step for the following

iteration is

n’1

’1 T

(7.39) sn = ’Bn F (xn ) = ’ (I ’ wj vj )F (xn ).

j=0

Since the product

n’2 T

(I ’ wj vj )F (xn )

j=0

must also be computed as part of the computation of wn’1 we can combine

the computation of wn’1 and sn as follows:

n’2 T

w = (I ’ wj vj )F (xn )

j=0

(7.40) + vn’1 w)’1

T

wn’1 = Cw w where Cw = ( sn’1 2

T

sn = ’(I ’ wn’1 vn’1 )w.

One may carry this one important step further [67] and eliminate the need

to store the sequence {wn }. Note that (7.40) implies that for n ≥ 1

T ’1

sn = ’w + Cw w(vn’1 w) = w(’1 + Cw (Cw ’ sn’1 2 ))

= ’Cw w sn’1 = ’ sn’1 2 wn’1 .

2

Hence, for n ≥ 0,

(7.41) wn = ’sn+1 / sn 2 .

Therefore, one need only store the steps {sn } and their norms to construct

the sequence {wn }. In fact, we can write (7.38) as

sj+1 sT

n’1 j

’1

(7.42) Bn = I+ .

sj 2

j=0

2

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.

126 ITERATIVE METHODS FOR LINEAR AND NONLINEAR EQUATIONS

We cannot use (7.42) and (7.39) directly to compute sn+1 because sn+1 appears

on both sides of the equation

sj+1 sT

sn+1 sT n’1 j

n

sn+1 =’ I+ I+ F (xn+1 )

2 sj 2

sn 2 j=0

2

(7.43)

sn+1 sT

n ’1

=’ I+ Bn F (xn+1 ).

2

sn 2

Instead, we solve (7.43) for sn+1 to obtain

’1

Bn F (xn+1 )

(7.44) sn+1 =’ .

’1 2

1 + sT Bn F (xn+1 )/ sn

n 2

By Proposition 7.3.1, the denominator in (7.44)

1 + sT Bn F (xn+1 )/ sn

’1 2 T ’1

= 1 + vn Bn un

n 2

is nonzero unless Bn+1 is singular.

The storage requirement, of n + O(1) vectors for the nth iterate, is of the

same order as GMRES, making Broyden™s method competitive in storage with

GMRES as a linear solver [67]. Our implementation is similar to GMRES(m),

restarting with B0 when storage is exhausted.

We can use (7.42) and (7.44) to describe an algorithm. The termination

criteria are the same as for Algorithm nsol in § 5.6. We use a restarting

approach in the algorithm, adding to the input list a limit nmax on the number

of Broyden iterations taken before a restart and a limit maxit on the number

of nonlinear iterations. Note that B0 = I is implicit in the algorithm. Our

looping convention is that the loop in step 2(e)ii is skipped if n = 0; this is

consistent with setting the empty matrix product to the identity. Our notation

and algorithmic framework are based on [67].

Algorithm 7.3.1. brsol(x, F, „, maxit, nmax)

1. r0 = F (x) 2 , n = ’1,

s0 = ’F (x), itc = 0.

2. Do while itc < maxit

(a) n = n + 1; itc = itc + 1

(b) x = x + sn

(c) Evaluate F (x)

(d) If F (x) ¤ „r r0 + „a exit.

2

(e) if n < nmax then

i. z = ’F (x)

ii. for j = 0, n ’ 1

z = z + sj+1 sT z/ sj 2

2

j

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.

127

BROYDEN™S METHOD

iii. sn+1 = z/(1 ’ sT z/ sn 2 )

n 2

(f) if n = nmax then

n = ’1; s = ’F (x);

If n < nmax, then the nth iteration of brsol, which computes xn and

sn+1 , requires O(nN ) ¬‚oating-point operations and storage of n + 3 vectors,

the steps {sj }n , x, z, and F (x), where z and F (x) may occupy the same

j=0

storage location. This requirement of one more vector than the algorithm

in [67] needed is a result of the nonlinearity. If n = nmax, the iteration

will restart with xnmax , not compute snmax+1 , and therefore need storage for

nmax + 2 vectors. Our implementation of brsol in the collection of MATLAB

codes stores z and F separately in the interest of clarity and therefore requires

nmax + 3 vectors of storage.

The Sherman“Morrison approach in Algorithm brsol is more e¬cient, in

terms of both time and storage, than the dense matrix approach proposed

in [85]. However the dense matrix approach makes it possible to detect

ill conditioning in the approximate Jacobians. If the data is su¬ciently

good, bounded deterioration implies that the Broyden matrices will be well

conditioned and superlinear convergence implies that only a few iterates will

be needed. In view of all this, we feel that the approach of Algorithm brsol

is the best approach, especially for large problems.

A MATLAB implementation of brsol is provided in the collection of

MATLAB codes.

7.4. Examples for Broyden™s method

√

In all the examples in this section we use the l2 norm multiplied by 1/ N .

We use the MATLAB code brsol.

7.4.1. Linear problems. The theory in this chapter indicates that Broy-

den™s method can be viewed as an alternative to GMRES as an iterative method

for nonsymmetric linear equations. This point has been explored in some detail

in [67], where more elaborate numerical experiments that those given here are

presented, and in several papers on in¬nite-dimensional problems as well [91],

[104], [120]. As an example we consider the linear PDE (3.33) from § 3.7. We

use the same mesh, right-hand side, coe¬cients, and solution as the example

in § 3.7.

We compare Broyden™s method without restarts (solid line) to Broyden™s

method with restarts every three iterates (dashed line). For linear problems,

we allow for increases in the nonlinear residual. We set „a = „r = h2 where

h = 1/32. In Fig. 7.1 we present results for (3.33) preconditioned with the

Poisson solver fish2d.

Broyden™s method terminated after 9 iterations at a cost of roughly 1.5

million ¬‚oating point operations, a slightly higher cost than the GMRES

solution. The restarted iteration required 24 iterations and 3.5 million ¬‚oating

point operations. One can also see highly irregular performance from the

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.

128 ITERATIVE METHODS FOR LINEAR AND NONLINEAR EQUATIONS

restarted method.

In the linear case [67] one can store only the residuals and recover the

terminal iterate upon exit. Storage of the residual, the vector z, and the

sequence of steps saves a vector over the nonlinear case.

1

10

0

10

Relative Residual

-1

10

-2

10

-3

10

0 5 10 15 20 25

Iterations

Fig. 7.1. Broyden™s method for the linear PDE.

7.4.2. Nonlinear problems. We consider the H-equation from Chapters 5

and 6 and a nonlinear elliptic partial di¬erential equation. Broyden™s method

and Newton-GMRES both require more storage as iterations progress, Newton-

GMRES as the inner iteration progresses and Broyden™s method as the outer

iteration progresses. The cost in function evaluations for Broyden™s method is

one for each nonlinear iteration and for Newton-GMRES one for each outer

iteration and one for each inner iteration. Hence, if function evaluations are

very expensive, the cost of a solution by Broyden™s method could be much less

than that of one by Newton-GMRES. If storage is at a premium, however,

and the number of nonlinear iterations is large, Broyden™s method could be at

a disadvantage as it might need to be restarted many times. The examples

in this section, as well as those in § 8.4 illustrate the di¬erences in the two

methods. Our general recommendation is to try both.

Chandrasekhar H-equation. As before we solve the H-equation on a 100-

point mesh with c = .9 (Fig. 7.2) and c = .9999 (Fig. 7.3). We compare

Broyden™s method without restarts (solid line) to Broyden™s method with

restarts every three iterates (dashed line). For c = .9 both methods required six

iterations for convergence. The implementation without restarts required 153

thousand ¬‚oating-point operations and the restarted method 150, a substantial

improvement over the Newton-GMRES implementation.

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.

129

BROYDEN™S METHOD

For c = .9999, the restarted iteration has much more trouble. The un-

restarted iteration converges in 10 iterations at a cost of roughly 249 thousand

¬‚oating-point operations. The restarted iteration requires 18 iterations and

407 ¬‚oating-point operations. Even so, both implementations performed bet-

ter than Newton-GMRES.

0

10

-1

10

-2

10

Relative Nonlinear Residual

-3

10

-4

10

-5

10

-6

10

-7

10

-8

10

0 1 2 3 4 5 6

Iterations

Fig. 7.2. Broyden™s method for the H-equation, c = .9.

0

10

-1

10

Relative Nonlinear Residual

-2

10

-3

10

-4

10

-5

10

-6

10

-7

10

0 2 4 6 8 10 12 14 16 18

Iterations

Fig. 7.3. Broyden™s method for the H-equation, c = .9999.

One should take some care in drawing general conclusions. While Broyden™s

method on this example performed better than Newton-GMRES, the storage

in Newton-GMRES is used in the inner iteration, while that in Broyden in the

nonlinear iteration. If many nonlinear iterations are needed, it may be the case

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.

130 ITERATIVE METHODS FOR LINEAR AND NONLINEAR EQUATIONS

that Broyden™s method may need restarts while the GMRES iteration for the

linear problem for the Newton step may not.

Convection-di¬usion equation. We consider the nonlinear convection-

di¬usion equation from § 6.4.2. We use the same mesh width h = 1/32 on

a uniform square grid. We set C = 20, u0 = 0, and „r = „a = h2 . We

preconditioned with the Poisson solver fish2d.

We allowed for an increase in the residual, which happened on the second

iterate. This is not supported by the local convergence theory for nonlinear

problems; however it can be quite successful in practice. One can do this in

brsol by setting the third entry in the parameter list to 1, as one would if the

problem were truly linear. In § 8.4.3 we will return to this example.

In Fig. 7.4 we compare Broyden™s method without restarts (solid line) to

Broyden™s method with restarts every 8 iterates (dashed line). The unrestarted

Broyden iteration converges in 12 iterations at a cost of roughly 2.2 million

¬‚oating-point operations. This is a modest improvement over the Newton-

GMRES cost (reported in § 6.4.2) of 2.6 million ¬‚oating-point operations.

However, the residual norms do not monotonically decrease in the Broyden

iteration (we ¬x this in § 8.3.2) and more storage was used.

1

10

0

10

Relative Nonlinear Residual

-1

10

-2

10

-3

10

-4

10

0 2 4 6 8 10 12 14 16

Iterations

Fig. 7.4. Broyden™s method for nonlinear PDE.

At most ¬ve inner GMRES iterations were required. Therefore the GMRES

inner iteration for the Newton-GMRES approach needed storage for at most

10 vectors (the Krylov basis, s, xc , F (xc ), F (xc + hv)) to accommodate the

Newton-GMRES iteration. The Broyden iteration when restarted every n

iterates, requires storage of n + 2 vectors {sj }n’1 , x, and z (when stored in the

j=0

same place as F (x)). Hence restarting the Broyden iteration every 8 iterations

most closely corresponds to the Newton-GMRES requirements. This approach

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.

131

BROYDEN™S METHOD

took 15 iterations to terminate at a cost of 2.4 million ¬‚oating-point operations,

but the convergence was extremely irregular after the restart.

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.

132 ITERATIVE METHODS FOR LINEAR AND NONLINEAR EQUATIONS

7.5. Exercises on Broyden™s method

7.5.1. Prove Proposition 7.2.1.

7.5.2. Prove Proposition 7.3.1.

7.5.3. Extend Proposition 7.3.1 by showing that if A is a nonsingular N — N

matrix, U is N — M , V is M — N , then A + U V is nonsingular if and

only if I + V A’1 U is a nonsingular M — M matrix. If this is the case,

then

(A + U V )’1 = A’1 ’ A’1 U (I + V A’1 U )’1 V A’1 .

This is the Sherman“Morrison“Woodbury formula [68], [199]. See [102]

for further generalizations.

7.5.4. Duplicate the results reported in § 7.4. Vary the parameters in the

equations and the number of vectors stored by Broyden™s method and

report on their e¬ects on performance. What happens in the di¬erential

equation examples if preconditioning is omitted?

7.5.5. Solve the H-equation with Broyden™s method and c = 1. Set „a = „r =

10’8 . What q-factor for linear convergence do you see? Can you account

for this by applying the secant method to the equation x2 = 0? See [53]

for a complete explanation.

7.5.6. Try to solve the nonlinear convection-di¬usion equation in § 6.4 with

Broyden™s method using various values for the parameter C. How does

the performance of Broyden™s method di¬er from Newton-GMRES?

7.5.7. Modify nsol to use Broyden™s method instead of the chord method after

the initial Jacobian evaluation. How does this compare with nsol on the

examples in Chapter 5.

7.5.8. Compare the performance of Broyden™s method and Newton-GMRES on

the Modi¬ed Bratu Problem

’∇2 u + dux + eu = 0

on (0, 1) — (0, 1) with homogeneous Dirichlet boundary conditions.

Precondition with the Poisson solver fish2d. Experiment with various

mesh widths, initial iterates, and values for the convection coe¬cient d.

7.5.9. The bad Broyden method [26], so named because of its inferior perfor-

mance in practice [63], updates an approximation to the inverse of the

Jacobian at the root so that B ’1 ≈ F (x— )’1 satis¬es the inverse secant

equation

’1

(7.45) B+ y = s

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.

133

BROYDEN™S METHOD

with the rank-one update

(s ’ Bc y)y T

’1

’1 ’1

B+ = Bc + .

y2

2

Show that the bad Broyden method is locally superlinearly convergent

if the standard assumptions hold. Try to implement the bad Broyden

method in a storage-e¬cient manner using the ideas from [67] and § 7.3.

Does anything go wrong? See [67] for more about this.

7.5.10. The Schubert or sparse Broyden algorithm [172], [27] is a quasi-Newton

update for sparse matrices that enforces both the secant equation and the

sparsity pattern. For problems with tridiagonal Jacobian, for example,

the update is

(y ’ Bc s)i sj

(B+ )ij = (Bc )ij + i+1 2

k=i’1 sk

for 1 ¤ i, j ¤ N and |i ’ j| ¤ 1. Note that only the subdiagonal, su-

perdiagonal, and main diagonal are updated. Read about the algorithm

in [172] and prove local superlinear convergence under the standard as-

sumptions and the additional assumption that the sparsity pattern of B0

is the same as that of F (x— ). How would the Schubert algorithm be

a¬ected by preconditioning? Compare your analysis with those in [158],

[125], and [189].

7.5.11. Implement the bad Broyden method, apply it to the examples in § 7.4,

and compare the performance to that of Broyden™s method. Discuss the

di¬erences in implementation from Broyden™s method.

7.5.12. Implement the Schubert algorithm on an unpreconditioned discretized

partial di¬erential equation (such as the Bratu problem from Exer-

cise 7.5.8) and compare it to Newton-GMRES and Broyden™s method.

Does the relative performance of the methods change as h is decreased?

This is interesting even in one space dimension where all matrices are

tridiagonal. References [101], [93], [92], and [112], are relevant to this

exercise.

7.5.13. Use your result from Exercise 5.7.14 in Chapter 5 to numerically estimate

the q-order of Broyden™s method for some of the examples in this

section (both linear and nonlinear). What can you conclude from your

observations?

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.

134 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 8

Global Convergence

By a globally convergent algorithm we mean an algorithm with the property

that for any initial iterate the iteration either converges to a root of F or fails

to do so in one of a small number of ways. Of the many such algorithms we

will focus on the class of line search methods and the Armijo rule [3], [88],

implemented inexactly [24], [25], [70].

Other methods, such as trust region [153], [154], [181], [129], [63], [24],

[25], [173], [70], and continuation/homotopy methods [1], [2], [109], [159],

[198], can be used to accomplish the same objective. We select the line search

paradigm because of its simplicity, and because it is trivial to add a line search

to an existing locally convergent implementation of Newton™s method. The

implementation and analysis of the line search method we develop in this

chapter do not depend on whether iterative or direct methods are used to

compute the Newton step or upon how or if the Jacobian is computed and

stored. We begin with single equations, where the algorithms can be motivated

and intuition developed with a very simple example.

8.1. Single equations

If we apply Newton™s method to ¬nd the root x— = 0 of the function

f (x) = arctan(x) with initial iterate x0 = 10 we ¬nd that the initial iterate

is too far from the root for the local convergence theory in Chapter 5 to be

applicable. The reason for this is that f (x) = (1 + x2 )’1 is small for large

x; f (x) = arctan(x) ≈ ±π/2, and hence the magnitude of the Newton step

is much larger than that of the iterate itself. We can see this e¬ect in the

sequence of iterates:

10, ’138, 2.9 — 104 , ’1.5 — 109 , 9.9 — 1017 ,

a failure of Newton™s method that results from an inaccurate initial iterate.

If we look more closely, we see that the Newton step s = ’f (xc )/f (xc ) is

pointing toward the root in the sense that sxc < 0, but the length of s is too

large. This observation motivates the simple ¬x of reducing the size of the step

until the size of the nonlinear residual is decreased. A prototype algorithm is

given below.

135

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.

136 ITERATIVE METHODS FOR LINEAR AND NONLINEAR EQUATIONS

Algorithm 8.1.1. lines1(x, f, „ )

1. r0 = |f (x)|

2. Do while |f (x)| > „r r0 + „a

(a) If f (x) = 0 terminate with failure.

(b) s = ’f (x)/f (x) (search direction)

(c) xt = x + s (trial point)

(d) If |f (xt )| < |f (x)| then x = xt (accept the step)

else

s = s/2 goto 2c (reject the step)

When we apply Algorithm lines1 to our simple example, we obtain the

sequence

10, ’8.5, 4.9, ’3.8, 1.4, ’1.3, 1.2, ’.99, .56, ’0.1, 9 — 10’4 , ’6 — 10’10 .

The quadratic convergence behavior of Newton™s method is only apparent

very late in the iteration. Iterates 1, 2, 3, and 4 required 3, 3, 2, and 2

steplength reductions. After the fourth iterate, decrease in the size of the

nonlinear residual was obtained with a full Newton step. This is typical of the

performance of the algorithms we discuss in this chapter.

We plot the progress of the iteration in Fig. 8.1. We plot

|arctan(x)/arctan(x0 )|

as a function of the number of function evaluations with the solid line. The

outer iterations are identi¬ed with circles as in § 6.4. One can see that most of

the work is in the initial phase of the iteration. In the terminal phase, where

the local theory is valid, the minimum number of two function evaluations per

iterate is required. However, even in this terminal phase rapid convergence

takes place only in the ¬nal three outer iterations.

Note the di¬erences between Algorithm lines1 and Algorithm nsol. One

does not set x+ = x + s without testing the step to see if the absolute value of

the nonlinear residual has been decreased. We call this method a line search

because we search along the line segment

(xc , xc ’ f (xc )/f (xc ))

to ¬nd a decrease in |f (xc )|. As we move from the right endpoint to the left,

some authors [63] refer to this as backtracking.

Algorithm lines1 almost always works well. In theory, however, there is

the possibility of oscillation about a solution without convergence [63]. To

remedy this and to make analysis possible we replace the test for simple

decrease with one for su¬cient decrease. The algorithm is to compute a search

direction d which for us will be the Newton direction

d = ’f (xc )/f (xc )

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.

137

GLOBAL CONVERGENCE

0

10

-1

10

-2

10

Relative Nonlinear Residual -3

10

-4

10

-5

10

-6

10

-7

10

-8

10

-9

10

-10

10

0 5 10 15 20 25 30 35

Function Evaluations

Fig. 8.1. Newton“Armijo iteration for inverse tangent.

and then test steps of the form s = »d, with » = 2’j for some j ≥ 0, until

f (x + s) satis¬es

(8.1) |f (xc + »d)| < (1 ’ ±»)|f (xc )|.

The condition in (8.1) is called su¬cient decrease of |f |. The parameter

± ∈ (0, 1) is a small, but positive, number intended to make (8.1) as easy

as possible to satisfy. We follow the recent optimization literature [63] and

set ± = 10’4 . Once su¬cient decrease has been obtained we accept the step

s = »d. This strategy, from [3] (see also [88]) is called the Armijo rule. Note

that we must now distinguish between the step and the Newton direction,

something we did not have to do in the earlier chapters.

Algorithm 8.1.2. nsola1(x, f, „ )

1. r0 = |f (x)|

2. Do while |f (x)| > „r r0 + „a

(a) If f (x) = 0 terminate with failure.

(b) d = ’f (x)/f (x) (search direction)

(c) » = 1

i. xt = x + »d (trial point)

ii. If |f (xt )| < (1 ’ ±»)|f (x)| then x = xt (accept the step)

else

» = »/2 goto 2(c)i (reject the step)

This is essentially the algorithm implemented in the MATLAB code nsola.

We will analyze Algorithm nsola in § 8.2. We remark here that there is

nothing critical about the reduction factor of 2 in the line search. A factor of

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.

138 ITERATIVE METHODS FOR LINEAR AND NONLINEAR EQUATIONS

10 could well be better in situations in which small values of » are needed for

several consecutive steps (such as our example with the arctan function). In

that event a reduction of a factor of 8 would require three passes through the

loop in nsola1 but only a single reduction by a factor of 10. This could be

quite important if the search direction were determined by an inexact Newton

method or an approximation to Newton™s method such as the chord method

or Broyden™s method.

On the other hand, reducing » by too much can be costly as well. Taking

full Newton steps ensures fast local convergence. Taking as large a fraction as

possible helps move the iteration into the terminal phase in which full steps

may be taken and fast convergence expected. We will return to the issue of

reduction in » in § 8.3.

8.2. Analysis of the Armijo rule

In this section we extend the one-dimensional algorithm in two ways before

proceeding with the analysis. First, we accept any direction that satis¬es the

inexact Newton condition (6.1). We write this for the Armijo sequence as

(8.2) F (xn )dn + F (xn ) ¤ ·n F (xn ) .

Our second extension is to allow more ¬‚exibility in the reduction of ». We

allow for any choice that produces a reduction that satis¬es

σ0 »old ¤ »new ¤ σ1 »old ,

where 0 < σ0 < σ1 < 1. One such method, developed in § 8.3, is to minimize

an interpolating polynomial based on the previous trial steps. The danger is

that the minimum may be too near zero to be of much use and, in fact, the

iteration may stagnate as a result. The parameter σ0 safeguards against that.

Safeguarding is an important aspect of many globally convergent methods and

has been used for many years in the context of the polynomial interpolation

algorithms we discuss in § 8.3.1 [54], [63], [86], [133], [87]. Typical values of σ0

and σ1 are .1 and .5. Our algorithm nsola incorporates these ideas.

While (8.2) is motivated by the Newton-iterative paradigm, the analysis

here applies equally to methods that solve the equations for the Newton step

directly (so ·n = 0), approximate the Newton step by a quasi-Newton method,

or even use the chord method provided that the resulting direction dn satis¬es

(8.2). In Exercise 8.5.9 you are asked to implement the Armijo rule in the

context of Algorithm nsol.

Algorithm 8.2.1. nsola(x, F, „, ·).

1. r0 = F (x)

2. Do while F (x) > „r r0 + „a

(a) Find d such that F (x)d + F (x) ¤ · F (x)

If no such d can be found, terminate with failure.

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.

139

GLOBAL CONVERGENCE

(b) » = 1

i. xt = x + »d

ii. If F (xt ) < (1 ’ ±») F (x) then x = xt

else

Choose σ ∈ [σ0 , σ1 ]

» = σ»

goto 2(b)i

Note that step 2a must allow for the possibility that F is ill-conditioned

and that no direction can be found that satis¬es (8.2). If, for example, step

2a were implemented with a direct factorization method, ill-conditioning of F

might be detected as part of the factorization.

Let {xn } be the iteration produced by Algorithm nsola with initial iterate

x0 . The algorithm can fail in some obvious ways:

1. F (xn ) is singular for some n. The inner iteration could terminate with

a failure.

2. xn ’ x, a local minimum of F which is not a root.

¯

3. xn ’ ∞.

Clearly if F has no roots, the algorithm must fail. We will see that if F has no

roots then either F (xn ) has a limit point which is singular or {xn } becomes

unbounded.

The convergence theorem is the remarkable statement that if {xn } and

{ F (xn )’1 } are bounded (thereby eliminating premature failure and local

minima of F that are not roots) and F is Lipschitz continuous in a

neighborhood of the entire sequence {xn }, then the iteration converges to a

root of F at which the standard assumptions hold, full steps are taken in the

terminal phase, and the convergence is q-quadratic.

We begin with a formal statement of our assumption that F is uniformly

Lipschitz continuous and bounded away from zero on the sequence {xn }.

Assumption 8.2.1. There are r, γ, mf > 0 such that F is de¬ned, F is

Lipschitz continuous with Lipschitz constant γ, and F (x)’1 ¤ mf on the

set

∞

„¦({xn }, r) = {x | x ’ xn ¤ r}.

n=0

Assumption 8.2.1 implies that the line search phase (step 2b in Algorithm

nsola) will terminate in a number of cycles that is independent of n. We show

this, as do [25] and [70], by giving a uniform lower bound on »n (assuming

that F (x0 ) = 0). Note how the safeguard bound σ0 enters this result. Note

also the very modest conditions on ·n .

Lemma 8.2.1. Let x0 ∈ RN and ± ∈ (0, 1) be given. Assume that {xn } is

given by Algorithm nsola with

(8.3) {·n } ‚ (0, · ] ‚ (0, 1 ’ ±)

¯

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.

140 ITERATIVE METHODS FOR LINEAR AND NONLINEAR EQUATIONS

and that Assumption 8.2.1 holds. Then either F (x0 ) = 0 or

r 2(1 ’ ± ’ · )

¯

¯

(8.4) »n ≥ » = σ0 min ,2 .

mf F (x0 ) mf F (x0 ) (1 + · )2 γ

¯

Proof. Let n ≥ 0. By the fundamental theorem of calculus and (8.2) we

have, for any » ∈ [0, 1]

1

F (xn + »dn ) = F (xn ) + »F (xn )dn + (F (xn + t»dn ) ’ F (xn ))»dn dt

0

1

= (1 ’ »)F (xn ) + »ξn + (F (xn + t»dn ) ’ F (xn ))»dn dt,

0

where, by the bound ·n ¤ · ,

¯

ξn ¤ ·n F (xn ) ¤ · F (xn ) .

¯

The step acceptance condition in nsola implies that { F (xn ) } is a decreasing

sequence and therefore

dn = F (xn )’1 (ξn ’ F (xn )) ¤ mf (1 + · ) F (xn ) .

¯

Therefore, by Assumption 8.2.1 F is Lipschitz continuous on the line segment

[xn , xn + »dn ] whenever

¯

» ¤ »1 = r/(mf F (x0 ) ).

¯

Lipschitz continuity of F implies that if » ¤ »1 then

γ2 2

F (xn + »dn ) ¤ (1 ’ ») F (xn ) + »¯ F (xn ) +

· » dn

2

m2 γ(1 + · )2 »2 F (xn ) 2

¯

f

¤ (1 ’ ») F (xn )| + »¯ F (xn ) +

·

2

m2 γ»(1 + · )2 F (x0 )

¯

f

¤ (1 ’ » + · ») F (xn ) + » F (xn )

¯

2

< (1 ’ ±») F (xn ) ,

which will be implied by

2(1 ’ ± ’ · )

¯

¯ ¯

» ¤ »2 = min »1 , .

(1 + · )2 m2 γ F (x0 )

¯ f

¯

This shows that » can be no smaller than σ0 »2 , which completes the proof.

Now we can show directly that the sequence of Armijo iterates either is

unbounded or converges to a solution. Theorem 8.2.1 says even more. The

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.

141

GLOBAL CONVERGENCE

sequence converges to a root at which the standard assumptions hold, so in the

terminal phase of the iteration the step lengths are all 1 and the convergence

is q-quadratic.

Theorem 8.2.1. Let x0 ∈ RN and ± ∈ (0, 1) be given. Assume that {xn }

is given by Algorithm nsola, {·n } satis¬es (8.3), {xn } is bounded, and that

Assumption 8.2.1 holds. Then {xn } converges to a root x— of F at which the

standard assumptions hold, full steps are taken for n su¬ciently large, and the

convergence behavior in the ¬nal phase of the iteration is that given by the local

theory for inexact Newton methods (Theorem 6.1.2).

Proof. If F (xn ) = 0 for some n, then we are done because Assump-

tion 8.2.1 implies that the standard assumptions hold at x— = xn . Otherwise

Lemma 8.2.1 implies that F (xn ) converges q-linearly to zero with q-factor at

¯

most (1 ’ ±»).

The boundedness of the sequence {xn } implies that a subsequence {xnk }

converges, say to x— . Since F is continuous, F (x— ) = 0. Eventually

|xnk ’ x— | < r, where r is the radius in Assumption 8.2.1 and therefore the

standard assumptions hold at x— .

Since the standard assumptions hold at x— , there is δ such that if x ∈ B(δ),

the Newton iteration with x0 = x will remain in B(δ) and converge q-

quadratically to x— . Hence as soon as xnk ∈ B(δ), the entire sequence, not

just the subsequence, will remain in B(δ) and converge to x— . Moreover,

Theorem 6.1.2 applies and hence full steps can be taken.

At this stage we have shown that if the Armijo iteration fails to converge

to a root either the continuity properties of F or nonsingularity of F break

down as the iteration progresses (Assumption 8.2.1 is violated) or the iteration

becomes unbounded. This is all that one could ask for in terms of robustness

of such an algorithm.

Thus far we have used for d the inexact Newton direction. It could well be

advantageous to use a chord or Broyden direction, for example. All that one

needs to make the theorems and proofs in § 8.2 hold is (8.2), which is easy to

verify, in principle, via a forward di¬erence approximation to F (xn )d.

8.3. Implementation of the Armijo rule

The MATLAB code nsola is a modi¬cation of nsolgm which provides for a

choice of several Krylov methods for computing an inexact Newton direction

and globalizes Newton™s method with the Armijo rule. We use the l2 norm in

that code and in the numerical results reported in § 8.4. If GMRES is used

as the linear solver, then the storage requirements are roughly the same as for

Algorithm nsolgm if one reuses the storage location for the ¬nite di¬erence

directional derivative to test for su¬cient decrease.

In Exercise 8.5.9 you are asked to do this with nsol and, with each

iteration, numerically determine if the chord direction satis¬es the inexact

Newton condition. In this case the storage requirements are roughly the same

as for Algorithm nsol. The iteration can be very costly if the Jacobian must be

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.

142 ITERATIVE METHODS FOR LINEAR AND NONLINEAR EQUATIONS

evaluated and factored many times during the phase of the iteration in which

the approximate solution is far from the root.

We consider two topics in this section. The ¬rst, polynomial interpolation

as a way to choose », is applicable to all methods. The second, how Broyden™s

method must be implemented, is again based on [67] and is more specialized.

Our Broyden“Armijo code does not explicitly check the inexact Newton

condition and thereby saves an evaluation of F . Evaluation of F (xc )dc if the

full step is not accepted (here dc is the Broyden direction) would not only verify

the inexact Newton condition but provide enough information for a two-point

parabolic line search. The reader is asked to pursue this in Exercise 8.5.10.

8.3.1. Polynomial line searches. In this section we consider an elabora-

tion of the simple line search scheme of reduction of the steplength by a ¬xed

factor. The motivation for this is that some problems respond well to one or

two reductions in the steplength by modest amounts (such as 1/2) and others

require many such reductions, but might respond well to a more aggressive

steplength reduction (by factors of 1/10, say). If one can model the decrease

in F 2 as the steplength is reduced, one might expect to be able to better

estimate the appropriate reduction factor. In practice such methods usually

perform better than constant reduction factors.

If we have rejected k steps we have in hand the values

F (xn ) 2 , F (xn + »1 dn ) 2 , . . . F (xn + »k’1 dn ) 2 .

We can use this iteration history to model the scalar function

2

f (») = F (xn + »dn ) 2

with a polynomial and use the minimum of that polynomial as the next

steplength. We consider two ways of doing this that use second degree

polynomials which we compute using previously computed information.

After »c has been rejected and a model polynomial computed, we compute

the minimum »t of that polynomial analytically and set

±

σ0 »c if »t < σ0 »c ,

σ1 »c if »t > σ1 »c ,

(8.5) »+ =

»t otherwise.

Two-point parabolic model. In this approach we use values of f and f

at » = 0 and the value of f at the current value of » to construct a 2nd degree

interpolating polynomial for f .

f (0) = F (xn ) 2 is known. Clearly f (0) can be computed as

2

f (0) = 2(F (xn )T dn )T F (xn ) = 2F (xn )T (F (xn )dn ).

(8.6)

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.

143

GLOBAL CONVERGENCE

Use of (8.6) requires evaluation of F (xn )dn , which can be obtained by

examination of the ¬nal residual from GMRES or by expending an additional

function evaluation to compute F (xn )dn with a forward di¬erence. In any

case, our approximation to f (0) should be negative. If it is not, it may

be necessary to compute a new search direction. This is a possibility with

directions other than inexact Newton directions, such as the Broyden direction.

The polynomial p(») such that p, p agree with f, f at 0 and p(»c ) = f (»c )

is

p(») = f (0) + f (0)» + c»2 ,

where

f (»c ) ’ f (0) ’ f (0)»c

c= .

»2

c

Our approximation of f (0) < 0, so if f (»c ) > f (0), then c > 0 and p(») has a

minimum at

»t = ’f (0)/(2c) > 0.

We then compute »+ with (8.5).

Three-point parabolic model. An alternative to the two-point model that

avoids the need to approximate f (0) is a three-point model, which uses f (0)

and the two most recently rejected steps to create the parabola. The MATLAB

code parab3p implements the three-point parabolic line search and is called

by the two codes nsola and brsola.

In this approach one evaluates f (0) and f (1) as before. If the full step is

rejected, we set » = σ1 and try again. After the second step is rejected, we

have the values

f (0), f (»c ), and f (»’ ),

where »c and »’ are the most recently rejected values of ». The polynomial

that interpolates f at 0, »c , »’ is

» (» ’ »’ )(f (»c ) ’ f (0)) (»c ’ »)(f (»’ ) ’ f (0))

p(») = f (0) + + .

»c ’ »’ »c »’

We must consider two situations. If

2

p (0) = (»’ (f (»c ) ’ f (0)) ’ »c (f (»’ ) ’ f (0)))

»c »’ (»c ’ »’ )

is positive, then we set »t to the minimum of p

»t = ’p (0)/p (0)

and apply the safeguarding step (8.5) to compute »+ . If p (0) ¤ 0 one could

either set »+ to be the minimum of p on the interval [σ0 », σ1 »] or reject the

parabolic model and simply set »+ to σ0 »c or σ1 »c . We take the latter approach

and set »+ = σ1 »c . In the MATLAB code nsola from the collection, we

implement this method with σ0 = .1, σ1 = .5.

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.

144 ITERATIVE METHODS FOR LINEAR AND NONLINEAR EQUATIONS

Interpolation with a cubic polynomial has also been suggested [87], [86],

[63], [133]. Clearly, the polynomial modeling approach is heuristic. The

safeguarding and Theorem 8.2.1 provide insurance that progress will be made

in reducing F .

8.3.2. Broyden™s method. In Chapter 7 we implemented Broyden™s

method at a storage cost of one vector for each iteration. The relation

y ’ Bc s = F (x+ ) was not critical to this and we may also incorporate a line

search at a cost of a bit of complexity in the implementation. As in § 7.3, our

approach is a nonlinear version of the method in [67].

The di¬erence from the development in § 7.3 is that the simple relation

between the sequence {wn } and {vn } is changed. If a line search is used, then

’1

sn = ’»n Bn F (xn )

and hence

(8.7) yn ’ Bn sn = F (xn+1 ) ’ (1 ’ »n )F (xn ).

If, using the notation in § 7.3, we set

yn ’ Bn sn sn ’1 T ’1

un = , vn = , and wn = (Bn un )/(1 + vn Bn un ),

sn 2 sn 2

we can use (8.7) and (7.38) to obtain

’dn+1 sn

+ (»’1 ’ 1)

wn = »n n

sn 2 sn 2

(8.8)

’sn+1 /»n+1 sn

+ (»’1 ’ 1)

= »n ,

n

sn 2 sn 2

where

’1

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

is the search direction used to compute xn+1 . Note that (8.8) reduces to (7.41)

when »n = 1 and »n+1 = 1 (so dn+1 = sn+1 ).

We can use the ¬rst equality in (8.8) and the relation

wn sT

n ’1

dn+1 =’ I’ Bn F (xn+1 )

sn 2

to solve for dn+1 and obtain an analog of (7.44)

sn 2 Bn F (xn+1 ) ’ (1 ’ »n )sT Bn F (xn+1 )sn

’1 ’1

n

2

(8.9) dn+1 =’ ’1

sn 2 + »n sT Bn F (xn+1 )

n

2

’1

We then use the second equality in (8.8) to construct Bc F (x+ ) as we did in

Algorithm brsol. Veri¬cation of (8.8) and (8.9) is left to Exercise 8.5.5.

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.

145

GLOBAL CONVERGENCE

Algorthm brsola uses these ideas. We do not include details of the line

search or provide for restarts or a limited memory implementation. The

MATLAB implementation, provided in the collection of MATLAB codes, uses

a three-point parabolic line search and performs a restart when storage is

exhausted like Algorithm brsol does.

It is possible that the Broyden direction fails to satisfy (8.2). In fact, there

is no guarantee that the Broyden direction is a direction of decrease for F 2 .

The MATLAB code returns with an error message if more than a ¬xed number

of reductions in the step length are needed. Another approach would be to

compute F (xc )d numerically if a full step is rejected and then test (8.2) before

beginning the line search. As a side e¬ect, an approximation of f (0) can be

obtained at no additional cost and a parabolic line search based on f (0), f (0),

and f (1) can be used. In Exercise 8.5.10 the reader is asked to fully develop

and implement these ideas.

Algorithm 8.3.1. brsola(x, F, „, maxit, nmax)

1. Initialize:

Fc = F (x) r0 = F (x) 2 , n = ’1,

d = ’F (x), compute »0 , s0 = »0 d

2. Do while n ¤ maxit

(a) n = n + 1

(b) x = x + sn

(c) Evaluate F (x)

(d) If F (x) ¤ „r r0 + „a exit.

2

(e) i. z = ’F (x)

ii. for j = 0, n ’ 1

a = ’»j /»j+1 , b = 1 ’ »j

z = z + (asj+1 + bsj )sT z/ sj 2

2

j

(f) d = ( sn 2 z + (1 ’ »n )sT zsn )/( sn 2 ’ »n sT z)

n n

2 2

(g) Compute »n+1 with a line search.

(h) Set sn+1 = »n+1 d, store sn+1 and »n+1 .

2

Since dn+1 and sn+1 can occupy the same location, the storage require-

ments of Algorithm brsola would appear to be essentially the same as those

of brsol. However, computation of the step length requires storage of both the

current point x and the trial point x + »s before a step can be accepted. So,

the storage of the globalized methods exceeds that of the local method by one

vector. The MATLAB implementation of brsola in the collection illustrates

this point.

8.4. Examples for Newton“Armijo

The MATLAB code nsola is a modi¬cation of nsolgm that incorporates the

three-point parabolic line search from § 8.3.1 and also changes · using (6.20)

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.

146 ITERATIVE METHODS FOR LINEAR AND NONLINEAR EQUATIONS

once the iteration is near the root. We compare this strategy, using GMRES

as the linear solver, with the constant · method as we did in § 6.4 with γ = .9.

As before, in all the ¬gures we plot F (xn ) 2 / F (x0 ) 2 against the number

of function evaluations required by all line searches and inner and outer

iterations to compute xn . Counts of function evaluations corresponding to

outer iterations are√ indicated by circles. We base our absolute error criterion

on the norm · 2 / N as we did in § 6.4.2 and 7.4.

In § 8.4.3 we compare the Broyden“Armijo method with Newton-GMRES

for those problems for which both methods were successful. We caution the

reader now, and will repeat the caution later, that if the initial iterate is far

from the solution, an inexact Newton method such as Newton-GMRES can

succeed in many case in which a quasi-Newton method can fail because the