Numerical Linear Algebra
where the algebra stops being the arithmetic
Field

Iterating, instead of factorising

At scale nobody factorises, and the reason is not accuracy — it is that the factors of a sparse matrix are not sparse. What replaces elimination is a sequence of approximations, and the question changes shape: not what the residual of the factorisation is, but how fast the error falls and whether the thing you can measure tracks the thing you cannot. On the model problem every rate is known in closed form before anything runs.
0408012016020024010⁻¹⁶10⁻¹³10⁻¹⁰10⁻⁷10⁻⁴10⁻¹iteration‖e‖_A / ‖e₀‖_Ameasuredκ bound119 steps40×40, spectrum spread evenly in logbound permits 1417

The rate the condition number predicts

Conjugate gradients converge at a rate governed by the square root of the condition number. That is a bound rather than an estimate, it is provable, and it is loose enough that provisioning iterations from it wastes nine out of ten.

 05101520253010⁻¹⁷10⁻¹⁴10⁻¹¹10⁻⁸10⁻⁵10⁻²10¹iteration‖RᵀR − I‖ and ‖r‖/‖b‖step n‖RᵀR − I‖residual30×30, run for exactly n stepsexact arithmetic would end here

An orthogonalisation nobody calls one

Conjugate gradients are derived as a minimisation and behave as an orthogonalisation, which is why the finite-termination property in every textbook is not a property the method has in floating point.

 051015202530354010⁻¹¹10⁻⁹10⁻⁷10⁻⁵10⁻³10⁻¹iteration‖r‖ / ‖b‖plain CGIC(0) CGwhat the preconditioner didκ(A)48κ(L⁻¹AL⁻ᵀ)5.1‖A − LLᵀ‖/‖A‖0.0832D Laplacian, n = 100√κ ratio predicts 3.07×

Changing the condition number on purpose

Preconditioning is usually introduced as a trick that makes an iteration converge faster. It is not a trick. It is solving a different system with the same solution and a condition number chosen rather than inherited, and the new condition number is computable.

 02468101210⁻¹⁴10⁻¹²10⁻¹⁰10⁻⁸10⁻⁶10⁻⁴10⁻²1iteration‖r‖ / ‖b‖Laplaciancyclic shiftno progress at allevery eigenvalue of the shift is on the unit circleand it predicts nothing

The spectrum that predicts nothing

For a symmetric matrix the eigenvalues govern how fast an iteration converges. Drop symmetry and they stop governing anything — there is a matrix whose eigenvalues are as evenly spread as eigenvalues can be, on which GMRES makes no progress at all until the last possible step.

 015030045060075090010⁻¹³10⁻¹¹10⁻⁹10⁻⁷10⁻⁵10⁻³10⁻¹iteration‖r‖ / ‖b‖JacobiGauss–SeidelSOR ω=1.78closed form vs measuredρ Jacobi, exact0.99ρ measured0.99ρ Gauss–Seidel, exact0.981D Laplacian, n = 24ω optimal = 1.777

A rate that is known in advance

On the model problem, Jacobi contracts by cos(π/(n+1)) per step, Gauss–Seidel by its square, and optimally relaxed SOR by a number given in closed form. Three rates, all known before anything runs, and all measurable against what runs.

All essays