Linear solution techniques; CG algorithm

Linear system $[A]\{x\}=\{b\}$ is solved by:

1. Direct methods (finite number of steps, exact apart from round-off): Gauss elimination (reduce to upper triangular form, back-substitute), Gauss-Jordan, LU decomposition, Cholesky decomposition $A=LL^T$ for symmetric positive definite $A$, and the matrix inverse. In FEM they use banded or skyline storage.
2. Iterative methods (start from a guess and improve): Jacobi, Gauss-Seidel, successive over-relaxation, and the conjugate gradient method. They suit very large sparse systems since the matrix is never changed and no fill-in occurs. Jacobi and Gauss-Seidel need diagonal dominance for convergence; CG needs a symmetric positive definite matrix.

Algorithm of the 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$.

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.