Conjugate Gradient

Worked Example

Consider the symmetric positive-definite system Ax = b with A and b defined as follows, and an initial guess x₀ = 0:

1. Initial State: x₀ = [0, 0]. The initial residual is r₀ = b - Ax₀ = [2, 1]. The initial search direction is p₀ = r₀ = [2, 1].
2. Step 1:
   • Compute Ap₀ = [{3, 1}, {1, 2}] [2, 1] = [7, 4].
   • Compute α₀ = (r₀ · r₀) / (p₀ · Ap₀) = (4 + 1) / (14 + 4) = 5 / 18.
   • Update x₁ = x₀ + α₀ p₀ = [0, 0] + (5/18) [2, 1] = [5/9, 5/18].
   • Update r₁ = r₀ - α₀ Ap₀ = [2, 1] - (5/18) [7, 4] = [1/18, -1/9].
3. Direction Update:
   • Compute β₀ = (r₁ · r₁) / (r₀ · r₀) = (1/324 + 1/81) / 5 = (5/324) / 5 = 1/324.
   • Update p₁ = r₁ + β₀ p₀ = [1/18, -1/9] + (1/324) [2, 1] = [5/81, -35/324].
4. Step 2:
   • Compute Ap₁. After exact math, we get α₁ = 18/25.
   • Update x₂ = x₁ + α₁ p₁ = [5/9, 5/18] + (18/25) [5/81, -35/324] = [3/5, 1/5].
   • Update r₂ = r₁ - α₁ Ap₁ = [0, 0].

The exact solution is x = [3/5, 1/5], which we obtained in exactly 2 steps, a characteristic of CG on a 2x2 SPD matrix.

Experiment

Open the Conjugate Gradient module to interact with this concept:

1. Set A11 to 3, A12 to 1, and A22 to 2. Set b1 to 2 and b2 to 1.
2. Set Updates to 1. The status will show "Finite iterate; not yet converged". The iterate is [5/9, 5/18] and the residual is [1/18, -1/9].
3. Increase Updates to 2. The algorithm converges exactly to [3/5, 1/5] with a zero residual.
4. Change A11 to 1, A12 to 2, and A22 to 1. This matrix is not positive-definite because its determinant is 1*1 - 2*2 = -3. The module will reject it and display "Outside SPD domain".
5. Set b1 to 0 and b2 to 0 with an SPD matrix. The initial residual is zero, so the algorithm stops immediately without dividing by zero.

Limitations

This explorer is designed for 2x2 real symmetric positive-definite (SPD) matrices. The Conjugate Gradient method requires the matrix to be SPD to guarantee convergence and avoid division by zero during the computation of the step sizes (α and β). If a non-SPD matrix is provided, it is rejected using exact principal-minor tests. The visualization only shows standard unpreconditioned linear CG. In a real-world scenario with large matrices, floating-point rounding errors might prevent the residual from becoming exactly zero, requiring a stopping criterion based on a tolerance (e.g., 1e-9).

References

The standard Conjugate Gradient algorithm:

• SciPy sparse.linalg.cg documentation describes the SPD system requirement and the actual residual convergence criterion. The exact rational CG coefficients were independently checked against the 2x2 inverse solution.