On Methods For Solving Symmetric Systems Of Linear

Transcript

On Methods for Solving Symmetric Systems of Linear Equations Arising in Optimization TOVE ODLAND Doctoral Thesis Stockholm, Sweden 2015 TRITA-MAT-A 2015:05 ISRN KTH/MAT/A-15/05-SE ISBN 978-91-7595-514-8 Department of Mathematics Royal Institute of Technology 100 44 Stockholm, Sweden Akademisk avhandling som med tillstånd av Kungl Tekniska Högskolan framlägges till offentlig granskning för avläggande av teknologie doktorsexamen, fredagen den 12 juni 2015 klockan 14.00 i sal F3, Lindstedtsvägen 26, Kungl Tekniska Högskolan, Stockholm. © Tove Odland, 2015 Print: Universitetsservice US AB, Stockholm, 2015 To F and S.B. "You are the sky; everything else is just weather." -Ani Pema Chödrön Abstract In this thesis we present research on mathematical properties of methods for solving symmetric systems of linear equations that arise in various optimization problem formulations and in methods for solving such problems. In the first and third paper (Paper A and Paper C), we consider the connection between the method of conjugate gradients and quasi-Newton methods on strictly convex quadratic optimization problems or equivalently on a symmetric system of linear equations with a positive definite matrix. We state conditions on the quasi-Newton matrix and the update matrix such that the search directions generated by the corresponding quasi-Newton method and the method of conjugate gradients respectively are parallel. In paper A, we derive such conditions on the update matrix based on a sufficient condition to obtain mutually conjugate search directions. These conditions are shown to be equivalent to the one-parameter Broyden family. Further, we derive a one-to-one correspondence between the Broyden parameter and the scaling between the search directions from the method of conjugate gradients and a quasi-Newton method employing some well-defined update scheme in the one-parameter Broyden family. In paper C, we give necessary and sufficient conditions on the quasi-Newton matrix and on the update matrix such that equivalence with the method of conjugate gradients hold for the corresponding quasi-Newton method. We show that the set of quasiNewton schemes admitted by these necessary and sufficient conditions is strictly larger than the one-parameter Broyden family. In addition, we show that this set of quasiNewton schemes includes an infinite number of symmetric rank-one update schemes. In the second paper (Paper B), we utilize an unnormalized Krylov subspace framework for solving symmetric systems of linear equations. These systems may be incompatible and the matrix may be indefinite/singular. Such systems of symmetric linear equations arise in constrained optimization. In the case of an incompatible symmetric system of linear equations we give a certificate of incompatibility based on a projection on the null space of the symmetric matrix and characterize a minimum-residual solution. Further we derive a minimum-residual method, give explicit recursions for the minimum-residual iterates and characterize a minimum-residual solution of minimum Euclidean norm. Keywords: symmetric system of linear equations, method of conjugate gradients, quasi-Newton method, unconstrained optimization, unconstrained quadratic optimization, Krylov subspace method, unnormalized Lanczos vectors, minimum-residual method vii Sammanfattning I denna avhandling betraktar vi matematiska egenskaper hos metoder för att lösa symmetriska linjära ekvationssystem som uppkommer i formuleringar och metoder för en mängd olika optimeringsproblem. I första och tredje artikeln (Paper A och Paper C), undersöks kopplingen mellan konjugerade gradientmetoden och kvasi-Newtonmetoder när dessa appliceras på strikt konvexa kvadratiska optimeringsproblem utan bivillkor eller ekvivalent på ett symmetrisk linjärt ekvationssystem med en positivt definit symmetrisk matris. Vi ställer upp villkor på kvasi-Newtonmatrisen och uppdateringsmatrisen så att sökriktningen som fås från motsvarande kvasi-Newtonmetod blir parallell med den sökriktning som fås från konjugerade gradientmetoden. I den första artikeln (Paper A), härleds villkor på uppdateringsmatrisen baserade på ett tillräckligt villkor för att få ömsesidigt konjugerade sökriktningar. Dessa villkor på kvasi-Newtonmetoden visas vara ekvivalenta med att uppdateringsstrategin tillhör Broydens enparameterfamilj. Vi tar också fram en ett-till-ett överensstämmelse mellan Broydenparametern och skalningen mellan sökriktningarna från konjugerade gradientmetoden och en kvasi-Newtonmetod som använder någon väldefinierad uppdateringsstrategi från Broydens enparameterfamilj. I den tredje artikeln (Paper C), ger vi tillräckliga och nödvändiga villkor på en kvasi-Newtonmetod så att nämnda ekvivalens med konjugerade gradientmetoden erhålls. Mängden kvasi-Newtonstrategier som uppfyller dessa villkor är strikt större än Broydens enparameterfamilj. Vi visar också att denna mängd kvasi-Newtonstrategier innehåller ett oändligt antal uppdateringsstrategier där uppdateringsmatrisen är en symmetrisk matris av rang ett. I den andra artikeln (Paper B), används ett ramverk för icke-normaliserade Krylovunderrumsmetoder för att lösa symmetriska linjära ekvationssystem. Dessa ekvationssystem kan sakna lösning och matrisen kan vara indefinit/singulär. Denna typ av symmetriska linjära ekvationssystem uppkommer i en mängd formuleringar och metoder för optimeringsproblem med bivillkor. I fallet då det symmetriska linjära ekvationssystemet saknar lösning ger vi ett certifikat för detta baserat på en projektion på nollrummet för den symmetriska matrisen och karaktäriserar en minimum-residuallösning. Vi härleder även en minimum-residualmetod i detta ramverk samt ger explicita rekursionsformler för denna metod. I fallet då det symmetriska linjära ekvationssystemet saknar lösning så karaktäriserar vi en minimum-residuallösning av minsta euklidiska norm. Nyckelord: symmetriska linjära ekvationssystem, konjugerade gradientmetoden, kvasiNewtonmetoder, optimering utan bivillkor, kvadratisk optimering utan bivillkor, Krylovunderrumsmetoder, icke-normaliserade Lanczosvektorer, minimum-residualmetoden viii Acknowledgments First of all I would like to express my deepest gratitude to my academic advisor and coauthor Professor Anders Forsgren. Your mathematical expertise, attention to detail and patience have made it a privilege to work with you these last five years. It has been an inspiration to work with a researcher who sees all of the possibilities and is fearless in the face of a sometimes conservative community. Secondly, I would like to thank my academic co-advisor Professor Krister Svanberg for supporting me and for giving me valuable pointers when writing this thesis. My officemate Hildur Æsa Oddsdóttir, I have really appreciated sharing a room with you. You are an expert, not only in your academic field, but also when it comes to the Stockholm housing situation, were to buy the best electronic equipment and so much more. We have had so many great discussions on everything from feminism to LaTeX. Thank you for being a great co-worker and friend. I would also like to thank the faculty members and postdocs at the Division of Optimization and Systems Theory at KTH and my fellow students, past and present, in the Applied and Computational Mathematics doctoral program. I gratefully acknowledge that the research described in this thesis has been partly financially supported by the Swedish Research Council (VR). Thank you to my family, my partner Jonas and my friends for your love and support. Stockholm, May 2015 Tove Odland ix Table of Contents Abstract vii Acknowledgments ix Table of Contents xi I 1 1 2 3 4 6 Introduction I.1 Introduction to optimization . . . . . . . . . . . . . . . . . . . . . . . . I.2 Unconstrained optimization . . . . . . . . . . . . . . . . . . . . . . . . . I.2.1 Linesearch methods for unconstrained optimization . . . . . . . . I.2.2 Quasi-Newton methods . . . . . . . . . . . . . . . . . . . . . . . I.2.3 The method of conjugate gradients . . . . . . . . . . . . . . . . . I.2.4 Connection between the method of conjugate gradients and quasiNewton methods on quadratic problems . . . . . . . . . . . . . . I.3 Symmetric systems of linear equations arising in constrained optimization I.3.1 Solving symmetric systems of linear equations . . . . . . . . . . I.3.2 Methods based on the symmetric Lanczos process . . . . . . . . I.4 Summary of papers . . . . . . . . . . . . . . . . . . . . . . . . . . . . . I.4.1 Contribution by co-author . . . . . . . . . . . . . . . . . . . . . I.5 References . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 8 8 10 12 13 15 15 A Paper A A.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . A.2 Background . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . A.2.1 Conjugate gradient method . . . . . . . . . . . . . . . . . . . . . A.2.2 Quasi-Newton methods . . . . . . . . . . . . . . . . . . . . . . . A.2.3 Background results . . . . . . . . . . . . . . . . . . . . . . . . . A.3 Results . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . A.3.1 Update matrices defined by Proposition A.3.3 . . . . . . . . . . . A.3.2 Relation between δk and φk . . . . . . . . . . . . . . . . . . . . A.3.3 Remarks on when the one-parameter Broyden family is well-defined A.4 Conclusions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . A.5 References . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 19 22 22 23 24 26 26 29 31 33 33 35 xi B Paper B B.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . B.1.1 Notation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . B.2 Background . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . B.2.1 An extended representation of the unnormalized Lanczos vectors . B.3 Properties of the unnormalized Krylov framework . . . . . . . . . . . . . B.3.1 Convergence in the unnormalized Krylov framework . . . . . . . B.3.2 An unnormalized Krylov algorithm . . . . . . . . . . . . . . . . B.3.3 On the case when normalization is well defined . . . . . . . . . . B.4 Connection to the minimum-residual method . . . . . . . . . . . . . . . B.4.1 Convergence of the minimum-residual method . . . . . . . . . . B.4.2 A minimum-residual algorithm based on the unnormalized Krylov method . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . B.5 Summary and conclusion . . . . . . . . . . . . . . . . . . . . . . . . . . B.6 Appendix . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . B.6.1 A result on the Lanczos vectors . . . . . . . . . . . . . . . . . . B.6.2 Properties of the sequence {γk } . . . . . . . . . . . . . . . . . . B.6.3 The method of conjugate gradients . . . . . . . . . . . . . . . . . B.7 References . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . C Paper C C.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . C.2 Background . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . C.2.1 Previous results on a sufficient conditions for equivalence for CG and QN . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . C.3 Results on the equivalence of CG and QN on (QP) . . . . . . . . . . . . . C.3.1 Necessary and sufficient conditions for equivalence of CG and QN C.3.2 Symmetric rank-one update schemes . . . . . . . . . . . . . . . . C.3.3 Review of known results for the one-parameter Broyden family . C.4 Conclusion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . C.5 Appendix . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . C.6 References . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . xii 37 39 41 41 44 45 46 46 48 49 49 53 54 55 55 56 58 59 63 65 66 67 69 69 73 76 78 78 80 Introduction I.1 Introduction to optimization Optimization is used in a great number of applications such as telecommunications, intensity-modulated radiation therapy treatment planning, crew-scheduling, modeling of metabolic reaction networks, structure/topology design, portfolio theory and traffic modeling to mention just a few. Depending on the application the difficulty may lay in either the modeling and/or solving the posed optimization problem. An optimization problem is a decision problem on the form min f (x), x∈F (I.1) where x ∈ Rn is the vector of decision variables, f , the objective function, is some function from Rn → R and F is the feasible region for x. The field of optimization can roughly be divided into two major parts. The first part is posing the problem (I.1), that is finding a suitable f and F to represent the problem at hand. This requires consideration since the choices made highly affects how hard the problem will be to solve. We will consider formulations where f is some smooth function over F, which means that f will be continuously differentiable as many times as required. The second part is to solve the posed problem (I.1) and determining the quality of the solution obtained. The choices made in the modeling will determine which procedures for solving (I.1) that will be adequate. In this thesis we focus on the mathematics behind the second part and consider in great detail methods for solving (I.1) and the symmetric systems of linear equations that arise in various formulations of optimization problems and methods for solving these. When solving a problem (I.1), the goal is to find a vector x� that is a global minimizer, i.e. x� ∈ F and f (x� ) ≤ f (x) for all x ∈ F. In some cases, a globally optimal solution can not be attained and then a locally optimal solution is sought, such a solution x� satisfy f (x� ) ≤ f (x) in some neighborhood of x� . In order to find x� one can use some method that systematically searches F without evaluating all feasible points. Depending on the characteristics of f and F there are different methods that are suitable. In general these methods work the following way: given a starting guess x0 , the method generates a sequence of iterates x1 , x2 , . . . , where each new iterate is in some way a better estimation than the previous of x� . Note that in some methods the intermediate iterates may not be feasible to F. This process continues until 1 2 O N M ETHODS FOR S YMMETRIC L INEAR E QUATIONS ARISING IN O PTIMIZATION one reaches an iterate that satisfies some stopping criteria based on feasibility and some conditions that guarantee local or global optimality. All the work in this thesis is done assuming exact arithmetic. There are many aspects that needs to be considered when implementing an algorithm and solving a problem in finite arithmetic that are beyond the scope of the work presented in this thesis. However, understanding the mechanics of the inner-most parts of optimization methods is crucial and we have devoted this thesis to the mathematics behind this. Notation The letters i, j and k are used to denote integer indices, other lowercase letters such as g, x, p and c denote column vectors in Rn , possibly with super- and/or subscripts. For a symmetric matrix H ∈ Rn×n , H � 0 denotes that H is positive definite, i.e. that v T Hv > 0 for all 0 �= v ∈ Rn . Analogously, H � 0 is used to denote that H is positive semidefinite, i.e. that v T Hv ≥ 0 for all v ∈ Rn . The null space and range space of H are denoted by N (H) and R(H) respectively. I.2 Unconstrained optimization Consider an unconstrained optimization problem, i.e. F = Rn in (I.1), given by min f (x), x∈Rn (P) where the function f : Rn → R is smooth. For such a general unconstrained problem it holds that if, i) ∇f (x� ) = 0, ii) ∇2 f (x� ) � 0, then x� is locally optimal solution. The above conditions are known as second-order sufficient optimality conditions in that they are sufficient for local optimality, see, e.g., [26, Chapter 11] and [1, Chapter 4]. An optimization problem is convex if the feasible region is convex and the objective function f is a convex function over the feasible region. If an optimization problem is convex, then all locally optimal solutions are globally optimal solutions and therefore this is a desirable property to have for a problem one wants to solve. See, e.g., [26, Chapter 11] for an introduction on the topic of convexity. It holds that an unconstrained problem (P) is convex if and only if ∇2 f (x) � 0 for all x ∈ Rn . Note that if f is a convex function on Rn then x� is a globally optimal solution if and only if ∇f (x� ) = 0, see. e.g. [1, Chapter 4]. For this reason, many methods for solving (P) will search for an iterate where the gradient of the objective function is zero, a so-called stationary point to the problem. In the next section a certain class of such methods known as linesearch methods is introduced. I NTRODUCTION 3 I.2.1 Linesearch methods for unconstrained optimization For our purposes, methods for unconstrained optimization problems (P) are divided into two classes, methods where second-derivative information is available and used and methods which only use first-derivative information. It deserves to be mentioned that there are also so-called derivative-free methods that uses neither first- nor second-derivative information, see, e.g., [9] for an introduction to such methods. We consider so-called linesearch methods that takes us from xk to the next iterate xk+1 by taking some step θk along a non-zero search direction pk , (I.2) xk+1 = xk + θk pk . Different linesearch methods arise depending on how the search direction pk is generated at each step k. A general requirement on a linesearch method is that it should have the descent property, i.e. f (xk + θk pk ) < f (xk ) for all k. This property is for example guaranteed if ∇f (xk )T pk < 0, i.e. pk is a descent direction of f at xk . Then there always exists a positive θk such that f (xk + θk pk ) < f (xk ) is satisfied, see, e.g. [22, Chapter 3]. There are many strategies on how to choose the steplength θk . Ideally, one would like to find the optimal value of θk corresponding to θk = argmin f (xk + θpk ), θ∈R given xk and pk . Using such a value of θk is referred to the linesearch method employing exact linesearch. However, it is not always worthwhile to use exact linesearch, especially for iterates that are far from the optimal solution. See, e.g., [41, Chapter 3] for a discussion on other strategies for choosing θk and their benefits. A common stopping criteria for linesearch methods is ∇f (x) = 0, i.e. the current iterate is a stationary point. If a linesearch method is implemented and used in finite precision then this stopping criteria will be �∇f (x)� < �, where � should be chosen as some small number in an appropriate manner. In addition it should be pointed out that for finite precision computations there are more stable choices for stopping criteria and one can also consider using a combination of such choices, see, e.g. [1, Chapter 11] and [22, Chapter 8]. When second-derivative information is available the perhaps most famous line-search method is Newton’s method, for which the search direction pk is given by solving ∇2 f (xk )pk = −∇f (xk ). (I.3) Newton’s method makes a quadratic model of f at xk . If ∇ f (xk ) � 0 then the unit step along pk takes us to the minimum of this quadratic model. If ∇2 f (xk ) is not positive definite then the quadratic model may not have a minimum and in this case there are strategies in which one uses a related positive-definite matrix instead of ∇2 f (xk ) in (I.3), see, e.g., [22, Chapter 4]. If Newton’s method is started sufficiently close to x� , then the convergence rate of the method is at least quadratic, see, e.g. [37, Chapter 7]. Note that if f is a strictly convex quadratic function then Newton’s methods finds the optimal solution in one iteration. 2 4 O N M ETHODS FOR S YMMETRIC L INEAR E QUATIONS ARISING IN O PTIMIZATION In this thesis we are concerned with methods for problems on the form (P) where second-derivative information is not available, but for which the gradient, ∇f (x), is available for all x ∈ Rn . One of the oldest and simplest methods of this kind is the method of steepest descent, also known as gradient descent, see, e.g., [45] for references to the origin of this method. In the method of steepest descent, the search direction pk is chosen in the negative direction of the gradient at xk , i.e. pk = −∇f (xk ). (I.4) Even though the search direction pk will always be a descent direction for ∇f (xk ) �= 0, (∇f (xk )T pk = −∇f (xk )T ∇f (xk ) < 0) this method may have a very poor convergence rate for some starting guesses and displaying a behavior known as ’zig-zagging’. Further, with the norm of the gradient as a stopping criteria one might run into problems as this norm tends to oscillate. For a discussion on the magnitude of the oscillations of the gradient norm in the method of steepest descent see, e.g., [40]. Next a class of first-order methods that is studied in two of the appended papers (Paper A and Paper C), namely quasi-Newton methods, is introduced. I.2.2 Quasi-Newton methods The first quasi-Newton method was suggested by Davidon in 1959 [11], using the term variable metric method. Davidon’s work caught the attention of the optimization community when Fletcher and Powell [17] presented a slightly modified method that outdid all other existing methods of its time. The paper by Fletcher and Powell marked the start of a period of intense research of quasi-Newton methods. For a thorough treatment of quasiNewton methods see, e.g., [41, Chapter 8], [16, Chapter 3] and [37, Chapter 9] and the references therein. An interesting note is that Davidon’s original paper was not published until 1991. Quasi-Newton methods comprise a class of methods that mimic Newton’s method by using quasi-Newton matrix Bk that is an approximation of the true Hessian ∇2 f (xk ) in (I.3). In effect, in a quasi-Newton method the search direction at step k is generated by solving Bk pk = −∇f (xk ), (I.5) usually with B0 = I, i.e. p0 is the steepest descent step. Different quasi-Newton methods arise depending on the choice of Bk in each iteration. The matrix Bk is usually given by adding an update matrix Uk to the previous approximation Bk−1 , i.e. Bk = Bk−1 + Uk . (I.6) One often considers updating a Cholesky factorization of Bk , then (I.5) can be solved in order of n2 operation. Further if the update matrix is of low rank, then one does not need to recompute the Cholesky factorization at each step, see, e.g., [22, Chapter 4]. In many paper, especially the earlier ones, an approximation of the inverse Hessian, Mk , is used and (I.5) then takes the form pk = −Mk gk . From a theoretical point of view I NTRODUCTION 5 it makes little difference if one works with Bk or Mk . In the work presented in this thesis we choose to work with quasi-Newton methods on the form given by (I.5). A famous class update schemes for quasi-Newton methods is the one-parameter Broyden family [3]. This family of update schemes depends on φk , a free parameter, known as the Broyden parameter. The one-parameter Broyden family includes the DFP (DavidonFletcher-Powell) update scheme, for φk = 1, which was the first quasi-Newton scheme to be suggested in [11, 17], the BFGS update scheme, for φk = 0, discovered independently by Broyden [4], Fletcher [15], Goldfarb [23] and Shanno [47] and the symmetric rank-one update scheme, SR1. All update schemes in the one-parameter Broyden family are of rank two except SR1 which is of rank one. There is a subclass of the one-parameter Broyden family called the convex Broyden family which is given by all convex combinations of the BFGS and DFP update schemes. Other parametric formulations of quasi-Newton update schemes are possible, see, e.g., [47, 49]. The update schemes in the one-parameter Broyden family can become not well defined e.g. for certain values (degenerate values) of the Broyden parameter φ, see, e.g., [19]. In many treatments the matrix Bk is assumed to be symmetric, for example the oneT parameter Broyden family has the property of hereditary symmetry, i.e. if Bk−1 = Bk−1 T then Bk = Bk . It is also possible to consider nonsymmetric update schemes, e.g. the Huang family of updates which has three degrees of freedom, see [34]. The symmetric part of the Huang family of updates can be shown to be equivalent to a class of updates schemes with two degrees of freedom defined by the self-scaling quasi-Newton algorithm by Oren and Luenberger [42]. In this thesis quasi-Newton methods for which Bk = BkT are the main focus. An interesting characteristic of the SR1 update scheme is that on a quadratic problem given any set of linearly independent search directions, the optimal solution will be found in at most n steps without the use of linear searches in each step, see, e.g. [10, 15]. In [52], the class of methods presented in [34] is modified to obtain this characteristic. For very large problems it might be too expensive to compute and store the full n × n matrix Bk at each step. In this case one can consider a limited-memory quasi-Newton method that uses an implicit representation of Bk , see, e.g., [41, Chapter 9]. Quasi-Newton methods can be shown to have Q-superlinear convergence, see, [13]. Note that, for this convergence, it is not required the {Bk } tends to the true Hessian at x� , only that it does so along the search directions {pk }. In 1972, Dixon [14] showed that all well-defined update schemes in the one-parameter Broyden family generate parallel search directions on (P) for any smooth function if perfect linesearch is used. Perfect linesearch is a generalization of exact linesearch for smooth nonconvex functions. For a strictly convex quadratic functions these search directions will in addition be mutually conjugate with respect to the Hessian and hence linearly independent. Hence, a quasi-Newton method using any well-defined update scheme in the one-parameter Broyden family has the property of finite termination on the corresponding quadratic optimization problem, see, e.g., [34]. Further, on a quadratic problem the search directions generated by a quasi-Newton method, using any well-defined update scheme from the one-parameter Broyden family, 6 O N M ETHODS FOR S YMMETRIC L INEAR E QUATIONS ARISING IN O PTIMIZATION will be parallel with those generated by the method of conjugate gradients(introduced in the next section) provided that both methods employ exact linesearch. I.2.3 The method of conjugate gradients Consider an unconstrained quadratic optimization problem 1 min q(x) = minn xT Hx + cT x, x∈R 2 x∈Rn (QP) where H = H T � 0. Note that solving (QP) is equivalent to solving a symmetric system of linear equations on the form Hx + c = 0, (I.7) where H = H T � 0. For (QP), the second-order sufficient optimality conditions are given by i) Hx� + c = 0, ii) H � 0. If H � 0 then there exists a unique solution to (I.7) or equivalently a unique optimal solution to (QP) given by x� = −H −1 c. For the case H � 0 there exist an infinite number of solutions to (I.7) if c ∈ R(H) and no solutions if c ∈ / R(H). In this section the method of conjugate gradients by Hestenes and Stiefel [33] is considered. This is a method for solving a system of symmetric linear equations on the form (I.7) and which can also be expressed as a linesearch method for (QP). The method of conjugate gradients is a first-order method and hence the matrix H is only needed to be available through matrix-vector multiplication. In a linesearch description of the method of conjugate gradients for solving (QP) the gradient of the objective function at xk is denoted by gk , i.e. gk = Hxk + c. In accordance with (I.2) the gradient is updated in each step as gk+1 = gk + θk Hpk . (I.8) On (QP) it is natural use exact linesearch for θk . Consider θ to be unknown and solve for the optimal θ, 1 min q(xk + θpk ) = min θ2 pTk Hpk + θpTk gk + q(xk ). θ∈R θ∈R 2 Note that q(xk ) is a constant term, hence, 0= � d� q(xk + θpk ) = θpTk Hpk + pTk gk . dθ It follows that the exact linesearch step is given by θk = − pTk gk . pTk Hpk (I.9) I NTRODUCTION 7 As Hpk is already needed for computation of gk+1 in (I.8), calculating the exact linesearch step only entails two vector-vector multiplications. In a linesearch description of the method of conjugate gradients the search directions are given by pCG = −g0 , 0 pCG = −gk + k gkT gk pCG k−1 , T g gk−1 k−1 (I.10) k ≥ 1. For the methods of conjugate gradients it holds that giT gj = 0 for all i �= j and also that gkT pCG = 0, for all j = 0, . . . , k − 1. In fact the gradients {gk } in the method of j conjugate gradients form an orthogonal basis for the Krylov subspaces generated by H and c, K0 (c, H) = {0}, Kk (c, H) = span{c, Hc, H 2 c, . . . , H k−1 c}, k ≥ 1, (I.11) i.e., it holds that gk ∈ Kk+1 (c, H) ∩ Kk (c, H)⊥ . Let r be the minimum index k for which Kk+1 (c, H) = Kk (c, H), then after generating r mutually orthogonal, and hence linearly independent, gradients the stopping criteria is reached as gr = 0. Hence, on (QP) the method of conjugate gradients has the property of finite termination. Further, the search directions given by (I.10) form an H-orthogonal basis for the expanding Krylov subspaces given by (I.11), i.e. the search directions are mutually conjugate with respect to H. See, e.g., [30] for an introduction to Krylov subspaces and Krylov subspace methods for solving systems of linear equations. There are several other ways to derive the method of conjugate gradients. For example, as a conjugate directions method where the Gram-Schmidt conjugation is based on the orthogonal gradients, as minimization of q(x) where x belongs expanding Krylov subspaces, as minimization of �xk −x� �H or equivalently �Hxk +c�H −1 , and as a method for solving a symmetric system of linear equations based on a normalized symmetric Lanczos process. For a variation of treatments of the method of conjugate gradients see, e.g. [12, Chapter 6], [25, Chapter 10], [41, Chapter 5] and [50] and the references therein. The method of conjugate gradients was first extended to a nonlinear version for solving general unconstrained problems (P) by Fletcher and Reeves [18]. For a survey on the Fletcher-Reeves and other nonlinear versions of the methods of conjugate gradients see, e.g., [32]. In [48], the topic of the method of conjugate gradients with inexact line searches is considered. Next the connection between the method of conjugate gradients and quasi-Newton methods on (QP) or equivalently for solving a symmetric system of linear equations with a positive definite matrix is discussed. This is one of the central themes of the thesis. Besides being important in its own right, quadratic optimization problems are often at the core of algorithms for more challenging non-linear optimizations problems. Further, many problems tend to behave quadratically in the region of an optimal solution. In addition, some methods, for example the method of conjugate gradients and BFGS-methods, 8 O N M ETHODS FOR S YMMETRIC L INEAR E QUATIONS ARISING IN O PTIMIZATION are considered effective when one wants to obtain an approximate solution to a quadratic problem with desirable characteristics in a timely manner, see, e.g. [6] for an article on intensity-modulated radiation therapy treatment planning where this property is used. I.2.4 Connection between the method of conjugate gradients and quasi-Newton methods on quadratic problems Even though the method of conjugate gradients and quasi-Newton methods may seem to be very different kind of methods, it is well-known that, when solving (QP) or equivalently when solving (I.7) with H = H T � 0, the method of conjugate gradients and a quasiNewton methods using any well-defined update from the one-parameter Broyden family generate parallel search directions provided that exact linesearch is used in both methods, see, e.g., [16, Chapter 3]. This implies that the corresponding x-iterates will be identical in exact arithmetic. In [38], this result is shown for the DFP-update scheme and in [39], Nazareth shows that for the BFGS update scheme the search directions are, in addition to being parallel, of the same length. The connection between the method of conjugate gradients and the BFGS update scheme for a quasi-Newton method is further explored by Buckley [5] where one of the results is that the BFGS update scheme can be carried out using any of the previous approximation matrices and equivalence with the method of conjugate gradients still holds on (QP). The search direction, pk , in a quasi-Newton method, given by solving (I.5), is determined by the choice of Bk in each iteration k. One can consider what the precise conditions are on Bk and/or the update matrix Uk such that pk and pCG are parallel, i.e. pk = δk pCG k k for some non-zero scalar δk . In two of the appended papers in this thesis (Paper A and Paper C), this topic is investigated. Next symmetric systems of linear equations that arise in problem formulations and methods for constrained optimization problems are considered. I.3 Symmetric systems of linear equations arising in constrained optimization Consider a constrained optimization problem on the form minx∈Rn subject to f (x) h(x) = 0, g(x) ≥ 0, (I.12) where f , h and g are smooth functions. Depending on further characteristics of the functions f , h and g there are different methods that are suitable when one wants to solve (I.12). Some methods in which symmetric systems of linear equations arise are mentioned next. Consider a problem on the form of (I.12) with a quadratic objective function and linear equality constraints, 1 T T minx∈Rn 2 x Hx + c x, (I.13) subject to Ax = b, I NTRODUCTION 9 where A ∈ Rm×n is of rank m. It holds that, if x� is a local minimizer of (I.13), then there exists λ� ∈ Rm such that � �� � � � x� −c H AT = . (I.14) � −λ b A 0 The vector λ� is the Lagrange multipliers corresponding to the iterate x� . The symmetric system of linear equations with a symmetric and indefinite matrix given in (I.14) are the first-order necessary optimality conditions for (I.13). Let Z be a matrix whose columns form a basis for the null space of A. It holds that x� is a global minimizer to (I.13) if and only if there exists λ� ∈ Rm such that (I.14) holds and the reduced Hessian Z T HZ is positive semi-definite. Solving a symmetric system of linear equations on the form (I.14) is important, not only for finding an optimal solution to (I.13) itself, but also in methods for solving more complex problems. An example of this is methods that employ sequential quadratic programming (SQP). An SQP-method is a method for solving (I.12) where in each step a search direction is found by solving a quadratic subproblem. For simplicity consider (I.12) with only equality constraints, then the subproblem will be on the form (I.13) where the objective function is a quadratic model of the Lagrangian function L(x, λ) = f (x) − λT h(x) at the current iterates (xk , λk ) and the constraints are given by a linearization of h around the current iterate xk . The first-order necessary optimality conditions of such a quadratic subproblem is given by (I.14), a symmetric system of linear equations with a matrix that is symmetric and indefinite. Such an SQP-method is equivalent to applying Newton iterations to the first-order necessary optimality conditions of (I.12) itself. Note that there are SQP-methods for (I.12) with both inequality and equality constraints. The iterates of an SQP-method will in general not be feasible to the original problem. See, e.g., [26, Chapter 15] for an introduction to SQP. Consider a quadratic problem with linear inequality constraints, minx∈Rn subject to 1 T 2 x Hx Ax ≥ b. + cT x, (I.15) When it comes to handling inequality constraints we present two approches. The first one is to consider an active-set method, in which one keeps track of the set of active constraints and treat them as equality constraints while ignoring the inactive constraints. In each iteration one then gets a problem on the form (I.13) and obtains a search direction to a new iterate at which the set of active constraints is updated. This implies that in each iteration one solves a symmetric system of linear equations with an indefinite symmetric matrix. The iterates of an active-set method will always be feasible to the original problem. See, e.g., [41, Chapter 15] for an introduction to active-set methods. The second approach entails moving in the interior of the feasible region, and in a sense letting all constraints be inactive but eventually reaching an optimal solution at which some constraints are active. Such a method is called an interior-point method. The main idea of interior-point methods is to perturb the first-order necessary optimality conditions of the 10 O N M ETHODS FOR S YMMETRIC L INEAR E QUATIONS ARISING IN O PTIMIZATION given problem. The perturbed optimality conditions for (I.15) for some µ > 0 are given by Ax − s = b, Hx − AT λ = −c, SΛe = µe, where S and Λ are diagonal matrices with the elements of s and λ on their respective diagonals and e is a vector of ones. The inequalities s ≥ 0 and λ ≥ 0 are kept implicitly in the method. These perturbed optimality conditions are called the primal-dual nonlinear equations. An primal-dual interior-point method is based on applying Newton iterations to the primal-dual nonlinear equations. Such an iteration can be written on the form of a symmetric system of linear equations with a matrix that is symmetric and indefinite given by � �� � � � ∆x H AT Hx + c − AT λ =− , −1 −1 −∆λ A −SΛ Ax − b − µΛ e where ∆s has been eliminated. Interior-point methods are illustrated here for solving (I.15), however these methods can be applied directly to (I.12) where the inequality constraints can be handled by a barrier transformation. There are also interior-point methods for solving linear programming problems on standard form. Symmetric systems of linear equations with an indefinite symmetric matrix arise in all of the above mentioned cases. See, e.g., [20, 21, 41, 53] for introductions to interior-point methods. I.3.1 Solving symmetric systems of linear equations Symmetric systems of linear equations on the form Hx + c = 0, (I.16) with H = H T arise in many methods for constrained optimization as indicated in the previous section. In addition, finding a solution to (I.16), if such a solution exists, or in the incompatible case a minimum-residual solution, argminx∈IRn �Hx + c�22 , are both important problems in the field of numerical linear algebra. As the only assumption on (I.16) is that H = H T , it may be the case that H is indefinite/singular and the system of linear equations (I.16) may be incompatible. Note that H and c in (I.16) represents the whole matrix on the left-hand side and the whole right-hand side of a symmetric system of linear equations, not to be confused with the H and c that appears as a components in e.g. the matrix and right-hand side of (I.14). In this thesis the mathematics behind so-called Krylov subspace methods for solving (I.16) is considered. These methods are based on generating an orthogonal basis for the expanding Krylov subspaces generated by H and c given in (I.11). In particular, this orthogonal basis will be generated by a symmetric version of a process given by Lanczos in [35, 36] and this process is commonly referred to as a symmetric Lanczos process. There has been a significant amount of research devoted to the theory on using Lanczos I NTRODUCTION 11 processes to solve (I.16) both for H symmetric and H non symmetric, see, e.g., [24,29,43]. The discussion here on solving (I.16) is aimed towards the treatment of this problem in the appended papers, it should be pointed out that solving (I.16) can also be done by applying different matrix-factorization techniques to H, see, e.g., [25, Chapter 5]. The symmetric Lanczos process is given as follows. Choose q0 = c/�c�2 such that �q0 �2 = 1, and generate {qk }rk=1 as −β0 q1 = −Hq0 + α0 q0 −βk qk+1 = −Hqk + αk qk + βk−1 qk−1 , k = 1, . . . , r − 1, (I.17) where the coefficients α and β are chosen such that one gets orthonormal vectors {qk }r−1 k=1 . These choices are given by, α0 = and β0 = q0T Hq0 , q0T q0 q0T Hq1 , q0T q0 αk = βk = qkT Hqk , qkT qk (I.18) k = 1, . . . r − 1, T qk−1 Hqk , T q qk−1 k−1 (I.19) k = 1, . . . r − 1. ⊥ Then {qk }r−1 k=0 are such that qk ∈ Kk+1 (c, H) ∩ Kk (c, H) and �qk �2 = 1 for all k = 0, . . . , r−1 and qr = 0. The three-term recurrence in (I.17) is due to the fact that H = H T , see, e.g. [46, Chapter 6]. Note that if one only requires that the generated vectors should be orthogonal then there is one degree of freedom in each iteration and (I.17) may be expressed as q1 = θ0 (−Hq0 + α0 q0 ) qk+1 = θk (−Hqk + αk qk + βk−1 qk−1 ), (I.20) k = 1, . . . , r − 1, where the scalar θk represents the available scaling factor in each iteration k and α and β are as in (I.18) and (I.19) respectively. These vectors {qk }rk=0 given by (I.20) are referred to as Lanczos vectors. Let Qk be the matrix with the Lanczos vectors q0 , q1 , . . . , qk as columns vectors, then (I.20) may be written on matrix form as, HQk = Qk+1 T¯k = Qk Tk − 1 qk+1 eTk+1 , θk where  α0   − θ1 0 Tk =    β0 .. . .. . .. .. . . 1 − θk−1  βk−1 αk   ,   T¯k = � Tk − θ1k eTk+1 � . (I.21) 12 O N M ETHODS FOR S YMMETRIC L INEAR E QUATIONS ARISING IN O PTIMIZATION For the choice of θk = − β1k one obtains the recursion (I.17) and it holds that ||qk ||2 = ||q0 ||2 for all k. Further, for this choice of θk , Tk will be a symmetric matrix. Changing the set of {θk }r−1 k=0 can be seen as a similarity transform of Tk , see, e.g., [29]. Many methods for solving (I.16) use matrix-factorization techniques on Tk or T¯k given in (I.21). For an introduction to how Krylov subspace methods are formalized in this way, see, e.g., [8, 43]. In the work presented in this thesis we choose to work with the recursion in (I.20) directly. As the set of vectors {qk }r−1 k=0 form an orthogonal basis for the Krylov subspaces in (I.11), each vector may be expressed as qk = k � (j) γk H j c = H j=0 k �� j=1 � (j) (0) γk H j−1 c + γk c, k = 0, . . . , r, (j) for some scalars {γk }kj=0 . The following choice θ0 = 1 , α0 θk = 1 , αk + βk−1 (I.22) k = 1, . . . , r − 1, (0) of θk in (I.20) will give a vector qk+1 on the form of a gradient gk+1 , i.e. with γk = 1. (0) Requiring γk = 1 is usually referred to as the normalization/consistency condition. We will denote by a normalized symmetric Lanczos process the symmetric Lanczos process where the particular choice of scaling (I.22) is used and the generated vector are made to satisfy the normalization condition. Note that (I.22) may not always be a well defined expression, so a normalized symmetric Lanczos process may suffer a so-called pivot breakdown. The idea of utilizing an unnormalized symmetric Lanczos process in order to avoid such breakdowns is due to Gutknecht [27, 28].1 In an unnormalized setting (0) there is no restriction on the sign or value of γk . One can show that this quantity can not be zero in two consecutive iterations. This property is used in for example composite step biconjugate gradient methods and other methods employing look-ahead techniques to avoid breakdowns, see, e.g., [2, 7, 31]. I.3.2 Methods based on the symmetric Lanczos process A symmetric Lanczos process (normalized or not) is by itself not sufficient to obtain a solution to (I.16). Hence, a method for solving (I.16) needs to manage additional quantities such that a solution can be extracted when the stopping criteria is reached. We shall mention a few such Krylov subspace methods, but for a thorough treatment of such methods, both for the symmetric and non symmetric case, see, e.g., [46, Chapters 5, 6, 7] and [29]. The method of conjugate gradients may be viewed as a method for solving (I.16) based on the normalized symmetric Lanczos process. In the normalized case search directions 1 Gutknecht considers the non symmetric case and uses the unnormalized framework to avoid pivot breakdowns. When solving (I.16) for a nonsymmetric H there are several possible breakdowns for methods based on a nonsymmetric Lanczos process, see, e.g., [29, 46]. For the symmetric case the pivot breakdown is the only one that can happen. I NTRODUCTION 13 may be defined and the three-term recurrence in (I.20) may be reduced to a two-term recurrence using these search directions. Note that when viewing the method of conjugate gradients as a method based on a normalized symmetric Lanczos process it becomes clear that the search directions are not necessary for the formulation but rather may be viewed as a product of the normalization itself. The breakdown of the method of conjugate gradients when solving (I.16) corresponds exactly to the case when normalization is not possible. For the case H = H T � 0 no breakdown of the method of conjugate gradients will occur, and further it holds that if H = H T � 0 and c ∈ R(H) normalization will always be possible. As mentioned it is also possible to consider a method based on the unnormalized symmetric Lanczos process in order to avoid pivot-breakdowns. Such an unnormalized Krylov subspace method is described in Paper B and in terms of the survey [29] this method is called an inconsistent ORes version of the method of conjugate gradients. Further in the case when (I.16) is incompatible one can consider a minimum-residual method that delivers a minimum-residual solution to (I.16). Minimum residual methods was first suggested in Lanczos’ early paper [36] and by Stiefel in [51]. The name, minimum-residual method, is adopted from the implementation of the method, MINRES, by Paige and Saunders, see [44]. In this implementation, which is based on the symmetric Lanczos process, the optimal solution of minx∈IRn �Hx + c�22 obtained is in general not of minimum Euclidean norm. In [8], Choi, Paige and Saunders present a MINRESlike algorithm, MINRES-QLP, that does deliver a minimum-residual solution of minimum Euclidean norm. In the second paper of this thesis (Paper B) the mathematical properties of methods based on an unnormalized symmetric Lanczos process are considered and both a minimumresidual solution and a minimum-residual solution of minimum Euclidean norm in the case when (I.16) is incompatible are characterized. I.4 Summary of papers In this section a short summary of each of the appended papers is given, highlighting the main contributions. The common theme of the three papers is the in-depth study of the mathematical properties of methods for solving symmetric systems of linear equations arising in various optimization problem formulations and derivation of new theoretical results. Paper A A. Forsgren and T. Odland. On the connection between the conjugate gradient method and quasi-Newton methods on quadratic problems. Computational Optimization and Applications, 60 (2015), pp 377-392. We derive conditions on the update matrix Uk , and hence on the quasi-Newton matrix Bk such that the search directions generated by the corresponding quasi-Newton method and the method of conjugate gradients are parallel when applied to an unconstrained strictly 14 O N M ETHODS FOR S YMMETRIC L INEAR E QUATIONS ARISING IN O PTIMIZATION convex quadratic optimization problem or equivalently a symmetric system of linear equations with a symmetric positive definite matrix. The conditions state that the range of the update matrix must lie in the last two dimensions of the Krylov subspace generated by H and c and that in addition the so-called quasi-Newton condition must be satisfied. These conditions on Uk , that are based on a sufficient condition to obtain mutually conjugate search directions, are shown to be equivalent to the quasi-Newton scheme belonging to the one-parameter Broyden family. An explicit formula for the one-to-one correspondence between the Broyden parameter and the scaling of the search direction from the corresponding quasi-Newton method compared to the search direction in the method of conjugate gradients is stated. Paper B A. Forsgren and T. Odland. On solving symmetric systems of linear equations in an unnormalized Krylov subspace framework. arXiv:1409.4937 [math.OC], 2014. We study the mathematical properties of an unnormalized Krylov subspace framework for solving a symmetric system of linear equations Hx + c = 0. This framework is based on an unnormalized symmetric Lanczos process for generating an orthogonal basis for the expanding Krylov subspaces generated by H and c. Using only the three-term recurrence relations for a triple of quantities associated with each orthogonal vector from the unnormalized symmetric Lanczos process we give conditions on whether the symmetric system of linear equations is compatible or not. In the compatible case a solution is obtained and in the incompatible case we obtain a certificate of incompatibility based on a projection on the null space of H and characterize a minimum-residual solution to the symmetric system of linear equations. Further we derive a minimum-residual method in the unnormalized Krylov subspace framework and give explicit recursions of the minimum residual iterates. In the case of an incompatible symmetric system of linear equations we characterize a minimum-residual solution of minimum Euclidean norm. Paper C A. Forsgren and T. Odland. On the equivalence of the method of conjugate gradients and quasi-Newton methods on quadratic problems. arXiv:1503.01892 [math.OC], 2015. Submitted to Computational Optimization and Applications. We derive necessary and sufficient conditions on the quasi-Newton matrix Bk and the update matrix Uk such that the search direction in the corresponding quasi-Newton method is parallel to the search direction in the method of conjugate gradients on an unconstrained strictly convex quadratic problem or equivalently a symmetric system of linear equations with a symmetric positive definite matrix. We show that the set of quasi-Newton schemes admitted by these necessary and sufficient conditions are strictly larger than the one-parameter Broyden family. Further we show that this set of quasi-Newton schemes contains an infinite number of symmetric rank-one I NTRODUCTION 15 update schemes. The necessary and sufficient conditions derived in this paper admit nonsymmetric quasi-Newton schemes as well, however it is shown that symmetry of Bk can be achieved without effort. I.4.1 Contribution by co-author Throughout the work Anders Forsgren has acted as an advisor in the classical sense and is therefore co-author of all three papers. I.5 References [1] N. Andréasson, A. Evgrafov, and M. Patriksson. An Introduction to Continuous Optimization. Professional Publishing Svc., 2007. [2] R. E. Bank and T. F. Chan. An analysis of the composite step biconjugate gradient method. Numer. Math., 66(3):295–319, 1993. [3] C. G. Broyden. Quasi-Newton methods and their application to function minimisation. Math. Comp., 21:368–381, 1967. [4] C. G. Broyden. The convergence of a class of double-rank minimization algorithms. II. The new algorithm. J. Inst. Math. Appl., 6:222–231, 1970. [5] A. Buckley. Termination and equivalence results for conjugate gradient algorithms. Math. Programming, 29(1):64–76, 1984. [6] F. Carlsson and A. Forsgren. Iterative regularization in intensity-modulated radiation therapy optimization. Medical physics, 33(1):225–234, 2006. [7] T. F. Chan and T. Szeto. Composite step product methods for solving nonsymmetric linear systems. SIAM J. Sci. Comput., 17(6):1491–1508, 1996. [8] S.-C. T. Choi, C. C. Paige, and M. A. Saunders. MINRES-QLP: a Krylov subspace method for indefinite or singular symmetric systems. SIAM J. Sci. Comput., 33(4):1810–1836, 2011. [9] A. R. Conn, K. Scheinberg, and L. N. Vicente. Introduction to derivative-free optimization, volume 8 of MPS/SIAM Series on Optimization. Society for Industrial and Applied Mathematics (SIAM), Philadelphia, PA; Mathematical Programming Society (MPS), Philadelphia, PA, 2009. [10] J. Cullum and R. K. Brayton. Some remarks on the symmetric rank-one update. J. Optim. Theory Appl., 29(4):493–519, 1979. [11] W. C. Davidon. Variable metric method for minimization. SIAM J. Optim., 1(1):1–17, 1991. 16 O N M ETHODS FOR S YMMETRIC L INEAR E QUATIONS ARISING IN O PTIMIZATION [12] J. W. Demmel. Applied numerical linear algebra. Society for Industrial and Applied Mathematics (SIAM), Philadelphia, PA, 1997. [13] J. E. Dennis, Jr. and J. J. Moré. A characterization of superlinear convergence and its application to quasi-Newton methods. Math. Comp., 28:549–560, 1974. [14] L. C. W. Dixon. Quasi-Newton algorithms generate identical points. Math. Programming, 2:383–387, 1972. [15] R. Fletcher. A new approach to variable metric algorithms. The Computer Journal, 13(3):317–322, Jan. 1970. [16] R. Fletcher. Practical methods of optimization. A Wiley-Interscience Publication. John Wiley & Sons Ltd., Chichester, second edition, 1987. [17] R. Fletcher and M. J. D. Powell. A rapidly convergent descent method for minimization. Comput. J., 6:163–168, 1963/1964. [18] R. Fletcher and C. M. Reeves. Function minimization by conjugate gradients. Comput. J., 7:149–154, 1964. [19] R. Fletcher and J. W. Sinclair. Degenerate values for Broyden methods. J. Optim. Theory Appl., 33(3):311–324, 1981. [20] A. Forsgren, P. E. Gill, and J. R. Shinnerl. Stability of symmetric ill-conditioned systems arising in interior methods for constrained optimization. SIAM J. Matrix Anal. Appl., 17:187–211, 1996. [21] A. Forsgren, P. E. Gill, and M. H. Wright. Interior methods for nonlinear optimization. SIAM Rev., 44(4):525–597 (2003), 2002. [22] P. E. Gill, W. Murray, and M. H. Wright. Practical optimization. Academic Press Inc. [Harcourt Brace Jovanovich Publishers], London, 1981. [23] D. Goldfarb. A family of variable-metric methods derived by variational means. Math. Comp., 24:23–26, 1970. [24] G. H. Golub and D. P. O’Leary. Some history of the conjugate gradient and Lanczos algorithms: 1948–1976. SIAM Rev., 31(1):50–102, 1989. [25] G. H. Golub and C. F. Van Loan. Matrix computations, volume 3 of Johns Hopkins Series in the Mathematical Sciences. Johns Hopkins University Press, Baltimore, MD, 1983. [26] I. Griva, S. G. Nash, and A. Sofer. Linear and nonlinear optimization. Society for Industrial and Applied Mathematics (SIAM), Philadelphia, PA, second edition, 2009. [27] M. H. Gutknecht. The unsymmetric Lanczos algorithms and their relations to Pade approximation, continued fractions, and the QD algorithm. Preliminary Proceedings of the Copper Mountain Conference on Iterative methods, April, 1990. I NTRODUCTION 17 [28] M. H. Gutknecht. A completed theory of the unsymmetric Lanczos process and related algorithms, part I. SIAM J. Matrix Anal. Appl., 13(2):594–639, 1992. [29] M. H. Gutknecht. Lanczos-type solvers for nonsymmetric linear systems of equations. In Acta numerica, 1997, volume 6 of Acta Numer., pages 271–397. Cambridge Univ. Press, Cambridge, 1997. [30] M. H. Gutknecht. A brief introduction to Krylov space methods for solving linear systems. In Frontiers of Computational Science, pages 53–62. Springer-Verlag Berlin, 2007. International Symposium on Frontiers of Computational Science, Nagoya Univ, Nagoya, Japan, Dec 12-13, 2005. [31] M. H. Gutknecht and K. J. Ressel. Look-ahead procedures for Lanczos-type product methods based on three-term Lanczos recurrences. SIAM J. Matrix Anal. Appl., 21(4):1051–1078, 2000. [32] W. W. Hager and H. Zhang. A survey of nonlinear conjugate gradient methods. Pac. J. Optim., 2(1):35–58, 2006. [33] M. R. Hestenes and E. Stiefel. Methods of conjugate gradients for solving linear systems. J. Research Nat. Bur. Standards, 49:409–436 (1953), 1952. [34] H. Y. Huang. Unified approach to quadratically convergent algorithms for function minimization. J. Optimization Theory Appl., 5:405–423, 1970. [35] C. Lanczos. An iteration method for the solution of the eigenvalue problem of linear differential and integral operators. J. Research Nat. Bur. Standards, 45:255–282, 1950. [36] C. Lanczos. Solution of systems of linear equations by minimized-iterations. J. Research Nat. Bur. Standards, 49:33–53, 1952. [37] D. G. Luenberger. Linear and nonlinear programming. Addison-Wesley Pub Co, Boston, MA, second edition, 1984. [38] G. E. Myers. Properties of the conjugate-gradient and Davidon methods. J. Optimization Theory Appl., 2:209–219, 1968. [39] L. Nazareth. A relationship between the BFGS and conjugate gradient algorithms and its implications for new algorithms. SIAM J. Numer. Anal., 16(5):794–800, 1979. [40] J. Nocedal, A. Sartenaer, and C. Zhu. On the behavior of the gradient norm in the steepest descent method. Computational Optimization and Applications, 22(1):5–35, 2002. [41] J. Nocedal and S. J. Wright. Numerical optimization. Springer Series in Operations Research. Springer-Verlag, New York, 1999. 18 O N M ETHODS FOR S YMMETRIC L INEAR E QUATIONS ARISING IN O PTIMIZATION [42] S. S. Oren and D. G. Luenberger. Self-scaling variable metric (SSVM) algorithms. I. Criteria and sufficient conditions for scaling a class of algorithms. Management Sci., 20:845–862, 1973/74. Mathematical programming. [43] C. C. Paige. Krylov subspace processes, Krylov subspace methods, and iteration polynomials. In Proceedings of the Cornelius Lanczos International Centenary Conference (Raleigh, NC, 1993), pages 83–92, Philadelphia, PA, 1994. SIAM. [44] C. C. Paige and M. A. Saunders. Solutions of sparse indefinite systems of linear equations. SIAM J. Numer. Anal., 12(4):617–629, 1975. [45] S. S. Petrova and A. D. Solov’ev. The origin of the method of steepest descent. Historia Mathematica, 24(4):361 – 375, 1997. [46] Y. Saad. Iterative methods for sparse linear systems. Society for Industrial and Applied Mathematics, Philadelphia, PA, second edition, 2003. [47] D. F. Shanno. Conditioning of quasi-Newton methods for function minimization. Math. Comp., 24:647–656, 1970. [48] D. F. Shanno. Conjugate gradient methods with inexact searches. Math. Oper. Res., 3(3):244–256, 1978. [49] D. F. Shanno and P. C. Kettler. Optimal conditioning of quasi-Newton methods. Math. Comp., 24:657–664, 1970. [50] J. R. Shewchuk. An introduction to the conjugate gradient method without the agonizing pain. Technical report, Pittsburgh, PA, USA, 1994. [51] E. Stiefel. Relaxationsmethoden bester Strategie zur Lösung linearer Gleichungssysteme. Comment. Math. Helv., 29:157–179, 1955. [52] A. Tits. Some investigations about a unified approach to quadratically convergent algorithms for function minimization. J. Optimization Theory Appl., 20(4):489–496, 1976. [53] S. J. Wright. Primal-dual interior-point methods. Society for Industrial and Applied Mathematics (SIAM), Philadelphia, PA, 1997.