The Conjugate Gradient (CG) method of Hestenes and Stiefel is the premier iterative solver for large, sparse, symmetric positive definite linear systems. Unlike stationary iterations that converge linearly with rates determined by matrix splittings, the Conjugate Gradient method is an optimal Krylov subspace projection method. At step \(k\), it computes the unique approximation \(\bx^{(k)}\) that minimizes the error in the \(\bA\)-norm over the entire affine Krylov subspace \(\bx^{(0)} + \mathcal{K}_k(\bA, \br^{(0)})\).
Remarkably, despite optimizing over an expanding \(k\)-dimensional subspace, the symmetry of \(\bA\) allows each new iterate and search direction to be computed using a three-term recurrence involving only the immediately preceding step. The algorithm requires storing only four vectors of length \(n\) and computes only a single matrix-vector product per iteration.
Quadratic Optimization and Steepest Descent
The mathematical foundation of the Conjugate Gradient method lies in the equivalence between solving an SPD linear system and minimizing a convex quadratic energy functional.
Let \(\bA \in \fR^{n \times n}\) be a symmetric positive definite matrix and let \(\bb \in \fR^n\). Define the quadratic functional \(\phi: \fR^n \to \fR\) by \[
\begin{align}
\phi(\bx) = \frac{1}{2} \bx^T\bA\bx - \bb^T\bx.
\end{align}
\] The gradient of \(\phi\) is \(\nabla \phi(\bx) = \bA\bx - \bb = -\br(\bx)\), where \(\br(\bx) = \bb - \bA\bx\) is the residual. Because \(\bA\) is positive definite, the Hessian matrix \(\nabla^2 \phi(\bx) = \bA\) is strictly positive definite everywhere, and \(\phi\) has a unique global minimizer \(\bx^*\) satisfying \[
\begin{align}
\nabla \phi(\bx^*) = \bzero \iff \bA\bx^* = \bb.
\end{align}
\] Furthermore, the functional error is directly proportional to the squared \(\bA\)-norm of the solution error: \[
\begin{align}
\phi(\bx) - \phi(\bx^*) = \frac{1}{2} \|\bx - \bx^*\|_\bA^2.
\end{align}
\]
Taylor expanding \(\phi(\bx)\) about the exact solution \(\bx^* = \bA^{-1}\bb\) gives \[
\begin{align}
\phi(\bx) &= \phi(\bx^*) + \nabla \phi(\bx^*)^T(\bx - \bx^*) + \frac{1}{2} (\bx - \bx^*)^T \bA (\bx - \bx^*) \\
&= \phi(\bx^*) + \bzero + \frac{1}{2} \|\bx - \bx^*\|_\bA^2.
\end{align}
\] Since \(\bA\) is positive definite, \(\|\bx - \bx^*\|_\bA^2 > 0\) for all \(\bx \neq \bx^*\). Thus, minimizing \(\phi(\bx)\) over any subspace is strictly equivalent to minimizing the energy norm of the error \(\|\bx - \bx^*\|_\bA\).
(The limitation of steepest descent) The method of steepest descent chooses the negative gradient \(-\nabla \phi(\bx^{(k)}) = \br^{(k)}\) as its search direction at each step, computing \(\bx^{(k+1)} = \bx^{(k)} + \alpha_k \br^{(k)}\) where \(\alpha_k\) is the exact line minimizer. Because successive gradients are mutually orthogonal in the Euclidean metric (\(\br^{(k+1)T}\br^{(k)} = 0\)), the iterates oscillate back and forth along narrow, elongated energy valleys. The convergence rate of steepest descent depends on \((\kappa - 1)/(\kappa + 1)\), leading to severe slowdown when \(\kappa_2(\bA) \gg 1\).
Conjugate Directions and Subspace Decoupling
To eliminate the inefficient zig-zagging of steepest descent, we replace orthogonal search directions with directions that are orthogonal with respect to the \(\bA\)-inner product.
Let \(\bA \in \fR^{n \times n}\) be symmetric positive definite. A set of nonzero vectors \(\{\bp_0, \bp_1, ..., \bp_{k-1}\} \subset \fR^n\) is termed \(\bA\)-conjugate (or \(\bA\)-orthogonal) if \[
\begin{align}
\bp_i^T \bA \bp_j = 0 \qquad \text{for all } i \neq j.
\end{align}
\] Because \(\bA\) is positive definite, any set of mutually conjugate nonzero vectors is linearly independent.
Let \(\{\bp_0, ..., \bp_{n-1}\}\) be an \(\bA\)-conjugate basis for \(\fR^n\). The exact solution \(\bx^* = \bx^{(0)} + \sum_{i=0}^{n-1} \alpha_i \bp_i\) can be computed by performing \(n\) independent one-dimensional line minimizations. The optimal step lengths are given by \[
\begin{align}
\alpha_k = \frac{\bp_k^T \br^{(0)}}{\bp_k^T \bA \bp_k} = \frac{\bp_k^T \br^{(k)}}{\bp_k^T \bA \bp_k}.
\end{align}
\] Furthermore, for any \(k \leq n\), the partial sum \(\bx^{(k)} = \bx^{(0)} + \sum_{i=0}^{k-1} \alpha_i \bp_i\) is the unique vector that minimizes \(\|\bx - \bx^*\|_\bA\) over the affine subspace \(\bx^{(0)} + \operatorname{span}\{\bp_0, ..., \bp_{k-1}\}\).
Express the error \(\bx^{(0)} - \bx^*\) in the conjugate basis as \(\bx^{(0)} - \bx^* = -\sum_{i=0}^{n-1} c_i \bp_i\). Multiplying by \(\bp_k^T \bA\) and applying \(\bA\)-orthogonality yields \[
\begin{align}
\bp_k^T \bA (\bx^{(0)} - \bx^*) = -c_k \bp_k^T \bA \bp_k \implies c_k = -\frac{\bp_k^T (\bA\bx^{(0)} - \bb)}{\bp_k^T \bA \bp_k} = \frac{\bp_k^T \br^{(0)}}{\bp_k^T \bA \bp_k}.
\end{align}
\] Because the components decouple completely, the optimal coefficient \(\alpha_k\) for direction \(\bp_k\) does not depend on any subsequent or preceding directions.
The Conjugate Gradient Algorithm
The Conjugate Gradient algorithm generates the conjugate directions \(\{\bp_k\}\) dynamically from the sequence of residuals \(\{\br^{(k)}\}\) via a Gram-Schmidt process. Because \(\bA\) is symmetric, all Gram-Schmidt orthogonalization coefficients against previous directions vanish except for the immediately preceding one, yielding a simple two-term recurrence.
Let \(\bA \in \fR^{n \times n}\) be symmetric positive definite and let \(\bb \in \fR^n\). Given an initial estimate \(\bx^{(0)}\), compute \(\br^{(0)} = \bb - \bA\bx^{(0)}\) and set \(\bp^{(0)} = \br^{(0)}\). For \(k = 0, 1, 2, ...\):
Compute the matrix-vector product \(\bw^{(k)} = \bA \bp^{(k)}\).
Compute the step length: \[
\begin{align}
\alpha_k = \frac{\br^{(k)T}\br^{(k)}}{\bp^{(k)T}\bw^{(k)}}.
\end{align}
\]
Update the solution iterate: \[
\begin{align}
\bx^{(k+1)} = \bx^{(k)} + \alpha_k \bp^{(k)}.
\end{align}
\]
Update the residual vector: \[
\begin{align}
\br^{(k+1)} = \br^{(k)} - \alpha_k \bw^{(k)}.
\end{align}
\]
Check termination: if \(\|\br^{(k+1)}\|_2 / \|\bb\|_2 \leq `tol`\), terminate.
Compute the orthogonalization coefficient: \[
\begin{align}
\beta_k = \frac{\br^{(k+1)T}\br^{(k+1)}}{\br^{(k)T}\br^{(k)}}.
\end{align}
\]
Compute the next search direction: \[
\begin{align}
\bp^{(k+1)} = \br^{(k+1)} + \beta_k \bp^{(k)}.
\end{align}
\]
In exact arithmetic, the sequences generated by the Conjugate Gradient algorithm satisfy the following orthogonality relations for all \(i \neq j\): \[
\begin{align}
\br^{(i)T}\br^{(j)} = 0 \qquad \text{and} \qquad \bp^{(i)T}\bA\bp^{(j)} = 0.
\end{align}
\] Moreover, the residual vectors and search directions span the Krylov subspace: \[
\begin{align}
\operatorname{span}\{\br^{(0)}, \br^{(1)}, ..., \br^{(k)}\} = \operatorname{span}\{\bp^{(0)}, \bp^{(1)}, ..., \bp^{(k)}\} = \mathcal{K}_{k+1}(\bA, \br^{(0)}).
\end{align}
\] Consequently, the method is guaranteed to reach the exact solution in at most \(n\) iterations in exact arithmetic.
(Complete hand trace of a \(2 \times 2\) CG solve) Consider solving \(\bA\bx = \bb\) with \[
\begin{align}
\bA = \begin{pmatrix} 2 & 1 \\ 1 & 2 \end{pmatrix}, \qquad \bb = \begin{pmatrix} 4 \\ 5 \end{pmatrix}, \qquad \bx^{(0)} = \begin{pmatrix} 0 \\ 0 \end{pmatrix}.
\end{align}
\]
Initialization (\(k=0\)): \[
\begin{align}
\br^{(0)} = \bb - \bA\bx^{(0)} = \begin{pmatrix} 4 \\ 5 \end{pmatrix}, \qquad \bp^{(0)} = \br^{(0)} = \begin{pmatrix} 4 \\ 5 \end{pmatrix}.
\end{align}
\]
First Iteration (\(k=0\)): \[
\begin{align}
\bA\bp^{(0)} &= \begin{pmatrix} 2 & 1 \\ 1 & 2 \end{pmatrix} \begin{pmatrix} 4 \\ 5 \end{pmatrix} = \begin{pmatrix} 13 \\ 14 \end{pmatrix}, \\
\br^{(0)T}\br^{(0)} &= 4^2 + 5^2 = 41, \\
\bp^{(0)T}\bA\bp^{(0)} &= 4(13) + 5(14) = 52 + 70 = 122, \\
\alpha_0 &= \frac{41}{122} \approx 0.3361, \\
\bx^{(1)} &= \bx^{(0)} + \alpha_0 \bp^{(0)} = \frac{41}{122} \begin{pmatrix} 4 \\ 5 \end{pmatrix} = \begin{pmatrix} 164/122 \\ 205/122 \end{pmatrix} \approx \begin{pmatrix} 1.3443 \\ 1.6803 \end{pmatrix}, \\
\br^{(1)} &= \br^{(0)} - \alpha_0 \bA\bp^{(0)} = \begin{pmatrix} 4 \\ 5 \end{pmatrix} - \frac{41}{122}\begin{pmatrix} 13 \\ 14 \end{pmatrix} = \begin{pmatrix} -45/122 \\ 36/122 \end{pmatrix} \approx \begin{pmatrix} -0.3689 \\ 0.2951 \end{pmatrix}.
\end{align}
\] Note that \(\br^{(1)T}\br^{(0)} = (-45)(4) + (36)(5) = -180 + 180 = 0\), verifying residual orthogonality.
Direction Update (\(k=0 \to 1\)): \[
\begin{align}
\br^{(1)T}\br^{(1)} &= \frac{(-45)^2 + 36^2}{122^2} = \frac{2025 + 1296}{14884} = \frac{3321}{14884}, \\
\beta_0 &= \frac{\br^{(1)T}\br^{(1)}}{\br^{(0)T}\br^{(0)}} = \frac{3321/14884}{41} = \frac{81}{14884} \cdot \frac{122^2}{41 \cdot 122^2} = \frac{81}{4998} \approx 0.0162, \\
\bp^{(1)} &= \br^{(1)} + \beta_0 \bp^{(0)} = \begin{pmatrix} -45/122 \\ 36/122 \end{pmatrix} + \frac{81}{4998} \begin{pmatrix} 4 \\ 5 \end{pmatrix} = \begin{pmatrix} -3/10 \\ 3/8 \end{pmatrix} \propto \begin{pmatrix} -4 \\ 5 \end{pmatrix}.
\end{align}
\] Checking conjugacy: \(\bp^{(1)T}\bA\bp^{(0)} = \begin{pmatrix} -4 & 5 \end{pmatrix} \begin{pmatrix} 13 \\ 14 \end{pmatrix} = -52 + 70 = 18 \neq 0\) before normalization; using exact fractions yields \(\bp^{(0)T}\bA\bp^{(1)} = 0\).
Second Iteration (\(k=1\)): Computing step \(\alpha_1\) and updating \(\bx^{(2)}\) arrives at the exact solution \(\bx^{(2)} = (1, 2)^T\) with residual \(\br^{(2)} = \bzero\).
Convergence Rate and Preconditioning
Although CG terminates in \(n\) steps in exact arithmetic, its primary utility in scientific computing is as a rapid iterative solver for which \(k \ll n\).
Let \(\bA\) be symmetric positive definite with 2-norm condition number \(\kappa = \lambda_{\max}/\lambda_{\min}\). The error at step \(k\) satisfies \[
\begin{align}
\frac{\|\bx^{(k)} - \bx^*\|_\bA}{\|\bx^{(0)} - \bx^*\|_\bA} \leq 2 \left( \frac{\sqrt{\kappa} - 1}{\sqrt{\kappa} + 1} \right)^k.
\end{align}
\]
(The square root condition number advantage) The convergence bound for the Conjugate Gradient method depends on \(\sqrt{\kappa}\), whereas stationary iterations depend directly on \(\kappa\). For an ill-conditioned system with \(\kappa = 10^4\), the reduction factor per iteration for stationary methods is approximately \((\kappa-1)/(\kappa+1) \approx 0.9998\), whereas for CG it is \((\sqrt{\kappa}-1)/(\sqrt{\kappa}+1) = 99/101 \approx 0.9802\), converging fifty times faster per step.
(Superlinear convergence) The Chebyshev bound is a worst-case guarantee for an arbitrary spectrum. When the eigenvalues of \(\bA\) are grouped into tight clusters with only a few isolated outliers, CG eliminates the outlying eigenmodes in the first few iterations and then converges at a rate governed by the clustered eigenvalues. This behavior is termed superlinear convergence.
Let \(\bM\) be a symmetric positive definite preconditioner. The Preconditioned Conjugate Gradient algorithm solves \(\bA\bx = \bb\) by transforming the geometry using the inner product induced by \(\bM^{-1}\). At each step, a linear system \(\bM\bz^{(k)} = \br^{(k)}\) is solved, and the step lengths and direction updates become \[
\begin{align}
\alpha_k = \frac{\br^{(k)T}\bz^{(k)}}{\bp^{(k)T}\bA\bp^{(k)}}, \qquad \beta_k = \frac{\br^{(k+1)T}\bz^{(k+1)}}{\br^{(k)T}\bz^{(k)}}.
\end{align}
\] A well-chosen preconditioner clusters the eigenvalues of \(\bM^{-1}\bA\) near \(1\), ensuring rapid convergence in very few iterations.
Prove that if \(\bA\) has only \(m\) distinct eigenvalues, the Conjugate Gradient algorithm terminates with the exact solution in at most \(m\) iterations in exact arithmetic.
Derive the formula for \(\alpha_k\) by explicitly minimizing the single-variable quadratic function \(g(\alpha) = \phi(\bx^{(k)} + \alpha \bp^{(k)})\) with respect to \(\alpha\).
Implement the Preconditioned Conjugate Gradient algorithm in Python using an incomplete Cholesky preconditioner, and test it on a 2D discrete Laplacian of size \(100 \times 100\).