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 positive-definite. The conjugate gradient method is often implemented as an iterative algorithm, applicable to sparse systems that are too large to be handled by a direct implementation or other direct methods such as the Cholesky decomposition. Large sparse systems often arise when numerically solving partial differential equations or optimization problems. The conjugate gradient method can also be used to solve unconstrained optimization problems such as energy minimization. It is commonly attributed to  Magnus Hestenes and Eduard Stiefel, who programmed it on the Z4, and extensively researched. The biconjugate gradient method provides a generalization to non-symmetric matrices. Various nonlinear conjugate gradient methods seek minima of nonlinear optimization problems. Description of the problem addressed by conjugate gradients. Suppose we want to solve the system of linear equations for the vector formula_2, where the known formula_3 matrix formula_4 is symmetric (i.e., A = A), positive-definite (i.e. xAx > 0 for all non-zero vectors formula_2 in R), and real, and formula_6 is known as well. We denote the unique solution of this system  by formula_7. Derivation as a direct method. The 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 well-known succinct formulation of the method. We say that two non-zero vectors u and v are conjugate (with respect to formula_4) if Since formula_4 is symmetric and positive-definite, the  left-hand side defines an inner product Two vectors are conjugate if and only if they are orthogonal with respect to this inner product. Being conjugate is a symmetric relation: if formula_12 is conjugate to formula_13, then formula_13 is conjugate to formula_12. Suppose that is a set of formula_17 mutually conjugate vectors with respect to formula_4, i.e. formula_19 for all formula_20. Then formula_21 forms a basis for formula_22, and we may express the solution formula_7 of formula_24 in this basis: Left-multiplying by formula_26 yields and so This gives the following method for solving the equation : find a sequence of  formula_17 conjugate directions, and then compute the coefficients formula_30. As an iterative method. If we choose the conjugate vectors formula_31 carefully, then we may not need all of them to obtain a good approximation to the solution formula_7. So, we want to regard the conjugate gradient method as an iterative method. This also allows us to approximately solve systems where "n" is so large that the direct method would take too much time. We denote the initial guess for by (we can assume without loss of generality that , otherwise consider the system Az = b − Ax instead). Starting with  x we search for the solution and in each iteration we need a metric to tell us whether we are closer to the solution (that is unknown to us). This metric comes from the fact that the solution is also the unique minimizer of the following quadratic function The existence of a unique minimizer is apparent as its Hessian matrix of second derivatives is symmetric positive-definite and that the minimizer (use D"f"(x)=0) solves the initial problem is obvious from its first derivative This suggests taking the first basis vector p to be the negative of the gradient of "f  at x = x. The gradient of f" equals . Starting with an initial guess x, this means we take p = b − Ax. The other vectors in the basis will be conjugate to the gradient, hence the name "conjugate gradient method". Note that p is also the residual provided by this initial step of the algorithm. Let r be the residual at the "k"th step: As observed above, formula_37 is the negative gradient of formula_38 at formula_39, so the gradient descent method would require to move in the direction r. Here, however, we insist that the directions formula_31 be  conjugate to each other. A practical way to enforce this is by requiring that the next search direction be built out of the current residual and all previous search directions. The conjugation constraint is an orthonormal-type constraint and hence the algorithm can be viewed as an example of Gram-Schmidt 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 follows from the definition of formula_37 . The expression for formula_45  can be derived if one substitutes the expression for x into "f" and minimizing it w.r.t. formula_45 The resulting algorithm. The above algorithm gives the most straightforward explanation of the conjugate gradient method. Seemingly, the algorithm as stated requires storage of all previous searching directions and residue vectors, as well as many matrix-vector multiplications, and thus can be computationally expensive. However, a closer analysis of the algorithm shows that formula_48 is orthogonal to formula_49, i.e. formula_50 , for i ≠ j. And formula_51is formula_4-orthogonal to formula_53 , i.e. formula_54 , for formula_20. This can be regarded that as the algorithm progresses , formula_51 and formula_48 span the same Krylov subspace. Where formula_48 form the orthogonal basis with respect to the standard inner product, and formula_51 form the orthogonal basis with respect to the inner product induced by formula_4. Therefore, formula_39 can be regarded as the projection of formula_2 on the Krylov subspace. The algorithm is detailed below for solving Ax = b where formula_4 is a real, symmetric, positive-definite matrix. The input vector formula_64 can be an approximate initial solution or 0. It is a different formulation of the exact procedure described above. This is the most commonly used algorithm. The  same formula for is also used in the Fletcher–Reeves nonlinear conjugate gradient method. Restarts. We note that formula_66 is computed by the gradient descent method applied to formula_67. Setting formula_68 would similarly make formula_69 computed by the gradient descent method from formula_70, i.e., can be used as a simple implementation of a restart of the conjugate gradient iterations. Restarts could slow down convergence, but may improve stability if the conjugate gradient method misbehaves, e.g., due to round-off error. Explicit residual calculation. The formulas formula_71 and formula_72, which both hold in exact arithmetic, make the formulas formula_73 and formula_74  mathematically equivalent. The former is used in the algorithm to avoid an extra multiplication by formula_4 since the vector formula_76 is already computed to evaluate formula_30. The latter may be more accurate, substituting the explicit calculation formula_74 for the implicit one by the recursion subject to round-off error accumulation, and is thus recommended for an occasional evaluation. A norm of the residual is typically used for stopping criteria. The norm of the explicit residual formula_74 provides a guaranteed level of accuracy both in exact arithmetic and in the presence of the rounding errors, where convergence naturally stagnates. In contrast , the implicit residual formula_73 is known to keep getting smaller in amplitude well below the level of rounding errors and thus cannot be used to determine the stagnation of convergence. Computation of alpha and beta. In the algorithm, is chosen such that formula_81 is orthogonal to formula_82. The denominator is simplified from since formula_84. The is chosen such that formula_85 is conjugate to formula_86. Initially, is using and equivalently formula_89 the numerator of is rewritten as because formula_81 and formula_82 are orthogonal by design. The denominator is rewritten as using that the search directions p are conjugated and again  that the residuals are orthogonal. This gives the in the algorithm after cancelling . Example code in MATLAB / GNU Octave. function x = conjgrad(A, b, x) end Numerical example. Consider 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. For reference, the exact solution is Our first step is to calculate the residual vector r associated with x. This residual is computed from the formula r = b - Ax, and in our case is equal to  Since this is the first iteration, we will use the residual vector r as our initial search direction p; the method of selecting p will change in further iterations. We now compute the scalar using the relationship We can now compute x using the formula This result completes the first iteration, the result being an "improved" approximate solution to the system, x. We may now move on and compute the next residual vector r using the formula Our next step in the process is to compute the scalar that will eventually be used to determine the next search direction p.  Now, using this scalar , we can compute the next search direction p using the relationship We now compute the scalar using our newly acquired p using the same method as that used for . Finally, we find x using the same method as that used to find x. The result, x, is a "better" approximation to the system's solution than x and x. If exact arithmetic were to be used in this example instead of limited-precision, then the exact solution would theoretically have been reached after "n" = 2 iterations ("n" being the order of the system). Convergence properties . The conjugate gradient method can theoretically be viewed as a direct method, as in the absence of round-off error it produces the exact solution after a finite number of iterations, which is not larger than the size of the matrix. In practice, the exact solution is never obtained since the conjugate gradient method is unstable with respect to even small perturbations, e.g., most directions are not in practice conjugate, due to a degenerative nature of generating the Krylov subspaces. As an iterative method, the conjugate gradient method monotonically (in the energy norm) improves approximations formula_70 to the exact solution  and 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 formula_106 of the system matrix formula_107: the larger formula_106 is, the slower the improvement. If formula_106 is large, preconditioning is commonly used to replace the original system formula_110 with formula_111 such that formula_112 is smaller than formula_113, see below. Convergence theorem. Define a subset of polynomials as where formula_115 is the set of polynomials of maximal degree formula_116. Let formula_117 be the iterative approximations of the exact solution  formula_118, and define the errors as formula_119. Now, the rate of convergence can be approximated as where formula_121 denotes the spectrum, and formula_122 denotes the condition number. Note, the important limit when formula_122 tends to formula_124 This limit shows a faster convergence rate compared to the iterative methods of Jacobi or Gauss–Seidel which scale as formula_126. No round-off error is assumed in the convergence theorem, but the convergence bound is commonly valid in practice as theoretically explained by Anne Greenbaum. Practical convergence. If initialized randomly, the first stage of iterations is often the fastest, as the error is  eliminated within the Krylov subspace that initially reflects a smaller effective condition number. The second stage of convergence is typically well defined by the theoretical convergence bound with formula_127, but may be super-linear, depending on a distribution of the spectrum of the matrix formula_107 and the spectral distribution of the error. In the last stage, the smallest attainable accuracy is reached and the convergence stalls or the method may even start diverging. In typical scientific computing applications in double-precision floating-point format for matrices of large sizes, the conjugate gradient method uses a stopping criteria with a tolerance  that terminates the iterations during the first or second stage. The preconditioned conjugate gradient method. In most cases, preconditioning is necessary to ensure fast convergence of the conjugate gradient method. The preconditioned conjugate gradient method takes the following form: The above formulation is equivalent to applying the conjugate gradient method without preconditioning to the system where The preconditioner matrix M has to be symmetric positive-definite 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. An example of  a commonly used preconditioner is the incomplete Cholesky factorization. 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 positive-definite 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. This version of the preconditioned conjugate gradient method can  be called flexible, as it allows for variable preconditioning. The flexible version is also shown to be robust even if the preconditioner is not symmetric positive definite (SPD). The implementation of the flexible version requires storing an extra vector. For a fixed SPD preconditioner, formula_144 so both formulas for are equivalent in exact arithmetic, i.e., without the round-off 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 does not converge slower than the locally optimal steepest descent method . Vs. the locally optimal steepest descent method. In both the original and the preconditioned conjugate gradient methods one only needs to set formula_145 in order to make them 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 and/or non-SPD preconditioner is used, see above. Conjugate gradient method as optimal feedback controller for  double integrator. The conjugate gradient method can also be derived using optimal control theory. In this approach, the conjugate gradient method falls out as an optimal feedback controller,formula_146 for the double integrator system,formula_147 The quantities formula_148 and formula_149 are variable feedback gains. Conjugate gradient on the normal equations. The conjugate gradient method can be applied to an arbitrary "n"-by-"m" matrix by applying it to normal equations AA and right-hand side vector Ab, since AA is a symmetric positive-semidefinite matrix for any A. The result is conjugate gradient on the normal equations (CGNR). As an iterative method , it is not necessary to form AA explicitly in memory but only to perform the matrix-vector and transpose matrix-vector 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 κ(AA) is equal to κ(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 ill-conditioned, i.e., A has a large condition number. Conjugate gradient method for complex Hermitian matrices. The conjugate gradient method with a trivial modification is extendable to solving, given complex-valued matrix A and vector b, the system of linear equations formula_150 for the complex-valued vector x, where A is Hermitian (i.e., A' = A) and positive-definite matrix, and the symbol ' denotes the conjugate transpose using the MATLAB/GNU Octave style. The trivial modification is simply substituting the conjugate transpose for the  real transpose everywhere. This substitution is backward compatible, since conjugate transpose turns into real transpose on real-valued vectors and matrices. The provided above Example code in MATLAB/GNU Octave thus already works for complex Hermitian matrices needed no modification.