Conjugate gradient method and algorithm

Conjugate gradient method

The conjugate gradient (CG) method solves $Ax = b$ for a symmetric, positive definite matrix $A$. It is equivalent to minimising the quadratic energy

$$\Pi(x) = \tfrac{1}{2}x^T A x - b^T x$$

whose minimum is exactly the solution of $Ax=b$. The minimum is searched along directions $d_k$ that are A-conjugate, i.e. $d_i^T A d_j = 0$ for $i \ne j$. Each step minimises $\Pi$ along one direction without spoiling the earlier ones, so in exact arithmetic the solution is reached in at most $n$ steps for an $n \times n$ system. The residual $r = b - Ax$ is the negative gradient of $\Pi$.

Compared with direct elimination, CG needs only matrix-vector products, so it stores only the non-zero terms of $A$. It is good for large sparse systems, and rounding errors do not accumulate in a long elimination.

Algorithm

1. Choose a starting vector $x_0$ (usually $0$). Compute the residual $r_0 = b - Ax_0$ and set the first search direction $d_0 = r_0$.
2. For $k = 0, 1, 2, \dots$ repeat:
   • Step length: $\alpha_k = \dfrac{r_k^T r_k}{d_k^T A d_k}$
   • Update the solution: $x_{k+1} = x_k + \alpha_k d_k$
   • Update the residual: $r_{k+1} = r_k - \alpha_k A d_k$
   • Check convergence: stop if $\|r_{k+1}\| <$ tolerance.
   • Direction factor: $\beta_k = \dfrac{r_{k+1}^T r_{k+1}}{r_k^T r_k}$
   • New direction: $d_{k+1} = r_{k+1} + \beta_k d_k$
3. The last $x_{k+1}$ is the solution.

Convergence is faster if $A$ has clustered eigenvalues, and preconditioning is used to improve this.