Newton Inexact Methods Applied to UnderDetermined Systems by Joseph P. Simonis A Dissertation Submitted to the Faculty of WORCESTER POLYTECHNIC INSTITUTE in Partial Fulfillment of the Requirements for the Degree of Doctor of Philosophy in Mathematics
May 2006
APPROVED: Professor Homer F. Walker, Dissertation Advisor Mathematical Sciences Department Professor Joseph D. Fehribach, Committee Member Mathematical Sciences Department Professor Arthur C. Heinricher, Committee Member Mathematical Sciences Department Dr. John N. Shadid, Committee Member Sandia National Laboratories Professor Suzanne L. Weekes, Committee Member Mathematical Sciences Department
Contents
1 Introduction 2 Overview 2.1 Preliminaries . . . . . . . . . . . . . . . . . 2.2 Newton-like Methods . . . . . . . . . . . . . 2.2.1 Newton's Method . . . . . . . . . . . 2.2.2 Inexact Newton Methods . . . . . . 2.2.3 Globalized Inexact Newton Methods 2.2.4 Trust Region Methods . . . . . . . . 2.3 Normal Flow Method . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4 6 6 9 9 10 12 13 14 16 16 22 26 32 34 41 45 45 45 46 47 49 51 51 53 56 56 57 57 59 60 63 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 78 78 79 80 82 87
3 Methods and Theories 3.1 Inexact Newton Methods for UnderDetermined Systems . . . . . . . 3.2 A Globalized Inexact Newton Method for UnderDetermined Systems 3.3 Backtracking Methods . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.3.1 Choosing the Scaling Factor . . . . . . . . . . . . . . . . . . . . 3.4 TrustRegion Methods for UnderDetermined Systems . . . . . . . . . 3.4.1 UnderDetermined Dogleg Method . . . . . . . . . . . . . . . . 4 Numerical Experiments 4.1 Test Problems . . . . . . . . . . . . . . . 4.1.1 The Bratu Problem . . . . . . . 4.1.2 The Chan Problem . . . . . . . . 4.1.3 The LidDriven Cavity Problem 4.1.4 The 1D Brusselator Problem . . 4.1.5 The 2D Brusselator Problem . . 4.2 The Solution Algorithms . . . . . . . . . 4.3 Results . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
5 Summary 5.1 Summary . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 5.2 Additional Applications . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 5.3 Future Work . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . Appendices A Adaptation of GMRES B Matlab Files C Raw Data C.1 Chan Problem Results . . . . . . C.2 Bratu Problem Results . . . . . . C.3 1D Brusselator Problem Results . C.4 2D Brusselator Problem Results . C.5 Driven Cavity Problem Results .
1
List of Figures
4.1 4.2 4.3 4.4 4.5 4.6 4.7 4.8 The Bratu Solution Space w/ Solution . . . . . . . . . . . . . . . . . Chan Solution . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . Lid Driven Cavity . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1D Brusselator Solution Space w/ Solution . . . . . . . . . . . . . . . 2D Brusselator Solution . . . . . . . . . . . . . . . . . . . . . . . . . Table of Compute Times . . . . . . . . . . . . . . . . . . . . . . . . . Table of Normalized Compute Times . . . . . . . . . . . . . . . . . . Residual Norms of Nonlinear functions . . . . . . . . . . . . . . . . . 47 48 49 50 52 54 54 55
2
Abstract
Consider an under-determined system of nonlinear equations F (x) = 0, F : IR m IR n , where F is continuously differentiable and m > n. This system appears in a variety of applications, including parameterdependent systems, dynamical systems with periodic solutions, and nonlinear eigenvalue problems. Robust, efficient numerical methods are often required for the solution of this system. Newton's method is an iterative scheme for solving the nonlinear system of equations F (x) = 0, F : IR n IR n . Simple to implement and theoretically sound, it is not, however, often practical in its pure form. Inexact Newton methods and globalized inexact Newton methods are computationally efficient variations of Newton's method commonly used on large-scale problems. Frequently, these variations are more robust than Newton's method. Trust region methods, thought of here as globalized exact Newton methods, are not as computationally efficient in the largescale case, yet notably more robust than Newton's method in practice. The normal flow method is a generalization of Newton's method for solving the system F : IR m IR n , m > n. Easy to implement, this method has a simple and useful local convergence theory; however, in its pure form, it is not well suited for solving large-scale problems. This dissertation presents new methods that improve the efficiency and robustness of the normal flow method in the largescale case. These are developed in direct analogy with inexactNewton, globalized inexactNewton, and trustregion methods, with particular consideration of the associated convergence theory. Included are selected problems of interest simulated in MATLAB.
3
Chapter 1 Introduction
This dissertation presents newly developed methods for solving the under-determined nonlinear system of equations F (x) = 0, F : IR m IR n with m > n 1. These methods are shown to be globally robust, locally fast, and computationally efficient on large-scale systems of equations. We define an underdetermined system of nonlinear equations to be any system of nonlinear equations with more unknowns than equations, regardless of the uniqueness or existence of its solutions. These systems appear in a variety of applications. After discretization, certain nonlinear partial differential equation (PDE) eigenvalue problems (e.g. the Bratu problem) and some parameterdependent PDE's (e.g. the driven cavity problem) take the form F (x, ) = 0 with x IR n and IR. Under-determined systems also sometimes appear when calculating periodic orbits of dynamical systems. Here, one seeks x(0) and T satisfying f (x(t), t) =
dx(t) . dt T 0
f (x(t), t)dt = 0, where
The function f (x(t), t) is assumed to be nonlinear.
This dissertation is divided into three main sections following the introduction. The first of these presents background material. Here, general notation is discussed along with useful lemmas and definitions. This section also includes descriptions of Newton's method and relevant Newton-like methods for solving the nonlinear sys4
tem F (x) = 0, F : IR n IR n ; these are important since they will be used later as models for the new methods. Relevant Newton-like methods include inexact Newton methods, globalized inexact Newton methods and trustregion methods. Inexact Newton methods only approximately solve the Newton system at each iteration. Details about all of these methods can be found in [2, 21, 31]. Globalized inexact Newton methods impose additional requirements on the generated iterates to improve robustness of the algorithms. See [5, 6, 22] and included references. Trustregion methods for nonlinear systems of equations stem from methods in unconstrained optimization. For an understanding of both the methods and their subsequent adaptation to nonlinear systems see [4, 18, 32, 20] and the references therein. The normal flow method, an adaptation of Newton's method for solving under-determined systems, is also presented here. Further material about the adaptation of Newton's method to underdetermined systems can be found in [16, 35] and the included references. The second section presents new methods for solving the under-determined system of nonlinear equations. Here, the methods from the first section are used to motivate generalizations of the normal flow method. The discussion and development of the new methods closely parallel the development of methods in [2, 5]. These generalizations are shown to have fast local convergence properties and to be globally robust. Additionally, if suitably implemented, they are computationally efficient for solving large-scale systems of equations. The final section discusses numerical experiments. Several specific methods are coded in MATLAB and applied to model problems of interest. The test problems include nonlinear eigenvalue problems (the Bratu problem [13] and the Chan problem [1]), a parameterdependent fluid flow problem (the driven cavity problem [7, 30]), and periodic orbit calculations (the Brusselator problem in one and two dimensions [15, 10, 11]). 5
Chapter 2 Overview
2.1 Preliminaries
We begin with a brief overview of assumptions, notation, and a few definitions and lemmas used throughout the text. The norm is assumed to be the Euclidean norm on vectors or the induced
norm on matrices throughout. Most results can be extended to the case of an arbitrary innerproduct vector norm. The function F is from IR m to IR n and is continuously differentiable. It has a Ja cobian matrix denoted by F , an nm matrix with Fij Fi (x) xj 1in,1jm
.
The -neighborhood of a point x IR n is the set N (x) {y IR n | x-y < }. A stationary point of F is a point x IR m for which there does not exist an s IR m such that F (x) + F (x)s < F (x) . The stationary points include local minimizers of F .
6
Definition 1 ([4]).
Let x IR m and xk IR m , k = 1, 2, . . . Then {xk } is said to converge to x if limk xk - x = 0. ^ If there exists a constant c [0, 1) and an integer k 0 such that for all ^ k k, xk+1 -x c xk -x , then {xk } is said to be q-linearly convergent to x . If for some sequence {ck } that converges to 0, xk+1 - x ck xk - x for each k, then {xk } is said to converge q-superlinearly to x . If {xk } converges to x and there exist constants p > 1 and c 0 such that xk+1 - x c xk - x
p
for each k, then {xk } is said to converge to x with
q-order at least p. If p = 2 or p = 3, then the convergence is said to be q-quadratic or q-cubic, respectively. Throughout, we use q-convergence as opposed to r-convergence. The q stands for "quotient", and r stands for "root". A sequence {xk } converges to xk with r-order p if { xk+1 - x } is bounded above by a sequence in IR that converges to zero with q-order p. Definition 2 ([4]). A function g is Hlder continuous with exponent p (0, 1] o and constant in a set IR m if, for every x, y , g(x) - g(y) x - y p . Definition 3 ([4]). A function g is Lipschitz continuous with constant in a set IR m , written g Lip (), if for every x, y , g(x) - g(y) x - y . Definition 4 ([3]). Given F : IR m IR n continuously differentiable and x IR m , F (x)+ is the pseudoinverse of F (x), if, given b IR n , F (x)+ b IR m is the 7
solution of F (x)s = b having minimal Euclidean norm. When F (x) is of full rank, the pseudoinverse has the form F (x)+ = F (x)T (F (x)F (x)T )-1 and is called the MoorePenrose pseudoinverse. Lemma 1 ([21]). Let F : IR m IR n be a continuously differentiable function. For any x IR m and > 0, there exists a > 0 such that F (z) - F (y) - F (y)(z - y) z - y whenever y,z N (x). Lemma 2 ([21]). Let F : IR m IR n be continuously differentiable in the open convex set IR m , let x , and let F be Hlder continuous with exponent p and o constant at x in the neighborhood . Then, for any x + s , F (x + s) - F (x) - F (x)s Proof. By Lemma (4.1.9) in [4], s 1+p
1+p
(2.1)
.
(2.2)
1
F (x + s) - F (x) - F (x)s = =
0
0 1
F (x + ts)s dt - F (x)s
[F (x + ts) - F (x)]s dt.
We then obtain
1
F (x + s) - F (x) - F (x)s
0 1
F (x + ts) - F (x) s dt ts
p
s dt
1
0 1+p 0 1+p
= s
tp dt .
= s 1+p
8
Lemma 3. Assume F (x) IR nm is a continuous function of x and is of full rank. Then F (x)+ is a continuous function of x. Proof. Because F (x) is of full rank, the MoorePenrose pseudoinverse can be written F (x)+ = F (x)T (F (x)F (x)T )-1 . Since F (x)T is a continuous function of x, we have that F (x)F (x)T is a continuous function of x, and it follows that (F (x)F (x)T )-1 and, hence, F (x)+ are continuous functions of x.
2.2
Newton-like Methods
This section presents three classes of Newton-like methods designed to solve F (x) = 0 with F : IR n IR n . We begin with a description of Newton's method and follow with descriptions of inexact Newton methods, globalized inexact Newton methods, and trust region methods.
2.2.1
Newton's Method
Consider the problem: find x IR n such that F (x) = 0, (2.3)
where F : IR n IR n is continuously differentiable. Newton's method begins by assuming an initial guess, x0 , and generates a sequence of iterates via xk+1 = xk - F (xk )-1 F (xk ). In practice, this involves solving the linear system F (xk )sk = -F (xk ) 9 (2.4)
for the Newton step sk , and then defining xk+1 = xk + sk . Algorithm NM: Newton's Method Let x0 be given. For k = 0 step 1 until do: Set xk+1 = xk + sk . Solve F (xk )sk = -F (xk )
Under mild assumptions the sequence will approach a root of F provided x0 is sufficiently near the root. Theorem 1 ([34, 4]). Suppose F is Lipschitz continuously differentiable at x , F (x ) = 0 and F (x ) is nonsingular. Then for x0 sufficiently near x , {xk } produced by Newton's method is well-defined and converges to x with xk+1 - x c xk - x for a constant c independent of k. The method is simple to implement and theoretically sound, but, in its pure form, not often used to solve large-scale problems. The exact linear solve at each iteration makes the method computationally inefficient.
2
2.2.2
Inexact Newton Methods
Inexact Newton methods [2] are variations of Newton's method designed to be computationally efficient on largescale problems, and are commonly used in the large-scale case. Recall, that the general idea of Newton's method is to linearize F around a current guess, xk , in hope that the root of the linear model, xk+1 , is a better approximation of the root of the nonlinear problem than was xk . There are two drawbacks 10
with this method. First, solving for a root of the linear model, in practice, may be computationally timeconsuming. Second, when far from a solution, the root of the linear model may not be a good approximation of the root of the nonlinear problem. We replace the Newton step with an "inexact" Newton step. We no longer require sk to exactly solve F (xk )sk = -F (xk ), rather only that sk be a point where the norm of the local linear model has been reduced. Precisely, we find some k [0, 1) and require sk to satisfy F (xk ) + F (xk )sk k F (xk ) . (2.5)
Notice that as k approaches zero, sk approaches the Newton step. This replacement allows sk to be calculated "cheaply." Often, an efficient iterative linear solver such as the generalized minimal residual method (GMRES)1 is used for this calculation. The following algorithm is the inexact Newton method (INM). Algorithm INM: Let x0 be given. Find some k [0, 1) and sk that satisfy F (xk ) + F (xk )sk k F (xk ) Set xk+1 = xk + sk . The scalar k is called the forcing term and its choice affects both local convergence properties and the robustness of the method [2, 6]. Assume x is a solution of (2.3) at which the Jacobian is of full rank. If x0 is sufficiently close to x and 0 k max < 1 for each k, then {xk } converges to x q-linearly in some norm. Furthermore, q-superlinear convergence is obtained if limk k = 0. Finally, if k = O( F (xk ) ), then the convergence is q-quadratic. See [2] for further details.
1
For k = 0 step 1 until do:
See Appendix A
11
2.2.3
Globalized Inexact Newton Methods
The robustness of an inexact Newton method often is enhanced by "globalizations," i.e., augmentations of the basic method that test and modify steps to ensure adequate progress toward a solution [5, 22]. A step satisfying the inexact Newton condition (2.5) yields a decrease in the local linear model norm, yet this decrease is not always reflected in the nonlinear residual norm. In other words, the chosen step may not actually reduce F . To ensure a reduction of F , an additional step selection criterion is added. The step should reduce F at least some fraction of the reduction predicted by the local linear model of F . More precisely, given a t (0, 1), sk should be chosen to satisfy (2.5) and a sufficient decrease condition: F (xk + sk ) [1 - t(1 - k )] F (xk ) . The resulting algorithm is a globalized inexact Newton method (GINM). Algorithm GINM: For k = 0 step 1 until do: Let x0 and t (0, 1) be given. (2.6)
Find some k [0, 1) and sk that satisfy F (xk ) + F (xk )sk k F (xk ) and F (xk + sk ) [1 - t(1 - k ) F (xk ) Set xk+1 = xk + sk . The following is a global convergence theorem for algorithm GINM. Theorem 2 ([5]). Assume that algorithm GINM does not break down. If k0 (1- k ) is divergent, then F (xk ) 0. If, in addition, x is a limit point of {xk } such that F (x ) is invertible, then F (x ) = 0 and xk x . 12
2.2.4
Trust Region Methods
A general trust region method produces a sequence of iterates using the following procedure: at each iteration we assume the local linear model is an accurate representation of the nonlinear function within some closed -ball around the current iterate. We choose a step, s, to minimize F (x) + F (x)s over all s satisfying s . Then, we check to see if s is acceptable. If it is not acceptable, this indicates the local linear model is not a good representation of F in the -ball; is decreased and a new s is chosen. This is repeated until an acceptable step is found. The value is a measure of our "trust" of the local linear model. One commonly used test for step acceptability is the ared/pred condition[5]. Let ared be the actual reduction in function norm obtained by taking a step s: ared(s) F (x) - F (x + s) . (2.7)
The predicted reduction, pred, is the reduction predicted by the local linear model; pred(s) F (x) - F (x) + F (x)s . (2.8)
The ared/pred condition for step acceptability requires the actual reduction in function norm to be at least some fraction of the predicted reduction; ared(s) t pred(s), t [0, 1). (2.9)
The following general trust region method (TR) from [5] is similar in spirit to the method in [18]. Algorithm TR: Trust Region Method Let x0 , 0 > 0, 0 < t u < 1, and 0 < min < max < 1 be given For k = 0 step 1 until do: Set k = k and 13
choose sk arg min
s k
F (xk ) + F (xk )s
While aredk (sk ) < t predk (sk ) do: Choose [min , max ] sk arg min
s k
Update k k and choose
F (xk ) + F (xk )s
Set xk+1 = xk + sk If aredk (sk ) u predk (sk ) choose k+1 k ; Else choose k+1 min k
Theorem 3 ([5]). Assume that Algorithm TR does not break down. Then every limit point of {xk } is a stationary point of F . If x is a limit point of {xk } such that F (x ) is invertible, then F (x ) = 0 and xk x ; furthermore, sk = -F (xk )-1 F (xk ), the full Newton step, whenever k is sufficiently large.
2.3
Normal Flow Method
Now, consider an under-determined rootfinding problem: find x IR m such that F (x) = 0, (2.10)
where F : IR m IR n is a continuously differentiable function with m > n. Here, the linear system (2.4) is under-determined, i.e., it may have an infinite number of solutions. In order to develop a welldefined algorithm, an additional constraint must be imposed so that a unique step, sk , can be defined. Choosing sk to be the solution of the linear system (2.4) with minimum Euclidean norm gives the normal flow method [35]. The pseudoinverse solution of the linear system is a natural choice for a "Newton" step because it is the shortest step from the current iterate to a root of the linear problem and, therefore, the linear model is likely to be a better representation of the nonlinear function at that step than at other solutions of (2.4). Hereafter, the 14
normal flow algorithm will be called Newton's method for under-determined systems (NMU). Algorithm NMU: Let x0 be given. For k = 0 step 1 until do: Set xk+1 = xk + sk . Let sk = -F (xk )+ F (xk )
Mathematically, we have F (xk )sk + F (xk ) = 0 and sk Null(F (xk )). When F (xk ) is of full rank, sk will hereafter be referred to as the Moore-Penrose step. A local convergence theory for Algorithm NMU is given in [35] and generalized in [16]. The central result from [35] with respect to this method follows. Hypothesis 1. F is differentiable and F is of full rank n in an open convex set , and the following hold: (i) There exist 0 and p (0, 1] such that F (y) - F (x) y - x for all x, y . (ii) There is a constant for which F (x)+ for all x . Definition 5. For > 0, let = {x : y - x < y }. Theorem 4 ([35]). Let F satisfy Hypothesis 1 and suppose is given by Definition (5) for some > 0. Then there is an > 0 which depends only on , p, , and such that if x0 and F (x0 ) < , then the iterates {xk }k=0,1,... determined by Algorithm NMU are well defined and converge to a point x such that F (x ) = 0. Furthermore, there is a constant for which xk+1 - x xk - x
p+1 p
, k = 0, 1, . . .
(2.11)
If F (x) is Lipschitz continuous in then p = 1 and the iterates produced by Algorithm NMU converge q-quadratically. 15
Chapter 3 Methods and Theories
3.1 Inexact Newton Methods for UnderDetermined Systems
The previous subsection introduced a variation of Newton's method for solving the under-determined system F (x) = 0, F : IR m IR n , and briefly discussed its local convergence theory. This section presents a class of inexact Newton methods for application to the under-determined system. A convergence theory is developed for these new methods. Each iteration of NMU requires the step, sk , satisfying F (xk ) + F (xk )sk = 0 with sk NullF (xk ). Calculation of sk requires solving a linear system of equations. When n is large, this may be computationally expensive. To improve the computational efficiency, we allow for an approximate solution of the linear system. We seek an sk satisfying F (xk ) + F (xk )sk k F (xk ) , where k [0, 1). 16 (3.1)
and sk Null(F (xk )). (3.2)
Constraint (3.1) will henceforth be called the inexact Newton condition. The inexact Newton method for under-determined systems (INMU) follows: Algorithm INMU: Let x0 be given. Find some k [0, 1) and sk that satisfy F (xk ) + F (xk )sk k F (xk ) sk Null(F (xk )) For k = 0 step 1 until do:
Set xk+1 = xk + sk . The remainder of this section presents a theoretical foundation for this algorithm. Lemma 4. Assume s Null(F (x)), then s = F (x)+ F (x)s . Proof. Define s = F (x)+ F (x)s. Then s is the pseudoinverse solution of F (x) = F (x)s. s (3.3)
Additionally, s Null(F (x)) because s is the minimum norm solution of the linear problem. Rearranging equation (3.3) gives F (x)( - s) = 0; s therefore, the vector ( - s) Null(F (x)). However, because both s and s are s orthogonal to the null space of the Jacobian, it is also true that (-s) Null(F (x)) . s Therefore, (-s) Null(F (x))Null(F (x)) = {0}. We conclude that s = s. Then, s s = s = F (x)+ F (x)s . 17
Theorem 5. Let F satisfy Hypothesis 1 and suppose > 0. Assume that k max < 1 for k = 0, 1, . . . . Then there is an > 0 depending only on , p, , and max such that if x0 and F (x0 ) , then the iterates {xk } determined by Algorithm INMU are welldefined and converge to a point x such that F (x ) = 0.
Proof. Suppose x and s is such that s Null(F (x)) and F (x) + F (x)s max F (x) . Define x+ x + s and suppose x+ . F (x)+ F (x)+ F (x)s - F (x) + F (x) + F (x)s We can write s = F (x)+ F (x)s because s Null(F (x)). Then s =
F (x)+ ( F (x) + F (x) + F (x)s )
( F (x) + max F (x) ) = (1 + max ) F (x)
and F (x+ ) F (x+ ) - F (x) - F (x)s + F (x) + F (x)s + max F (x) max )]1+p F (x)
1+p
s 1+p 1+p [(1 + 1+p
+ max F (x) .
Choose > 0 sufficiently small that
[(1 1+p
+ max )]1+p p + max < 1 and
(1+max ) 1-
< .
If F (x) , then F (x+ ) F (x) . We argue by induction that if x0 and F (x0 ) then F (xk+1 ) F (xk ) < and xk for all k. First we show that F (x1 ) F (x0 ) < and x1 . We have that s0 Null(F (x0 )) and F (x0 ) + F (x0 )s0 max F (x0 ) , so
18
s0
(1 + max ) F (x0 ) (1 + max ) <
(1+max ) 1-
< . Therefore we have x1 because x0 . By the above argument, F (x1 ) F (x0 ) < . Now assume F (xj+1 ) F (xj ) < and xj for all j k. We show that this is true for j = k + 1. As before xj , sj satisfies sj Null(F (xj )), and F (xj ) + F (xj )sj max F (xj ) . Then xj+1 - x0
j l=0 j l=0
sl (1 + max ) F (xl )
j l=0 l=0
(1 + max ) =
(1+max ) 1-
l l
< (1 + max ) < ,
so x0 implies xj+1 . Again, using the earlier argument, F (xj+1 ) F (xj ) j F (x0 ) < . Now, F (xk+1 ) F (xk ) implies the sequence { F (xk ) } converges to zero, and since sk (1 + max ) F (xk ) , it must be that sk 0 as k . Note xk+l - xk
l-1 j=0 l-1 j=0
sk+j (1 + max ) F (xk+j )
l-1 j=0 j=0
(1 + max ) < (1 + max )
(1+max ) 1- (1+max ) k 1-
j F (xk ) j F (xk )
=
F (xk ) F (x0 )
(1+max ) k . 1-
19
Therefore, {xk } is a Cauchy sequence. It has a limit x . The continuity of F yields F (x ) = 0. Theorem 6. Let F satisfy Hypothesis 1 and suppose > 0. Assume that 0 k max < 1 for k = 0, 1, . . . and that {xk } produced by Algorithm INMU converges to x such that F (x ) = 0. Let M F (x ) and assume F (xk ) 2M for all xk . If > 0 and max in the proof above are chosen sufficiently small such that
((1 1+p
+ max ))1+p p + max <
1 (1+max )2M +1
and
(1+max ) 1-
< ,
then {xk } converges to x q-linearly. Proof. Start with xk+1 - x
j=k+1 j=k+1
sj (1 + max ) F (xj )
j=k+1 j j=0
= (1 + max ) =
(1+max ) 1-
F (xj ) F (xk+1 )
(1 + max ) = By the choice of , the term convergent.
(1+max ) 1- (1+max ) 1-
F (xk+1 ) F (xk ) F (xk ) - F (x ) xk - x .
2M (1+max ) 1-
2M (1+max ) 1-
is less than 1. Therefore {xk } is q-linearly
Theorem 7. Let F satisfy Hypothesis 1 and suppose > 0. Assume that 0 k max < 1 for k = 0, 1, . . . and that {xk } produced by Algorithm INMU converges to x such that F (x ) = 0. If k 0, then {xk } converges to x q-superlinearly. If k = O( F (xk ) p ), then {xk } x with q-order 1 + p. 20
Proof. Let M
F (x ) . There exists a > 0 such that F (x)
1+p
2M and
F (x + s) - F (x) - F (x)s <
s
1+p
whenever x - x . Assume that k is
sufficiently large that xk - x . We first show superlinear convergence. We have F (xk+1 ) F (xk+1 ) - F (xk ) - F (xk )sk + F (xk ) + F (xk )sk sk
1+p
1+p 1+p 1+p
+ k F (xk )
[(1 + k ) F (xk ) ]1+p + k F (xk ) [2M (1 + k ) xk - x ]1+p + k 2M xk - x .
Previous calculations give xk+1 - x Now, let ck
1+p (1+max ) 1- (1+max ) 1- (1+max ) 1-
F (xk+1 )
[2M (1 1+p [2M (1 1+p
+ k ) xk - x ]1+p + k 2M xk - x + k )]1+p xk - x
p
+ k 2M
xk - x .
[2M (1 + k )]1+p xk -x p +k 2M . Combined, limk xk -x =
0 and limk k = 0 imply limk ck = 0. Thus, {xk } is q-superlinearly convergent. Now assume k = O( F (xk ) p ). Because k is on the order of F (xk ) p , there exists a constant C independent of k such that k C F (xk ) all sufficiently large k. Then xk+1 - x
(1+max ) 1- (1+max ) 1- 2M (1 1+p 2M (1 1+p p
C(2M )p xk - x
p p
p
for
+ k )1+p xk - x + k )1+p xk - x
p
+ k
xk - x
+C(2M )p xk - x
(1+max ) 1- (1+max ) 1- 2M (1 1+p 2M (1 1+p
] xk - x
+ k )1+p + C(2M )p + max )1+p + C(2M )p
xk - x
1+p 1+p
xk - x
,
which gives q-order 1 + p convergence. It follows that, if F is Lipschitz continuous, then p = 1, and we have q-quadratic convergence. 21
3.2
A Globalized Inexact Newton Method for Under Determined Systems
Inexact Newton methods for under-determined systems can achieve fast local convergence rates. Under mild assumptions, including an x0 such that F (x0 ) is sufficiently small, the sequence generated by Algorithm INMU converges to a solution of problem (2.10). However, if an acceptable x0 cannot be found, the sequence may fail to converge. Here, the goal is to augment the step selection criteria of Algorithm INMU with a sufficient decrease condition. In analogy with GINM, the additional requirement on the chosen steps is meant to increase the likelihood that the iterates converge to a solution, given an arbitrary x0 . Additionally, we seek a modification that retains the fast rates of convergence. Each step, sk , must still satisfy the inexact Newton conditions (3.1) and (3.2). We now also require that the step reduce the norm of F at least some fraction of the reduction predicted by the local linear model. Given some t (0, 1), sk should be chosen such that F (xk + sk ) [1 - t(1 - k )] F (xk ) , (3.4)
the same criterion chosen by Eisenstat and Walker in [5]. They note that step criteria (3.1) and (3.4) are similar to acceptability tests used in certain minimization algorithms [18, 32] and methods for solving nonlinear equations [12, 26]. Imposing this additional constraint yields our globalized inexact Newton method for underdetermined systems (GINU) Algorithm GINMU: Global Inexact Newton Method For UnderDetermined Systems For k = 0 step 1 until do: Let x0 and t (0, 1) be given.
22
Find some k [0, 1) and sk that satisfy F (xk ) + F (xk )sk k F (xk ) sk Null(F (xk )) and F (xk + sk ) [1 - t(1 - k )] F (xk ) Set xk+1 = xk + sk . The first lemma shows that an sk satisfying (3.1), (3.2) and (3.4) can be found for each k so long as F (xk ) = 0 and xk is not a stationary point of F . Lemma 5. Let x and t (0, 1) be given and assume that there exists an s satis fying F (x) + F (x) < F (x) and s Null(F (x)). Then there exists min [0, 1) s such that, for any [min , 1), there is an s satisfying F (x) + F (x)s F (x) F (x + s) [1 - t(1 - )] F (x) s Null(F (x)). Proof. Clearly F (x) = 0 and s = 0. Set
F (x)+F (x) s F (x)
, , ,
(1-t)(1-) F (x) s
min max , 1 - where > 0 is sufficiently small that
(1-) s
F (x + s) - F (x) - F (x)s s whenever s . For any [min , 1), let s
1- s. 1-
Then 23
F (x) + F (x)s
= =
1- - (F (x)) + 1- (F (x) + F (x)) s 1- 1- - F (x) + 1- F (x) + F (x) s 1- - 1- F (x) + 1- F (x) 1-
= F (x) , and, since s = it follows that F (x + s) F (x + s) - F (x) - F (x)s + F (x) + F (x)s
1- 1- 1- 1-
s
1-min 1-
s ,
s + F (x) =
s + F (x)
= (1 - t)(1 - ) F (x) + F (x) = [1 - t(1 - )] F (x) . Assume y is an element of the null-space of F (x). Then
1- sT y = ( 1- s)T y
=
1- ()T y s 1-
= 0. Thus s Null(F (x)). Theorem 8. Assume that {xk } is generated by Algorithm GINMU. If F (x ) is of full rank, and there exists a independent of k for which sk (1 - k ) F (xk ) (3.5)
k0 (1-k )
is divergent, then F (xk ) 0. If, in addition, x is a limit point of {xk } such that
whenever xk is sufficiently near x and k is sufficiently large, then F (x ) = 0 and xk x .
24
Proof. From equation (3.4), F (xk ) [1 - t(1 - k-1 )] F (xk-1 ) F (x0 )
0j<k
[1 - t(1 - j )] (1 - j ) .
F (x0 ) exp -t
0j<k
Since t > 0 and 1 - j > 0, the divergence of
k0 (1
- k ) implies F (xk ) 0.
Suppose that x is a limit point of {xk } such that F (x ) is of fullrank and that {xk } does not converge to x . Let > 0 be such that there exist infinitely many k for which xk N (x ) and sufficiently small that (3.5) holds whenever xk N (x ) / and k is sufficiently large. Since x is a limit point of {xk }, there exist {kj } and {lj } such that, for each j, xkj N/j (x ),
kj + lj < kj+1 . Then for j sufficiently large, /2 = cannot hold for large j.
t t
xkj +lj N (x ), /
xkj +i N (x ), i = 0, . . . , lj - 1
xkj +lj - xkj
kj +lj -1 k=kj
sk (1 - k ) F (xk ) { F (xk ) - F (xk+1 ) } .
kj +lj -1 k=kj
kj +lj -1 k=kj t
F (xkj ) - F (xkj +lj ) F (xkj ) - F (xkj+1 )
But the last right-hand side converges to zero since xkj x ; hence, this inequality
An alternate proof of the first half of the theorem follows: 25
Proof. (Walker, private communication) From equation (3.4), t(1 - k-1 ) F (xk-1 ) and
k
F (xk-1 ) - F (xk )
F (x0 ) - F (xk )
=
i=1 k-1
( F (xi-1 ) - F (xi ) ) ( F (xi-1 ) - F (xi ) ) + F (xk-1 ) - F (xk ) .
=
i=1
This implies
k-1
F (x0 )
=
i=1 k-1
( F (xi-1 ) - F (xi ) ) + F (xk-1 ) ( F (xi-1 ) - F (xi ) ) (t(1 - i-1 ) F (xi-1 ) ) (1 - i-1 ) F (xi-1 ) .
k0 (1
i=1 k-1
i=1
k-1
= t
i=1
Since t > 0 and 1 - j > 0, the divergence of
- k ) implies F (xk ) 0.
3.3
Backtracking Methods
The global inexact Newton method for an under-determined system presented above generates a sequence of steps satisfying the inexact Newton and sufficient decrease conditions. This section discusses methods for determining satisfactory steps. Assume that an initial step satisfying (3.1) and (3.2) can be found, i.e., an sk approximating the MoorePenrose step and a forcing term, k , are computed. Furthermore, assume this step does not satisfy the sufficient decrease condition. Backtracking methods 26
systematically scale sk and k to find an sk and k satisfying all three conditions (3.1), (3.2), and (3.4). This leads to the under-determined backtracking method (BINMU): Algorithm BINMU: Let x0 and t (0, 1), max [0, 1), and 0 < min < max < 1 be given. For k = 0 step 1 until do:
Find some k [0, max ] and sk that satisfy F (xk ) + F (xk )k k F (xk ) , s sk Null(F (xk )).
Evaluate F (xk + sk ). Set k = k and sk = sk . While F (xk + sk ) > [1 - t(1 - k )] F (xk ) , do Update sk sk and k 1 - (1 - k ). Choose [min , max ].
Evaluate F (xk + sk ). Set xk+1 = xk + sk .
If a step satisfying the original inexact Newton condition is found, then properly scaling the step will yield a step satisfying both a modified inexact Newton condition and the associated sufficient decrease condition. Theorem 9. Assume that at the k th step of Algorithm BINMU there exists an k [0, max ] and sk satisfying F (xk ) + F (xk )k k F (xk ) s sk Null(F (xk )). Also, assume F (xk ) is of full rank. Then the whileloop will terminate in a finite
27
number of steps with an sk and k satisfying F (xk ) + F (xk )sk k F (xk ) sk Null(F (xk )) F (xk+1 ) [1 - t(1 - k )] F (xk ) .
Proof. The ith iteration of the whileloop scales sk by some i [min , max ]. At the mth step of the whileloop sk =
m i=1 i m i=1 i m i=1 i sk
and k = 1 -
m i=1 i (1
- k ). Notice
< .
m i=1 max
m = max . Given any > 0, an m can be found such that m
Choose m large enough such that F (xk + m sk ) - F (xk ) - F (xk )m sk C m sk , where C =
1 F (xk )+ ( (1-t)(1-max ) ). (1+max )
We claim that sk = m sk and k = 1 - m (1 - k ) are satisfactory. Indeed, F (xk ) + F (xk )sk = (1 - m )F (xk ) + m F (xk ) + m F (xk )k s
(1 - m ) F (xk ) + m F (xk ) + F (xk )k s (1 - m ) F (xk ) + m k F (xk ) = [1 - m + m k ] F (xk ) = [1 - m (1 - k )] F (xk ) = k F (xk ) . Further, assume y is an element of the null-space of F (xk ); 28
sT y = (m sk )T y k = m (k )T y s = 0. Finally, F (xk + sk ) F (xk ) + F (xk )sk + F (xk + sk ) - F (xk ) - F (xk )sk
k F (xk ) + F (xk + sk ) - F (xk ) - F (xk )sk k F (xk ) + C sk k F (xk ) + Cm F (xk )+ F (xk )k s
k F (xk ) + Cm F (xk )+ ( F (xk ) + F (xk ) + F (xk )k ) s k F (xk ) + Cm F (xk )+ (1 + k ) F (xk ) = = k + Cm F (xk )+ (1 + k ) k + k + F (xk ) F (xk ) F (xk )
m (1 - t)(1 - max )(1 + k ) 1 + max m (1 - t)(1 - max )(1 + max ) 1 + max
[k + m (1 - t)(1 - k )] F (xk ) = [k + 1 - 1 + m (1 - k ) - t + t - tm (1 - k )] F (xk ) = [k + 1 - (1 - m (1 - k )) - t + t(1 - m (1 - k ))] F (xk ) = [k + 1 - k - t + tk ] F (xk ) = [1 - t + tk ] F (xk ) = [1 - t(1 - k )] F (xk ) . Therefore, satisfactory sk and k are found in at most m steps, so the whileloop always terminates in a finite number of steps. 29
Theorem 10. Assume that {xk } is generated by Algorithm BINMU. Assume x is a limit point of {xk } such that F (x ) is of full rank. Then F (x ) = 0 and xk x . Furthermore, k = k and sk = sk for all sufficiently large k. Proof. Set K = F (x )+ and let > 0 be sufficiently small that F (x)+ 2K whenever x N (x ). Let =
1 (1-t)(1-max ) ( (1+max ) ). 2K
^ There exists a > 0 such that
F (z) - F (y) - F (y)(z - y) z - y ^ for all z, y N (x ). Let = min{, /2}. Suppose that xk N (x ). Let m be the ^
m smallest integer such that max 2K(1 + max ) F (x0 ) < . Then, with m as in the
proof of Theorem 9, m sk
m max sk
m = max sk m max 2K(1 + k ) F (xk ) m max 2K(1 + max ) F (xk ) m max 2K(1 + max ) F (x0 )
< . This implies F (xk + m sk ) - F (xk ) - F (xk )m sk m sk ,
which, as in the proof of Theorem 9, guarantees that the whileloop terminates in at most m iterations. Therefore, when xk N (x ) we have 1 - k = m (1 - k ) m (1 - max ) 30
m min (1 - max ).
m It is clear that min (1 - max ) > 0. It is given that x is a limit point of {xk },
so there are an infinite number of xk N (x ). Therefore, the sum diverges. when xk N (x ). Indeed, sk F (xk )+ F (xk )sk 2K ( F (xk ) + F (xk ) + F (xk )sk ) 2K(1 + k ) F (xk ) with = 2K(1 + max ) . m min (1 - max )
2K(1+max )(1-k ) m min (1-max )
k0 (1
- k )
We claim that sk (1 - k ) F (xk ) for some independent of k
2K(1 + max ) F (xk ) F (xk ) = (1 - k ) F (xk )
Therefore, by Theorem 8 we have that F (x ) = 0 and xk x . for all sufficiently large k, it is sufficient to show that F (xk + sk ) - F (xk ) - F (xk )k sk s
To show k = k
(3.6)
for all sufficiently large k. Equation (3.6) is true if sk < . Note that xk N (x ) for all sufficiently large k, and, therefore F (xk )+ 2K. Now sk F (xk )+ F (xk )k s
2K( F (xk ) + F (xk ) + F (xk )k ) s 2K(1 + k ) F (xk ) 2K(1 + max ) F (xk ) . It is clear that F (xk ) 0, so there exists some k such that for all k > k we have 2K(1 + max ) F (xk ) < . Therefore, for k > k, sk < , which implies (3.6). 31
3.3.1
Choosing the Scaling Factor
This section does not contribute to the mathematical literature, but it is included here for completeness. In-depth descriptions of the subsequent methods can be found in [4] and [34]. A scaling factor [min , max ] must be chosen at each iteration of the whileloop in Algorithm BINMU. The goal is to choose such that xk+1 = xk + sk is an acceptable next iterate. Ideally, minimizes F (xk + sk ) , or equivalently F (xk + sk ) 2 . However, this 1-dimensional minimization problem may be computationally expensive. An alternative is to find an easy-to-minimize approximation to F (xk + sk ) 2 . Two of the most popular schemes for choosing are described below. The quadratic backtracking method chooses to be the minimizer of a quadratic polynomial approximating the function g() = F (xk + sk ) 2 . Let p() denote this quadratic polynomial. The polynomial can be defined using three pieces of information; g(0) = F (xk ) 2 , g(1) = F (xk + sk )
2
and g (0) = 2F (xk )T F (xk )sk . Notice
the first two are already known, and the third is relatively inexpensive to calculate. Using these three values, p() is determined, and its minimizer can be calculated. The quadratic polynomial is given by p() = [g(1) - g(0) - g (0)]2 + g (0) + g(0). The derivatives are: p () = 2[g(1) - g(0) - g (0)] + g (0) and p () = 2[g(1) - g(0) - g (0)]. (3.7)
32
If p () 0, then the quadratic function is concave down, so choose = max . If p () > 0, then find such that p () = 0: 0 = 2[g(1) - g(0) - g (0)] + g (0) =
-g (0) . 2[g(1)-g(0)-g (0)]
Correcting to ensure [min , max ], one obtains . The method updates sk sk and 1 - (1 - k ) k and then checks to see if the updated sk and k satisfy the stepselection criterion. The cubic backtracking method uses a cubic polynomial to approximate g() = F (xk + sk )
2
after the first step-length reduction. The cubic polynomial, p(), is
constructed using four interpolation values. Again, is chosen to be the minimizer of p() over [min , max ]. On the first step reduction there is no clear way to choose a fourth value, so just three are chosen, and a quadratic polynomial is minimized, yielding 1 . If subsequent reductions are necessary four values are available. The two values g(0) and g (0) are used, along with values of g at the two previous values. For example, 2 is found using g(0), g (0), g(1 ), and g(1). Generalizing, i uses g(0), g (0), g(i-1 ), and g(i-2 ). As in [4], denote the two previous values as prev and 2prev . The cubic polynomial approximation of the function F (x + sk )
2
is
p() = a3 + b2 + g (0) + g(0) with a b =
1 prev -2prev 1 2 prev -2prev 2 prev -1 2 2prev prev 2 2prev
g(prev ) - g(0) - g (0)prev g(2prev ) - g(0) - g (0)2prev
-b+
.
The local minimizer of the model is given by + =
b2 -3ag (0) . 3a
As before, update
the step and forcing term by sk sk and 1 - (1 - k ) k . 33
3.4
TrustRegion Methods for UnderDetermined Systems
Trustregion methods for an under-determined system of equations are very similar to the methods for the fully determined system. That is, we define a region in which the local linear model is expected to be an accurate representation of the nonlinear function. A step is chosen to minimize the local linear model norm within this region and tested to see whether it satisfies the ared/pred condition. If it does not, the trust region is shrunk and a new minimum is calculated. To ensure locally fast convergence, the steps must approach the Moore-Penrose steps as {xk } approaches a root of F (x). This condition suggests that each step be chosen such that it is orthogonal to the null space of the Jacobian. The trustregion method for under-determined systems (UTR) becomes: Algorithm UTR: Let x0 , 0 > 0, 0 < t u < 1, and 0 < min < max < 1 be given. For k = 0 step 1 until do: Set k = k and choose sk arg min
s k
F (xk ) + F (xk )s
While aredk (sk ) < t predk (sk ) do: Update k k , and choose sk arg min
s k
with sk Null(F (xk )). Choose [min , max ].
F (xk ) + F (xk )s
Set xk+1 = xk + sk .
with sk Null(F (xk )).
If aredk (sk ) u predk (sk ) choose k+1 k ; Else choose k+1 min k .
The analogy between UTR and TR is manifest in this section. Parallels between 34
the theorems presented here and in [5] reflect this close relationship. The following Lemma was shown in [5] in the fully determined case and is extended here to the underdetermined case. Lemma 6. Assume that {xk } is such that predk (sk ) (1 - k ) F (xk ) and aredk (sk ) t (1 - k ) F (xk ) for each k, where t (0, 1) is independent of k. If x is a limit point of {xk } such that there exists a independent of k for which sk predk (sk ) whenever xk is sufficiently near x and k is sufficiently large, then xk x . Proof. Assume that {xk } does not converge to x . Let > 0 be such that there exist infinitely many k for which xk N (x ) and sufficiently small that sk predk (sk ) / holds whenever xk N (x ) and k is sufficiently large. Since x is a limit point of {xk }, there exist {kj } and {lj } such that, for each j, x kj xkj +i xkj +lj N (x ), i = 0, . . . , lj - 1 N (x ), / N/j (x ),
kj + lj < kj+1 . Then for j sufficiently large, /2 =
xkj +lj - xkj
kj +lj -1 k=kj
sk predk (sk )
kj +lj -1 k=kj
kj +lj -1 aredk (sk ) k=kj t { t { t
F (xkj ) - F (xkj +lj ) } F (xkj ) - F (xkj+1 ) }.
But the last right-hand side converges to zero because F is continuous and xkj x ; hence, this inequality cannot hold for large j.
35
Lemma 7. Assume that {xk } is a sequence generated by Algorithm UTR. Suppose that x is a limit point of {xk } such that there exists a independent of k for which sk predk (sk ) (3.8)
whenever xk is sufficiently near x and k is sufficiently large. Then xk x and lim inf k k > 0. Proof. It is clear that {xk } satisfies the hypotheses of Lemma 6, and it follows immediately that xk x . Choose > 0 such that (3.8) holds whenever xk N (x ) and k is sufficiently large and also such that F (y) - F (x) - F (x)(y - x) 1-u y-x (3.9)
whenever x, y N (x ). Let k0 be such that if k k0 , then xk N/2 (x ) and (3.8) holds. We claim if k k0 , then the while-loop in Algorithm UTR terminates with k min{k0 , min /2}. Note that if k k0 and if sk is a trial step for which sk /2, then (3.8) and (3.9) give aredk (sk ) F (xk ) - F (xk + sk )
1-u
predk (sk ) - (1 - u)predk (sk ) u predk (sk ) From this, we see that the whileloop terminates when k /2. So, if the whileloop never reduces the radius, we have k = k . If reductions occur, the loop terminates
predk (sk ) -
F (xk ) - F (xk ) + F (xk )sk - F (xk + sk ) - F (xk ) - F (xk )sk sk
36
on or before the iteration that first brings k /2, which implies k min /2. So k min{k , min /2}. (3.10)
Furthermore, if k /2 on termination, then k+1 k ; whereas, if k > /2 on termination, then k+1 min k > min /2. Thus k+1 min{k , min /2}. Induction on (3.10) and (3.11) gives that the whileloop terminates with k min{k0 , min /2}. Finally, k min{k0 , min /2} implies lim inf k k min{k0 , min /2} > 0. Lemma 8. If x is such that F (x ) is of full rank, then there exist and > 0 such that for any > 0, s arg min and s Null(F (xk )) satisfies s { F (x) - F (x) + F (x)s } whenever x N (x ). Proof. Set K F (x )+ and let > 0 be sufficiently small that F (x) is of full rank and F (x)+ 2K whenever x N (x ). Suppose that x N (x ) and s is given by (3.12) and (3.14) for an arbitrary > 0. Denote sM P -F (x)+ F (x). sM P is the unique global minimizer of the norm of the local linear model orthogonal to the null space of F (x). We claim that s sM P . Indeed, if sM P then 37 (3.14) (3.13)
s
(3.11)
F (x) + F (x) s
(3.12)
s = sM P , otherwise the radius is too small to include sM P and s < sM P . If sM P = 0, then s = 0 and (3.14) holds trivially for any . If sM P = 0 then we know that, since s minimizes F (x) + F (x) over all s with s Null(F (x)), it s must be that F (x) + F (x)s F (x) + F (x) F (x) - F (x) + F (x)s = = = =
s sM P
sM P ; therefore,
s sM P s sM P
F (x) - F (x) + F (x) F (x) - F (x) + F (x) F (x) - F (x) - F (x) +
s sM P s - sM P sM P s sM P
sM P (-F (x)+ F (x))
F (x)
F (x)
F (x) .
Since sM P = -F (x)+ F (x) and therefore, sM P F (x)+ F (x) 2K F (x)
we have F (x) 1 MP . 2K s Hence, F (x) - F (x) + F (x)s and (3.14) holds with 2K for all x N (x ). Lemma 9. If x is not a stationary point of F , then there exist , > 0 and > 0 such that s given by (3.12) satisfies (3.14) whenever x N (x ) and 0 < < . Proof. Let > 0 be such that if x N (x ), then F (x) that F (x ) + F (x )s < F (x ) . Choose such that F (x ) + F (x )s / F (x ) < < 1.
1 2
s , 2K
F (x ) . Let s be such
38
Since F and F are continuous, there exists (0, ] such that F (x) + F (x)s F (x) whenever x N (x ). Choose (0, s ). Suppose that x N (x ) and 0 < . For s given by (3.12), we have F (x) - F (x) + F (x)s and (3.14) holds with 2 s . (1 - ) F (x ) F (x) - F (x) + F (x) F (x) - 1 - -
s s (1- ) F (x) s (1- ) F (x ) 2 s s s s s
s
F (x)
F (x) + F (x)s s s ,
Theorem 11. Assume that {xk } is a sequence produced by Algorithm UTR. Then every limit point of {xk } is a stationary point of F . If x is a limit point of {xk } such that F (x ) is of full rank, then F (x ) = 0 and xk x ; furthermore, sk = -F (xk )+ F (xk ) whenever k is sufficiently large. Proof. Assume that x is a limit point of {xk } that is not a stationary point of F . We claim that, for any > 0, there exists an > 0 such that if xk N (x ) and k is sufficiently large, then k . If this were not true, then there would exist some
39
> 0 and xkj xk such that xkj x and kj > for each j. Then 0 = =
j
lim { F (xkj ) - F (xkj+1 ) }
j
lim { F (xkj ) - F (xkj +1 ) } lim aredkj (skj )
j
j
t lim predkj (skj ) = t lim { F (xkj ) - F (xkj ) + F (xkj )skj }
j
= t lim { F (xkj ) - min
j j s
s kj
F (xkj ) + F (xkj )s }
t lim { F (xkj ) - min F (xkj ) + F (xkj )s } = t { F (x ) - min F (x ) + F (x )s }.
s
However, the last righthand side must be positive since x is not a stationary point. Now let , and be as in Lemma 9. By the above claim, there exists (0, ] such that if xk N (x ) and k is sufficiently large, then k . By Lemma 9, we have sk predk (sk ) for independent of k whenever xk N (x ) and k is sufficiently large. Then Lemma 7 implies that xk x and lim inf k k > 0. But since the claim implies that k 0, this is a contradiction. Hence, x must be a stationary point. Suppose that x is a limit point of {xk } such that F (x ) is of full rank. Since x must be a stationary point, we must have F (x ) = 0. It follows from Lemma 8 that there exists a independent of k for which sk predk (sk ) whenever xk is sufficiently near x . Then Lemma 7 implies that xk x , and there exists a > 0 such that k for sufficiently large k. Since xk x and F (xk ) F (x ) = 0, we have that
k
lim F (xk )+ F (xk ) lim F (xk )+
k
F (xk ) = F (xk )+
F (xk ) = 0
40
implies that, for sufficiently large k, the Moore-Penrose step will be accepted. Thus sk = -F (xk )+ F (xk ) whenever k is sufficiently large.
3.4.1
UnderDetermined Dogleg Method
At each iteration, calculating the trustregion step requires the minimization of F (x) + F (x)s over all s such that s . Calculating the trustregion step to a high degree of accuracy is usually prohibitively expensive. Many methods of approximating the step have been developed. Examples include the locally constrained optimal "hook" step method, the dogleg step method and the double dogleg step method, see [4] and [34]. The dogleg method ([25, 24]) is the focus of this section. The method builds a piecewise linear curve, DL , which approximates the curve minimizing the local linear model in the trustregion. In the fullydetermined case, the dogleg curve connects the current point to the Cauchy point and subsequently the Newton point. The Cauchy point is defined to be "the minimizer of l(s)
1 2
F (x) + F (x)s
2 2
in the steepest
descent direction -l(0) = -F (x)T F (x), the steepest descent point,"[34] sCP = - k - F (x)T F (x) 2 T 2 F (x) F (x). F (x)F (x)T F (x) 2 2
Here, we replace the Newton point with the MoorePenrose point. We denote the Cauchy point by sCP and the MoorePenrose point by sMP . In the fullydetermined k k case, the dogleg curve, DL , intersects the trust region boundary at a single point [4]. This result is easily extended to the underdetermined case. The dogleg step is the step from the current point to the intersection point. We introduce the underdetermined Dogleg method (UDL):
41
Algorithm UDL: For k = 0 step 1 until do: If sMP , sk = sMP . k k If sMP > then do: k Compute sCP = - k Let x0 , 0 < min < max < 1 and 0 < min be given. Calculate sMP = -F (xk )+ F (xk ). k
F (x
-F (xk )T F (xk ) 2 2 T k )F (xk ) F (xk ) sCP k
2 2
F (xk )T F (xk ).
If sCP , then sk = k
sCP . k
If sCP < , then sk = sCP + (sMP - sCP ), k k k k While aredk (sk ) < t predk (sk ) Do: Update max{, min } Choose [min , max ]
where is uniquely determined by sCP + (sMP - sCP ) . k k k
Set xk+1 = xk + sk and update . The dogleg method chooses a step to minimize the linear model norm along the piecewise linear curve within the trustregion. The point where DL intersects the boundary is the minimizer within the trustregion, and can be computed analytically [20]. The following argument verifies this last statement. Let s be a step from the current point, x0 , to a point along DL . We claim that the length of s increases as DL is traced from xc to sCP to sM P . Indeed, s as we move from x0 to sCP . To show that s
2 2
2 2
Redetermine sk DL k
2 2
increases
increases from sCP to sM P we define
s()
the function s() = sCP + (sM P - sCP ) for [0, 1]. We have s()
2 2
0 because,
= =
sCP + (sM P - sCP ) sCP
2 2
2 2 2 2
+ 2 sM P - sCP
+ 2(sCP )T (sM P - sCP )
implies
s()
2 2
= 2 sM P - sCP 42
2 2
+ 2(sCP )T (sM P - sCP ),
and therefore
s()
2 2
0 if and only if (sCP )T (sM P - sCP ) 0. We introduce
J F (x) and F F (x), and this notation will be used in subsequent equations. Now (sCP )T (sM P - sCP ) = - = = = yet JT F
4 2 JT F JJ T F
2 2 2 2
2 2 2 2 JT F 2 2 JJ T F 2 2 JT F 2 F 2 2 2 JJ T F 2 2
JT F JJ T F
(J T F )T -J + F (x) +
JT F JJ T F
2 2 2 2
JT F JJ T F
2 2 2 2
(J T F )
F T JJ + F - F
2 2
F T JJ T F
-
1-
4 2 2 2 JT F 4 2 JJ T F 2 F 2 2 2
JT F JJ T F
,
= (J T F )T (J T F )(J T F )T (J T F ) = F T JJ T F F T JJ T F = F T JJ T F = JJ T F
2 2 2 2F
F
2 2,
which implies JT F 4 2 JJ T F 2 F 2 and (sCP )T (sM P - sCP ) = 0. Therefore
s()
2 2
2 2
= 1,
0. This proves the claim that the length of s increases as DL is
traversed from x0 to sCP to sM P . We also claim F (x)+J(x)s J(x)s
2 2 2 2
is decreasing along DL . It has been shown F (x)+
decreases along the dogleg curve from x0 to sCP [4]. We must also show that
2 2
F (x) + J(x)s
decreases from sCP to sM P along DL . Let s() = sCP + (sM P - sCP )
Then F + Js()
2 2
= F T F + 2[J T F ]T (sCP + (sM P - sCP )) 43
+(sCP + (sM P - sCP ))T J T J(sCP + (sM P - sCP )).
Taking the derivative with respect to ,
F +Js()
2 2
= 2(J T F )T (sM P - sCP ) + (sCP )T J T J(sM P - sCP )+ = 2[(J T F )T + (sCP )T J T J](sM P - sCP )+ (sM P - sCP )T J T J(sM P - sCP ) J(sM P - sCP ) ,
F +Js()
2 2
(sM P - sCP )T J T JsCP + 2(sM P - sCP )T J T J(sM P - sCP )
= 2[(J T F )T + (sCP )T J T J](sM P - sCP )+
from which it follows that
is an increasing function of . So if the right
2 2
hand side is not positive at = 1, then F + Js() know F + Js()
2 2
is decreasing. However, we
= 0 at sM P , and so the righthand side is negative for [0, 1].
2 2
Hence F (x) + J(x)s
is decreasing along DL .
Finally, we must verify that steps along the dogleg curve are orthogonal to the null space of the Jacobian. It is already known that the MoorePenrose step is orthogonal to the null space. Assume t is an element of Null(F (x)). The Cauchy step is sCP = -
-F (x)T F (x) 2 2 F (x)F (x)T F (x)
2 2
F (x)T F (x). Then - F (x)T F (x) 2 T 2 F (x) F (x))T t (x)F (x)T F (x) 2 F 2 - F (x)T F (x) 2 2 F (x)T F (x)t F (x)F (x)T F (x) 2 2 - F (x)T F (x) 2 2 F (x)T 0 F (x)F (x)T F (x) 2 2
(sCP )T t = (- = - = - = 0
The dogleg curve consists of linear combinations of zero, the Cauchy step, and the MoorePenrose step; therefore, we can conclude that steps along the dogleg curve are orthogonal to the null space of the Jacobian.
44
Chapter 4 Numerical Experiments
The numerical experiments do not provide an extensive comparison of the methods but are meant to highlight three things. First, these new methods are able to solve under-determined systems of nonlinear equations. Second, the methods achieve fast rates of convergence when near a solution of the problem. Finally, the inexact methods are computationally more efficient than the exact methods.
4.1
Test Problems
The test problems arise from typical test problems for nonlinear system solvers. Here, they are modified to be underdetermined problems.
4.1.1
The Bratu Problem
The Bratu (or Gelfand) problem is a nonlinear eigenvalue problem of the form u + eu = 0, in , u = 0 on . (4.1)
A detailed description of the Bratu problem can be found in [8] and [19], with additional information on its solution in [6] and [33]. Figure 4.1, provided by H. Walker, 45
it shows the solution space for the Bratu problem along with a typical solution. In practice, an initial u is calculated by fixing and then applying a nonlinear solver to the arising system. For our tests, we treat as an additional unknown, and solve the corresponding underdetermined system. For our tests, we assume that = [0, 1] [0, 1]. We discretized using centered differences on a 50 50 uniform grid. This leads to 2501 total unknowns in 2500 equations. For our tests we used an initial guess of u = 2 sin(x) sin(y), = 7.0. Previous work [8] shows that the Newton equation is difficult to solve for fine grids. Therefore, preconditioning the linear system is often necessary. As done in [33], we rightprecondition using a Poisson solver.
4.1.2
The Chan Problem
The Chan problem is a nonlinear eigenvalue problem similar to the Bratu problem. A description can be found in [1] with solutions in [34] and [33]. The problem is u + 1 + u + u2 /2 1 + u2 /100 = 0, in , u = 0 on . (4.2)
Figure 4.2 shows a solution of the Chan problem. For our tests we assume that = [0, 1] [0, 1] and is treated as an unknown. We discretized using centered differences on a 50 50 uniform grid. This leads to 2501 total unknowns in 2500 equations. The initial guess for the Chan problem is u = 1.0, = 0.0
46
12
10 12 8 10
6
8
6 4 4 2 1
2 0.5
0 0
0 2 4 6 0 0.5 1 0
Figure 4.1: The left panel depicts the solution space of the discretized Bratu problem. The x-axis is the value and the y-axis is u . The right panel is a representative solution of the problem.
4.1.3
The LidDriven Cavity Problem
The liddriven cavity problem involves a confined fluid flow in a square box. Circulation of the fluid is driven by a moving top boundary, or "lid". The two sides and bottom are held fixed while the top is moving from left to right. We use the stream function formulation of this problem. The Reynolds number is a nondimensional parameter; as the Reynolds number is increased, areas of counter circulation appear in the corners of the domain, and the problem becomes increasingly difficult to solve [22]. We denote the stream function by u, and the Reynolds number by Re. Usually, this problem is solved by fixing the Reynolds number and solving for the stream 47
Chan Problem: Solution
2
1.5
1
0.5
0 60 40 20 0 10 0 30 20 50 40
Figure 4.2: A plot of a solution of the Chan problem, at = 7.740455
function using a standard nonlinear solver. Here, we treat the Reynolds number as an additional unknown, and solve the corresponding underdetermined system. The linear problem is preconditioned using a biharmonic solver as in [33]. The discretization of the domain is a 40 40 equallyspaced grid. Treating the Reynolds number as an unknown leads to 1601 unknowns in 1600 equations. The code for this problem was provided H. Walker. The problem is a fourth order system 1 2 u - (uy ux + ux uy ) = 0 on Re with u = 0 on ,
u n u n
(4.3)
= 0 on the sides and bottom, = 1 on the top
Here, is the Laplacian operator, and 2 is the biharmonic operator,i.e., the Laplacian applied twice. The initial data for the liddriven cavity problem are u = 0, Re = 1000. Figure 4.3 contains a contour plot of a solution to the liddriven problem.
48
Lid cavity Driven Cavity Problem: Solution 1 0.9 0.8 0.7 0.6 0.5 0.4 0.3 0.2 0.1 0 0 0.2 0.4 0.6 0.8 1 0.1 0.08 0.06 0.04 0.02 0 -0.02 -0.04 -0.06 -0.08 -0.1
Figure 4.3: A plot of the stream function for the liddriven cavity with Re = 1000.00043784100.
4.1.4
The 1D Brusselator Problem
The onedimensional Brusselator problem [27] involves a coupled nonlinear system of equations derived from a hypothetical set of chemical reactions. The system of partial differential equations is u Du 2 u = 2 2 + u2 v - (B + 1)u + A t L x Dv 2 v v = 2 2 - u2 v + Bu t L x (4.4) (4.5)
where u and v represent the concentrations of different chemical species; the parameters A, B, Du and Dv are fixed at values of 2, 5.45, .008 and .004 respectively. The characteristic length L is the bifurcation parameter of the problem. The Dirichlet boundary conditions are u(t, x = 0) = u(t, x = 1) = A u(t, x = 0) = v(t, x = 1) = B/A. (4.6) (4.7)
A steady-state solution of this problem is u = A and v = B/A. From this we may form a trivial steady-state branch of the corresponding continuation problem. Hopf 49
bifurcations occur at the parameter values of Lk = 0.5130k, k = 1, 2, . . . . A Hopf bifurcation is characterized by "the appearance, from equilibrium state, of small amplitude periodic oscillations."[14] Branches of periodic solutions extend from these bifurcation points. Our goal is to solve for a periodic solution somewhere along one of these branches. We denote the period of a solution by T . The domain is divided into 32 evenly spaced points. Solving for u, v, L, and T yields 64 total unknowns in 62 equations.The initial guess for this problem is the same as used by K. Lust et al. in [17]: u = (sin(3x) + sin(x))/100 + 2, v = .3 sin(x) + 5.45/2, T = 3.017, L = .55551. Figure 4.4 depicts part of the solution space for this problem along with an example solution.
3.4
3.2
3
2.8
2.6
2.4
2.2
2
0
10
20
30
40
50
60
70
Figure 4.4: The left panel depicts the solution space for the 1D Brusselator problem. The right panel shows a plot of the initial concentration vector [u, v] for a periodic solution with period T = 3.020249 at parameter L = 0.557423.
50
4.1.5
The 2D Brusselator Problem
The twodimensional Brusselator problem is a version of the above chemical reaction occurring on a square grid [10, 11]. As for the previous problem, the goal is to find a periodic solution. However, here I fix the characteristic length and solve only for the two initial concentrations of reactants, u and v, and the period T . The partial differential equation system for this reactiondiffusion problem is:
u = t
2u 2u + x2 y 2
+ 1 + u2 v - 4.4u + 3.4u - u2 v,
(4.8) (4.9)
v = t
2v 2v + x2 y 2
with = .002. On the boundary, we have Neumann boundary conditions: u = 0, n v = 0. n (4.10) (4.11)
We discretized on a 21 21 uniform grid. With the addition of T as an unknown, there are 883 unknowns in 882 equations. The initial conditions are given by u = 0.5 + y, v = 1 + 5x, T = 7.5.
4.2
The Solution Algorithms
For the numerical tests, I coded four methods in MATLAB. They are NMU of section 2.3, INMU of section 3.1, the backtracking method of section 3.3 using quadratic backtracking (QINMU), and UDL of section 3.4.1. The methods are applied to each 51
2D Brusselator Problem: u
2D Brusselator Problem: v
6 4 2 0 30 20 10 0 0 10 20
6 4 2 0 30 20 10 0 0 10 20
30
30
Figure 4.5: The initial concentrations u (left panel) and v (right panel) for a periodic solution of the 2D Brusselator problem with period T = 7.47997.
of the five test problems from the previous section. Remember, that NMU is not a new method. It is the model method for comparisons with the the subsequent methods and is used here as a control. Initial inexact Newton steps were calculated using the GMRES method described in Appendix A. The methods discussed in the previous sections have many parameters that must be set, e.g., the GMRES restart value, maximum forcing term, etc. These parameters affect the performance of the methods. I chose parameters commonly used in the literature. This section lists some of the more important parameters and the values chosen. GMRES: For the Bratu Problem, the Chan Problem and the 1D Brusselator Problem, I used GMRES with a restart value of 20; the maximum allowed number of iterations is 100. For the 2D Brusselator and the LidDriven Cavity problem, the restart value is 50, with a maximum of 500 allowed iterations. 52
Forcing Terms: For the forcing terms k used in both INMU and QINMU, I used Choice 1 from [6], with 0 = .9, max = .9, min = .1. Backtracking Parameters: I used min = .1, max = .5; the maximum number of backtracks allowed is 20. Both inexact methods (INMU,QINMU) require an orthonormal basis of the null space of the Jacobian at each iteration. To accomplish this, we first build and store the full Jacobian at x0 , and use MATLAB's null command. The call returns an orthonormal basis for the null space of the Jacobian. It is calculated by a singular value decomposition. At each subsequent iteration, a corrector to the null space is calculated using the adapted GMRES method discussed in the appendix. Mathematically, if B is an orthonormal basis for the null space of F , with columns bi , then, at each iteration we calculate correctors, bi , by solving F (bi + bi ) = 0, or F bi = -F bi , and orthonormalizing the resulting updated vectors. For further discussion about the computational expense in the largescale case, see the Summary chapter.
4.3
Results
NMU successfully found a solution for each of the Bratu, Chan, and 1D Brusselator problems. However, the method failed to solve the 2D Brusselator problem and the Driven Cavity problem. In the case of the former, NMU returned a trivial solution; a solution with zero period. INMU, QINMU, and UDL successfully found solutions to 53
Method/Problem Bratu Chan Driven NMU 2403 3130 25253 INMU 722 1150 372 QINMU 720 1174 541 UDL 2347 3152 4068
1D Bruss 2D Bruss 573 13661 828 1746 874 1545 572 21349
Figure 4.6: The total compute times, in seconds, for the five test problems. Method/Problem Bratu Chan Driven NMU 801.00 782.50 252.53 INMU 90.25 115.00 10.05 QINMU 90.00 117.40 16.39 UDL 782.33 788.00 339.0 1D Bruss 2D Bruss 191 525.42 118.29 47.19 124.86 30.90 190.67 533.73
Figure 4.7: The compute times per nonlinear iteration, in seconds, for the five test problems. every test problem given. They did not always find the same solution. For example, on the Bratu problem INMU converged to a solution with = 6.569 while UDL converged to a solution with = 6.349. Table 4.6 shows the amount of time, in seconds, the methods needed to solve each problem. Table 4.7 takes the same times and scales them by the number of iterations. We see INMU and QINMU are indeed more computationally efficient than NMU. UDL has times comparable to NMU. This is expected because UDL calculates the exact MoorePenrose step at each nonlinear iteration. Figure 4.8 plots the iteration number against the log of the nonlinear residual norm for each method and problem. The Bratu, Chan, and 1D Brusselator problems did not require the globalizations of QINMU and UDL. In these cases, the methods INMU and QINMU produce the same nonlinear residual norms. The methods NMU and UDL, also, produce the same values at each iteration. Because of the overlap in the plots, circles and diamonds are used for plotting the data from QINMU and UDL.
54
Bratu Problem 4 NMU INMU QINMU UDL 5
Chan Problem NMU INMU QINMU UDL
2
0 0 log10(||F||) -2 log10(||F||) -5 -6 0 1 2 3 4 Iteration 5 6 7 8 -10 0
-4
-8
-10
2
4 Iteration
6
8
10
1D Brusselator Problem 0 NMU INMU QINMU UDL 4
2D Brusselator Problem NMU INMU QINMU UDL
-1
2
-2
0
log10(||F||)
-3
log10(||F||) 1 2 3 4 Iteration 5 6 7
-2
-4
-4
-5
-6
-6
-8
-7
-10
0
10
20
30 Iteration
40
50
60
Driven Cavity Problem 15
10
5 log10(||F||) NMU INMU QINMU UDL
0
-5
-10
-15
0
20
40 Iteration
60
80
100
Figure 4.8: For each problem we plot the log10 ( F ) at each nonlinear iteration. Where the plots of different methods overlap, we use circles and diamonds to label the iterations.
55
Chapter 5 Summary
5.1 Summary
This work introduces three classes of methods for solving underdetermined systems. Chapter 2 presents background material, including Newton's method for under determined systems (NMU). Section 3.1 introduces inexact Newton methods for underdetermined systems (INMU). Included in 3.1 is a local convergence theory for these new methods. We show, theoretically, that the methods have fast local convergence rates under appropriate assumptions on the forcing terms. In 3.2 we seek to improve the robustness of INMU by imposing a sufficientdecrease condition. This leads to the globalized inexact Newton methods for underdetermined systems (GINMU) and, in 3.3, the backtracking method (BINMU). Important here is that BINMU becomes INMU close to a solution; therefore, BINMU can achieve fast local rates of convergence under appropriate assumptions on the forcing terms. Finally, in 3.4, we adapt a general trustregion method to solve an underdetermined system. This section includes a discussion of the underdetermined dogleg method (UDL), a specific underdetermined trustregion method (UTR). A general convergence theory for these methods is presented. 56
Chapter 4 presents the results from preliminary numerical experiments. Five test problems were coded in MATLAB, and solved with each of four methods: NMU, INMU, QINMU (BINMU with quadratic steplength reduction), and UDL. We found that the three new methods all produced solutions of the problems; additionally, they often exhibited the predicted fast local rates of convergence. The two inexact methods, INMU and QINMU, were computationally more efficient than NMU and UDL, probably because each of the latter required a full linear solve for each nonlinear iteration.
5.2
Additional Applications
The methods presented in this dissertation were originally developed for solving under-determined systems arising from the discretization of parameterdependent partial differential equations. Often, the dimension of the null space of the Jacobian in these problems is only of dimension one or two. However, these methods are applicable to problems with a null space of much larger dimension. Another application is solving continuation, homotopy and bifurcationtracking problems. The new methods may be utilized in two ways. First, one may use an underdetermined system method to find an initial point in the solution set of one of these problems, as done above in the Bratu, Chan, and liddriven cavity problems. Additionally, an underdetermined system method may be used for the corrector iterations in a predictor/corrector method for tracing the solution set.
5.3
Future Work
The next logical step, for further study of these methods, is to code the methods in C + +. This would allow for application to much larger problems. (The current 57
MATLAB code is limited in the number of total unknowns allowed.) Additionally, a parallel implementation would further increase the maximum allowable size of the problems. An area which still must be explored is the efficient calculation of the vectors spanning the null space of the Jacobian. The k th nonlinear iteration of our inexact Newton methods requires an orthonormal basis of the null space of F (xk ). Section 4.2 discusses the method used in our numerical tests. To recap, we build and store the full Jacobian at the initial guess, x0 and use MATLAB's null command to calculate orthonormal spanning vectors. Each subsequent iteration updates the spanning vectors individually. In order to make these methods viable for largescale simulations, we must find a more computationally efficient method of obtaining the initial nullspace basis. A possibility is to begin with a random initial guess for each vector and use the update procedure outlined in 4.2, perhaps repeatedly. Another possible avenue of research is inexact dogleg methods for underdetermined systems. Previously, a general dogleg method was extended to the inexact Newton context in the case of a fullydetermined system [23]. In [23], the second point of the dogleg curve (the Newton step) is replaced by an approximation. A similar idea may be applied to the underdetermined case; the MoorePenrose step would be replaced with an inexactNewton step. A global convergence analysis similar to that in [23] could then be performed.
58
Appendices
59
Appendix A Adaptation of GMRES
GMRES is a Krylov subspace method. These methods are designed to solve iteratively the linear problem: find x Rn such that Ax = b, A Rnn , b Rn . A Krylov subspace method begins with an initial x0 and at the k th step, determines an iterate xk through a correction in the k th Krylov subspace Kk span{r0 , Ar0 , . . . , An-1 r0 }, (A.1)
where r0 b - Ax0 is the initial residual. A strong attraction of these methods is that implementations only require products Av and sometimes AT v. Thus, within the methods, no direct access to or manipulation of the entries of A is required. GMRES uses only Av products. Detailed descriptions of GMRES and other Krylov subspace methods can be found in [3, 29, 28]. In [33], Walker adapts Krylov subspace methods to solve the problem: find x IR n+1 such that Ax = b, A IR n(n+1) , b IR n . This involves imposing the additional constraint that x Null(A) on the solution to guarantee uniqueness. It is our goal to extend his method to handle the more general case: find x IR n+p such that Ax = b, A IR n(n+p) , b IR n . First note that this linear system is underdetermined; an additional constraint must be imposed to have a unique solution. 60
Following the method in [33], assume the constraint is of the form T T x = 0, T IR (n+p)p , (A.2)
where the columns of T form an orthonormal basis of the null space of A. This condition is equivalent to choosing the solution of Ax = b with minimum Euclidean norm. The adapted Krylov methods work by imposing the constraint (A.2) directly on the iterations of standard Krylov methods. This is done as follows: 1. Find Q IR (n+p)n such that Range(Q) = {colspaceT } and Qy all y IR n . Then AQ IR nn . 2. Apply the Krylov subspace method to solve approximately AQy = b for y IR n . Then set x = Qy. Just as in [33], x = Qy satisfies (A.2) regardless of how well it approximately satisfies Ax = b. Use p Householder transformations to transform T into a triangular matrix. 0 0 0 . . ... 0 . . . . 0 0 x Pp P2 P1 T = 0 x x . . . . . . . . . x x x x The product of these transformations with the (n + p) n identity matrix yields an acceptable Q. Algorithm HH forms our Householder transformation vectors [9]. = y for
2
2
Algorithm HH: Let T IR (n+p)p be given. 61
for j = 1 : p
(j) = T (1 : n - j + 1, j) ; for k = j + 1 : p
T (n - j + 1, j) = T (n - j + 1, j) + sign(T (n - j + 1, j)) (j); T (1 : n - j + 1, k) = T (1 : n - j + 1, k)... T (1 : n - j + 1, j) T (1 : n - j + 1, k); -2/ T (1 : n - j + 1, j)
2
T (1 : n - j + 1, j)...
end end
This method produces the Householder vectors ui in the columns of T . The Householder transformation matrices are then formed by Pi = I -
2 ui
2 2
ui uT . Let Ip i
denote the (n + p) (n + p) identity matrix with the final p columns deleted. The matrix Q is then formed as Q = P1 . . . Pp-1 Pp Ip . Once Q has been created, it is straightforward to apply a Krylov method.
62
Appendix B Matlab Files
NMU.m
0001 0002 0003 0004 0005 0006 0007 0008 0009 0010 0011 0012 0013 0014 0015 0016 0017 0018 0019 0020 0021 0022 0023 0024 0025 0026 0027 0028 0029 0030 0031 function [x,resids]=NMU(x,f,jac,tol,maxits) % Newton's Method for Under-determined systems % Author: Joseph Simonis % Latest update: 03-01-06 % % x initial guess of solution % f function to compute F(x) % jac function to compute J(x) % tol solution tolerance % maxits maximum number of nonlinear iterations % F=feval(f,x); residual=norm(F); its=1; resids(its,1)=residual; fprintf('\nIt.No. ||F(u)|| fprintf(' %d %e %c
GMRES Its. Lin Mod Norm Eta \n'); %e %c\n', 0,residual,'*',0,'*');
while(residual > tol & its<maxits) J=feval(jac,x); s=-pinv(full(J))*F; % Solve the under-determined lin. sys. x=x+s; % Update x linnorm=norm(J*s+F); F=feval(f,x); residual=norm(F); fprintf(' %d %e %c %e %c\n',... its,residual,'*',linnorm,'*'); its=its+1; resids(its,1)=residual; end 63
QINMU.m
0001 0002 0003 0004 0005 0006 0007 0008 0009 0010 0011 0012 0013 0014 0015 0016 0017 0018 0019 0020 0021 0022 0023 0024 0025 0026 0027 0028 0029 0030 0031 0032 0033 0034 0035 0036 0037 0038 0039 0040 0041 0042 0043 0044 0045 0046 0047 0048 function [x,resids,fail count]=QINMU(x,f,jac,jacv,tol,maxits,... use precond,pcond fun,eta choice,eta0) % This function was written as a simple implementation % of the QINMU method from my dissertation. % Inexact Newton Method for Under-determined systems % using quadratic backtracking. % % Author: Joseph Simonis % Latest update: 03-01-06 % % x initial guess of solution % f function to compute F(x) % jac function to compute J(x) % jacv function to compute J(x)*v % tol solution tolerance % maxits maximum number of iterations % use precond flag for turning on preconditioning % pcond fun function to perform preconditioning % eta choice flag for choosing eta choice % eta 0 the initial eta value % If eta choice==0 eta=constant % If eta choice==1 use eta=|F-linear residual|/Fprev % If eta choice==2 use eta=gamma*(F/Fprev)^alpha % % % etamin=.1; etamax=.9; maxbtsteps=20; thetamin=.1; thetamax=.9; gmres restart=20; gmres max=100; % %We will keep track of the number of backtracking failures % in with the fail count fail count=0; % The first step is to compute the Null vectors of the % Jacobian. A=feval(jac,x); t=null(full(A)); for j=1:size(t,2) t(:,j)=t(:,j)./norm(t(:,j)); % Scale the null vectors end
64
0049 0050 0051 0052 0053 0054 0055 0056 0057 0058 0059 0060 0061 0062 0063 0064 0065 0066 0067 0068 0069 0070 0071 0072 0073 0074 0075 0076 0077 0078 0079 0080 0081 0082 0083 0084 0085 0086 0087 0088 0089 0090 0091 0092 0093 0094 0095 0096 0097 0098 0099 0100
% Begin Method F=feval(f,x); % Evaluate F and the residual r=||F||. residual=norm(F) n=length(F); its=1; % Nonlinear iterations count. resids(its,1)=residual; % For output eta=eta0; % The initial forcing term. fprintf('\nIt.No. ||F(u)|| fprintf(' %d %e %d GMRES Its. Lin Mod Norm Eta \n'); %e %e\n', 0,residual,0,0,eta0);
while(residual > tol & its<maxits) % Here I want to update the null vectors t. % To do this I solve J(deltat)=-Jt for a correction % to t for each direction in the null space. for k=1:size(t,2) temp = feval(jacv,x,t(:,k)); [deltat(:,k),error,dummy]=GMRES House(jacv,temp,x,gmres restart,... (1.0e-3)/norm(temp),gmres max,t,use precond,pcond fun); t(:,k)=t(:,k)+deltat(:,k); end for j=1:size(t,2) t(:,j)=t(:,j)./norm(t(:,j));% Scale t end % Scale t % Now call the under-determined GMRES method. [s,rho,ftjs,succ,linits]=GMRES House(jacv,F,x,gmres restart,eta,... gmres max,t,use precond,pcond fun); % Catching errors... if ftjs >= 0, error('IN step is not a descent direction.'); end if (eta choice==0) %constant eta [s,F,residual,Failure]=quadbt(x,residual,s,eta,ftjs,thetamin,... thetamax,f,n+1,maxbtsteps); x=x+s; residual=norm(F); if (Failure==1) % If the linear solve failed to meet solve tol. fail count=fail count+1; end else Fprev=F; fnrmprev=residual; etaprv=eta; [s,F,residual,Failure,etaused]=quadbt(x,residual,s,eta,ftjs,... thetamin,thetamax,f,n+1,maxbtsteps); % Recalculate rho here. The step has been % backtracked, so I need to give etaupdate \|F+J\theta*s\|, not % \|F+J*s\|. %-----------------FpJS=Fprev+feval(jacv,x,s); 65
0101 0102 0103 0104 0105 0106 0107 0108 0109 0110 0111 0112 0113 0114 0115 end
fpjsnorm=norm(FpJS); %------------------x=x+s; [eta]=etaupdate(eta choice,fnrmprev,residual,fpjsnorm,... etamax,etamin,eta); if (Failure==1) fail count=fail count+1; end end fprintf(' %d %e %d linits,rho,etaprv); its=its+1; resids(its,1)=residual; %e %e\n',its,residual,...
66
INMU.m
0001 0002 0003 0004 0005 0006 0007 0008 0009 0010 0011 0012 0013 0014 0015 0016 0017 0018 0019 0020 0021 0022 0023 0024 0025 0026 0027 0028 0029 0030 0031 0032 0033 0034 0035 0036 0037 0038 0039 0040 0041 0042 0043 0044 0045 0046 0047 0048 function [x,resids]=INMU(x,f,jac,jacv,tol,maxits,... use precond,pcond fun,eta choice,eta0) % Inexact Newton Method for Under-determined systems % Author: Joseph Simonis % Latest update: 03-01-06 % % x initial guess of solution % f function to compute F(x) % jac function to compute J(x) % jacv function to compute J(x)*v % tol solution tolerance % maxits maximum number of iterations % use precond flag for turning on preconditioning % pcond fun function to perform preconditioning % eta choice flag for choosing eta choice % eta 0 the initial eta value % If eta choice==0 eta=constant % If eta choice==1 use eta=|F-linear residual|/Fprev % If eta choice==2 use eta=gamma*(F/Fprev)^alpha % etamin=.1; etamax=.9; maxbtsteps=20; thetamin=.1; thetamax=.9; gmres restart=20; gmres max=100; % % The first step is to compute the Null vector of the % Jacobian. A=feval(jac,x); t=null(full(A)); for j=1:size(t,2) t(:,j)=t(:,j)./norm(t(:,j));% Scale t end % Begin Method F=feval(f,x); % Evaluate F and the residual r=||F||. residual=norm(F) n=length(F); its=1; resids(its,1)=residual; eta=eta0; % The initial forcing term. fprintf('\nIt.No. ||F(u)|| fprintf(' %d %e %d GMRES Its. Lin Mod Norm Eta \n'); %e %e\n', 0,residual,0,0,eta0);
67
0049 while(residual > tol & its<maxits) 0050 % Here I want to update the null vectors t. 0051 % To do this I solve J(deltat)=-Jt for a correction 0052 % to t for each direction in the null space. 0053 for k=1:size(t,2) 0054 0055 temp = feval(jacv,x,t(:,k)); 0056 [deltat(:,k),error,dummy]=GMRES House(jacv,temp,x,... 0057 gmres restart,(1.0e-3)/norm(temp),gmres max,t,... 0058 use precond,pcond fun); 0059 t(:,k)=t(:,k)+deltat(:,k); end 0060 0061 for j=1:size(t,2) 0062 t(:,j)=t(:,j)./norm(t(:,j)); % Scale t 0063 end 0064 0065 % Now call the under-determined GMRES method. 0066 [s,rho,ftjs,succ,linits]=GMRES House(jacv,F,x,gmres restart,eta,... 0067 gmres max,t,use precond,pcond fun); 0068 fnrmprev = residual; 0069 x=x+s; 0070 F=feval(f,x); 0071 residual=norm(F); 0072 etaprv=eta; 0073 [eta]=etaupdate(eta choice,fnrmprev,residual,rho,... 0074 etamax,etamin,eta); %e %d %e %e\n', its,residual,... 0075 fprintf(' %d 0076 linits,rho,etaprv); 0077 its=its+1; 0078 resids(its,1)=residual; 0079 end
68
UDL.m
0001 0002 0003 0004 0005 0006 0007 0008 0009 0010 0011 0012 0013 0014 0015 0016 0017 0018 0019 0020 0021 0022 0023 0024 0025 0026 0027 0028 0029 0030 0031 0032 0033 0034 0035 0036 0037 0038 0039 0040 0041 0042 0043 0044 0045 0046 0047 0048 function [x,resids]=UDL(x,f,jac,tol,maxits) % Under-Determined Dogleg method % Author: Joseph Simonis % Latest update: 03-01-06 % % x initial guess of solution % f function to compute F(x) % jac function to compute J(x) % tol solution tolerance % maxits maximum number of iterations % % t=10^-4; thetamin = 0.1; thetamax = 0.5; u=.75; v=0.1; inneritsmax=20; % % Algorithm F=feval(f,x); residual = norm(F); fprintf('\nIt.No. ||F(u)|| GMRES Its. Lin Mod Norm Delta \n'); %e %c %e %e\n', 0,residual,'*',0,0); fprintf(' %d its=1; innerits=0; resids(its,1)=residual; while(residual > tol & its<maxits) J=feval(jac,x); % Calculate the Moore-Penrose step. snewt=-pinv(full(J))*F; snewtnorm=norm(snewt); if (its==1) delta=snewtnorm; end % Calculate the Dogleg step. dogleg step = Dogleg(F,J,snewt,snewtnorm,delta); Fpls = feval(f,x+dogleg step); Fplsn = norm(Fpls); Js = J*dogleg step; lin res = norm(F+Js); ared = residual-Fplsn; pred = residual-lin res; % Inner Dogleg loop. while (ared<t*pred & innerits < inneritsmax); if (snewtnorm < delta) delta = snewtnorm; 69
0049 0050 0051 0052 0053 0054 0055 0056 0057 0058 0059 0060 0061 0062 0063 0064 0065 0066 0067 0068 0069 0070 0071 0072 0073 0074 0075 0076 0077 0078 0079 0080 0081 0082 0083 0084 0085 0086 0087 end
end d = Fplsn^2-residual^2-2*F'*Js; if (d <= 0) theta = thetamax; else theta = -(F'*Js)./d; if (theta > thetamax) theta = thetamax; end if (theta < thetamin) theta = thetamin; end end delta = theta*delta; % Update Delta % Recalculate dogleg step. dogleg step = Dogleg(F,J,snewt,snewtnorm,delta); Fpls = feval(f,x+dogleg step); Fplsn = norm(Fpls); Js = J*dogleg step; lin res = norm(F+Js); ared = residual-Fplsn; pred = residual-lin res; innerits=innerits+1 end innerits = 0; x=x+dogleg step; F=Fpls; residual=Fplsn; %e %c %e %e\n', its,residual,... fprintf(' %d '*',lin res,delta); its = its+1; resids(its,1) = residual; % Update delta if (ared > u.*pred & snewtnorm > delta) delta = 2.*delta; elseif (ared < v.*pred) delta = .5.*delta; end
70
Dogleg.m
0001 0002 0003 0004 0005 0006 0007 0008 0009 0010 0011 0012 0013 0014 0015 0016 0017 0018 0019 0020 0021 0022 0023 0024 0025 0026 0027 0028 0029 function [dogleg step]=Dogleg(F,J,snewt,snewtnorm,delta) % Computes the dogleg step % Inputs: % F = F(xcurrent) % J = J(xcurrent) % snewt = Moore-Penrose Step % snewtnorm = norm of Step % delta = current trust region radius if (snewtnorm <= delta) dogleg step = snewt; else JTF=J'*F; JTFnorm=norm(JTF); JJTFnorm=norm(J*JTF); sdescent = -(JTFnorm./JJTFnorm)^2.*JTF; sdescentnorm = norm(sdescent); if (sdescentnorm >= delta) dogleg step = (delta./sdescentnorm).*sdescent; else sdiff = snewt-sdescent; a = norm(sdiff).^2; b = sdescent'*sdiff; c = sdescentnorm.^2-delta.^2; tao = -c./(b+sqrt(b.^2-a*c)); dogleg step = sdescent + tao.*sdiff; end end
71
GMRES_House.m
0001 0002 0003 0004 0005 0006 0007 0008 0009 0010 0011 0012 0013 0014 0015 0016 0017 0018 0019 0020 0021 0022 0023 0024 0025 0026 0027 0028 0029 0030 0031 0032 0033 0034 0035 0036 0037 0038 0039 0040 0041 0042 0043 0044 0045 0046 0047 0048 function [step,rho,ftjs,success,ittot] = GMRES House(jacv,fval,... u,m,eta,itmax,T,preconflag,pcond fun) % This function is a GMRES routine for under-determined % systems. It solves J(u)*step=-fval for step. We assume % J(u):R(n+d)->R(n). T is a normalized basis for the Null(J(u)). % It contains d vectors. This function uses a starting guess % of zeros. % % INPUTS: % jacv = routine for computing jacobian-vector products. % fval = current function value F(u) % u = current approx. solution of F(u) = 0 % m = GMRES restart value % eta = forcing term % itmax = maximum number of GMRES iterations % T = orthonormalized basis of Null(J(u)) % preconflag = 0-> no preconditioning used % = 1-> preconditioning used % pcond fun = preconditioning function % % OUTPUTS: % step = approx. solution of J(u)*step=-fval % rho = \|-fval-J(u)*step\| % ftjs = fval'*J(u)*step % success = 1 if tol was acheived % ittot = total number of gmres iterations % Compute d Householder reflections and store them in T. % Reference "An Adaptation of Krylov subspace methods to path % following problems" H. Walker 1999 d=size(T,2); n=size(u,1); for j=1:d beta(j)=norm(T(1:n-j+1,j)); T(n-j+1,j)=T(n-j+1,j)+sign(T(n-j+1,j))*beta(j); for k=j+1:d T(1:n-j+1,k)=T(1:n-j+1,k)-2/norm(T(1:n-j+1,j))^2*T(1:n-j+1,j)... *T(1:n-j+1,j)'*T(1:n-j+1,k); end end % The GMRES routine % SETUP n = size(u,1)-d; %the length of vectors in the range of J. D = zeros(d,1); % used for appending zeros onto n-vectors. r = -fval; rho = norm(r); tol = eta*rho; 72
0049 0050 0051 0052 0053 0054 0055 0056 0057 0058 0059 0060 0061 0062 0063 0064 0065 0066 0067 0068 0069 0070 0071 0072 0073 0074 0075 0076 0077 0078 0079 0080 0081 0082 0083 0084 0085 0086 0087 0088 0089 0090 0091 0092 0093 0094 0095 0096 0097 0098 0099 0100
step = zeros(size(u,1),1); % set the initial guess to zero ftjs=0; V = zeros(n,m+1); % Holds the Krylov subspace basis R = zeros(m,m); c = zeros(m,1); % s and c are used in computing the Givens rotations s = c; w = zeros(m+1,1); ittot = 0; % END SETUP % OUTER LOOP while (rho > tol & ittot < itmax) V(:,1) = r/rho; w(1) = rho; % INNER LOOP for k = 1:m ittot = ittot + 1; if ittot > itmax, break, end % % % % % % % The first step is to calculate V (k+1)=A*V k Multiplying a vector by A first requires (maybe) multiplying by a preconditioner M^-1. The second step is applying Q=P1...Pd*I which means applying the Householder reflections to the vector appended with zeros. The final step is to send it to the Jacv routine.
% Step 1 apply preconditioner and append zeros % or just append zeros. if (preconflag == 1) Hv = feval(pcond fun, V(:,k)); Hv = [Hv;D]; else Hv = [V(:,k);D]; end % Step 2 apply Householder reflections for j=d:-1:1 Hv(1:n+d-j+1,1)=Hv(1:n+d-j+1,1)-2/norm(T(1:n+d-j+1,j))^2 ... *T(1:n+d-j+1,j)*T(1:n+d-j+1,j)'*Hv(1:n+d-j+1,1); end % Step 3 Multiply by the Jacobian Matrix. V(:,k+1) = feval(jacv,u,Hv); % With V (k+1) calculated it is time to continue with GMRES % as normal. for i = 1:k R(i,k) = V(:,k+1)'*V(:,i); V(:,k+1) = V(:,k+1) - R(i,k)*V(:,i); end for i = 1:k-1 73
0101 0102 0103 0104 0105 0106 0107 0108 0109 0110 0111 0112 0113 0114 0115 0116 0117 0118 0119 0120 0121 0122 0123 0124 0125 0126 0127 0128 0129 0130 0131 0132 0133 0134 0135 0136 0137 0138 0139 0140 0141 0142 0143 0144 0145 0146 0147 0148 0149 0150 0151 0152
temp = R(i,k); R(i,k) = c(i)*temp + s(i)*R(i+1,k); R(i+1,k) = -s(i)*temp + c(i)*R(i+1,k); end tempnorm = norm(V(:,k+1)); temp = sqrt(R(k,k)^2 + tempnorm^2); c(k) = R(k,k)/temp; s(k) = tempnorm/temp; R(k,k) = temp; w(k+1) = -s(k)*w(k); w(k) = c(k)*w(k); rho = abs(w(k+1)); if rho <= tol, break, end V(:,k+1) = V(:,k+1)/tempnorm; end % END INNER LOOP % Solve for step using back substitution. for i = k:-1:1 w(i) = w(i)/R(i,i); if i>1 w(1:i-1) = w(1:i-1) - w(i)*R(1:i-1,i); end end % the step y which solves JQM^-1y=-F % is a linear combination of the V's. % Here we put y in the temp vector. tempvec = V(:,1:k)*w(1:k); % To now calculate our step correction we % apply M^-1 if neccessary, and then apply % Q. if (preconflag == 1) Hv = feval(pcond fun,tempvec); Hv=[Hv;D]; else Hv=[tempvec;D]; end for j=d:-1:1 Hv(1:n+d-j+1,1)=Hv(1:n+d-j+1,1)-2/norm(T(1:n+d-j+1,j))^2 ... *T(1:n+d-j+1,j)*T(1:n+d-j+1,j)'*Hv(1:n+d-j+1,1); end % Update the step. step = step + Hv; V(:,m+1) = feval(jacv,u,step); ftjs = fval'*V(:,m+1); if ittot > itmax, break, end r = - fval - V(:,m+1); rho = norm(r); 74
0153 0154 0155 0156 0157 0158 0159
end if (rho<tol) success = 1; else success = 0; end % END OUTER LOOP
75
quadbt.m
0001 0002 0003 0004 0005 0006 0007 0008 0009 0010 0011 0012 0013 0014 0015 0016 0017 0018 0019 0020 0021 0022 0023 0024 0025 0026 0027 0028 0029 0030 0031 0032 0033 0034 0035 0036 0037 0038 0039 0040 0041 0042 0043 0044 0045 0046 0047 0048 function [step,trialf,trialn,Fail,eta]=quadbt(xcur,fcnrm,step,eta,... oftjs,thetamin,thetamax,fh,meshsize,maxbtsteps) % Quadratic Backtracking Method % Author: Joseph Simonis % Latest update: 03-01-06 %Find a suitable step through bactracking, also return new eta. % % % % % % % % % % % % % % % % % % INPUT xcur fcnrm step eta oftjs thetamin thetamax fh meshsize maxbtsteps OUTPUT step trialf trialn Fail current value of x norm of F at xcur initial trial step the forcing term F'JS the minimum scaling factor per iteration the maximum scaling factor per iteration handle for function evaluations the size of the mesh maximum allowable backtracking steps final step F at xcur+step ||trialf|| 1 if backtracking failed to produce an acceptable step 0 otherwise
Fail=0; t=10^-4; accept='no'; redfac=1.0; ibt=0; %Bactracking iterations. while (strcmp(accept,'no') & ibt<maxbtsteps) trials=xcur+step; %Take the step trialf=feval(fh,trials); %Determine f at the new value. trialn=norm(trialf); %Find the norm % Uncomment the following for printing % fprintf('trialn='); % fprintf(' %e\n',trialn); % fprintf('(1-t*(1-eta))*fcnrm'); % fprintf(' %e\n',(1-t*(1-eta))*fcnrm); if trialn<=(1-t*(1-eta))*fcnrm %Test our condition residual reduction accept='yes'; else ibt=ibt+1; %Find theta to reduce our step size. phi=trialn^2-fcnrm^2-2*oftjs*redfac; 76
0049 if phi <= 0 0050 theta=thetamax; 0051 else 0052 theta=-(oftjs*redfac)/phi; 0053 end 0054 % Keep theta within bounds. 0055 if theta<thetamin 0056 theta=thetamin; 0057 end 0058 if theta>thetamax 0059 theta=thetamax; end 0060 0061 step=theta*step; 0062 %Increase eta. 0063 eta=1-theta*(1-eta); 0064 %Update the reduction factor 0065 redfac=theta*redfac; 0066 stepnorm=norm(step); if (stepnorm)<sqrt(eps) break;end 0067 0068 end 0069 end 0070 if (strcmp(accept,'no')) 0071 % We have Backtracking failure, so we take the full step. 0072 step = (1/redfac)*step; 0073 disp('Backtracking Failure: Taking full step'); 0074 Fail=1; 0075 trials=xcur+step; %Take the step 0076 trialf=feval(fh,trials); %Determine f at the new value. 0077 trialn=norm(trialf); %Find the norm 0078 0079 end
77
Appendix C Raw Data
C.1
NMU It.No. 0 1 2 3 4 INMU It.No. 0 1 2 3 4 5 6 7 8 9 10 F (u) GMRES Its. Lin Mod Norm Eta 3.751216e + 004 0 0.000000e + 000 9.000000e - 001 1.700332e + 004 1 1.700252e + 004 9.000000e - 001 1.186035e + 004 1 1.185984e + 004 8.432626e - 001 6.940770e + 003 2 6.938583e + 003 7.589363e - 001 3.719799e + 003 3 3.714728e + 003 6.399826e - 001 1.608581e + 003 6 1.595137e + 003 4.857060e - 001 4.998264e + 002 11 4.810233e + 002 3.108434e - 001 6.449217e + 001 18 6.253675e + 001 1.509785e - 001 2.355349e - 001 38 2.413941e - 001 3.912209e - 003 2.137906e - 005 64 2.049753e - 005 9.085024e - 005 8.904233e - 011 99 7.830998e - 011 3.742667e - 006 F (u) GMRES Its. Lin Mod Norm Eta 3.751216e + 004 0.000000e + 000 3.318422e + 002 2.601388e - 009 1.627407e + 000 2.896797e - 010 9.151679e - 005 1.796905e - 012 4.070212e - 011 1.151032e - 016
Chan Problem Results
78
QINMU It.No. 0 1 2 3 4 5 6 7 8 9 10 UDL It.No. 0 1 2 3 4 F (u) Inner Its. Lin Mod Norm Delta 3.751216e + 004 0.000000e + 000 0.000000e + 000 3.318422e + 002 0 2.601388e - 009 2.613503e + 001 1.627407e + 000 0 2.896797e - 010 2.613503e + 001 9.151679e - 005 0 1.796905e - 012 2.613503e + 001 4.070212e - 011 0 1.151032e - 016 2.613503e + 001 F (u) GMRES Its. Lin Mod Norm Eta 3.751216e + 004 0 0.000000e + 000 9.000000e - 001 1.700332e + 004 1 1.700252e + 004 9.000000e - 001 1.186035e + 004 1 1.185984e + 004 8.432626e - 001 6.940770e + 003 2 6.938583e + 003 7.589363e - 001 3.719799e + 003 3 3.714728e + 003 6.399826e - 001 1.608581e + 003 6 1.595137e + 003 4.857060e - 001 4.998264e + 002 11 4.810233e + 002 3.108434e - 001 6.449217e + 001 18 6.253675e + 001 1.509785e - 001 2.355349e - 001 38 2.413941e - 001 3.912209e - 003 2.137906e - 005 64 2.049753e - 005 9.085024e - 005 8.904233e - 011 99 7.830998e - 011 3.742667e - 006
C.2
NMU
Bratu Problem Results
It.No. 0 1 2 3 F (u) GMRES Its. Lin Mod Norm Eta 3.391596e + 003 0.000000e + 000 2.984006e + 000 2.878279e - 010 1.141729e - 003 3.152350e - 012 9.795232e - 011 1.518800e - 015
INMU It.No. 0 1 2 3 4 5 6 7 8 F (u) GMRES Its. Lin Mod Norm Eta 3.391596e + 003 0 0.000000e + 000 9.000000e - 001 9.236813e + 002 1 9.219854e + 002 9.000000e - 001 2.311795e + 002 1 2.251650e + 002 8.432626e - 001 2.090284e + 001 2 2.156339e + 001 7.589363e - 001 2.293760e + 000 2 2.313649e + 000 6.399826e - 001 5.251453e - 001 1 5.251154e - 001 4.857060e - 001 1.297609e - 002 2 1.298191e - 002 3.108434e - 001 1.556908e - 003 2 1.556908e - 003 1.509785e - 001 2.685922e - 010 7 1.496858e - 012 1.614431e - 008
79
QINMU It.No. 0 1 2 3 4 5 6 7 8 UDL It.No. 0 1 2 3 F (u) Inner Its. Lin Mod Norm Delta 3.391596e + 003 0.000000e + 000 0.000000e + 000 2.984006e + 000 0 2.878279e - 010 2.883900e + 000 1.141729e - 003 0 3.152350e - 012 2.883900e + 000 9.795232e - 011 0 1.518800e - 015 2.883900e + 000 F (u) GMRES Its. Lin Mod Norm Eta 3.391596e + 003 0 0.000000e + 000 9.000000e - 001 9.236813e + 002 1 9.219854e + 002 9.000000e - 001 2.311795e + 002 1 2.251650e + 002 8.432626e - 001 2.090284e + 001 2 2.156339e + 001 7.589363e - 001 2.293760e + 000 2 2.313649e + 000 6.399826e - 0...