Conjugate gradient method

For a positive-definite matrix $A$, the conjugate gradient method chooses mutually $A$-conjugate search directions. Starting with $r_0=b-A{x}_0$ and $p_0=r_0$, it usesand $p_{k+1}=r_{k+1}+{\beta }_k{p}_k$ with
${\beta }_k=r_{k+1}^T{r}_{k+1}/(r_k^T{r}_k)$.

For a symmetric positive-definite matrix, choose an orthonormal eigenbasis and write the exact solution as $x={\sum }_i{\alpha }_i{v}_i$. Starting at zero, the conjugate gradient method chooses an admissible polynomial $p(0)=1$ of degree at most the iteration count that minimizes the displayed error. The weights depend on both the eigenvalues and the solution coefficients. For two eigenvalues $1,L$, the first polynomial is $1-\eta t$, with $\eta =({\alpha }^2+{\beta }^2{L}^2)/({\alpha }^2+{\beta }^2{L}^3)$. Comparable nonzero coefficients and large $L$ give $\eta \approx 1/L$, so the high-eigenvalue error is removed first. A sufficiently tiny high-mode coefficient changes this behavior; the spectrum alone does not determine the first step.

Nonzero vectors are conjugate directions with respect to a symmetric positive-definite matrix $A$ when they are pairwise orthogonal in the inner product $⟨u,v{⟩}_A=u^{T}Av$. They need not be orthogonal in the ordinary Euclidean inner product. They are linearly independent: multiplying a vanishing linear combination by $p_j^{T}A$ isolates its $j$th coefficient times $p_j^{T}A{p}_j>0$.

Let $p_1,\dots ,p_n$ be conjugate directions for a symmetric positive-definite matrix $A$. Starting from any $x_1$, minimize $x^{T}Ax/2-b^{T}x$ successively along those full lines, stopping early if the linear-system residual $r_j=b-A{x}_j$ vanishes. The quadratic minimization along a line uses ${\alpha }_j=p_j^T{r}_j/(p_j^{T}A{p}_j)$ and gives $r_{j+1}=r_j-{\alpha }_{j}A{p}_j$.

This makes $p_j^T{r}_{j+1}=0$. For $i<j$, conjugacy preserves $p_i^T{r}_{j+1}=p_i^T{r}_j$. By mathematical induction, $r_{j+1}$ is orthogonal to all preceding search vectors. The $n$ conjugate directions are linearly independent and hence form a basis, so $r_{n+1}=0$. Thus the method reaches the unique solution of the linear system in at most $n$ steps in exact arithmetic. The conjugate gradient method generates such directions without storing a complete prescribed basis in advance.

For a symmetric positive-definite matrix $H$, the Conjugate gradient method minimizes the error in the induced weighted norm over its affine Krylov subspace. Writing the error as $p(H)e_0$ givesFor ${\lambda }_{\min }<{\lambda }_{\max }$, choose the scaled Chebyshev polynomial $p(\lambda )=T_k(({\lambda }_{\max }+{\lambda }_{\min }-2\lambda )/({\lambda }_{\max }-{\lambda }_{\min }))/T_k(({\lambda }_{\max }+{\lambda }_{\min })/({\lambda }_{\max }-{\lambda }_{\min }))$. On the spectral interval its numerator has absolute value at most one. With $\rho =(\sqrt{\kappa }+1)/(\sqrt{\kappa }-1)$, the denominator equals $({\rho }^k+{\rho }^{-k})/2\ge {\rho }^k/2$. This proves the displayed bound. If the two endpoints coincide, $H$ is a positive scalar multiple of the identity and one iteration solves the system exactly.

The residual of an approximate solution $x$ of $Ax=b$ is $r=b-Ax$. It is the unbalanced right-hand side, while the solution error is $x_{\ast }-x=A^{-1}r$. The conjugate gradient method makes residuals orthogonal to successive Krylov subspaces.

The norm of a linear-system residual, used to measure how closely an approximate solution satisfies its linear system. It differs from the solution-error norm $\Vert x_{\ast }-x\Vert$. Their relationship depends on the matrix inverse and the condition number, so a small residual alone is not a scale-free relative-error guarantee.

At iteration $m$ of the conjugate gradient method,and $r_m$ itself lies in ${\mathcal{K}}_{m+1}(A,r_0)$. These Galerkin conditions imply finite termination once the Krylov subspace becomes invariant under $A$.

A symmetric preconditioner transforms $Ax=b$ intoThe aim is to make the transformed positive-definite matrix have clustered eigenvalues while keeping application of $P$ inexpensive.

For a positive-definite matrix $A$, choosing $P=A^{-1/2}$ gives $PA{P}^T=I$, so conjugate gradients terminates in one step. Constructing or applying the exact inverse square root is generally as expensive as solving the original system.

The order-$k$ Krylov subspace generated by $A$ and $b$ is

With zero initial guess, the generalized minimal residual method minimizes the Euclidean norm of the linear-system residual over each Krylov subspace. The Arnoldi decomposition reduces this to ${\min }_y\Vert {\tilde{H}}_{k}y-\Vert b{\Vert }_2{e}_1{\Vert }_2$, followed by $x_k=Q_{k}y$. Givens rotations solve the small problem in $O(k^2)$ operations once the decomposition is available. For an invertible matrix $A$, the same iterate minimizes the error in the weighted norm induced by $A^{T}A$, because $\Vert Ax-b{\Vert }_2=\Vert x-x_{\ast }{\Vert }_{A^{T}A}$.

One step of restarted GMRES searches along the current residual. In the displayed nonsingular system, the first residual minimization is ${\min }_{\alpha }[(1-\alpha {)}^2+(1+\alpha {)}^2]$, whose minimizer is zero. Restarting then repeats the same problem forever and never moves. Unrestarted GMRES instead solves exactly at step two, since $b,Ab$ span the entire space. Thus a smaller search-space cost can come at the price of stagnation; nonsingularity alone does not guarantee progress for one-step restarting.

Starting with $q_1=b/\Vert b{\Vert }_2$, the Arnoldi iteration constructs an orthonormal basis of successive Krylov subspaces. At step $j$, orthogonalize $A{q}_j$ against $q_1,\dots ,q_j$, and normalize the remainder as $q_{j+1}$ if it is nonzero. The coefficients form an upper Hessenberg matrix ${\tilde{H}}_k$ and satisfy $A{Q}_k=Q_{k+1}{\tilde{H}}_k$. A zero remainder means the current Krylov subspace is invariant under $A$.

For a real symmetric matrix $A$, the leading block $Q_k^{T}A{Q}_k$ of the Arnoldi decomposition is both symmetric and an upper Hessenberg matrix, hence a tridiagonal matrix. The resulting three-term recurrence is the Lanczos algorithm. If $A$ is a symmetric positive-definite matrix, its compression is too, since $y^T{Q}_k^{T}A{Q}_{k}y=(Q_{k}y{)}^{T}A(Q_{k}y)>0$ for nonzero $y$.

The factorization produced by the Arnoldi iteration, with orthonormal $Q_{k+1}$ and an upper Hessenberg matrix ${\tilde{H}}_k$. Multiplying by $Q_k^T$ identifies its leading square block as $H_k=Q_k^{T}A{Q}_k$.

For a diagonalizable finite-dimensional operator, write $v={\sum }_j{u}_j$ with nonzero $u_j$ in distinct eigenspaces of eigenvalues ${\lambda }_j$. Then $A^{r}v={\sum }_j{\lambda }_j^r{u}_j$. The Vandermonde matrix of the first $k$ powers is invertible for $k$ distinct eigenvalues, so these powers span exactly the $k$ independent components $u_j$. Thus the stabilized Krylov subspace dimension is $k$. It counts distinct spectral components, not the number of nonzero coefficients in an arbitrary eigenbasis: $A=I$, $v=e_1+e_2$ has Krylov dimension one. Equivalently it is the minimum number of eigenvectors needed when those vectors may be chosen freely.

In exact arithmetic, conjugate gradients reaches the exact solution in at most the number of distinct eigenvalues of $A$, and therefore in at most $n$ steps for an $n$-dimensional system.