EXPLOITING SPARSITY
IN NEWTON-LIKE METHODS
by
Earl Stanley Marwil
TR78--H335
EXPLOITING SPARSITY IN NEWTON-LIKE METHODS
A Thesis
Presented to the Faculty of the Graduate School
of Cornell University
in Partial Fulfillinent for the Degree of
Doctor of Philosophy
by
Earl Stanley Marwil
January 1978
Biographical Sketch
The author was born in lartleaville, Oklahoma on July 31, 1950.
He spent most of his youth there, attended the public schools, and
graduated from Sooner Nigh School in Bartlesville in Nay 1968. In
May 1971. he graduated from the University of Oklahoma with high
honors, having earned a Bachelor of Science degree in mathewatica.
Ne then began graduate study in the Center for Applied Nathematics
at Cornell University. In Nay 1974, he was awarded a Master of
ScIence, and he completed his Doctor of Philosophy at Cornell in
January 1978.
--H Ii --H
To my friend Denise
iii			--H
Acknowledge?nents
The author is very grateful to Cornell University for the financial
support it has given and for the use of its facilities. My most sincere
thanks go to my parents for their love and generosity during?my graduate
study. I would like to thank the members of my special committee for
their patience: Larry fayne, minor member, Al Schat:, who served as
chairman of the special committee, and John Dennis. who advised me on
this thesis. I am very thankful to John for the valuahle comments,
for the time, and for the encouragement and understanding. which he haa
given so freely. Also, a special thanks to Jorge More', who served on
the special committee and advised me on parts of this thesis when he
was at Cornell. finally I wish to thank Denise, to whom this work is
dedicated, for "the eggs".
--H iv -
Table of Contents
Bin anhical sketch			--H
ir.
ii
i'edicatton--H--H			----- iii
Ac?now1ed8ements--H------H--H			tv
Table of contents--H--H			v
Chapter 1 --H introduction--H--H--H--H----H-			I
1.1			General			I
1.2			?xamples			4
1.3			Literature--H--H
1.4			Organization--H-			8
1.5			Frequently quoted theorerns			9
Chapter 2 --H Local and Superlinear Convergence of			Srhuhcr('a Nethod			II
2.1			IreIiminarles			II
2.2			Schubert's update						12
2?3			Local and superlinear convergence --H----			15
2.4			A ?Cantorovich analysis			21
Chapter 3 --H Sysmetric Schubert update 33
3.1			IreIi?inarlas--H			33
3.2			Notation and technical lemmas ----			34
3.3			Symmetric Schubert update			42
3.4			Properties of the update			43
3.3			An algorithm--H------H--			49
3.6			Convergence						31
Chapter			4 --H Natrix Factorizations ----------H--			33
4.1			Basic considerations--H--H---			33
 --------H--						36
4.2 Some algorithms
?.3 Continuity of the
factorization --H--H
4.4			Convergence--H----H
Cbapter 5 --H Further Applications and Conclusion
1.1			Further appIications --H--H--H--H--?--H
5.2			Conclusion--H
fliblfography--H
`9
66
26
76
27
81
--H vi --H
1. Introduction
1.1 General.
Nonlinear problems occur in physics, engineerins, economics,
ecology. and many other fields. Let F: ??. 1i? be a nonlinear
function on fl--Hdimensional real linear space. the problem of solving
a system of nonlinear equations is to find a solution to the equa-
tion
(1.1.1)			f(x)			0.
Numerical analysts seek solutions by iterative techniques.
The classical Newton's method for solving (1.1.1) is a linear--H
ization of the Taylor expansion of F. Given a starting point or
initial approximation, x0 to the solution, the iteration can be
defined by
(1.1.2)			xk+i ??5k --H F?<xk)?1F(xk)
Theorem 1.1.1 [27] If there is a point xSD0 an open, convex set
such that F(x*) 0, F is differentiabie in D0, F' is Lipschttz
continuous in D0, and F'(xt) is nonsingular, then the Newton iter--H
stion, (1.1.2), converges quadratically to x*, a solution of (1.1.1)
for any x0 cD0.
Quadratic convergence means for a positive constant a, the
sequence (x?l behaves as
(1.1.3)			115k+l --H x*li <a ixk --H x*l 12.
A Newton--Hlike method replaces the Jacobian, F' , in (1.1.2)
--H1--H
--H2
by an invertible operator. Given x0  IRn and a sequence of invert--H
ible operators (?l,
(1.1.4)
Xk+l = xk - Ak1F(xk)
defines an iterative process which approximates Newton?s method. The
invertible operator Ak is meant to approximate the Jacobian at xk
F?(xk)
There are numerous articles that have been written about iter--H
attons of the form (1.1.4) and other generalizations of Newton's
method. See [33? for an excellent bibliography including those topics.
Recently a subciasa of iterations of the form (1.1.4) has received a
lot of attention. This was begun by flavidon [7], lietcher and Powell
[18], and Rroyden [2]. The idea is to update the approximate Jacoblan at
each iteration to obtain the next approximate Jacobian.
Given x0 and B0 an initial approximation to the Jacobtan
at define ?xt'8k1 iteratfely by
(1.1.5)			xk+l = Xk --H Bk11(xk)
and
(1.1.6) llk+1 = a function of Bk xk+i?xk and
This class of methods are frequently called quasi--HNewton methods. The
reason for thia is the requirement that llk+l satisfy the matrix--Hvec-
tor equation:
(1.1.7)
Bk+l(x?+i --H x?) = F(x?+1) -
--H3--H
This equation is called the quasi--HNewton equation. N'.?t?'n?s method is
based on the truncated Taylor expansion:
F?(x?+i ??5k+i --H xk) = i(x?+i) --H 1?5k?
For Xk+l close to the solution, F(x ) 0, and i??5k? ii place ci
k+l
defines Newton's method. Equation (1.1.7) uses lk.i in place
ci E'(xk+l) and in one dimension It reduces to tbe secant equation. The
quasi-Newton methods are also called secant methods.
The quasi--HNewton methods exhibit good local convergence proper--H
ties [4], [15] and behave well in practice for a large variety of
problems [is]. The type of convergence usually exhiblied is super--H
linear. That is, for (x? I
(1.1.8)			lIXk+i --H x*li < mkllxk --H x*ii
with 0k 0. Many authors have discussed various techniques to en--H
large tbe region where an initial guess will lead to convergence.
Problems of the form (1.1.1) which have a special structure will
be considered here. Intuitively. a system ci nonlineat equstlons
is sparse if each equation depends on only a iew variables In the
domain apace. A linear problem is considered to be sparse when the
matrix has fewer than lOt noneero entries [16].
The quasi--HNewton methods are usually implemented replacing (1.1.3)
by:
(1.1.0)
(1.1.10)
Solve 8k5k = ?1(x?) for
then set 1k?i = 5k + 5k'
-4--H
This is done for numerical stability.
Since there is a linear system to be solved at each iteration,
it is very natural to define sparse nonlinear equations as those
which have fewer than lOZ nonzeroes in the Jacobian matrix. In the
`ethods discussed here linear systems with the same zero--Hnonzero
structure or sparseness structure occur at each iteration.
1.2 Examples.
Specific types-?f nonlinear problems are abundant. They include
power flow systems [3], integrating stiff ordinary differential e
quations [5], nonlinear electrical networks [5], flow networks in
physiolo8ical models, finite analog of nonlinear boundary value
problems [37], certain econometric models, and fixed point problems [34].
Esample 1.2.1 A t?o point boundary value problem, Rheinboldt [34].
Consider the differential equation
(1.2.1)			--H d			du 2
= f(t) for tc (0,1)
and where
(1.2.2)
u(O) = u(l) = 0,
c(a) = I + a , and
f(t) = 600(1 --H 2t)2 + 20.
flheinboldt notes that this problem exhibits properties typical of
elliptical boundary value problems in higher dimensions. The exact
solution to (1.2.1) is
u(t) = lOt(l--Ht).
--H5--H
On a uniform mesh tj jh. j 0,1n+l, h
a standard approximation of (1.2.1) is of the for.
(1.2.3)			A(x) = b,			?? 1n
where			A(x) = (a???5??n?n' wiLh
--H ??1)1h			for J			1,
(l.2?4)			a (x)			e?!h			for i--Hil = 1, k			o+in(i,J),
0 otherwise,
for f,j			= 1,2n, and where
x			=x			-O
0			n+1 -
and mi = c[xi?i --H x1)!h)
Here, xi is the ith component of the vector x.
Example 1.2.2 A problem of heated networks, Porachina (28). which
generarca a sparse nonlinear system.
Althou8h in the context that is considered in [28] the end result
is in the form of the fixed point iteration, the underlying nonlinear
system is sparse. The eq?ations governing the network flow are:
(1.2.5)			Aw = b
B(w)h = q
w = DATp +
The following definitions identify the matrices and vectors occurring
in (1.2.5).
A? ?n ?m is the incidence matrix of positive flow in arcs k at nodes i.
depending on direction;
aik			ot connected.
wci is the vector of unknown flow rates.
be 1tn is a known vector cf the boundary flow rates.
n+l
B(w) r 1n+l			o+ with
b (w) -			IwkI			if arc k connects nodes I and j and flow is toward I,
ii --H
for I 6 1. and			hii			--H ? bij
j6i
h e 1tn+l is the vector of unknown eothalpy (a measure of thermal energy)
q is a known vector of the thermal energies of the boundary flows.
p r IL is a vector of pressures.
T
r<h)			(..., ckP((h --Hh )12),...)
ii
with V as the fluid density.
D = diag(d?)
The constants ck dk are determined by the Laminar Flow Law, which is
an equation of the form
wk = dk(pi --H Pj) + cko((h +h )f2) k = 1,2
I j
with the flow in arc k from node i to node j.
The nonlinear system (1.2.5) is sparse if A is sparse.
The preceding examples give an indication of the problems which
lead to a system of sparse nonlinear equations. In example 1.2.1,
when the mesh size becomes small, the dimension of the nonlinear differ-
ence equations is large. In many large scale problems the system is
--H7
likely to be sparse. The choice of the difference formula and the
ordering of the nodes of the mesh in higher order problems can be made
so that each equation involves only a few variables.
1.3 Literature.
Tewarson and Stephenson ?7? state some considerations for design
and selection of methods to solve large systems of nonlinear equations.
These considerations include storage requirements, rate of convergen-
ce. the domain of attraction and choosing the initial guess, computa--H
tional instability and large programming costs. Algorithms are de-
veloped which require only a few of the equations at a time. These
are based on a discretired Newton method and have as a particular-
ly useful application flow network problems.
porsching [28? presents two source problems, tncluding Example
1.2.2, and discusses some methods for solution of sparse nonlinear
equations. For the network flow problem, a continuation method is
applied to the system posed as a f?sed point problem. tie also gives
an elliptic boundary value problem arising in nonlinear diffusion
and discusses solution by accelerating the Gauss--HSeidel method by
overrelaxation.
In lievorkian and Snoek 2? large scale steady state design pro-
blems are posed as a system of nonlinear equations. Graph theoretic
arguments are used to develop a tearing procedure. ileid ?2] gives
FORTRAN subroutines, adapting the Levenbera--HMarquardt aigorithm for
least squares, to solve nonlinear equations. Rheinboldt [343 dis-
cusses the underlying reasons for some of the difficulties encounter-
ed in large dimensional problems. Another aspect of that paper is
--H8--H
the discussion of some techniques for reducing the dimension of tbe
nonlinear system.
The idea of adapting quasi--HNewton methods to Lake advantage of
sparsity was first published in 1970 by Schubert [36]. Ne modifled
Broyden's method [2? to preserve sparsity in the a?proximate Jacobian.
Broyden [3] later gave a local and linear convergence result for
Schubert's method.
Powell [29] suggested using a varlational method to develop
algorithms to minimize functions which have many zero second deriv-
atives. This is reiterated in his survey paper [3I]. Toint [3?] has
applied the variational techniques as in Creenstadt [23]. This is
further discussed in Chapter 3.
For general surveys. Rheinholdt [33] contains good background
material, an overview of many different types of available methods,
and an extensive bibliography. For a discussion of quasi--HNewton
methods. the survey by Dennis and More [15] reviews the development
of those methods, gives convergence results and suggests many open
problems. The recent paper by Duff [16] has a comprehensive bibliog-
raphy for sparse linear equations.
I.' Organization.
in Chapter 2 superlinear convergence is shown for Schubert's
sethod; a Kantorovich analysis is also given. Chapter 3 contains
the developaent of a symmetric version of Schubert's update. An
aigorithm is given and a linear convergence result is proved. in
Chapter 4 a matrix factorization is used for the solution of nonlinear
equations. The Schubert update formula is used to update a trian-
--H9-
gular factor while the other factor Ia not changed for a finite
number of iterations. Continuity of the factorization is shown be
fore proving the local and superlinear convergence of the algorithia.
Chapter 3 suggests some further applicationa of these sparse tech
niques to other problems. Also the findi??s are autinarized.
All of the algorithms described here are considered in the con--H
text of being used locally. New notation wilt be introduced as re-
quired. Some frequently quoted theorems are stated in the next sec-
tion for reference.
1.3 Irequently quoted theorems.
Theorem 1.3.1 Bounded Deterioration Theorem Broyden Dennis, and
?ore [4J. (The reference for the theorem as stated here is Ll5?.)
Let F: ?n ?n be continuously differenttable in an open convex
set D0. Let xtcD0 be such that F(x*) 0 and F'(xt) is nonsingu--H
lar. Suppose F' is Lipschitz continuous at 5* there extsts k > 0
such that
IIF'(x) --H F'(x')II < elix --H 5*11 xcD0.
If there exist constants O1 and o2?> 0 for which
(1.5.1) IlBkel?F'(S*)II < ?I+a1o(x?x?+?)]iIBkF'(x*)II + m2O(s?x?+i)
with o(x?`5k+l? max(llx??x*11 15k+l?5*11?' then there exist ?,? >0
such that if ? D0 and 11x0--Ha*ll <  and IIB0?l?(xt)II < I then
the iteratton (1.1.5) and (1.1.6) is well defined and converges linear--H
ly to x'.
- 10 -
Theorem l.?.2 Banach Perturbation Lemma [27?.
--HI
Let ACr L(1n) and assume A is Invertible, with IA II < m.
If IA --H CII < 8 and 8? < I then C is also invertible and
IIC?1I I < --H1f--H58
Theor?m 1.5.3 Superlinear convergencei4].
Let ?: iRn ?? he continuously differentiable in D0, open. convex.
Suppose there exists x?cDo such that F(x*) 0 and F'(x*) is non--H
singular. Let ?8k? be a sequence of nonsingular matricex in
If there exists an x0 such that (xkl, defined by
xk+l xk --H -IE(xk) k
8k
remains in D0 and Xk I x* for any k, and xk?xt Then (x I conver--H
k
ges superlinearly to 5* if and only if
(1.5.2)			tim 118k--HF'(x*)J(xk+l--H5k?11
15k+l			--H			5k1
2. Local and Superlinear Convergence ci Schubert's Method
2.1 Preliminaries.
for continuously difierentiable 1, Broyden's method [2? for solv-
ing (1.1.1) generates a sequence ci matrices which approximate the
Jacobian. Civen x0 and B0 nonsingular, define the iteration:
(2.1.1)			Xk+l = xk --H Bk1F(xk)			and
--H Bksk)aTk
(2.1.2)			Bk+l = 8k +			T
kk
where			= F(x			) --H			and 5k = xk4l --H xk			The matrix Bk is
k+l			F(xk)
corrected or updated to obtain the next approximate Jacoblan Bk+l.
The formula (2.1.2) is called the Broyden update; it replaces (1.1.6)
in ?he general formulation ci the quasl:Newton methods.
Definition 2.1.1 For vectors u,v ?n define the subspace ?v,u by
)v,u iA L(lln) : Au = v). This is the quotient ci v by u.
Broyden's update is a projection of Bk onto Q This oh-
servation is due to the fact (l? that ?k+l is the unique solution to
(2.1.3) min I IS --H BkIIF : SCQy?e? I
where jj?j1? is the Frobenius norm. Another way ci viewing the
Frobenius norm of a matrix in L(lln) is as the L2 norm ci a vector in
n5
The convergence properties of Broyden's method are given in the
--H Ii --H
--H 12 --H
following theorea. Let J(x) = be the Jacobian of F at the
point x.
Theorem 2.1.2 [l5? Suppose F: ?n ?n satisfies (2.1.4) and (2.1.5):
(2.l.4)(a) I is continuously differentiable in an open convex set D0.
(b) There exists an x*cD0 such that I(xt) = 0 and J(x*)
is nonsingular.
(2.1.5) There exists a constant, % > 0 such that
iIJ(x) -			< Kix --H xII? for all x,xED0.
Then there exist c,i >0 such that if x0 rO0 and ?o nonaingular
aatisfy 1ix0 --H a*ll2 < E and 1180 --H J(x?)tl? < 1 then Broyden's
method generates 18 1 with Ilk nonsingular and (a ) which converges
k k
to a* and convergence is superlinear in tbe sense of (1.1.8).
The update (2.1.2) does not take advantage of any sero--Hnon:ero
structure. The second term on the right in (2.1.2) is a vector outer--H
product, and is in general a full matrix.
2.2 Schubert's update.
Schubert [36] proposed a modification of Broyden's method to accom--H
o+odate systems of nonlinear equationa with sparse Jacobian. Following
the definition for linear problems [16], a Jacobian will be considered
"sparse" if it has less than l0? nonrero entries. Recall that the
iteration (2.1.1) and (2.1.2) would be implemented as (1.1.9), (1.1.10),
and (2.1.2). The solution of a linear system with 8k as the coeffi-
cient matrix is required. The pattern of zero and nonzero entries, or
--H 13
sparsity pattern, in the Jacobian is assu?ed to be fixed for all x.
Definition 2.2.1 For j = 1,2n, define zjC3t? by
(2.2.1) Zj = (v ?n : etv = 0 for all i such that [J(x)%1 = 0
for all x? 1:n1
th
Z is the subspace determined by the sparsity pattern of the J row
of i(s).
Also. define `be set of matrices Z that preserve the sparseness
structure of J as
(2.2.2) Z (A c L(lin) : Ate?r Zj for j = 1,2
A matrix, 11, is said to preserve the sparseness structure of i if B?Z.
Corresponding to each Zj there is a projection operator that
lin
`ups			onto Z
Definition 2.2.2 For j = 1,2n, define v tL(? ) by
(2.2.3)			Dj = diag (dj1?dj2djn?' where
djj =			?, ot1??????se????			= 0 for all ??
It is clear that Dj is a projection.
In
Lemesa 2.2.3 Let J(x) be the Jacobian of FcC (I ) . Then for j = 1,2,
(2.2.4)			eTJ(x) = eTJ(x)D			for all xc 1in
j j			j
(The vector Cj is the ?th column of the nen identity matrix.)
Proof: By Definition 2.2.1, ?eTJ(x)?T t Zj For any vc 2, DjV = V
--H 14
Equ?tlon (2.2.4) follows.
To simplify notation for a given iteration, no subscript will de-
note the current iterate, and the subscript "plus", +, will indicate
the next iterate. For esaeple, (2.1.2) would be
T
(2.2.3)			B			= B + ?Bs)s
+			T
The idea Schubert had was that the update should be asade row--Hwise
to preserve the sparsity pattern. The Schubert update can be susmar--H
ized as
(2.2.6)
tn=1 ejeT (VT?BS) ?T
5 ?o			? 5i5J
where 5 = D s. Note that for = 0, the zero correction is added to
j			j			si
th
the j row of B. ?or a scalar ar ? the notation of the pseudo--Hin
verse is a suterscr? "plus":
if a 0.
The update (2.2.6) can be rewritten as
+			j=l			j
(2.2.7)			B = B + 1tn			(5T5jj			--H Bs) ejaT.
Lesma 2.2.4 Let F r c1(ii0) Let Bc Z; x+,xc ?n Let a = x+ --H x and
y = F(x ) --H F(x). If 8+ is defined by (2.2.?), then 8+ tQ nZ.
+ y,s
Proof: For I = 1,2n, multiply both sides of (2.2.?) by 5T1 to
obtain
(2.2.8)			T T			??i?i? ?T1(y?85) 5T
e1B+ = e1B +			T+
--H 15
z . Lin--H
Since B , Z it follows that BTe1 t Z ? By definiti?n. a
earity implies
B+e1 rzi. and thus 5 cZ.
To prove 5+ C Q			exa?ine B 5 and y compnnent he cosponent.
y.s			+
If s? i 0 then
(2.2.9)			T5 ?			eTss + (5T5)+			T(			B)			T
i +			I ii			s1y			--H 5			5j5
? . For 5 6 0
Since ?i is a projection of a, s?s = ?i?i I
(aT			+ T			--H1
i?i? = ??i ?i?			by definition, and equation (2.2.9) reduces to
(2.2.10)			o+T1B+s = 5T155 + eT(y --H sa) = eTy.
If St 0 then by the Mean Value Theorem (273, there exists
t1t(o,i) such that
(2.2.Ji)
T			T
CiY			e1?(x+) --H 1(x))
= eTJ(x + --H
i			ti 5)5
=5T??5 + t1s)s1, because J C Z,
= 0, since ?i = 0.
Thus, for I = 1,2			n, eT1s+5 = eTy. which implies that B+s =
Therefore B c Q
+ y,s
The Schubert update is the unique solution to 5 minimirstion prob-
1cm. The update B+ can be viewed as the projection of B onto Q flZ.
y,s
Theore? 2.2.3 Let BC z; y,sC			with a # 0. Define B+ by (2.2.7).
Then B			is the unique solution to			min(ljS --H B1i1: Sc Q			flZl.
+			y,s
--H It --H
Proofs Let Sq nZ. Then ?s = y end
y,e
T			+T			T
(2.2.12) IIB+ --H B11g = II ? ?5i5i? ei(S?B)si cici II?
1=1
? n? (sTisi )+eT(??B)sil2 15i1122??
1=1
? n? IeTi(S?B) I 12214 I Is - ?l F
i=l
Uniqueness follows aince f: L(lin) lin defined by f(A) = IIA-BII? is
stricLly convex and ??5nz is convex.
2.3 Local and superlinear convergence.
In l9?l, Broyden gave a local convergence proof for Schubert's
o+ethod [3?. He conjectured that the method was superlinearly conver-
gent based on numerical results. A different, simple proof of local
convergence is by an application of the Bounded Deterioration Theorem
of Broyden, Dennis and More [4?. Although superlinear convergence was
first proved for Broyden's method in [4], the techniques for the proof
in [133 are extended to Schubert's method here.
The Bounded Deterioration Theorem gets its nawe from the hypothe-
sis on the behavior of the approximate Jacobians (1.5.1). Define
= B --H J(x*) and E = B --H J(x*) for B,B+ c L(? ) . Then (1.5.1)
may be written as
(2.3.1) lIE_ II < [1 + a?0(s?x+)]IIfII + m2o(x,x+).
The norm ??.l is unspecified. The error IIE+II may be worse than
IlEll, but only ins controlled way. The difficulty in applying the
--H 17 --H
Bounded Deterioration Theorem is to show for a particular norm that the
estimate (2.3.1) is satisfied for the update method generating the 5--H
quence of ??proxinate Jacobians.
For Broyden?s method. the Frohenius norm was one for which (2.3.1)
was proved in 113] That is the same norm in which Broyden?s update was
a minimum norm correction. See (2.1.3). Theorem 2.2.1 shows that
Sthuhert?s update is also a minimum norm correction in the Frobenius
norm. The following tteorem provides an estimate for the Schubert up--H
date in the Frobenius norm. This is the crux of the subsequent conver-
gence results.
Theorem 2.3.1 Let B.JtZ and ???? ?n with a ? 0. If B+ is de-
fined by (2.2.?). then
2 2			Il(B--HJ)5lI22			?T)+eT(?Js)]2.
(2.3.2) llB+--HJIl? < lIB?JIIF --H			2 i?1 ?1?i			1?
115112			+
proof: By definition of the Frobenius norm,
IIB+?JlI2F =			t (eT1(B+?J)ejJ2
i.i=l
Substituting for B+ using (2.2.?) gives
SiCjJ
(2.3.3) lIB+?JllF2 =			??			[eT1(B?J)ej + (aT1si)+eT(y?Bs) T			2
n			T			2 T			T
? i[e?(B--HJ)e + 2e1(B--HJ)e?(5T151)+eT1(y?Bseja1
i,i=i
+ (Ts)+2eT(y?Bs)]2(eTs )2j
5iji
--H 18 -
18??11E2 + 2 n? (sTisi)?eT(B?J)sl eTj(y?ls) ???1(sT1s?)+?eT(y?Bs)]2
The second equality follows by expanding the quadratic; the third is by
summing each term over J 1,2n. ?ow, add and subtract Js in both
the second and third terms. Then,
+ 2 n? (sTisi)+eT(B?J)si [eT?(y?Js) + 5i?8?3?53
(2.3.4)			18+??112F = I IB--HJI 2			T
+ ???j(sTs?)+[eT(y?Js) + eTi(J?B)s?2
= 18??112F + fl? (sT?s?)+i2eT(B?J)a?[eT(y?Js)] + 2eTi(B?J)si[eT(B?J)s)
+ ?eT(y?Js)?2 + 2eT?(y?J5)[eT?(J?B)s? + [eTj(J?B)s]2)
T			T(8?j)5 = eT(B-J)a. The first and
Since 8,3  Z, ?8?3?5i r Zi and e? ?
fourth terms inside the summation cancel; the second and fifth combine
to yield
(235)			BJI2 = 118?3112 + n? (T)+(eT(?J5)]2			T(g?J)52)
-.			+?F			F			___ s?s1??			ei.
Finally observe that 115i112 < 115112 implies that (5T5_)+<?IIs1I?22
Therefore
(2.3.6) l18+?3112F <			2 --H 11(83)5112			(s?s?)+eT?y?s 2
18?311F			2			+ n
115112
Theorem 2.3.2 Suppose F: ?n??n satisfies (2.1.4), and (2.3.7).
(2.3.7) There exists K = (k1,k2,... kn) , 1n with k > 0 for i = I,
- 19 --H
2n, such that
lleTJJ(x) - J(x)?lI2 ?<killx --H xll? for all xxcD0.
Then there exist ,i ) 0			such that if			"o cD0 and			80 nonsingular
satisfy 11x0 --H x*l 12 < f			and 1180 --H			11F < ?			then Schubert's
method generates 18k1 with 8k nonsingular and 15k1 converges locally
and superlinearly to x+.
Proof Let K = 11K112; ? is a Lipschitz constant for J in the Fro--H
benius norm. From Theorem 2.3.1 with J replaced by J+
E = 8 - J+, and Ii = 8 -
2 2			11EslI22			n
(2.3.8) llE+ lIF < jIElIF - 1151122 + ?=?1(sT?s?)*[eT(y?J*s)]2.
For each St ? 0 estimate the term in the sum.
For j = 1,2n, define			1i?I by
(2.3.9)			o+ (t) = eT??F(x) + F(x+ts) --H tJ*sJ.
I			I
Then (0) = 0, o+J(l) = eT1[?F(x) + F(x+) --H j*s) = T --H J*s), and
e1(y
= eTj[J(x+ts)s --H J+s? = eT[J(x+ts) --H J+Jsj since J , Z. 8y the
Mean Value Theorem l2??, there exists t ? (0,1) auch that o+ (1) = s'(t )
I			I			ii
Thus,
let (y--HJ+s)l2			leTt[J(x+tts)--Hjt)a12
2			=			115112
115i112			t2
lleTt[J(x+t?)?J+]ll2
--H< , by the Cauchy--HSchwartx inequality,
2 --H 2
< killx+tis --H x*l 12 by hypothesis.
(2.3.10)
- 20 -
<2			2
k1c(x,x+)			since tie (0,1),
and "here o(x,x+) = `sax (lix --H x*l12,llx+ --H x*1i2).
If any St = 0 then 0 =			T			+T			2 2			2
??1?1? [e(y?J+s)] < ko(x,x+) -
Therefore,
(2.3.11)			________
2 2 I1Es1122
IlE+llF < 111?F --H _______
+ eu(x,x+)2
Disregarding the negative term it lollows that
(2.3.12)			iIE+IlF < llEllF +			n(x
2 2 2 2
since if a,b,c > 0 and a < (b + c ) < (b + c) then a < b + c.
Inequality (2.3.1) is satisfied in the Frobenius norm with a1 0
and a2 = . The hypotheses cf the Bounded Deterioration Theorem 1.3.1
(iSJ are satisfied and this proves local and linear convergence of Schu-
bert's method.
Theorem 1.3.3 gives a necessary and sufficient condition for
superlinear convergence. To show that (1.5.2) is satisfied for Schu-
bert's method requires 5 more careful estimate of the right hand side
of inequality (2.3.10). Observe that
T			2
IIE(I --H 55			2 _______
T?11F =
This and (2.3.1)) give
(2.3.13) ilE+llF ? [IlEll2?I1Esl122 ?? +
I			1151122
--H 21 --H
The proof now proceeds as in Dennis and More [p. 38,153. For a,b > 0
(a2 --H b2)4 <a --H b2
--H			--H			ao that
2a
(2.3.14) IIE+ llF<llEllF?(2llEllF1) IlEsli2 2
115112			+
Let ?			lIEkIlF lIB			--H J?Jl			I IEkskI 12
k			F'			0k = 0?5k'?t+1?
k2
Inequality (2.3.14) is now
(2.3.15)
< [1 - (2r?2?1 *k23'k +
This gives ?+l < rk + Kck The sequence ?5k1 is linearly convergent
which implies 1nk? is bounded. Let n be an upper bound for 1nk1. Then
(2.3.16)			(2??l??k <			--H `k+l + ?0k' and
? o+k2 < + __ 0k
k=o			0
forcing lim ?k			?			Condit4on (1.5.2) is satisfied and superlinear
convergence of Schubert's method follows.
2.4 A lantorovich analysis.
A general Newton--Hlike method is of the form
(2.4.1)			Xk+l = `k --H Ak1F(xk)
with x0 given and (Ak1 a sequence of invertible operators. Dennis [8)
gave a Kantorovich analysis for the iteration (2.4.1). In [?) he studied
--H 22 --H
Broyden's method with a Karrorovich analysis. The Kantorovich analysis
consists of constructing a scalar iteration to majorire tit' ?ccror
iteratton. The acalar sequence is used to study the properties of the
wector sequence.
The first lemma is a corollary to a theorem in ?7?. The bound will
be used in proving the Kantorovich theorem.
Lemesa 2.4.1			Let			FeC1?to? fori0clI . Let F satisfy
(2.3.7); let			11K112. Then for any			x,yr D0,
(2.?.2) 111(y) --H 1(x) --H F'(x)(y--Hx)?i2 < ?e?iy-xI(2.
Proof: In vector notation, with the integration occurring cosponent--Hwise.
(2.4.3) F(y)--HF(x)--HF'(x)(y--Hx) = j01[F'((l--Ht)x+ty)--HF'(x)?(y--Hx)dt.
Thus,
(2.4.4)			11(y)-lt(x)-F'(x)(y-x)l12 --H< lli01[F'((l--Ht)x?ty)--H1'(x)](yx)dtJl2
<!011lF'((i--Ht)x+ty)--HF'(x)lI2lly--Hxll2dt.
Now ll?ll2 ? ll?111 By definition of 11?llF and (2.3.7),
lF'(x)--HF'(y)l11 < eIIx?y?I?
Thus
(2.4.5)			111(y) --H 1(x) --H F'(x)(y--Hx)ii2 ? lol etiiy-xii? dt
= ?eIIy--HxII22.
--H 23 --H
The following is another technical lemma that is essential to the
Kantorovich analysis. The proof depends on the estimate obtained in
Theorem 2.3.1. The estimate in this lemna shows the relationship of
the error of 8k+l as an approximation to J(xk+l) to the error of
the initial approximation and the total length of the steps taken.
Lenuna 2.4.2 Let 5: 3R? ain satisfy (2.l.4)(a) and (2.3.7). Set
lIKIl2. Let 1Xk? and 18k1 be generated by Schubert's method. If
(xi)t=+11CD0 then
(2.4.6) llB?+1--HJ(x?+1)ll? < 11B?--HJ(50)l1, + 2?? ?? lIxk+l?s1t1l2
Proof: An application of Theorem 2.3.1 yields
(2.4.7) 1lB?+1--HJ(x?+1)ll2 < 118 --HJ(x )jl2 --H ll(B?--Hj(51+1))s.112
5 --H			in			in+l F			2
115.112
+ n? (sTjs1)+[eT?(y??J(x?+1)s )J2
< 1lim??5in+i)II2 + ???1(sT?s?)+[eT?(y?J(x?+1??m???
liere s? Dj5? The techniques in Lemma 2.4.1 can be applied to obtain
leTiym?3(xin+l)a?1 < kkillStll2llSmll2 Thus,
(2.4.8)			llBm+i?J(Xm+i )ll 25< ll8,,J(x?+i)lls2 + t1??1k2il15m
2 %llKll2llXm+l?Xa112 and
l1B??J(x?<1)lls +			2			2
(2.4.9) 1l8?+1--HJ(5?+1)ll? < 118m??5m+i )ll? + ?ells.+i?smll2
< 118_--HJ(x_)ll? + 11??5m+l??3?Sm?1lF + ?ells?+1--Hx1ll2
--H 24 -
< IIB1--HJ(x?)II, + ???IIxm+i?x.II?
since IIJ(x1+1)?J(xm)I 1F < e11xm+l?Mm112 Therefore
IIBm?i?J(xm?i)11F --H< IIBo?J(xo)II? +			? IIxk+l?xkII2
The following theorem asserts the existence of a zero of F. The
technique for proving the result is the ?antorovich analysts. The details
are almost identical to the analysis of Broyden's method in 9].
Theorem 2.4.3 Let F: lin lin satisfy (2.l.4)(a) and (2.3.?). Let
= IlK?I2 and let B0c L(lln) be nonsingular and suppose there exist
constants f.i,n >0 such that IIBD1112 ? f IIJ(x0)--HB011, <i
and IlBo F(x0)112 < n.
(a)			Define h' =			--H and r'			lt(l--H2h')? (i?Bl). If gi < I
and h' <4 and B(x0,rf)C00 then F has a root x+cB(x0.r') and if
h < 4, ?t is unique in D0?B(x0,r+'). The sequence (x?') generated by
--H --Hl
(2.4.10)			Xk+l = k			?0 ??5k'?
converges to x* from any x0c000B(x0,r+').
gel			1--H(l--Hah)4
(b) Define h = (l?3ii)2 and r = ,ge (l--H31i). If 3f1 < 1,
b <?8 and B(x0,r)C00, then Schubert's method converges to x+.
x'cD			and			e			? +
Proof: (a) Let			suppos			x0 --H x0112 <r < r'. Consider
the iteration (2.4.10), and note that the approximste Jacobian B0 is
used throughout the iteration. In tbe first part of the proof, (2.4.10)
--H 25 --H
will be used to establish the existence and uniqueness of tbe root X?
of F. The convergence of the sequence (xk) is established at the same
time.
Define a scalar function f1: It It by
(2.4.11)			f1(t)			4?t2 --H (b0 --H 6)t + b011B01F(x0)ll2
where b0			?--H1			Define the scalar iteration
(2.4.12)			ts+1			tk' + Sfl(tk')			for t0'c[llx0--Hx0112,r0).
The next step in the proof is to show that (tk') inajorizes the vector
sequence (x?). Let t ? [rf,r+'). The proof proceeds by tnduction:
lx? --H xO'112			IIB01F(x0)112 < tj --H t0,
Now let x?rB(x0?t?)flDO and suppose that tk'C (c,t') for k = 0,1
I, and
(2.4.13)			llxk'.l --H xk'll2 --H< tk'+l --H t' for k = 0,1`1.
Then ? llxk'+l --H xk'112 < te' --H t0,. From the vector iteration,
k
xk'+l - xk' = -B?1E(x'), and
0			k
o = F(x?'?1) --H Bo(x?' --H x?'?1) so that
--H Xm' = Bo11E(xm)?F(xm?i)?Bo(xm'?x?i)) Thus,
(2.4.14) llxm'+i --H xm'll? = llBo1[t(xm)?l(x?1)?Bo(xm?xm?1)]ll2
--H 26 -
< B[llF(xin)?F(xin?1 )?J(X??1)(X?Xin'?1)ll2 + ( lJ(Xin?1 )-J(x0')l12
+ 1lJ(x0)--HJ(x0)ll2 + IlJ(Xo)?BoIl2)1lXin?Xin?1Il2]
2			+			+ Kin?2
--H 6[??llXin'?Xin?1ll2 + (6 Kllxo?xo'll2 ? llxk'+1?xk'll2)
I1Xin??Xmi II 2J
< B[??tin?t?I)2 + (6 + Kt' )(f?--Ht' )],
in-I			in in-I
by (2.4.13), and llx0--Hx0ll2?< to' ?ow from (2.4.12) and (2.4.11),
fi(tm?'I)
(2.4.15)			+ -f'(r' ) --H b0 + Kt'			--H (b --H 6)
in in-I			I in-I			?-I			0
6+ t'
in-I
SubstItuting (2.4.15) 6n (2.4.14) yIelds
(2.4.16) 1Xm+l?Xm'112 < 1?I?(tm?tin?1)2 + f1(tm?'i) + f?(tin'?l)(tin'?tm'?1)]
< If (r') --H t'			--H
--H			I in			m+I			in
Therefore ? llXk+I?kll2 < ?m'?l?to' < t' for all in. Thus the se--H
k=o
quenre (xk'l Is Cauchy and must converge to some x+r lI(x?,t?). By def--H
0
inition of r+, r is the only solutIon of f1(t) 0 in the interval
[o,t'?; ?, + ??I(t') < t' since f1(t) < 0 for r ? ?r,r;). with
equaliry only for t' = r'. By a theorem in Kantorovich and Akilov ?4]
and quoted 6n ?(?, x+ is unique In 1I(x r') and x converges to
0' -			k
x* from any x0'r1I(x0,r').
--H 2? --H
if h' < ?, then r' and r' are distinct, an t?? s?ae theorem implies
x+ is unique In D0r'B(x0,r+') and ixk') converge' t s, t'?a any
x0, r D0r'B(xo,r:).
(b) The first step is to show that ?Bk1 are invercible and to pr?
duce bounds on IIBk1 12 Assume for ? >0 that
(2.4.1?)			rn? llxk+l --H xkiI2 < r.
By hypothesis, 11801112 < 8 so that
(2.4.18)			ll8?l8			--Hill < 8118			--H
0 m+l F --H m+l 8O11F
< 8?ll8m+1?J(Xm?i)llg + 1j?5rn+1)--HJ(x0) 11F + 1??'O??8Ol1F?'
using the triangle inequality after adding and subtta.tIng both
and J(x0). Now Le?a 2.4.2 can be applied to the term 18m+l???5m+l?11F
yielding
a			m
(2.4.l?) ll8?l8			3
0 vn+l 111F --H< 8?211??5O??8O11F + 2'_? 11,k112
+ IIJ(x,,,,1)-J(x0)ll?3
< 8[2i + 32Kk?OIlSkll2 + `Il?m.t?XOII2]
since llJ(xO)?8ollF < I by hypothesis, and
IIJ(Sm+1)??(xO)llF < ( ?? lleTi(J(xe,,l)?J(xm))ll2)?
< ( ? k1 llx.+l?xOIl2)
1=1
--H 28
I Ixm+i?xoI 12
Applying the t?iang1e ineqoality to IJXm?l?XoII2 gives
(2.4.20) II8o1Bm.i?IIIi ? 281 + ?8e w? 11Sk112
<28i + ?f?r			by (2.4.17)
--H 281 + 8?t?l--H(i--H8h)4](1--H38i)
by deflni?ion of r. Now, 0 < h < 8 implies that			0 < (l-8h) -<
and [i--H(l--H8h)43 < 1. Then c --H85[1--H(1-8h)43 < 1.			Therefore,
(2.4.21)			--Hl
18o 8m+1?1 j? < 381			c(i--H3?i)
(l--Hc)3ii			c
< (I--Hc) ? c			1,
since 388 < I by hypothesis. Since l.?J2 < I1?11F' (2.4.21) 1w--H
--H1
piles 1IB0 8m+i?1112 < I and by the Banach Lem,na 1.5.2,
--Hl			--Hi --Hl			--H1
(B0 8m?i? 80 = 8m+i exists and
(2.4.22)			IIB?1 II			(1--H28i)i$f 15k1 2)
?+I 2
k=0
The next step in the proof is to define the scalar iteration which
wajorizes the sequence (Xk) produced by Schubert's method. Let b0
h0= s and set to = 0 i'efine the scalar function f: Ii Ii by
--H 29
(2.4.23)			f(t)			2?? 2 --H (b0 --H 3i)t			h&?n.
Define hk			i			and b --H b?1			Dci Inc tIe ,.?l'r iteration
(1--H2ii?8ct?)			k --H k
(2.4.24) tk,i = tk + bkf<tk)
Next, the se?uence ?`k1 is shown to be strL?tl? in ,easi,? and
bounded above, hencergent. An induction is u?,j t `,Cove itkl
remains bounded. Observe tbat
t =t? +b0f(t0)			b0b0n
and ? > I?bclF(xO)II2 1151112 Since r is the miii s-uot of `be
quadratic (2.4.23) or equivalently 26cr --H (i--H3iii)i `I. elementary
analysis shows that r > ?. Thus
(2.4.25)			r > 1151112
Assume that			tkc (O,r) for k = 1,2m.			Then by (2.4.24)
(2.4.26)			r --H tm+i = r --H tm + h??f(r)			--H f(tm)
(r?tm)[l+hmf?(i)]. for ? (im,r),
since f(r) = 0 and by an application of the Ne?n Value iheorem [27?.
Next, the derivative is shown to be negative. Differentiating (2.4.23),
and evaluating f' at t gives
(2.4.27)			f'(?)			4 i --H (b0 --H 3i)
< 4cr --H (b --H 3ij
0
- 30
?l--H(1--H8b) ](b0--H3i) --H (b0--H3i), by definition of r,
= --H(l--H8h) (b0--H3i) < 0.
This implies I + hmf'(i) < I and with (2.4.26) yields t?tm4l <t?tm
Also, hm? 8 = h0, so that
(2.4.28) h,,,f'( i) > --Hhm(bo?3i) ? --Hh0(b0--H38) = --Hl?3fi > --HI.
Thus I + hmf'(i) > 0, and 0 < r --H tm?i < r0 --H tm. Therefore itk1 is
strictly increasing in ?0,r]. Hence itk1 converges to t* < r. if
ft < r, then t* = t* + (1 --H 2ii --H --H5get*)--Hlf(t*) and f(t*) must be
zero. Therefore t* = r since r is the smallest root of f in [O,r1.
it remains to show that k=?o''5k''2 < r for all m. The first
tep of the induction is satisfied by (2.4.25). Suppose that x0,x1
Xm satisfy
(2.4.29)			115k+i --H XkII2 < tk+l --H tk k			0,1m--Hl.
m-l
Then ? 115k1 12 < tm < r and the first part of the proof ot (b) shows
1im1 exists with I 8m1 I 12 < hm To complete the induction, it remains
to show that lix			--H x			< t			t			F
m+l m 2 --H m+l - m or Schubert's method,
Xm+l Xm = 8m1F(xm) and 0 F(Xm?i??8m--Hiam?i.
Now, Xm+l Xm = Bm1EF(xm) - E(xm?i) - ?.s-l?mi?`
--H 31 --H
and taking norrns yields
(2.4.30)			llSm=1?Xmll2 =			I
1?m [F(X?)?F(Xm?i)?B?iS?i]ll2
< h_[llI(xm)?F(xrn?1)?J(Xm?i???--Hl 112+ llJ(x???1)--HB1?1 121 l5?i 1121
< h_??<?fr_ Xm?l ii2 + (6 + 23??2llskll2)lfr.xIIH1ll2l
by an application of Le?as 2.4.1 and 2.4.2. The induction hypothesis
(2.4.29) with (2.4.30) gives
(2.4.31) llXm+l?Xmll2<hm?4?(tm?tm?l)2 + (6 + 22?t?i)(tmt?i)]
The scalar function f satisfies
f(tm?1)
(2.4.32) (tm?tm?i) + f'(tm?i) bm?i + f'(t?1), by (2.4.24)
bo(i?2h06??ho?tm?1) + 4et?1 --H (b0--H36).
by definition of h??1 bm?i,h0 and b0, and evaluating f at t?1. This
sieplifies to
_________			5
(2.4.33) (tm?tm?i) + f'(tm?i) =			b0 + 26 --H --H2?t?1 + 4?t?1 --H b0			36
=			6 + 2?tm?l
Substituting this expression in			(2.4.31) yields
(2.4.34) IlXm+l?X?ll2 < hm???(tm?t?1)2 + f(t?1) + f'(t?1)(t1--Ht?1)1
< h_[2e(t --Ht			?2			)(t
+ f'(t			--Ht			) + f(t
__			m ?l			?l			m ?1
--H 32 --H
b f(t ) by the Taylor expansion of f at
t1+1 L1 by (2.4.18).
Tbe induction is complete.
Therefore ? llXk+l?xkll2 < t1 < r for all I. The sequence (xk)
k=o
is Cauchy and must converge to some x+*tB(x0,r). From Lemma 2.4.2
(2.4.35)			IBm --H ?(51)lI2 < a + 31
2 ? llxk+l?xklI2
ko
2m
and			1i8m112 < 1lB1-J(x1)1i2			IIJ(x1)J12
< 11J(xm)I12			+ a			+ --H2?t1
< IIJ(x0)lI2			+ a			+ --H2?t1
< IIJ(x0)112			+ a			+ --H2er c, a constant.
flence
(2.4.36) IIF(xt*)I12 = lim IiF(xk)I12 = lim llBkll2llxk+l?xkll2
c lim llxk+l?xkll2 = 0
NoM, r < r' and the
+ uniqueness of x* implies ?*A =
3. Symmetric Scbubert Update
3.1 Preliminaries.
Let f: 1f0 ? be given. Consider the prohiem t finding a local
minimum of f on some open set DDC? . Let ?t fl?? b? the minimum.
so that
min (f(x): at D0I.
in the following, assume that f is twice cont1nu?sly differentiable.
When solving this problem using a Newton--Hlike m?th?J tu find a
?eru of g(x) Vf(x). it is usually the case that an a?proaimation to
J (a) v f(x), the Jacobian of g or the Nessian of I, luks as much
like J as possible. Since v f is symmetric, it I :?simjted by a
8
sy'nmetric matrix.
Definition 3.1.1 Define the set of symmetric matri?e? 5 by
5 - ?At L(1f0) : A AT).
A Newton--Hlike method, (1.1.4), for the solution of ti? minimization
problem can he written as: given xc D0 and Bc 5.singular,
(3.1.1)			x+			x --H B?1g(x)
8 cSflQ
+ y,s
where y g(x+) --H g(x) and 5 x+ --H x. Specific m?ii,?ds for choos--H
ing a 8			in S?Q			have the property that near the solution, 8
+			y,s			+
is nonsingular.
Suppose that v2f is sparse. The quasi--HNewt?. m?thods, (3.1.1)
require the solution of a linear system, namely
--H 33 --H
34 -
(3.1.2)			11= =
The objective Is to define a method of updating B to preserve the
sparseness structure of the Nessian. Some quasi--HNewton methods for
solving the minImi?ation problem are designed to update B?1 = H [15].
However, there is no guarantee that II will be sparse if B is. Here,
only methods for updating the Hessian are considered. As in section
2.2, V f will he considered sparse if it has less than 1i37. nonzeroes.
The sparsity pattern is again assumed to be fixed for all x.
Assume further that
(3.1.3)			[s2f311 ? 0 for i = 1,2n.
Otherwise the 1th			component of St is linear in			x1 , so that			x1
can be written as a function of the remaining Xj			i = Ii--HI
n, and the dimension of the problem can be reduced.
3.2 Notation and technical iem,mas.
Definition 3.2.1 Let ? r ?n and define the components of 5 by
(3.2.1)			= eis for i = 1,2
The subspaces Zj? j = 1,2n, have the same meaning as in
Definition 2.2.2 with J(x) now being replaced by J (x). Instead of using
the projection operators defined in Chapter 2, it will be more conven--H
ient, because of the complexity, to work with the individual components.
Definition 3.2.2			For			J = 1,2n define the t2 projection of
5 c Ii onto Zj			by			s? with			x?			defined component--Hwise by
--H 35 -
(3.2.2)			eTisj			0:			?? [V2f]ji=?00
if [V2f]jj
Lemma 3.2.3			Let			ac?n and suppose f:I?.? t			?J [V2f?ii #0
for I = 1,2			n. Then
(3.2.3)(a)
T			T
ea=es,
i			ii
eTs eTs			= eTs			ets , anA
J ii			ii			ii
T
w = ? C16?a e r Z1 for a? c Ii , j - I
i=i
Proof: By hypothesis, [V2fijj # 0, which implie? eT? - = eta.
(b): if [V2fJj? # 0, then by syrr?netry [521 ?ii ?			a that
T			T			T
e s			Pj =			5			?ence (3.2.3)(b) holds. If [v'			- 0 then ejSj = 0
ii			1'
and both sides of (3.2.3)(b) are zero.
T T			2
(c): CjW = e1a1a? = 0 if [VfJ?j = 0 The,?I..,e a's Z
Definition 3.2.4 Let ? ?n			Define PrL(l1?) Iv
(3.2.4)			P = ?[I --H j?=?1(sTjsj)+eIjs			T
aici'.
Define A,DsL(IIn) by
(3.2.5)			A = ? PjSjCj
i=1
= [p1s1 I ?2?2 I . I snan I,
and
- 36 --H
(3.2.6)			D diag (5T151?5T252
T
Observe that A is the Frobenius projection of as onto Z.
Also, P can be written as
(3.2.7)			p			--H AD+J.
Detinition 3.2.5 For v  ? and Mc L(?n) and an integer 0 < m
< n, def?ne M c L(??) to be the lower right o+? aubmatris of
M			M			and v			,,:`			to he the la?t m coeponents of v.
n			m
The matrix f is a projectton operator. This is a consequence
of the follo?ing lemma on the elgenvalues of P.
Lema 3.2.6 The zigenvalues of f are in fo, max (?, I --H min ?T? ??
p ?0 s s
i			ii
and the Cigenvalues of (I--Hf) are in [mi?(4, min
? ), lJ.
ii
Proof: Without loss of generality, assume all the zero components
of 5 are ordered first. That is.
0, i = 1,2f, and
# 0, i = 1+1
Then An?f is the nonzero part of A, and 0n--Ht is invertibie by
(3.2.3)(a) and Definition 3.2.5. Let m = n--Ht.
From (3.2.7), partitioning th? matrix,
--H 37
(3.2.8)			I			--H
ol o
o A D?1
$1			4[lm0AmD?'J
It is clear that I has i elgenvalues equal to ?, a., t?,c rewaining
fl--HI are the etgenvalues of I[lm?Aml)m1J Similarly. I?l' ha? t cigen--H
values equal to 4, and its other fl-I are the eigenva?u?s
4?1m + AmDm13
It is now sufficient to consider I and (1-1) .1 tmeusfon m,
so the suhscripts are omitted temporarily. Observe ii..
P 4[i --H AD?1J and
(I--HI) = 4[l + AD?1j
are similar to
I			-?
- A)1)?? and 41) (1) + A)D
respectively.
Schnabel [351 observed that D + A is the sum of a `ifa?w.al pcsi--H
tive definite matrix and positive semi--Hdefinite ran? ,.,.. ?trices. Thus,
- ?8 -
(3.2.9) D?A			2			?22
2			??
J			ii
2			2
Fv ?3ij?0			p
I			?i
By the Interleaving Elgenvalo Lemma, see Wilkinson ?pp. 95--H98,41]
for exawple,
(3.2.10)
2			(t)
win e.v. (D + A) > win 2p > 0.
i			I
Siwilarly D --H A can be written as a sum of positive semi--Hdefi-
nite rank--Hone watrices
(3.2.11)			D --H A
i<J
v2f]
-p
-J			?P1?j -
?PjPi			p2ij
So win e.v. (0 --H A) > 0.
Therefore P is positive semi--Hdefinite and
win e.v. (1 -P) > wIn			> 0.
`s			T
I			as
Ii
In general, win e.v. P > 0 and
win e.v. (I--Hp) > min(?, win
ilo
ii
This gives the lower bounds.
The bounds on the largest Pigenvalues follow from
Toint [38] also showed that 0 + A is positive definite.
- 39 -
?Tx = 5T(1?p)5 xTpx.
This completes the proof of Lermsa 3.2.6.
Lemma 3.2.7 Let u,v,s ? 11 and P be given b? Then
(3.2.12)			(D+(l?P)??u.v? =
Proof: Suppose that the ?ero component? of a are ?.4e?eJ l(rst. Then
I --H P -
Now,
(3.2.13) (`D+(1--HP)?iu.) = ? [2(sTisi)+eTujeTv (?1tr ?)atum.vm)
-Note that (i?P)mDm = ??0m + Am) which is symsef?i p?sitive deft--H
nite by (3.2.10). Its inverse must be symesetri= a.'(?(ve definite
also, so that
(3.2.14)
I			-l
?1???m um?vm? = (`um?Dm1(I?P)m'vm?
From (3.2.13) and (3.2.14),
(3.2.15) ?D+(l?P)?1u.v,) = ? eTu[2(sTjsj)+eTvJ o+ 4'.m?I'm'(I?P)m'wm)
= (u.?+(I??)?1v?
and the result follows.
Le?a 3.2.e Let arlln and P be defined by (?.??1).
- 40 -
T
(a) If v r ?n ?s such that e1v 0 for all I such that 5? = 0,
then
(3.2.16)
- 4D+(l?f)?1v.v? --H< -
(b) If u ? 1tn then
(3.2.1?) (D+(I--HP)?1u,u? < max (2
n? (sT1s?)+(eTju)2
Proof: A8ain assume that the zero components of 5 are ordered first.
(a): Tten,
(3.2.18) --H?fl+(I-P)?1v,v? = --H ? 2(sT?s?)+(eT1v)2 - 4Dm1(l??)m1Vm?Vm)
2			v v
T2			IT
-			I-PI)
T			I			ma			mm
5 Ao
T			T			2
since ?t?i <5 5 and (Mx,x) > 11X112 for nonsingufar N.
I HI 12
By Le-a 3.2.6,			1?'?I?m11 < I so
2--H
(3.2.19) 1?1???m0a112 --H< 1?1?1?m11211?m112
< IIDmII2
T
< max 5 a
--H L<J<n ?
T
< 5 5.
therefore. (3.2.18) and (3.2.19) yield
(3.2.20)
--H< -			T2			T
sTs ?2 ? (c v) + v v
5 60
I			L
--H?			T2			T
6s1?1?Io+ m m
5 60
T
Tv o
since for			= 0			e1
(b): Direct ev?luation gives
(3.2.21)			(?+(i-?r1u.u) =			T +T 2			D.(i?P)?u??Um
2?5t5? (ejul +
L T			+T2			--Ht
? 2(515 j) (e jU) + 2(Da+Am)U??U?
1=1
<2 ? (sTjSj)+(eTju)2 + ll2Di'(Dm+Am)?1Dm4ll2 llDm?umll22
=2			?? (5T1s1)+(eTju)2			ll('--HP)?1ll2 ?Dm1U.uo)
L			T
< 2 ? (5T151)+ T 2
(e1u) + max ??(D11u1.u1?
+ 60 p
i			i
by Lenna (3.2.6) and
--H1			_			I
l(1???m 112 - o?in e.v.(1--HP)?
for positive definite (1?f?m			Therefore,
D(l--Hf)0,0? max (2, max -14			? (sTint)+(eTu)2
P ?O o
J			j
completing the proof.
Set
3.3 Symnetric Schubert update.
fo?ell f3o? derived the Powell symmetric Itroyden update from the
Broyden update (2.1.2) by iteratively projecting B+ into 5 and
then projecting back to Q . Using the same technique, Dennis [113
y,s
derived most of the syrrnet?ic updates satisfying the linear equation
B+s y from the rank--Hone updates
(3.3.1)			It			It + (y?Its)cT
+			T
with various choices of c r
?ere, the iterative double projection technique is applied to
Schubert's update, alternately projecting into 5 and back to ZflQ
y,s
to derive the synjnetric Schubert update.
Let Itt 5 be given. Define C1 by ?chubert's update (2.2.5),
C1			It + ? (5T?)+?T(y??5)			T
j=l			e?s?
C2 - ?(C1 + CT)
1'
C			T+T			T
2k+l C2k + ? ( es
j=l			jj
(3.3.2)
- 43 -
and
and
C2k+2			?(C2k+l + T
C2k+l)
st + a e
= C2k + 4j??1(sjsj)+e3(Y?C2?s)(cJ J i 3
Let Ck = C2kI tien
(3.3.3) Gk+l = Ck + 4?:??(sjsl)+ej(y?Cks)(e1s3 + sje3)
Ck+la = Gks + 4 ? (sTs)+eT(y?Cs)(eaTa + ajeTs).
j=jiJJkJJi			i
and Y?Gk+la = Y?Gks --H 4 n? (sTJsj)+e3(YG?a)(ejs3sj aje3s).
Set Wk = Y?Cka			Then
wk+i = wk --H 4??=?1(a3sj)+[ejs3s1 + ajeTJaJe3w?
= --H 4 ? (sTa )+(e sTa +seTs)eTJw
i=1 ii			ii jjjjk
= Pwk with P defined by (3.2.4)
0
Now, lim 0 -
k --H 0o + ? ?0k+i --H 0k? , and from (3.3.3),
n T			T
) ?Ck+1 --H 0k? = 4 ? (a a )+(e a + a cT)e w
k=o j=1 ii			ii			jj j k
= 4 ?? (a3aj)+(ejs3 + aje3) e3_I wk
--H 44
Nence.
4 f-(5T5j)+(ejaTj + sjeTj) eT_o+? pkw0
i=l
lim G			B +			+ 5 eT) eT ( ? 1k)(y=Bs).
k			? (mTjsj)+(ejsT			? ?			k?O
i-I
The update is given by
(3.3.4)			+ 4 n? (sTjsj)+e3A (ejaTj + 5 CT)
i=l			i i
k
where			A =			P )(y--HBs)
k=O
and P is defined by (3.2.4).
Lena 3.2.6 implies that the maximum elgenvalue of f is less than
1, since min >0. Therefore. by the Neumann Lemma
v			60			(5 5
i			ii
pk = (i--HP)?1
k=O
and A is the solution of
(3.3.5)			(1 --H p)A = (y-Bs)
Observe that the matrix (I--HI) has the same sparseness tructute
as V f. Although (1--H?) is not symmetric we can rewrite the update
so that the linear system is symmetric and positive definite. Again
assume that the zero components of s are ordered first, then
- 45 -
(3.3.6)
11+ = B +			(e 5T ? a CT) eT?
i i			i i			i
J=l
5 #0
where ? is the solution to
(3.3.7)			________
I			o			(y--HBs)
fi
0			0m			Am
and G			diag(?i??2			with ?i defined by
T1			s1Tisi;? if?=0SfiO
3.4 Properties of the update.
The first lemsa of this section will verify that the update has the
desired structure.
lemma 3.4.1 Let BrSflZ and y,sc?0 . Then B defined by (3.3.4)
satisfies Ii t SflZ(\Qy,5
Proof: Obviously B cS. To show B cZ. it is sufficient to checI:
+			+
that
e r Zi for each i. For I 1,2
(3.4.1) eTB+ = CT1B , 4(sT1s1)+eT1AsT1 + 4 ?? eTaj(535j)+eT??T
J=1			J
The transposes of the first and second terms on the right are in Z by
Lemma 3.2.3(c), the transpose oi the third term is too. Therefore, B s
To see that B tQ , form the vector B s and cIe?k co?p?ncnt--Hwise
+			y,s			+
--H 46 -
that eTiB+5			eT1y, i			1,2n. If St f
(3.4.2)			eTs s = eTBs + 4eTA			?T			T			+T T
I I I ? 41=1e1s?(s1Sj) e?ACjS
= eT1BS + ?T1(1?p??, from (3.2.4),
= ejY.
If s?= 0, then
T T
(3.4.3)			e1B s			e B & -? 0
I + i --
By the Mean Value Theorem [27?, there exist t1t (0,1) such that
eTj(g+ g) = eT1Jg(x +
ti5)5
Then eT + ?i???i = 0 because J rZ. Therefore B ?Q
1Y = eT1Jg(x			g			+			y,s
and the proof Is complete.
The following estImate will be used to prove the convergence
properties of the algorlthm given in section 3.5.
Theorem 3.4.2 Let B,J c Sr'Z and y,s r lIn with a ? 0, and let B
be defined by (3.3.4). Then
(3.4.4) IIB+?II2f <			2 II(f-?)+II2
IIB?JIIF --H			2
115112
55			n			2
+ max(2, max ?) ? (s1s1)+[e1(y=Js]
p ?0			p			1=1
J			i
Proof: Expand the expression for the Irobenius norm of B --H
- 47 --H
(3.4.5)
IfB+ --H Jil2			T
F			? le (? --H
t,i=1 t +
?? [eTj(B?J)e? + 4(5T5 )+eTAsT			T			T
I,j--H1			I I			I I J			I			j
= n? ([eTj(B?J)ej]2 + 2eT1(B?J)ej(5T15j) eT??T5
1,1=1			I			t
+ 4[(sT1s1)+(eT1?)(eTjsj) + 4(sTjsj)+cT1IcTj?1(%Tjsj)+cTjAeTjsj)
= IIB?JIIF2+ n? (sT1sj)+eTjAI2eT1(B?J)s+ 4eT?? ? ? (5T5 )+eTs eTs eTA).
i-I ?? ? ? ?
?ow observe that
y--HBs = (I--HP)A = 41 + 4? ?? (sTsJ)+eTjs sjeTjJA
J?1
and
(3.4.6)
T T nT +T T T
c1(y--HBs) = 4e11 + ? ?(sjSj) CjSCg%jCjA
i=1
= 4eT11 + 4 n? (sTisJ)+cTsLeTISJeTiA by (3.2.3)(b).
i=1
Therefore, (3.4.6) applied to (3.4.5) yields
(3.4.?) IIB+--HJ112 = IIBJII2F + n? (sT1s1)+eTl[2e3(?J)s ?T1(y??5)?
= IIB?JIIF2+ n? (sT1s1)+e34e3(B?J)s o+ ?T1(y??5)j
To complete the proof, examine the sum on the ?i?bt hand side of
--H 48 -
(3.4.7). Let u y - 35 and v = (B --H J)s. Then
A (l--HP)?1(y-Bs) = (l-P) (u--Hv),
and
(3.4.8)
n? (sTjs1)+eT1A[eT?(B?J)s + eTj?y?Js)]
= n? (5T15?)+et(1?P)?1(u=v)[eT1v etjui
t=1
= n? et1D+(l?P)?l(u?v) eT1(v t
t=l
(D+(l-P)?1u,v? ?D+(l-P)?1v,v) + ?D+(1--HP)?1u,u?
- ?D+(1?P)?lv.u)
Now, substituting (3.4.8) into (3.4.?) gives
(3.4.9)
1l8+??11F2 = 18?311F2 + 4D+(l--HP)?1u,v? --H ?I)+(i--HP)?1v,v?
+ ?D+(I-?)=1u,u? --H
= 18???'F2 --H ?D+(1?P)?iv,v) +
by Lemma 3.2.?,
2 vTv 5T5 ? T +T2
18?311F --H 5?5 + max(2, :a;oJAe2) ? ??i?i? (e1u)
i--H--HI
(B?J)2			?T5			?
2 ___________ + rsax(2, max A??) ? (sTisi)+[eTi(y?Js)J2.
- 18311F --H			1151122			plO
i			J
the inequality in (3.4.9) follows from Lemma 3.2.8 and since s?= 0
--H 49 -
implies			?			T T			i by
Civ			Cj(B --H J)S			ei(B --H J)61			0. The last?equality s
substituting for u and v, and the proof is complete.
Finally, 8+ defined by (3.3.4) is the solution to
min 1118 - 811			8rS?ZflQ			).
F			y,s
it is straightforward to apply the variational techniques [233. [223 to
this constrained problem. The solution by Toint [383 to the problem
posed in the variational form [313 is the same as the 8 in (3.3.t).
3.5 An algurithm.
Let x  ? , BC SflZ, positive definite be given.
I. Solve Bix --H8 for ?x, the step to be taken; g g(x)
2.			Set X+			x + 8x.
3. Evaluate g+ = g(x+); test for convergence.
4.			Set y			?+ --H
5. For each i such that leTiAx < Il???i4Il2 with ii > 2 fixed,
T
set Pj = 0; otherwise p			e ix . This defines 5, the step
i i
used to update 8.
6. Compute 8 from B,s, and y using (3.3.6) and (3.3.7).
The convergence theorem in the next section will show that 8
is positive definite in the region of local convergence. An implementa-
tion of the algorithm might avoid singularity of the liessian approxima-
tion by re--Hevaluating the llessian 8k= J (x ) or by adding a positive
g k
diagonal matrix to a singular or nearly singular 8k One possible
choice for starting the iteration would be to take 80 = 1.
- 50 -
The test in step 5 meJos that oc co?rposeot of As should tend to zeto
faster than its correspooding projection. Also, this ensures that the
cigenvalues of P, defined by (3.2.4), are not tending to I as k?--H, or
those of the matrix associated with the linear system (3.3.7) are not
tending to 0 as k?--H. This is critical for the `conditioning of the linear
systems involving A or ?, and it is also important in the convergence
proof.
The test is about what would happen anyway on the machine. When
the si?e of any component is less than machine precision times the
size of the corresponding projection, that component is insignificant.
I			I
In other words, for Ii = n?mach.eps. , the algorithm is roughly the one
carried out on the computer.
On the other hand, for small Ii, say Ii = 2, it Is likely that many
0 are set to zero. In that case, the correction to Ii doesn't take
?,ch work. This is close to the idea of keeping the same approximation
to che Hessian for seseral it=rations [151.
In step 6, the form of the update (3.3.6) is preferable to (3.3.4)
since the linear system (3.3.7) is symmetric. This, of course, requires
some bookkeeping to do the permutation on the vector s and the corre-
sponding permutations of ? and y. However, the important block in
that sym'setric system has the same sparseness structure aa the corre-
sponding block of the permuted Hessian. Also, this is a possibly small-
er system to solve.
3.6 Convergence.
Theorem 3.6.1 Let f: ?n It be C (00) for ?o open and convex;
assuTne there exists ?? t l3o such that Vf(x?) :.			??t?5el ta post.
tive definite a?d [V2f(x)]jj ? 0			for all ? ?			?t ts??e exists
K = (k1,k2,... ,k0) ? ?n with k			0 such that
(3.6.1) lIeT1(V2f(x) --H V2f(x+))ll2<kilfr - x?
for i = 1,2n and for all xc then there esist c,S > 0 such
that for (x0,B0), B0cZflS which satisfy
11B0 --H v2f(x+)II2 < I
and			11x0 --H x*112 < c
ten the syrmuetric Schuhert method generates (Bk) with Bk well-defined
and (xk? which converges linearly to x*.
Proof: Choose M > 2, and set c = IlK112. From Ihc?rern 3.4.2 with
J = = V2f(xt),
2			2			5T151 ,, [?T1(y?j*5)J2
(3.6.2) II11+?J+IIF?< liB-J*il + max(2, max			2			T
j 0 Pj I I
`I
For xi ? 0, an application of the Mean Value Theor? 12i\ gives
(3.6.3)			leTi(y?J?s)12			2			2
2			--H< k1o(x,x+s)
115112
where o(x,x+s) = max(llx--Hx*ll2, IIx+s --H x*l 12) Therefore,
T			T			2
n [e (y?J*5)?			2			2
(3.6.4)			max(2, max			2			?			< M `(x.x.s)
#0			tj			si=#lo			SiSi
52 --H
and this with (3.6.2) implies
2 2			2			2
(3.6.5)			18+?3*11F			--H< llB=J+llF + M?o(x,s+s)
and
(3.6.6)			18+??*11F			--H< llB?Jall? + aa(x,x+s)
where a = M4.
Now,
(3.6.7)
so that
(3.6.8)
I Ix+s --H xal 12<115112 + lix --H x+l12
< ?II?xII2 + IIxx*112
< 4L0H[IIx?x*iI2 + IIX+?XtIi2] + IIx?x*II2
--H< (l+FN0)IIx?xtII2+ t?nIIx+?xaII?
a(x?x+s) < (l+FN0)IIa?x*Ii2 + F0M1 lx+?x*I 12
<?(l+?) o(x,x+).
Thus,
(3.6.9)			llB+?J*ll, < 18??*11F + &o(x
where a = ?1?tN0? An application of Theorem 1.5.1 completes the proof.
In the preceding proof observe that since lxii < Ilixil, F?n could
he replaced by 1, regardless of ii. In general, though. by appropriate
choice of N, F;;  I.
4. Natrix Factorizationa
4.1 Basic considerations.
The linear system.
(4.1.1)			Ax			b
with nonsingular A? L(?n) and bc lin is solved by Gaussian elim-
ination. The algorithas most frequently used are the Crout and Dolittle
decompositions with partial pivoting. The matrix A is factored into
a product of lower and upper triangular matrices:
(4.1.2)			PA			LU
where P is the nonsingular pivot matrix. In the case where no pivot
ing is required, the linear system becomes
(4.1.3)			LUx = b.
Then the solution of (4.1.3) is accomplished by aolving two triangular
linear systems. First Lw = b is solved for w by forward substitu--H
tion, then Ux = v is solved by backward substitution for x.
For rank--Hone and rank-two quasi-Newton methods, techniques are
available [I]. [201 for updating factors of the approximate Jacobian.
Brodlie, Courlay, and Creensradt [11 show how to express rank--Hone and
rank--Htwo corrections to symmetric positive definite matrices which re--H
suIt in positive definite matrices in product form. This is applied
to quasi--HNewton methods for unconstrained minimization. The important
fact is that the numerical stahility of aigorlthms for those problems
is increased. This is because the approximate Jacobian is 8=aranteed
--H 53 --H
to be positive definite, prodded the iteration hegan with a positive
definite matrix; the effect of roundingon updating formulae
which theoretically maintain positive definiteness is minimized. Gill
and Hurray L20] propose application of their techniques to the solution
of nonlinear equations in addition to modifying the factorization of the
ilessian to preserve positive definiteness. In the former case, they
tecosmend using ao orthogooal factorization of the approximate Jacobian
on the basis of numerical stability.
Gill et al. [i?] describe in general, modifying Cholesky factors
and QR lactors of a matrix following a rank--Hone change. in a later
paper, Gill et al. ?i] give methods to compute the LDV factoriazation
of a matrix Ii and the modified factors after the deletion or addi-
tion of a column to f.
liroyden ?3? proposed dealing with just the factors of a sparse
Jacobian. However, the Schubert update formula adds 0 rank--Hone correc-
tions and the syxaxetric version adds n rank?two corrections. There
seems to be little hope for a simple method to update factors of the
Schubert or the sytmwtric Schubert methods.
Gill et al. ?? note that the factors of a sparse matrix are
more likely to be sparse than its inverse.' Iteeping in mind the idea
that the Schubert method involves so ving sparse linear systems at
each Iteration, and these systems have the same structure, new algor--H
ithms are suggested. These algorlthms are based on an LU factorization
of the approximate Jacobian. The Schubert update is applied to one of
the factors, and preserves the triangular structure along with any addi-
tional sparseness structure. Several algorithms are proposed using the
Crout and Dolittle factorlitations. While we have done no numerical
--H 55 --H
testing on these algorithms, More (26] has implemc -:e. `s?
on ?pdating ? of the QR decomposition, and has hta
preliminary numerical results.
Let A (a ). The following algorith gi?t c?? t'?to,i.'tt0n
ii
(4.1.2) where L is lower triangular and U is n=tC n?ge' t'ia?$ular.
The pivots are stored in the noosingular mains T.
Al&orithm 4.1.1 Crout factorization [42].
For r			1,2			n, do steps I through `.
r--Hl
I.			Compute t?			= a			--H ? ilkUke			i = r,r'1
r			ir			k=I
2. Set int > r equal to the smallest k 10' ht.h
ILke1 = max
rJ??<n
th			th
3. Interchange the r and tnt rowa of the a',aV. Itlet to these
rows hy their new positions.
r--Hl
4. Compute urj = (arj --H ??1trkukj)1ttr , In, I - ,.l
The convention ?			0 is used.
The algorithm without steps 2 and 3 is calledf?, totization
without pivoting. The pivoting strategy of step'. ? l I Is Ih class1--H
Cal one and considered closed, in that there is nu ll..?lhllIty in the
choice of pivots.
Instead of the closed pivoting of the standard I,r I).?lItrle
factorications, an open pivoting strategy should he			1. lulf [it]
calls this "threshold pivoting". As Duff notes. the las'.t'al selection
of pivots was designed to maintain numerical stahillty However that
--H 56 -
selection may destroy the sparsity in the factorization. Threshold
pivoting allows selection ol noneero pivots which are at least l!u
times the largest element in the row or coluwn. Erisman and ileid [l7?
give a formula to monitor the growth of the matrix elements in the
factori?ation. This can be done working only with the sparsity pattern
and the fa?tor u. If the growth is too large, a smaller value of u
can be chosen. Puff also note that u 10 yields good rerention of
sparsity and generally good numerical accuracy. lie also recomwends
iterative refinement to improve the solution sin?e solutions for systems
with sparse triangular factorizations are fairly cheap.
4.2 Some algorithms.
Consider again problem (1.1.1). The Jacobi?,, or approxis?te
Jacobian, is decomposed into triangular factors by either the Crout
or Dolittle tactorizations. The unit triangular factor is held fixed
for a finite number of iterations-; the other factor is updated using
a Schubert correction to preserve the triangular structure and any
other zero-noozero structure of that factor. (if that unit triangular
factor were fixed for the entire iterative process, a term involving
1ix0 --H x*l 12 would occur as a linear growth, denying the bounded de-
terioration estimate, and the convergence proof cannot be obtained
using this type of analysis.)
The algorithms are presented without pivoting for simplicity of
the exposition. The implementations would e?ploy an open pivoting
strategy as discussed in the previous section. The convergence results
in section 4.4 are local and will be based on using a fixed pivoting
strategy throughout the iteration. The initial choice of pivots is
--H 5? --H
further discussed at the end of section 4.4.
The first algorithm uses the Dolittle factorization to decompose
the JacobIan into a unit lower triangular matrix L, end an upper trian-
gular mjtrix, U.
Algorithm 4.2.1
I. Given x0, the initial guess to the solution of F(x) = 0
m a fixed positive integer and
? > 0, a fixed constant
2.			Evaluate F0= F(x0) and J0-- 3F?'o?
3.			?erform a Dolittle factorization of J0 to obtain J0 L0U0
where			L0 is unit lower trian???lar and			U0 is upper triangular.
Set k			0.
4.			Solve			Low? = ?Fk and Uksk = for			wk and 5k respectively.
5. If = = m-l set x0 = 5k + 5k and go to 2;
else set xk+l = Xk + 5k
Fk+l =
and			= Fk+l - 1k
6.			Solve L0vk =
7.			Update the jth row of Uk			if 115k112 --H< 811?5k?j112
Let Xf(Sk) e (J: 115k112			--H< 611?5k?j112?
(4.2.1) Uk+l = Uk +			I			F(s?)Tj(sk)j?+eTj(v??U?s?) Cj?5k)T
18?5k?			i
Replace k by k+l and go to step 4.
Another technique for getting started in step 2 would be to make
--H 58 -
a finicc, difference approxi?acion ?6] to i?(xo) i?stcad of evaluating
all of the partial derivatives.
The other basic algorith is to use the Cruut factoti?atlon in place
of the Dolittle factorization.
Algotith 4.2.2
1.			Civen			x0 c Ii			?			a fixed positive integer, and f > 0 fixed.
2.			Evaluate F0 =			F(x0) and J0 JF(Xo)
Set k			= 0.
3.			Factor			J0 = L0U0			using the Crout decowposition.
4.			Solve			L?v? =			ik			and 1i05k wk
5. If k in-i, set ? xk 5k and go to step 2;
else xk+l = Xk 5k'
Fk+l = F(x?+1)
and			?k = Fk+l
6. Update the 1th row of Lk ii IIWkII2 < lII(w?)iIi?
Let xe(wk) = (j: IIwkIi2 ? BII(W?)jII?)
T+T			T
(4.2.2) Lk.l			Lk ?			?c			[(w?)(w?)1e(y??L?Uos?) e (wk)
x (v )iii			j			J
1k
Replace k by k.l and go to step 4.
Ther? is one less triang?iar linear system to le solved at each
iteracion in this algorithin. Cit algocithm 4.2.1 at step 6.
The trade?off is an additional matrix-vector multiplication In updat-
ing the triangular factor. implicit in these algori thms is the def ini--H
tion of the projections (.) . in algorithm 4.2.1, the projection
th
is onto the zero--Hnozero structure of the j row of the upper trian--H
59
glat factor. in algorichum 4.2.2, the projection is onto the zero--H
th
no??zero strucof the J row of the lower triangular factor.
The constant i is used to prevent a correction to a row when
the projected step along that row is too small relative to the full
step. Co?pate with N used in the algorithm in section 3.5. Also note
that this technique would he used in implementing the ordinary Schubert
update (2.2.h) or (2.2.7)
Other algorithos are readily suggested. After the initialization
of J0 by either direct evaluation or formIng a finite difference
approximation, a Schubert (or Broyden) update can be performed every
rn iterations. That is, form the product
Bm?i = Lottm?i
and obtain B from B?1, Sm and by the Schubert update, then
???????
11 =LU
by the tolittle decomposition. Similarly, algorithm 4.2.2 could be
modified.
Another modification in the basic algorithms would be to monitor
the growth ot the elements in the updated triangular factor and to re-
evaluate the Jacohia, if it becomes poorly conditioned. This would
supplem??t the cciterlno to re-evaluate every m iterations.
4.3 Continuity of tie factotfeation.
The coovergene properties for the algorithms will depend on the
continuity of the factorization. The following theorem establtshes
60
when the Crout factoriation without pivoting exists.
Theorem 4.3.1 [133 Let A rL(?0) be nonsingular. Then the following
are equivalent:
(a) Crout reduction without pivoting can he carried out;
(b) There is a unique nonsingular lower triangular matrix, L,
and a unique unit upper triangular matrix, U, such that A LU.
(c) All the leading principal submatrices of A are nonsingular.
A similar theorem holds for t-he tolittle iactori?ation.
The next theorem gives sufficient conditions for tl?e LU decompo
sition to be continuous in a neighborhood of a point where the Crout
factorization without plvotio+? is known to exist.
Theorem 4.3.2 Let A: lin 110 and suppose for some x e 110 that
0
A(x0) is nonsingular and that there exists t0 > 0 and > 0 such
that
(4.3.1)			Iei[A(x) --H A(Y)JejI < ?ij1? --H y112
for i,J = 1,20, and for all x,y t B(x0,e0). If the Crout factor--H
ization without pivoting exists at x0
A(x0)			L(x0)U(x0)
then there exists t > 0 such that the Crout factorization without
pivoting exists for all xrll(x0.e)
Eurthermore, there exist constants c0 and d0 > 0 such that
- 61 -
(4.3.2)(a)			IIL(x)			--H L(x0) 1tF <c0Ifr			--H x0l12
and
(4.3.2)(b)			IIU(x)			- U(x0) 11F ?d0ifr			--H xo112
for all a ? B(x0.?)
Proof: Since A(x0) is nonsingular and the Crout factorization Ithout
pivoting exists at xo , it follo=s frow Theorem 4.3.1 that the leading
principal submatrices of A(x0) kA(x0) k 1,2n, must be non-
singular. By (4.3.1), A(x) and all of its submatrices are Lip?hltr
continuous in B(x0, t0). Therefore by the Inverse Function Theorew ?7J
I' 2, such that kA(x)
there exist t ,			n			is nonsingular in 1i(x0``k
Let t win r Then all the leading prin?Lpal sutwatrices ut A(x)
l<t<n k
are nonsingular in B(x0,t). Another application of Theorem 4.3.1
establishes the existence of the Crout factorization "ithout pivoting
for all x ? B(x0,t)
A(x) = L(x)U(x).
Bestrict t so that the preceding factorization holds on B(x0,t).
The prool of (4.3.2) is by induction on the dimension [4oj. Let
a tB(a0,t) and A(x) = L(x)U(x).
For L 1,
LA(x) = I
A(x) = a11(x)
a11(x) = a11(x)l 111(x)u11(a) = 1L(x)1U(x)
and
62
Thus,
I I 1L(x) --H L(x0) II?			Ia11(x) --H a11(x0)
< ?lii Ix --H Xo I 12
c1lIx --H x0112
by (4.3.1), and
II			U(x) --H 1U(x0) 11F			--H I = 0 Ix --H X0II2
By the Inverse ?unction Theoree, 1L?1 and 1U?1 are Llpschitz
continuous in B(x0,t). Let c and a be the respective l'ipschitz
I			I
constants.
Assurne for t			2,3k, k c n--Hl, that
(4.3.3)			J?L(x)			--H ?L(Ko)IIF <CfI Ix			--H x0112
and			IItu(x)			--H U(xo)IIF < dtIIx			--H x0II2
t --HI			t --Hi
for all x t B(x ,r). Furthereore, the inverses, L (x) and U (x) are
0
Lipschitz continuous on B(x0,r) with Llpschitz constants c and a
L			L
respectively.
For L = k+l, let k+lA(x) he written as
k			k
(4.3.4)			k+lA(x)			v(x)
kwT(x)			ak+l,k+l
k k			k			k+l
Since			A(x) = L(x) 0(x)-			A(x) can be factored
kSlA(x)			kelL(x)k+I0(x)
63
0			U(x)			L (x) v(x)
where ?k+l = 8k+l,k+l			--H k T (x)kU?l (x) kL?l (x)kv(x)
Now,
(4.3.5)			11k+lL(x) --H k+1
L(x0)111
kL(x)			I			o			kL(x )
-			0
I			k T			k --HI
kwT (x)kU?l (x)			Gk+l (x)			w (x0) u (x0)
11kL(x) --H kL(xO)IfF + IIkwT(x)kU?l (x) --H kwT(x0)kU?l(x0)II2
+ a			(x) --H a			(x )I
k+l			k?l 0
An application of the trianglc Inequality to the second term yields
(4.3.6) IIkwT(x)kU?l(x) --H kwT(x0)ku?1(x0)If2
< lI[kwT(x) --H kwT(x0)?kU?l(x)1I2 + IIkwT(x0)[ku?l(x) --H ku?l(x0)iII2
11kw(x) --H kw(x0) II2IIku?l(x)II2 + lIkw(x0)II2IIkU?l(x) --H ku?l (xo) Iii
< ? k2			?			ki			k
I=?l ?k,i? Ifr-x011211 ti (x)112 + II w(xO)II2?kIfr?xOII2
k 2			4			ku?l			k
--H< ?? ? ?k j) Sup (x)112 + II w(xo)II2???Ifr?xoII2
1=1 `			x f B(x0 ,?)
ik lix --H x01l2
- 64
Similarly the third term on the right hand side of (4.3.5) can be
estimated:
(4.3.7)
jCk+l (x) --H a?+i(xo)I
< Iak.l?k+l			--H ak+l.??l(xo)I
+ IkwT(x)ku?l (x)kL?l (x)kv(x) kwT(x0)ku?1 (x0)kL?l(x0)kV(x0)I
--H< ?k+l,k+l I Ix-x0t12 + J(kwT(x) kwT(xo))ku?l(x)kL?l(x)kv(x)I
+ IkwT(x0)[kU?i(x) --H ku?l(x0)]kL?l(x) k
+ I kwT(xo)?u?i(x0)[kL?l(x) --H kL?l(x0)?kv(x)I
+ IkwT(x0)kU?l (x )kL?i(x )( kv(x) --H k
0			0 v(xo))I
< ??k+l.k+l + ? ?2			)4			supIIku?i(x)II2IIkL?l(x)II2IIkv(x)II2
?? k+l,i			xe B(x0 ,r)
+ IIk?(x0)II2 ?? :u:n?xot?IIkLl(x)IzIIkv(x)II2
+ IIkw(xO)II2IIkU?l(xO)II2?ksup			IIkv(x)II
xc B(x0 ,c)			2
+ IIkw(xO)II2IIkU?l(xo)II2IIkL?l(x0)1I2( k
?2i k+l)?)IIx?x0II2
- ?k IIx--Hx01I2
Set Ck+l = Ck +			+
To complete the proof, observe that
(4.3.8)			11kV(x) --H k
U(xo)IIF
kU(x) f kL?l (x)kv(x)
0
ku(x0) j kL?l(x0)kv(x0)
11ku(s) --H ku(xo)IIF + llkL?l (x)kv(x) --H kL?i (x0)kv(xo)IIF
?<dk lIxx?II?+ llkL?l(x)[kv(x)?kv(x0)Jl1 + ll?kL?l(x)?kL?l(x0)]
v(s0) 12
< (dk + 5S????50 ,c) llkL?l(s)ll2 ?i?1?2I?k+l?? + c?llkv(xo)ll2llls?soIl2
+ d			lxx II
k+1			0 2
k?l --HI			k+l -l
By the Inverse Function Theorem, L (x) and u (x) are
Lipxchitz continuous On 11(50 .t) with Lipschit constants c and
k.l
??,1 respectively. The induction is complete. Define c			c			and
0			n
d0 dn to prove (4.3.2).
A similar theorem holds for the Dolittle decomposition.
Theorem 4.3.3 Let A:lt?. ?nxn ad suppose for some			 ?n			that A(s0)
is nonsingular and that there exists 0 > 0 and			0			such that
(4.3.9)			leTi[A(s) --H A(Y)]ejl < ?ijIls --H y1l2
for i,J 1,2n, and for all x,yB(s0,r0). If the Dolittle fec--H
torization without pivoting exists at ?o' A(s0) L(x0)u(x0) then
there exists  > 0 such that the folittle factorization without pivot-
ing exists for all x B(x0, ).
66 --H
furthermore, there exist constants c0 and d0 > 0 such that
(4.3.iO)(a) IIL(x) - L(x0)ll? <c0 lix --H x0ii2
and
llu(x) --H U(xo)ilF ? d0 lix - x0ll2
tor all xc B(x0,r)
Proof: The proof is very similar to the proof of Theorem 4.3.2.
4.4 Convergence.
Before stating the convergence theorems, several technical lemmas
will be given. The first two lemmas relate the updated triangnlar
matrix factors to the preceding ones.
Lemna 4.4.1 Let f:1f If and suppose there exists x* such that
0 and J* J(x*) is nonsingular. Suppose there exists c > 0
and y ) 0 such that
(4.4.1)			leT1[J(x) - J+]ejl < Yij lix --H x*li2
for i,j			1,2n, and for all xc
Then there exists co ) 0 and a pivoting strategy, F0, such that
the Crout factorization without pivoting of fo? exists at x0c B(x*,c0):
(4.4.2)			P0J(x0)			L(x0)iJ(x0)
m-i
Furthermore, (L?1?? defined by (4.2.2) satisfy
2 2			2			2
(4.4.3) I ILini--HL*ii ? 5 I iLt?L+ii ? + i n(cn + I iL*i I2d0I boll 121 Ix0--Hx+I 12)
- 67 -
where			jrJJF and F			(? ), and n?			max(iix			--Hx'I
ii			???			2,Itx,--Hs?II2)
Proof The proof Is very si?ilar to that of the following les=a.
Le?na 4.4.2 Let f:?n ?? and suppose there exists x* such that
= 0 and J* = J(x*) is nonsingular. Suppose there exists c>O
and ?ij > 0 such that
(4.4.4)			leTi [J(5) - J*?eiI<?iiIIx --H x?II?
for i,i = 1,2n, and for all xc
Then there exists CO > 0 and a pivoting strategy, Po such that
the tolittie factorization without pivoting of P0J exists at x0e
(4.4.5)			P0J(x0) = L(x0)U(x0)
Further?oce, (Ui1m--Hl			by (4.2.1) satisfy
defined
2			2
(4.4.6)			11u5+1 U*IIF ? 11u1 ??? F
ii?01ii? n12[eni + IIu+I?2c0I?x0x+II2
where I >0 is set In step I of the algorithm 4.2.1, e 0 jirli
r = (? ), and o? = max(IIx --H x+li2' IIXL --H x?II2?
ii tol
Proof: For the first assertion, let P0 be any pivoting strategy such
that the Dolittle factoriration of
P03(x*)			L+u*
exists. Apply Theorem 4.3.3 to A(x) = P 3(x) to obtain to such that
-o
the Dolittle factori?ation without pivoting of
A(x) = L(x)U(x)
--H 68 --H
exists for all x ? B(x0.?0?. That theorem also gives a c0 > 0 such
IIL(x) --H L*ilF <Coils --H xt 112
for all
in particular, if ?o E 8(s*,c0), then the Dolittle factorization
P0J(x0) L0U0 exists. To prove the second assertion, evaluate
1lU1+1 --H U*11F2 using the update (4.2.1). Set X8(St) = (1:1 Is Ii
1 2--H
o+lI(sL?i11 21 then
(4.4.?)
IIU1+1 --H U?Il2
- ilU1 --H U* + ? [(51)T1(55)1?+eTj(L0iy1 --H U1st? e1(51)T111?
X8(si)
T			T			+			T 2
--H			?(lle (U			U*)(i?[(5 ) (s)i(s)(s)ll
I t			1 iLitii2
i8(s1)
+ lIeTi(Loly? --H U*sl)[(sl)Tl(sl)iJ+(s1)T1121
+ ?lleTj(U1 --H			2
U*)112
i8(s1)
2 T +T I
--H< ??l???? I + ? ll[(s11??t?iJ e1L0 (y1--HL*U*s1
?8??1?
+ L*U*s1 --H LOU*sl)(si)TiIl2
< llU?U*1l12 + I [11(51)T1(51)1?+eTjL0l(y1--HJ*s1)(51)T1112
?i??
+ 1[(5T1)1(51)1?+eT1L0I(L*?L0)U*s1(s1)T1112J2
--H t9
IIULIJ*1IF2 + ? ?IIL01 iI?(II(s?)iII?IIy?Jta?II?
xli?5t?
+ 1?5t?i1 2IIL*--HL0 2 ??*??211stII2)]2
? IIUi?U*IIF2 + iLjlii228??Y1i?ot+ IIu*112c011x0 --H x*i12)2
< IIl1i?U*IIF2 + 1iL0111282n(roi. IIu?iI2c0IIx0--Hx*II2)2
giving (4.4.6).
There are several norws which can he used for the error estimates.
The mostient weighted norm is the one established at the start of
the iteration.
Definition 4.4.3 Define the weighted norm II?II?0l.f:L(? ) ? by
(4.4.8)			1?11L?1 F			IL01 .11,
0'
for L0 nonsingular lower triangular. This is just the left weighted
Frobenius norm.
Lemma 4.4.4 Let L.?n ?n he a continuous lower triangular operator
-I			n
Suppose L exists and is continuous on B(a*,t) for some x* ?
and t > 0. Then there exist constants n and n > 0 such that
(4.4.9)			11.1 L?1(x),F --H< ni 1F			for all			xc B(x*,t) and
(4.4.10)			I I. 111 < nI 1.1 L?1(x),f			for all			xc
Proof: Let n			sup			IIL?1(x)112 and =n = sup			IIL(x)112
x c B(x*,t)			x t 8(x*.t)
70
Then
L?1(x),F = IIL?1(x) `11F--H< IIL?1(x)II2II.I			<??11?11F
and
=			I
F			x--Hl			I
iL(a)L (x)			2o+L(x),F --H<
Let n >0 he such that the matrix norms satisfy 1,11F < ?II?II2
Theorem 4.4.5 Let F:IIn lin be continuously differentiable in an open
convex set 00 Suppose there exists x* ? l7o such that F(x*) = 0 and
J(x*) is nonsingular. Suppose there exists r = ??ij? with
?ii ?70
such that
(4.4.11)			leTi?J(x) --H J(?)ie?1 <			jx --H y;f2
tor i,j			1,2n, and for all x,y r
?hen there exist ?,i > 0 and a pivoting strategy, P0, such that
if 11x0--Hx+II2			? and			?J --HJ*j			< I then the algocithm 4.2.1 gener
0			1F
ates (xkl which converges locally and supcrlin?arlv t x+.
Proof: Let c1 >0 such that B(x*,?1)C170. Set y 7 jj??1 ? and
define x IltilF
By Lemma 4.4.2, there exists > 0 and a pivoting strategy, P0,
such that the Dolittie factorization of P0J(x) exists for all
xc B(at,c ). Furthermore, (4.3.ll)(a) can now he written as
0
(4.4.12) IlL(s) --H L(x+)llF <c0Ifr --H x*1l2
for all x e B(x*
--H JI
1'ow choose r in (0??o)			I > 0 and r r
or			<I
2			1--Hr
(4.4.13)
+ 2n??I3 ? r
where o2 = ol + (?l)o and a --H nIlu*ll c fnn8?
o			o			20
and n			sup			llL(x)1l2 , n			sup			I.
x r 1i(x*,c0)			xc
Further restrict c so that llJ(x) --H
llx			--H x*l			< r. Let x0cB(x*,t), then
(4.4.14)			llL(x ) --H L(x+
o ?11F --H< c0lIx - s?If
by (4.4.12).
If ji0 --H ?*11F < I then 1lj0 --H i*ll2 < nA .?. ud (4.4.13)
implies that 3(l+r)2n1 < r. Therefore the Dana,'
(4.4.19)			I li?1			(l+r)y
0			112<
The iteration of Algoricha 4.2.1 can be sfe?e.'			i.,ii???			Let
p = im, I = 0,1and kr (p,p+I(mIll
(4.4.16)			`(k+1 = xk --H UIc1 L?1 Fk
with IL 1 and (u 1 determined in step 5, and iu? I"' ??			d?termined
p			p			?-l,'
in step ?. Then,
--H1 --H1			--H1 --Hl
xk+l?x*			xk?x*?U? Lp (Fk?F*) + Uk Lp J*(x??X+)			?? L?, Jt(X??xt)
and
72 --H
(4.4.17) IIXk+l?x*II2 < IIUk1Lp1!J2?IIF??F*?Jt(x?=x*)II2
+ IILptlk=JtII2IIxk?x*II2]
For p = 0, k			0, (4.4.17) hecomes
(4.4.18) IIx1=x*112< IIJj1II2[(K? + 2?i)IIx0x*II2?
< (ltr)?(=r + 2?i) I Ix0--Ha*l 12
=< rIIxo?x*II2
applying (4.4.1?) and (4.4.13).
Assume for p			Ia, I			0,1J			that
(4.4.19)			Illim --H x*II2 < rIfrjm?i			--H X*I 12
Then x			?B(X*,C)			and it foliows that			--H ?*IIF < I, and
11?jm --H J?II2 < 2n1. From (4.4.13), 2y(l+r)nI < r, and the Banach
Leawa ?7? gives IIJj?II2?< (l+r)?.
For p = Ja, let k			p, and (4.4.17) Is
(4.4.20) IIXiw+1 --H X*II2 < IIJ;m1II2[(e? + 2fl1)IIXim --H X*I I2?
< (1+r)?(?r + 2n1)IIx --H X*II
< rIIXja --H X*I I2
Assume for L p,p+lk--HI, k < p+(?1) that
IIXt+I --H x*jJ2 ? rIIx1 --H X*I 12
Now,
73 --H
(4.4.21) IlL UL+l?J*llL?l ,? llLpU1+1?L U*+LpU*?L*U*llL?1 F
p			p
<			llut+lu*llF			+			llLp1ll2llLpL*llFllU*1l2
< llUt?U*llF + lILp1Il2 \FniE?o1
+ llu*ll2c0 llXp?X*ll 2?
+ llL?1 ll2llu*ll2llL?L*ll
p			p
by (4.4.6) and the fact that for a,h,c > 0, c2 < a2 + b2 ? (a + b)2
implies c <a + b. Since lILp1Il2 <??
(4.4.22) IlL Ut+l?J+llL?l I < IlL Ui?L*U+llL?i I
p			p
p p
+ 2nllu*ll2 llLp?L+llF + m1o, + fl?illU*Jl2CoIlXL??*ll2
< IlLPUL L*U*llL?i I + m1o1 + Gollxt			x*l 12
Therefore,
(4.4.23)			IL u			--H L*u*ll --Hl			< IlL u --H L*u?lI --H1			+
p 1+1			L			I --H			p p			L			I			ol			? Cj
p			p			ip
+ m0(t+l)llx --H x*l
p			12
<a2 i=?p01 + IlL u --H L*U*llL?l
p p			I
and for I			k--Hl,
(4.4.24) IlLpUk?L*?*llL?1 I --H< IILpUp?L*U*llL?l,F + a2 k-i
p p i=p
74
tn-I
l-r			p
IL u --H L*u*ll I			+
p p			L ,F			2			l-r
< ni+ i <			2ni
Thus			IILpUt+1?JtII2 < nIIL?ut+1?J*lIi < fl?IILpUt?l?3*11L?1,? and this
gives			llUk1Lp1ll2 < (1+r)y			which with (4.4.17) yields
IIxk.l?xt 112 < (l?r)?[cc + 2nn?i]I lxk --H x+1 12
< rllskx*Il2 by (4.4.13).
For t p+(ml), (4.4.17) is
lIx??? --H x+ll2 < clix			- x*ll
--H			p+(??l)			2
so x			r B(xt,r) and			lj			3*11 < 6. This establishes the local
p?m			p+rn			F
and linear convergence of Algoritha 4.2.1. Superlinear convergence
tnllo? fro Ticocen l.?.S sic?c for p tTn, i 0,1
iI(Jp--H3+)(5pMff??ffiffi2			___
lIx??1 --H Xp112			< lIJp --H 3*112 ???` 0
The convergence theorem for algocithrn 4.2.2 is proved in a similar
fashion.
-Thence 4.4.6 Let F:1f Ii be continuously differentiable in an
open convex set to Suppose there exists x+ c D0 such that F(x?) = 0
and 3(x*) is nonsingular. Suppose thece exists t = (y ) with
ii
0 such that
--H 75 --H
(4.4.25)
IeT1[J(a) --H l(y)]e?I < ?ii'?5 --H yI12
for i,J			1,2n, for all x,y D0.
Then there exist r,i > 0 and a pivoting strategy, P0. such
that if 1ix0 --H x?II2 t and II?0 - J*IIF < i then algorithm 4.2.2
generates 15kI which converges locally and superlinearly to x*.
For the variations discussed following algorithms 4.2.1 and 4.2.2,
convergence theorems can also be proved. in the case where a Schubert
update is made every a iterations, the bounded deterioration estimate
for the Schubert update would be used every a steps in the inductive
proof of local cgence. Superlinear convergence is conjectured.
in the convergence proofs, the pivoting strategy was fixed through--H
our the iterations. flue to the nature of the local analysis, it was
chosen in such a way that the factori?ation p0J(x*) was well--Hdefined.
in practice, this is not feasible since xt and J(xt) are not known,
The pivoting strategy, p0 , would be chosen so that p0J(x0) had the
appropriate LU factorization. For x0 sufficiently close to the
solution. P0A(x+) will have all principal leading subastrices nonsingu-
lar.
5. Further applications and conclusion
5.1 Further applications.
Curtis, Powell and Reid [t? examine the problem of estimating a
Jacobian, J. by finite differences when J is sparse with a known
n			n			n
sparsity pattern. Let F:1R			?			have component functions ?i It			It
so that
ii(5
f2(x:
F(x) =
f (x
Then ?F - ??xfji)
The finite difference Newton method uses the approximation
(5.1.1)			?f?			f1(x+h1ei?--H?i(x)
hj
with the step length h This can also be written as a vector approxi
`nation to the ?th column of J:
(5.1.2)			Je?			h\?F(x + hf1) --H
The disadvantage here is that n additional function evaluations are
required at each iteration.
Curtis et al. suggested grouping columns that are mutually orthogonal
to allow approximating several col?mns at the same time by one differ-
ence. The groupings are accompli?ied by inspecting the columns and group--H
ing together those which have no unknowns in common with any other col--H
--H 76 -
--H 7? --H
umn in the group, Sparsity is what makes this feasible. This tech-
nique of approsimating the Jacobian reduces the numbe? ?t tuction
evaluations compared to the finite difference Newton me:kd. fle
choice of step length is discussed in Curtis and Reid ?
Tv above problvm is con?idvred sp?tse n [b1
wri?ten in the form
F(x)			o+(x) + Nx
where N is a matrix of known derivatives and J0 is o+%pd??e.
Consider now a Jacohian of the form
(5.1.3)			J(x) = C(x)			A(x)
where C(x) is easily computed and A(x) is the part ot the Jatobien
to be approximated. Dennis 112] had suggested approximating A(s)
using its structure as follows:
Given 5k' xk+l and 1tk approximating J(xk)
set 5k =			--H
k+1			k
= ?(x?+1) --H
Compute Ck?l
Define ?+? approximating A(sk+l) to satisfy
(3.1.4)			Ak+15k =			--H Ck+lsk
This is a quasi--HNewton equation for Ak+l, Now tha structure ci A can
be preserved by an appropriate updating formula.
In particular, if A(x) is known to be sparse, Schubert's method
can be used to update Ak to obtain a sparse Ak+i satlsfyin? (5.1.4).
78
Based on the convergence properties of Schubert's method, local and
auperlinear convergence for such an algorithm is conjectured.
3.5 Conclusions.
Ciii and Hurray ?O? comment that advantage should be taken of
special structure. The structure of a problem can be important in de-
veloping efficient methods for solution.
Schubert's method takes advantage of sparsity in solving nonlin--H
ear equations. An approximate Jacohian can be updated row--Hwise to
preserve the sparseness. The advantage here is the reduced storage
requirement over a non--Hsparse algorithm. There is less work; the work
to compute the new approximate Jarobian is the order of th
ol the Jacobjan. At each iteration in the algorithm, a sparse linear
system must be solved. These systems all have the same 5re.
The preprocessing can be done just once U6] and then applied to each
linear system. Schubert's method exhibits the local and superllncar
convergence properties that Broyden's method has. The peoofs of t-hese
properties depended on finding the bounded deterioration estimate (2.3.2).
Also. a miniproperty is satisfied by Schubert's update.
?he symmetric Schubert update can be derived in two ways. One
is by an iterative projection technique Pil The other i by varia--H
tional principles [23? which also establish a minimum norm property for
this update. As with the Schubert method, the advantage of reduced
storage is present. At each iteration of the symmetric Schubert algor--H
Ithm, two sparse symmetric linear systems must be solved. All of these
systems have the same structure. Once a pivoting strategy is chosen ??,
the preprocessing can be applied to all subsequent linear systems.
- 79 -
The algorithm is shown to he locally and linearly converget
to a solution. Again the proof depended on obtaining a bounded detet
loration estimate (3.f.2). Superlinear convergence of the symnetric
Schubert method is conjectured.
Natrix factorizations are suggested as a means of solving sparse
nonlinear equations. The underlying update is Schubert's update which
is used to preserve the structure in a triangular factor of an approxi-
mate Jacohian. In the factored form, solution of the linear systems
are easily performed. The convergence of the algorlthms is contingent
on the continuity of the factorization. Local and suparltnear--H
gence are proved for factorirations based on the Crout and Dolittle
decomposittons.
A further application of the sparse techniques is suggested for
problems in which the Jacobian can be written as J(x) = C(a)
Part of the Jacohian is easily computed, C(x); part of the Jacobi?n is
to be approximated. When A(s) is sparse, an algorithm is suggested
and superlinear convergence is coniectured.
The convergence results depend on the founded Deterioration Ticorem
[i]. No nu?erical results are presented here. Schubert's method has
received some testing with good results [3?. Some testing has als teen
done recently on the sparse symmetric updating formula [?8],[39]. Tuint
claims a superlinear convergence property for that algorithm, however
it is not presented. More' has an indication from his testing of a
factoriration that the factored update methods perform well [??.
Several questions have been raised. The next step is to implement
and test the algorithms and verify the theoretical convergence proper-
ties discussed here. Sparse nonlinear algorithma hava widespread applica--H
- 80 -
btllty since most large scale nonlinear problems are likely to be
sparse. with the interest in sparse linear problems, codes are avail--H
able which will make it easier to handle the nonlinear problems using
tbe aigorlthrns presented and discussed here.
Bibliography
[1] K.w. Brodlie, A.R. Gourlay. and J, Creenstadt, ?a?i ,`?e ad rank--H
t,,o corrections to positive definite matricex expre?,c ?, ?ruduct form,
3. lost. Math. AppI. 11(1973), pp. 73--H82.
[2] C.C. Broyden, A class of methods for solving nat Itaneous
equations, Math. Compi 19(1965), pp. 577--H593.
[3] C.C. liroyden, The convergence of an algorithm 1g sparse
nonlinear systems, Math. Comp. 25(1971), pp. 2ii5-7??
[4] C.G. Broyden, 3.E. Dennis, Jr., and 3.3. More, ?n (I, i,.?al and super--H
linear cce of quasi--HNewton methods, 3. lnsr, N?r6. Alpi. 12(1973),
pp. 223--H245.
[5] A.R. Curtis and J.K. Reid, The choice of step lentrh' he' using
differences to approximate Jacohian matrices, 3. las . ?art'. Ai'pt. 13
(1974), pp. 121-126.
[6] A.R. Curtis, M.3.D. Powell, and J.K. Reid, On the e.tta?tton ol
sparse Jacohian matrices, 3. lost. M'th. AppI. 13(101;,). 117-119.
[7] w.C. Davidon, Variable metric method for minimi?at j?.',, A?C Rts. and
tev. Report ANL--H5990, 1959 (revised 1966).
[8] 3.E. Dennis, Jr., On convergence of Newton--Hlike ett,,,dx. in Numerical
Methods for Non--Hlinear Algebraic Equations, cd. P. Ral,in"?it', Cordon and
Breach, Londoo, 1970, pp. 163--H181.
[9] 3.f. tennis, Jr., On the convergence of Broyden' .,,tl,?d for nonlinear
systems of equations, Mail'. Camp. 25(1971), pp. l59?.'t''
[10] 3.i. tennis, Jr., Toward a unified convergence ii?'??y I?r Newton--Hlike
methods, in Nonlinear Functional Analysis and Appliat I?,,,s, ed. L.l1. RaIl,
Academic Press, New York, 1972, pp. 425--H472.
[it] 3.i. Dennis, Jr., On some methods based on Brny'l?"' s?rant approxi-
mation to the hessian, in Nometical Methods for N',',1 i,.? ipt imization,
ed. F.A. Lootsma, Academic Press, London, 1972, pi' I'
[12] 3.f. Dennis, 3r. , A brief survey of convergea?e ,,v',lt' for quasi--H
Newton methods, in Nonlinear Progran,ming. ed. R. Cart I, ad C. Lempke
SIAN-AMS Proceedings 9,1976, pp. 185-l9?.
[13] 3.E. tennis, Jr., and 3.3. Note', Lecture note? 1', ti.aerical Analy-
xis, Computer Science 421--H422, Cornell University. `i," t',,3.
[14] 3.E. tennis, Jr. and 3.3. Mare', A characterixati"' i superlinear
rouvergaud Its applicotlon to quasi--HNewton melL's, M,,i,. Camp. 28
(l?74), pp. 549-560.
--H 81 --H
--H 82
[153 J.E. Oconis, Jr. and J.j. More, Q?,asi-Newton methods, motivation and
theory, 53kM Review 19(1977), pp. 46-8?.
[16?. 1.5. Off, A sf siarx resear?b, Proceediogs of the ifEf
65(1977), pp. 500--H535.
[17] A.M. Erisman and J.K. Reid, Monitoring the stability ol the triangu-
lar factoriz?tion of a oplx, Nuner. Math. 22(1974), pp. 183--H186.
[183 R. fletcher and M.d.f. Powell, A rapidly convergent descent method
for minimizatiun, Computer J. 6(19t3), pp. 163-168.
[193 P.li. Ciii, C. Colub, W. Murray, and M.A. Saunders, Methods for modi-
fying matrix factorizations, Math. Comp. 28(1974), pp. 505--H535.
[203 P.F. Ciii and 1'. Murray, tluasl--HNrwton metlods for uncoostralned min--H
1mi?ation, J. lost. Math. Appi. 9(1972), pp. 91--H108.
[213 P.f. Ciii, w. Murray and M.A. Saunders, Methods for computing the
Ll)V facrors of a matrix, Math. Comp. 29(1975), pp. 1051--H1077.
[223 I). Coidfarb, A family of variable-metric methods derived by varia-
tioo?I means, Math. Comt. 24(1970), pp. 23--H26.
[233 J. Creenstadt, Variations on variable-metric methods, Math. Comp.
24(1970), pp. 1--H18.
[741 L.V. Ranrorovich and C.P. Akilov, Functional Analysis in Nc"mel
Spaces. Fittmatgiz, Moscow, 1959; EngI . transi., Petamon Press, Oxford,
1964.
[25] A. Revorkian and J. Snoek,ompositino in large scale system
theory and appl ica- in solving large sets of nool inert simo Itancous
equain Orconposi t ions of large Scale Prohiems , ed. 1). Nimmeihian,
North Holland, New Sock, 19i3, pp. 4?i?Slh.
[261 J .j . More, private ronuno cat ion, 1977.
[273 J.M. Ortega and w.C. Rheinholdt, iterative Solution of Nunlinear
Equations In Several Variables, Acadexit Press, New York, 1970.
[283 T.A. Porsehiog, On the origins and numerical solution of some xparse
nonliorat systems, in Sparse Matrix Computations, Academic Press, New
York, 1976, pp. 439-453.
[293 M.i.O. Powell, A survey ofical methods for unconstralned opt 1
mixa, SiAll Review 12(1970), pp. 79-97.
[301 M.J.O. Powell, A new algorithe for unconstrained optimizatlon, In
Nonlinear Prograrrming, ed. J.li. Rosen, O.L. Maugasarlan, K. Ritter, Aca-
demic Press, New York, 1970, pp. 31-65.
[313 M.j.D. Powell, A view of uncon?trained optimizatlon, AfRE Harweli
Teth. Rep. C.S.S. 14, 1975.
- 83 --H
[32] J.K. Iteid, Fortran subroutines lor the solution of par- s?;f??
of nonlinear e9nations, AE?? ilarwell Tech. ?ep. AfRf--HN77Q1. l??2.
[33] w.c. Rhelnholdr, Methods for Solving Systens of Noulinrar Eqoatfons,
SlAM, Philadelphia, 1974.
[34] W.C. Rhefnholdt , On the solution of large sparse scts of 1 in--H
ear eqoations, In ConTotat ion?l Me?hanics, ed. J.T. Oden, Spr1n? verlag,
Berlin, 1975, pp. l69?i94.
[35] R. Schoahef, private connunication, 1977.
[36] L.K. Schuhert, Modifi?atlon of a quasl--HNewtQn method in, nonlinear
equations wlch a sparse Ja?ohtan, Math. Comp. 24(1970), pp. 27-30.
[37] R.P. Tewarson and .1. Stephenson, Solution of large sparse systems
of nonlinear equations, College of Engr. Rep. No. 265, StiNY at 5y
Brook, 1975.
[38] Fh.l. Toint, On 5d syrmnettir matrix updating subjert to a
linear equation, flept. of Appi. Math. and Theor. Phys. Tech. Rep. i)AMTP
77iNA1, University of Canfrldge, 1977.
[39] fir. L. Toint, lone nunorical results using a sparse mattla ntrdatlng
fornula in uronstrained opt iniration, nanuscript. tept. of Mathe,?rt ir-s,
FNDf, Nanur, Belgiun, 1977.
[40] C. Van Lean, private corm-anleation, 1977.
[41] J.N. Wilkinson, The Algefraic Bigenvalue Problem, Clarendun Press,
Oxford, 1965.
[42] J.M. Wilkinson and C. Reinsch, Linear Algebra, Springer-Veriag, New
York, 1971.