<<

. 6
( 7)



>>

and the proof is complete.

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

<<

. 6
( 7)



>>