How to Learn Krylov-Subspace Algorithms: Conjugate Gradient, GMRES, Residuals and Preconditioning

Wait, What?

A matrix with a billion mostly-zero entries may be easier to solve than a much smaller dense matrix—if you stop trying to factor it.

That is the doorway into Krylov-subspace methods. For large sparse systems, the professional question is often not “Which factorisation should I compute?” but “How much can I learn about the solution by repeatedly multiplying a vector by the matrix?”

This article owns that narrower learning job. The existing Numerical Linear Algebra Algorithms article owns the broader foundation of elimination, LU, QR, SVD and conditioning. Here we move into iterative sparse solves: residuals, Krylov spaces, Conjugate Gradient, GMRES, preconditioning and convergence diagnostics.

Quick Answer

Learn Krylov methods through the route linear system → residual → matrix-vector product → Krylov subspace → projection idea → Conjugate Gradient for symmetric positive-definite systems → GMRES for general systems → restart → preconditioning → stopping criteria → finite-precision effects → sparse storage → production diagnostics. A beginner should be able to trace residual updates on a tiny system. A professional should be able to justify the solver class, choose and assess a preconditioner, interpret stagnation, and distinguish a small residual from a trustworthy solution.

1. Start With the Equation, Not the Solver Name

We want to solve Ax = b. In a direct method, we transform or factor A and then solve. In an iterative method, we begin with a guess x₀ and repeatedly improve it.

The central observable is the residual r = b − Ax. If r is zero, x satisfies the equation exactly. If r is nonzero, it tells us how far the current iterate fails to satisfy the system.

2. Separate Residual From Error

The residual is computable because A, b and the current x are known. The true error e = x* − x is usually not computable because the exact solution x* is precisely what we do not know.

This distinction matters. A small residual can coexist with a larger solution error when the system is ill-conditioned. Conditioning therefore remains part of the story even when the solver is iterative.

3. The Matrix-Vector Product Is the Primitive Operation

Krylov methods can work through repeated products Av without needing every element of A to be treated as a dense table. This is why they pair naturally with sparse matrices and with operator interfaces where A is represented by a function that computes Av.

Current SciPy sparse linear algebra and PETSc KSP both expose iterative solvers around this operator viewpoint.

4. Build the Krylov Space by Repeated Action

Starting from an initial residual r₀, the kth Krylov subspace is generated by vectors like r₀, Ar₀, A²r₀, and so on. The solver searches for an improved answer inside a space built from information obtained by repeatedly applying A.

The learning insight is simple: the algorithm is not wandering through every possible vector. It builds a structured sequence of increasingly expressive search spaces.

5. Projection Is the Unifying Idea

Many Krylov methods can be understood as choosing the best available approximation in the current Krylov space subject to a particular orthogonality or minimisation condition. Different methods differ in which geometry they exploit and which quantity they optimise.

6. Conjugate Gradient Has a Special Contract

Conjugate Gradient, usually abbreviated CG, is designed for systems whose matrix is symmetric positive definite. That assumption is not a technical footnote. It gives the method a geometry in which carefully chosen search directions do not undo earlier progress.

PETSc’s current KSPCG documentation explicitly requires the matrix and preconditioner to satisfy the relevant symmetry and definiteness conditions.

7. Trace CG on a 2×2 System Before Coding It

For a tiny symmetric positive-definite system, write down x₀, compute r₀, set the first search direction, compute the step length, update x, update r, and then generate the next conjugate direction. Do this with exact arithmetic first.

The learner should be able to explain why each scalar exists rather than memorising formulas. What is being minimised? Why is the new direction not merely the new residual? What information from the previous direction is retained?

8. Conjugacy Is Stronger Than Ordinary Orthogonality

CG search directions are conjugate with respect to A: different directions are orthogonal under the A-induced inner product. This is what allows progress along a new direction without destroying the optimum already obtained along previous conjugate directions in exact arithmetic.

9. GMRES Handles a Wider Class of Matrices

GMRES—Generalized Minimal Residual—works for general nonsymmetric linear systems. It builds an orthonormal Krylov basis, commonly through the Arnoldi process, and chooses the iterate that minimises the residual norm over the current Krylov space.

That broader applicability comes with a cost: the basis grows, orthogonalisation work grows, and memory use grows unless the method is restarted.

10. Arnoldi Is the Engine Behind Standard GMRES

Arnoldi takes a starting vector and repeatedly applies A while orthogonalising the new vector against the basis already built. The result is an orthonormal basis plus a small upper-Hessenberg representation of A on that basis.

This reduction is powerful because the enormous original problem is projected into a much smaller least-squares problem at each stage.

11. Restarted GMRES Is a Trade-Off, Not a Free Optimisation

To cap memory and orthogonalisation cost, practical GMRES often keeps only a fixed number of basis vectors before restarting from the latest iterate. Restarting controls resources but may slow or even stall convergence because useful subspace information is discarded.

A professional should therefore treat restart length as an algorithmic parameter tied to memory, spectrum and convergence—not as a cosmetic setting.

12. Preconditioning Changes the Problem You Present to the Solver

A preconditioner M is chosen so that solving with M is cheap and the transformed system is easier for the iterative method. Informally, a good preconditioner makes the hard directions less hard.

PETSc’s KSP manual treats Krylov solvers and preconditioners as a paired design decision. This is the correct professional mental model: solver choice and preconditioner choice are not independent.

13. Learn Preconditioners by Increasing Structural Knowledge

Jacobi/diagonal scaling: cheap, weak, easy to understand.
Block Jacobi: exploits small block structure.
Incomplete LU or Cholesky: approximates a direct factorisation while controlling fill.
Domain decomposition: breaks a large problem into interacting subproblems.
Multigrid: attacks errors at multiple spatial scales and can be exceptionally effective for PDE-derived systems.

The strongest preconditioner is not automatically the best. Setup time, parallelism, memory, repeated solves and robustness all matter.

14. Stopping Criteria Must Match the Question

Common stopping rules compare a residual norm with an absolute or relative tolerance. But a tolerance is meaningful only relative to scaling, data uncertainty and the downstream use of the solution.

Do not teach “residual below 10⁻⁸ means solved.” Teach “the stopping rule is evidence that a specified criterion has been met.”

15. Finite Precision Breaks the Perfect Textbook Story

In exact arithmetic, Krylov methods have elegant orthogonality and termination properties. In floating-point arithmetic, basis vectors lose orthogonality, recurrences accumulate error, and residuals maintained by recurrence may drift from explicitly recomputed residuals.

Professionals therefore monitor behaviour rather than assuming the algebraic theorem will appear unchanged on hardware.

16. Sparse Storage Is Part of the Algorithmic System

Compressed sparse row, compressed sparse column and matrix-free operators can change memory access, parallel efficiency and the practical cost of Av. Two mathematically identical solver setups may perform very differently because of representation.

17. Diagnose Convergence Curves, Not Just Final Status

Rapid early decrease followed by stagnation.
Steady geometric reduction.
Long plateaus broken by sudden improvement.
Oscillation or instability.
Fast residual convergence but unacceptable application error.
Iteration count exploding as the mesh or problem size increases.

Each pattern is evidence about conditioning, spectrum, scaling, preconditioning and stopping logic.

18. Test With Problems Whose Structure You Control

Create diagonal matrices with known eigenvalues, SPD matrices with adjustable condition number, nonsymmetric systems, nearly singular systems and sparse PDE matrices. Change one property at a time and observe iteration counts and residual histories.

19. Compare Against a Direct Solve on Small Cases

For small validation systems, compute a trusted direct solution and compare. This gives the learner a reference for actual error, not only residual. Once behaviour is understood, move to sizes where direct factorisation becomes expensive.

20. Common Learning Failure States

Using CG because it is fast without checking SPD assumptions.
Equating residual with solution error.
Ignoring matrix scaling.
Choosing GMRES restart size without measuring the effect.
Treating preconditioning as an optional afterthought.
Comparing solvers with different tolerances.
Timing only iteration time while ignoring preconditioner setup.
Assuming fewer iterations always means less wall-clock time.
Ignoring finite-precision loss of orthogonality.
Benchmarking one easy matrix and generalising the result.

21. A Beginner-to-Professional Learning Ladder

Level 1: compute Ax and a residual by hand.
Level 2: explain why sparse Av can be cheap.
Level 3: construct the first few Krylov vectors.
Level 4: trace CG on a 2×2 SPD system.
Level 5: explain Arnoldi and the GMRES least-squares step.
Level 6: implement small CG and compare with SciPy.
Level 7: vary conditioning and measure convergence.
Level 8: add and compare preconditioners.
Level 9: profile sparse storage, setup cost and solve cost.
Level 10: justify a production solver-preconditioner-stopping configuration from evidence.

22. Teach by Prediction Before Execution

Before running code, ask learners to predict which system CG should solve quickly, which matrix violates CG’s assumptions, and how a diagonal preconditioner might change convergence. Then run the experiment and reconcile prediction with evidence.

This code-reading and prediction-first progression aligns with the PRIMM approach described in current programming-education work, including the 2026 study Adapting PRIMM to a Primary Computing Setting and the original ACM work on Predict–Run–Investigate–Modify–Make.

23. Use Faded Worked Examples

First provide every CG quantity. In the next problem remove the step length. Then remove the update coefficient. Then give only A, b and x₀. Finally ask the learner to choose between CG and GMRES and justify the choice before coding.

24. Immediate, Delayed and Transfer Checks

Immediate: calculate one residual and one Krylov step.
Concept: state the assumptions that make CG appropriate.
Delayed: explain why restarted GMRES can stagnate.
Transfer: choose a solver for SPD, nonsymmetric and matrix-free cases.
Professional: interpret a convergence history and propose the next diagnostic experiment.

25. AI Assistance Boundary

AI can generate test matrices, explain trace steps, produce residual plots and suggest experiments. The learner should still be able to verify matrix assumptions, derive the meaning of the residual, explain the solver’s search space, justify preconditioning and independently interpret convergence.

Professional Direction

Advanced study includes MINRES, BiCG, BiCGSTAB, flexible GMRES, deflation, recycling Krylov spaces, block Krylov methods, multigrid preconditioning, domain decomposition, mixed precision, communication-avoiding variants, matrix-free Newton–Krylov methods and scalable distributed implementations.

Algorithm-learning rule: never ask only whether an iterative solver converged. Ask what space it searched, what assumptions made the geometry valid, what the preconditioner changed, what the stopping test actually measured, and whether the returned answer is accurate enough for the real problem.