Conjugate Gradient on a Small SPD System

Solve a symmetric positive-definite linear system with conjugate gradient and watch it converge in at most n steps.

Problem

Conjugate gradient (CG) solves Ax=b when A is symmetric positive definite, using only matrix-vector products and no factorization. In exact arithmetic it converges in at most n iterations for an n-by-n system, and often far fewer when eigenvalues cluster.

Algorithm

CG builds a sequence of search directions that are A-conjugate (mutually orthogonal in the A inner product), so each step minimizes the error along a new direction without undoing previous progress. The residual doubles as the steepest-descent direction, corrected to stay conjugate.

import numpy as np
A=np.array([[4.,1.],[1.,3.]]); b=np.array([1.,2.])
x=np.zeros(2); r=b-A@x; p=r.copy(); rs=r@r
for k in range(2):
    Ap=A@p; alpha=rs/(p@Ap)
    x=x+alpha*p; r=r-alpha*Ap
    rs_new=r@r
    print(k,'x',np.round(x,4),'res',round(np.sqrt(rs_new),6))
    if rs_new<1e-12: break
    p=r+(rs_new/rs)*p; rs=rs_new

Result

For this 2x2 system CG reaches the exact solution in two steps, driving the residual to machine zero. The step length alpha is chosen to exactly minimize the quadratic along the current direction. In large sparse problems, such as discretized PDEs, CG is preferred over direct solvers because it never forms a dense factor and needs only the action of A.

• Convergence rate depends on the condition number; preconditioning clusters eigenvalues to speed it up.
• CG requires symmetry and positive definiteness; use GMRES or BiCGSTAB for general matrices.
• Kronos finite-element and finite-difference solvers rely on preconditioned CG for large stiffness and Poisson systems.