Conjugate gradient method: Information from Answers.com

A comparison of the convergence of gradient descent with optimal step size (in green) and conjugate gradient (in red) for minimizing a quadratic function associated with a given linear system. Conjugate gradient, assuming exact arithmetics, converges in at most n steps where n is the size of the matrix of the system (here n=2).

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 positive-definite. 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.

The biconjugate gradient method provides a generalization to non-symmetric matrices. Various nonlinear conjugate gradient methods seek minima of nonlinear equations.

1 Description of the method
2 Conjugate gradient on the normal equations
3 See also
4 References
5 External links

Description of the method

Suppose we want to solve the following system of linear equations

Ax = b

where the n-by-n matrix A is symmetric (i.e., A^T = A), positive definite (i.e., x^TAx > 0 for all non-zero vectors x in Rⁿ), and real.

We denote the unique solution of this system by x_*.

The conjugate gradient method as a direct method

We say that two non-zero vectors u and v are conjugate (with respect to A) if

$\mathbf{u}^\top \mathbf{A} \mathbf{v} = 0.$

Since A is symmetric and positive definite, the left-hand side defines an inner product

$\langle \mathbf{u},\mathbf{v} \rangle_\mathbf{A} := \langle \mathbf{A}^\top \mathbf{u}, \mathbf{v}\rangle = \langle \mathbf{A} \mathbf{u}, \mathbf{v}\rangle = \langle \mathbf{u}, \mathbf{A}\mathbf{v} \rangle = \mathbf{u}^\top \mathbf{A} \mathbf{v}.$

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ⁿ, so we can expand the solution x_* of Ax = b in this basis:

$\mathbf{x}_* = \sum^{n}_{i=1} \alpha_i \mathbf{p}_i$

The coefficients are given by

$\mathbf{b}=\mathbf{A}\mathbf{x}_* = \sum^{n}_{i=1} \alpha_i \mathbf{A} \mathbf{p}_i.$

$\mathbf{p}_k^\top \mathbf{b}=\mathbf{p}_k^\top \mathbf{A}\mathbf{x}_* = \sum^{n}_{i=1} \alpha_i\mathbf{p}_k^\top \mathbf{A} \mathbf{p}_i=\alpha_k\mathbf{p}_k^\top \mathbf{A} \mathbf{p}_k.$ (because $\forall i\neq k, p_i, p_k$ are mutually conjugate)

$\alpha_k = \frac{\mathbf{p}_k^\top \mathbf{b}}{\mathbf{p}_k^\top \mathbf{A} \mathbf{p}_k} = \frac{\langle \mathbf{p}_k, \mathbf{b}\rangle}{\,\,\,\langle \mathbf{p}_k, \mathbf{p}_k\rangle_\mathbf{A}} = \frac{\langle \mathbf{p}_k, \mathbf{b}\rangle}{\,\,\,\|\mathbf{p}_k\|_\mathbf{A}^2}.$

This result is perhaps most transparent by considering the inner product defined above.

This gives the following method for solving the equation Ax = b. We first find a sequence of n conjugate directions and then we 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₀. We can assume without loss of generality that x₀ = 0 (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 have gotten 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_*.

$f(\mathbf{x}) = \frac12 \mathbf{x}^\top \mathbf{A}\mathbf{x} - \mathbf{b}^\top \mathbf{x} , \quad \mathbf{x}\in\mathbf{R}^n.$

This suggests taking the first basis vector p₁ to be the gradient of f at x = x₀, which equals Ax₀−b. Since x₀ = 0, this means we take p₁ = −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:

$\mathbf{r}_k = \mathbf{b} - \mathbf{Ax}_k. \,$

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 are conjugate to each other, so we take the direction closest to the gradient r_k under the conjugacy constraint. This gives the following expression:

$\mathbf{p}_{k+1} = \mathbf{r}_k - \sum_{i\leq k}\frac{\mathbf{p}_i^\top \mathbf{A} \mathbf{r}_k}{\mathbf{p}_i^\top\mathbf{A} \mathbf{p}_i} \mathbf{p}_i$

(see the picture at the top of the article for the effect of the conjugacy constraint on convergence).

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. 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.

The algorithm is detailed below for solving Ax = b where A is a real, symmetric, positive-definite matrix. The input vector x₀ can be an approximate initial solution or 0.

$\mathbf{r}_0 := \mathbf{b} - \mathbf{A x}_0 \,$

$\mathbf{p}_0 := \mathbf{r}_0 \,$

$k := 0 \,$

repeat

$\alpha_k := \frac{\mathbf{r}_k^\top \mathbf{r}_k}{\mathbf{p}_k^\top \mathbf{A p}_k} \,$

$\mathbf{x}_{k+1} := \mathbf{x}_k + \alpha_k \mathbf{p}_k \,$

$\mathbf{r}_{k+1} := \mathbf{r}_k - \alpha_k \mathbf{A p}_k \,$

if r_k+1 is sufficiently small then exit loop end if

$\beta_k := \frac{\mathbf{r}_{k+1}^\top \mathbf{r}_{k+1}}{\mathbf{r}_k^\top \mathbf{r}_k} \,$

$\mathbf{p}_{k+1} := \mathbf{r}_{k+1} + \beta_k \mathbf{p}_k \,$

$k := k + 1 \,$

end repeat

The result is x_k+1

Example code in Octave

function [x] = conjgrad(A,b,x0)
 
   r = b - A*x0;
   w = -r;
   z = A*w;
   a = (r'*w)/(w'*z);
   x = x0 + a*w;
   B = 0;
 
   for i = 1:size(A)(1);
      r = r - a*z;
      if( norm(r) < 1e-10 )
           break;
      end
      B = (r'*z)/(w'*z);
      w = -r + B*w;
      z = A*w;
      a = (r'*w)/(w'*z);
      x = x + a*w;
   end
 
end

Numerical example

To illustrate the conjugate gradient method, we will complete a simple example.

Considering the linear system Ax = b given by

$\mathbf{A} \mathbf{x}= \begin{bmatrix} 4 & 1 \\ 1 & 3 \end{bmatrix}\begin{bmatrix} x_1 \\ x_2 \end{bmatrix} = \begin{bmatrix} 1 \\ 2 \end{bmatrix} ,$

we will perform two steps of the conjugate gradient method beginning with the initial guess

$\mathbf{x}_0 = \begin{bmatrix} 2 \\ 1 \end{bmatrix}$

in order to find an approximate solution to the system.

Solution

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

$\begin{bmatrix} 1 \\ 2 \end{bmatrix} - \begin{bmatrix} 4 & 1 \\ 1 & 3 \end{bmatrix} \begin{bmatrix} 2 \\ 1 \end{bmatrix} = \begin{bmatrix}-8 \\ -3 \end{bmatrix}.$

Since this is the first iteration, we will use the residual vector r₀ as our initial search direction p₀; the method of selecting p_k will change in further iterations.

We now compute the scalar α₀ using the relationship

$\alpha_0 = \frac{\mathbf{r}_0^\top \mathbf{r}_0}{\mathbf{p}_0^\top \mathbf{A p}_0} = \frac{\begin{bmatrix} -8 & -3 \end{bmatrix} \begin{bmatrix} -8 \\ -3 \end{bmatrix}}{ \begin{bmatrix} -8 & -3 \end{bmatrix} \begin{bmatrix} 4 & 1 \\ 1 & 3 \end{bmatrix} \begin{bmatrix} -8 \\ -3 \end{bmatrix} } = \frac{73}{331}.$

We can now compute x₁ using the formula

$\mathbf{x}_1 = \mathbf{x}_0 + \alpha_0\mathbf{p}_0 = \begin{bmatrix} 2 \\ 1 \end{bmatrix} + \frac{73}{331} \begin{bmatrix} -8 \\ -3 \end{bmatrix} = \begin{bmatrix} 0.2356 \\ 0.3384 \end{bmatrix}.$

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

$\mathbf{r}_1 = \mathbf{r}_0 - \alpha_0 \mathbf{A} \mathbf{p}_0 = \begin{bmatrix} -8 \\ -3 \end{bmatrix} - \frac{73}{331} \begin{bmatrix} 4 & 1 \\ 1 & 3 \end{bmatrix} \begin{bmatrix} -8 \\ -3 \end{bmatrix} = \begin{bmatrix} -0.2810 \\ 0.7492 \end{bmatrix}.$

Our next step in the process is to compute the scalar β₀ that will eventually be used to determine the next search direction p₁.

$\beta_0 = \frac{\mathbf{r}_1^\top \mathbf{r}_1}{\mathbf{r}_0^\top \mathbf{r}_0} = \frac{\begin{bmatrix} -0.2810 & 0.7492 \end{bmatrix} \begin{bmatrix} -0.2810 \\ 0.7492 \end{bmatrix}}{\begin{bmatrix} -8 & -3 \end{bmatrix} \begin{bmatrix} -8 \\ -3 \end{bmatrix}} = 0.0088.$

Now, using this scalar β₀, we can compute the next search direction p₁ using the relationship

$\mathbf{p}_1 = \mathbf{r}_1 + \beta_0 \mathbf{p}_0 = \begin{bmatrix} -0.2810 \\ 0.7492 \end{bmatrix} + 0.0088 \begin{bmatrix} -8 \\ -3 \end{bmatrix} = \begin{bmatrix} -0.3511 \\ 0.7229 \end{bmatrix}.$

We now compute the scalar α₁ using our newly-acquired p₁ using the same method as that used for α₁.

$\alpha_1 = \frac{\mathbf{r}_1^\top \mathbf{r}_1}{\mathbf{p}_1^\top \mathbf{A p}_1} = \frac{\begin{bmatrix} -0.2810 & 0.7492 \end{bmatrix} \begin{bmatrix} -0.2810 \\ 0.7492 \end{bmatrix}}{ \begin{bmatrix} -0.3511 & 0.7229 \end{bmatrix} \begin{bmatrix} 4 & 1 \\ 1 & 3 \end{bmatrix} \begin{bmatrix} -0.3511 \\ 0.7229 \end{bmatrix} } = 0.4122.$

Finally, we find x₂ using the same method as that used to find x₁.

$\mathbf{x}_2 = \mathbf{x}_1 + \alpha_1 \mathbf{p}_1 = \begin{bmatrix} 0.2356 \\ 0.3384 \end{bmatrix} + 0.4122 \begin{bmatrix} -0.3511 \\ 0.7229 \end{bmatrix} = \begin{bmatrix} 0.0909 \\ 0.6364 \end{bmatrix}.$

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 in stead of limited-precision, then the exact solution would theoretically have been reached after n = 2 iterations (n being the order of the system).

The preconditioned conjugate gradient method

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 A^TA and right-hand side vector A^Tb, since A^TA is a symmetric positive-semidefinite matrix for any A. The result is conjugate gradient on the normal equations (CGNR).

A^TAx = A^Tb

As an iterative method, it is not necessary to form A^TA 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 κ(A^TA) 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.

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:

Kendell A. Atkinson (1988), An introduction to numerical analysis (2nd ed.), Section 8.9, John Wiley and Sons. ISBN 0-471-50023-2.
Mordecai Avriel (2003). Nonlinear Programming: Analysis and Methods. Dover Publishing. ISBN 0-486-43227-0.
Gene H. Golub and Charles F. Van Loan, Matrix computations (3rd ed.), Chapter 10, Johns Hopkins University Press. ISBN 0-8018-5414-8.
Yousef Saad, Iterative methods for sparse linear systems (2nd ed.), Chapter 6, ???. ISBN 978-0898715347.

External links

^{^} An Introduction to the Conjugate Gradient Method Without the Agonizing Pain by Jonathan Richard Shewchuk.
Conjugate Gradient Method by Nadir Soualem.
Preconditioned conjugate gradient method by Nadir Soualem.
Iterative methods for sparse linear systems by Yousef Saad
LSQR: Sparse Equations and Least Squares by Christopher Paige and Michael Saunders.

This entry is from Wikipedia, the leading user-contributed encyclopedia. It may not have been reviewed by professional editors (see full disclaimer)

Donate to Wikimedia

	What is a gradient? Read answer...
	What is the Conjugation of pedir? Read answer...
	What is verb conjugation? Read answer...

	Why are bacterial transformation conjugation and transduction not considered to be methods of reproduction?
	What is the answer using Cartesian method to prove that gradient of vector b dot vector r equals b?
	What is a conjugate?

Conjugate gradient method

Contents