 Conjugate gradient method

In mathematics, the conjugate gradient method is an algorithm for the numerical solution of particular systems of linear equations, namely those whose matrix is symmetric and positivedefinite. The conjugate gradient method is an iterative method, so it can be applied to sparse systems that are too large to be handled by direct methods such as the Cholesky decomposition. Such systems often arise when numerically solving partial differential equations.
The conjugate gradient method can also be used to solve unconstrained optimization problems such as energy minimization. It was developed by Hestenes and Stiefel.^{[1]}
The biconjugate gradient method provides a generalization to nonsymmetric matrices. Various nonlinear conjugate gradient methods seek minima of nonlinear equations.
Description of the method
Suppose we want to solve the following system of linear equations
 Ax = b
where the nbyn matrix A is symmetric (i.e., A^{T} = A), positive definite (i.e., x^{T}Ax > 0 for all nonzero vectors x in R^{n}), and real.
We denote the unique solution of this system by x_{*}.
The conjugate gradient method as a direct method
We say that two nonzero vectors u and v are conjugate (with respect to A) if
Since A is symmetric and positive definite, the lefthand side defines an inner product
So, two vectors are conjugate if they are orthogonal with respect to this inner product. Being conjugate is a symmetric relation: if u is conjugate to v, then v is conjugate to u. (Note: This notion of conjugate is not related to the notion of complex conjugate.)
Suppose that {p_{k}} is a sequence of n mutually conjugate directions. Then the p_{k} form a basis of R^{n}, so we can expand the solution x_{*} of Ax = b in this basis:
The coefficients are given by
 (because are mutually conjugate)
This result is perhaps most transparent by considering the inner product defined above.
This gives the following method for solving the equation Ax = b: find a sequence of n conjugate directions, and then compute the coefficients α_{k}.
The conjugate gradient method as an iterative method
If we choose the conjugate vectors p_{k} carefully, then we may not need all of them to obtain a good approximation to the solution x_{*}. So, we want to regard the conjugate gradient method as an iterative method. This also allows us to solve systems where n is so large that the direct method would take too much time.
We denote the initial guess for x_{*} by x_{0}. We can assume without loss of generality that x_{0} = 0 (otherwise, consider the system Az = b − Ax_{0} instead). Starting with x_{0} we search for the solution and in each iteration we need a metric to tell us whether we are closer to the solution x_{*} (that is unknown to us). This metric comes from the fact that the solution x_{*} is also the unique minimizer of the following quadratic function; so if f(x) becomes smaller in an iteration it means that we are closer to x_{*}.
This suggests taking the first basis vector p_{1} to be the negative of the gradient of f at x = x_{0}. This gradient equals Ax_{0}−b. Since x_{0} = 0, this means we take p_{1} = b. The other vectors in the basis will be conjugate to the gradient, hence the name conjugate gradient method.
Let r_{k} be the residual at the kth step:
Note that r_{k} is the negative gradient of f at x = x_{k}, so the gradient descent method would be to move in the direction r_{k}. Here, we insist that the directions p_{k} be conjugate to each other. We also require the next search direction is built out of the current residue and all previous search directions, which is reasonable enough in practice.
The conjugation constraint is an orthonormaltype constraint and hence the algorithm bears resemblance to GramSchmidt orthonormalization.
This gives the following expression:
(see the picture at the top of the article for the effect of the conjugacy constraint on convergence). Following this direction, the next optimal location is given by
with
where the last equality holds because p_{k+1} and x_{k} are conjugate.
The resulting algorithm
The above algorithm gives the most straightforward explanation towards the conjugate gradient method. However, it requires storage of all the previous searching directions and residue vectors, and many matrix vector multiplications, thus could be computationally expensive. In practice, one modifies slightly the condition obtaining the last residue vector, not to minimize the metric following the search direction, but instead to make it orthogonal to the previous residue. Minimization of the metric along the search direction will be obtained automatically in this case. One can then result in an algorithm which only requires storage of the last two residue vectors and the last search direction, and only one matrix vector multiplication. Note that the algorithm described below is equivalent to the previously discussed straightforward procedure.
The algorithm is detailed below for solving Ax = b where A is a real, symmetric, positivedefinite matrix. The input vector x_{0} can be an approximate initial solution or 0.
 repeat
 if r_{k+1} is sufficiently small then exit loop end if
 end repeat
 The result is x_{k+1}
This is the most commonly used algorithm. The same formula for β_{k} is also used in the Fletcher–Reeves nonlinear conjugate gradient method.
Example code in GNU Octave
function [x] = conjgrad(A,b,x) r=bA*x; p=r; rsold=r'*r; for i=1:size(A)(1) Ap=A*p; alpha=rsold/(p'*Ap); x=x+alpha*p; r=ralpha*Ap; rsnew=r'*r; if sqrt(rsnew)<1e10 break; end p=r+rsnew/rsold*p; rsold=rsnew; end end
Numerical example
To illustrate the conjugate gradient method, we will complete a simple example.
Considering the linear system Ax = b given by
we will perform two steps of the conjugate gradient method beginning with the initial guess
in order to find an approximate solution to the system.
Solution
Our first step is to calculate the residual vector r_{0} associated with x_{0}. This residual is computed from the formula r_{0} = b  Ax_{0}, and in our case is equal to
Since this is the first iteration, we will use the residual vector r_{0} as our initial search direction p_{0}; the method of selecting p_{k} will change in further iterations.
We now compute the scalar α_{0} using the relationship
We can now compute x_{1} using the formula
This result completes the first iteration, the result being an "improved" approximate solution to the system, x_{1}. We may now move on and compute the next residual vector r_{1} using the formula
Our next step in the process is to compute the scalar β_{0} that will eventually be used to determine the next search direction p_{1}.
Now, using this scalar β_{0}, we can compute the next search direction p_{1} using the relationship
We now compute the scalar α_{1} using our newlyacquired p_{1} using the same method as that used for α_{0}.
Finally, we find x_{2} using the same method as that used to find x_{1}.
The result, x_{2}, is a "better" approximation to the system's solution than x_{1} and x_{0}. If exact arithmetic were to be used in this example instead of limitedprecision, then the exact solution would theoretically have been reached after n = 2 iterations (n being the order of the system).
Convergence properties of the conjugate gradient method
The conjugate gradient method can theoretically be viewed as a direct method, as it produces the exact solution after a finite number of iterations, which is not larger than the size of the matrix, in the absence of roundoff error. However, the conjugate gradient method is unstable with respect to even small perturbations, e.g., most directions are not in practice conjugate, and the exact solution is never obtained. Fortunately, the conjugate gradient method can be used as an iterative method as it provides monotonically improving approximations to the exact solution, which may reach the required tolerance after a relatively small (compared to the problem size) number of iterations. The improvement is typically linear and its speed is determined by the condition number κ(A) of the system matrix A: the larger is κ(A), the slower the improvement ^{[2]} .
If κ(A) is large, preconditioning is used to replace the original system with so that gets smaller than , see below.
The preconditioned conjugate gradient method
See also: PreconditionerIn most cases, preconditioning is necessary to ensure fast convergence of the conjugate gradient method. The preconditioned conjugate gradient method takes the following form:
 repeat
 if r_{k+1} is sufficiently small then exit loop end if
 end repeat
 The result is x_{k+1}
The above formulation is equivalent to applying the conjugate gradient method without preconditioning to the system^{[1]}
where and .
The preconditioner matrix M has to be symmetric positivedefinite and fixed, i.e., cannot change from iteration to iteration. If any of these assumptions on the preconditioner is violated, the behavior of the preconditioned conjugate gradient method may become unpredictable.
The flexible preconditioned conjugate gradient method
In numerically challenging applications, sophisticated preconditioners are used, which may lead to variable preconditioning, changing between iterations. Even if the preconditioner is symmetric positivedefinite on every iteration, the fact that it may change makes the arguments above invalid, and in practical tests leads to a significant slow down of the convergence of the algorithm presented above. Using the Polak–Ribière formula
instead of the Fletcher–Reeves formula
may dramatically improve the convergence in this case.^{[3]} This version of the preconditioned conjugate gradient method can be called^{[4]} flexible, as it allows for variable preconditioning. The implementation of the flexible version requires storing an extra vector. For a fixed preconditioner, so both formulas for β_{k} are equivalent in exact arithmetic, i.e., without the roundoff error.
The mathematical explanation of the better convergence behavior of the method with the Polak–Ribière formula is that the method is locally optimal in this case, in particular, it converges not slower than the locally optimal steepest descent method.^{[5]}
The conjugate gradient method vs. the locally optimal steepest descent method
In both the original and the preconditioned conjugate gradient methods one only needs to always set β_{k}: = 0 in order to turn them into locally optimal, using the line search, steepest descent methods. With this substitution, vectors are always the same as vectors , so there is no need to store vectors . Thus, every iteration of these steepest descent methods is a bit cheaper compared to that for the conjugate gradient methods. However, the latter converge faster, unless a (highly) variable preconditioner is used, see above.
Derivation of the method
Main article: Derivation of the conjugate gradient methodThe conjugate gradient method can be derived from several different perspectives, including specialization of the conjugate direction method for optimization, and variation of the Arnoldi/Lanczos iteration for eigenvalue problems. Despite differences in their approaches, these derivations share a common topic—proving the orthogonality of the residuals and conjugacy of the search directions. These two properties are crucial to developing the wellknown succinct formulation of the method.
Conjugate gradient on the normal equations
The conjugate gradient method can be applied to an arbitrary nbym matrix by applying it to normal equations A^{T}A and righthand side vector A^{T}b, since A^{T}A is a symmetric positivesemidefinite matrix for any A. The result is conjugate gradient on the normal equations (CGNR).
 A^{T}Ax = A^{T}b
As an iterative method, it is not necessary to form A^{T}A explicitly in memory but only to perform the matrixvector and transpose matrixvector multiplications. Therefore CGNR is particularly useful when A is a sparse matrix since these operations are usually extremely efficient. However the downside of forming the normal equations is that the condition number κ(A^{T}A) is equal to κ^{2}(A) and so the rate of convergence of CGNR may be slow and the quality of the approximate solution may be sensitive to roundoff errors. Finding a good preconditioner is often an important part of using the CGNR method.
Several algorithms have been proposed (e.g., CGLS, LSQR). The LSQR algorithm purportedly has the best numerical stability when A is illconditioned, i.e., A has a large condition number.
See also
 Biconjugate gradient method (BiCG)
 Conjugate residual method
 Nonlinear conjugate gradient method
 Iterative method. Linear systems
 Preconditioning
 Gaussian Belief Propagation
References
The conjugate gradient method was originally proposed in
 Hestenes, Magnus R.; Stiefel, Eduard (December 1952). "Methods of Conjugate Gradients for Solving Linear Systems" (PDF). Journal of Research of the National Bureau of Standards 49 (6). http://nvl.nist.gov/pub/nistpubs/jres/049/6/V49.N06.A08.pdf.
Descriptions of the method can be found in the following text books:
 Atkinson, Kendell A. (1988). "Section 8.9". An introduction to numerical analysis (2nd ed.). John Wiley and Sons. ISBN 0471500232.
 Avriel, Mordecai (2003). Nonlinear Programming: Analysis and Methods. Dover Publishing. ISBN 0486432270.
 Golub, Gene H.; Van Loan, Charles F.. "Chapter 10". Matrix computations (3rd ed.). Johns Hopkins University Press. ISBN 0801854148.
 Saad, Yousef. "Chapter 6". Iterative methods for sparse linear systems (2nd ed.). SIAM. ISBN 9780898715347.
Notes
 ^ Straeter, T. A.. "On the Extension of the DavidonBroyden Class of Rank One, QuasiNewton Minimization Methods to an Infinite Dimensional Hilbert Space with Applications to Optimal Control Problems". NASA Technical Reports Server. NASA. http://hdl.handle.net/2060/19710026200. Retrieved 10 October 2011.
 ^ Saad, Yousef (2003). Iterative methods for sparse linear systems (2nd ed.  ed.). Philadelphia, Pa.: Society for Industrial and Applied Mathematics. pp. 195. ISBN 9780898715347.
 ^ Golub, Gene H.; Ye, Qiang (1999). "Inexact Preconditioned Conjugate Gradient Method with InnerOuter Iteration". SIAM Journal on Scientific Computing 21 (4): 1305. doi:10.1137/S1064827597323415.
 ^ Notay, Yvan (2000). "Flexible Conjugate Gradients". SIAM Journal on Scientific Computing 22 (4): 1444. doi:10.1137/S1064827599362314.
 ^ Knyazev, Andrew V.; Lashuk, Ilya (2008). "Steepest Descent and Conjugate Gradient Methods with Variable Preconditioning". SIAM Journal on Matrix Analysis and Applications 29 (4): 1267. doi:10.1137/060675290.
External links
 An Introduction to the Conjugate Gradient Method Without the Agonizing Pain by Jonathan Richard Shewchuk.
 Iterative methods for sparse linear systems by Yousef Saad
 LSQR: Sparse Equations and Least Squares by Christopher Paige and Michael Saunders.
Numerical linear algebra Key concepts Problems Hardware Software Categories: Numerical linear algebra
 Gradient methods
Wikimedia Foundation. 2010.
Look at other dictionaries:
Nonlinear conjugate gradient method — In numerical optimization, the nonlinear conjugate gradient method generalizes the conjugate gradient method to nonlinear optimization. For a quadratic function : The minimum of f is obtained when the gradient is 0: . Whereas linear conjugate… … Wikipedia
Derivation of the conjugate gradient method — In numerical linear algebra, the conjugate gradient method is an iterative method for numerically solving the linear system where is symmetric positive definite. The conjugate gradient method can be derived from several different perspectives,… … Wikipedia
Preconditioned conjugate gradient method — The conjugate gradient method is a numerical algorithm that solves a system of linear equations:A x= b.,where A is symmetric [positive definite] . If the matrix A is ill conditioned, i.e. it has a large condition number kappa(A), it is often… … Wikipedia
Conjugate residual method — The conjugate residual method is an iterative numeric method used for solving systems of linear equations. It s a Krylov subspace method very similar to the much more popular conjugate gradient method, with similar construction and convergence… … Wikipedia
Biconjugate gradient method — In mathematics, more specifically in numerical analysis, the biconjugate gradient method is an algorithm to solve systems of linear equations :A x= b.,Unlike the conjugate gradient method, this algorithm does not require the matrix A to be self… … Wikipedia
Gradient descent — For the analytical method called steepest descent see Method of steepest descent. Gradient descent is an optimization algorithm. To find a local minimum of a function using gradient descent, one takes steps proportional to the negative of the… … Wikipedia
Iterative method — In computational mathematics, an iterative method is a mathematical procedure that generates a sequence of improving approximate solutions for a class of problems. A specific implementation of an iterative method, including the termination… … Wikipedia
Nelder–Mead method — Nelder–Mead simplex search over the Rosenbrock banana function (above) and Himmelblau s function (below) See simplex algorithm for Dantzig s algorithm for the problem of linear opti … Wikipedia
Finite element method — The finite element method (FEM) (sometimes referred to as finite element analysis) is a numerical technique for finding approximate solutions of partial differential equations (PDE) as well as of integral equations. The solution approach is based … Wikipedia
Folded spectrum method — In mathematics, the folded spectrum method (FSM) is a iterative method for solving large eigenvalue problems.Here you always find a vector with an eigenvalue close to a search value varepsilon. This means you can get a vector Psi in the middle of … Wikipedia