Numerical Linear Algebra
where the algebra stops being the arithmetic
Concept

Condition number — where it appears

15 essays name this object, across 7 fields. What follows is each of them, and the objects they name alongside it.
the problem you posedA = H10b = A·(1, 2, …, 10)the problem it answered exactlyA + δA, b + δb‖δ‖ / ‖A‖ = 2.3·10⁻¹⁷the answer you wantedx = (1, 2, …, 10), exactlythe answer you gotx̂, wrong by 2.7·10⁻⁴ relativebackward error 2.3·10⁻¹⁷forward error 2.7·10⁻⁴κ = 1.6·10¹³κ · η = 3.6·10⁻⁴, and the measured forward error is 2.7·10⁻⁴.The algorithm is not at fault. The problem is.H10, LU with partial pivotingresidual and error differ

The exact answer to a nearby problem

A good algorithm does not give an approximate answer to your problem. It gives the exact answer to a problem very close to yours — and once that is the definition, a wrong result has two possible authors and they can be measured apart.

error · backward error A = H8 · κ = 1.5·10¹⁰ · both factorisations reconstruct A to 6·10⁻¹⁷the diagonal is 1 in both — every column is a unit vector either way1.000000000001.000000000001.000000000001.000000000001.00000.002-0.002000001.0000.125-0.13300000.0020.1251.000-1.0000000-0.002-0.133-1.0001.000classical Gram–Schmidt1.000000000001.000000000001.000000000001.000000000001.000000000001.000000000001.000000000001.000Householderclassical ‖QᵀQ − I‖1.4Householder ‖QᵀQ − I‖1.4·10⁻¹⁵largest off-diagonal 1 against 3.1·10⁻¹⁶length is not angle

Orthogonal is a number

"Q is orthogonal" is a claim about a measurable quantity, ‖QᵀQ − I‖, and on the eight-by-eight Hilbert matrix two standard algorithms return 10⁻¹⁵ and 1 for it. The one that returns 1 still reconstructs the matrix perfectly, which is why nothing warns you.

orthogonality · orthogonality 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.

iterative · krylov the matrix, lower triangle408 entriesits Cholesky factor1739 entries · 1331 created→‖A − LLᵀ‖/‖A‖1.4·10⁻¹⁶fill, symbolic1331fill, numeric1331n = 144 · density 3.2% · bandwidth 12same matrix, renumberedthe answer is identical to rounding

The factor is not sparse

A sparse matrix has a factor that is not sparse, and the gap between them is the entire reason iterative methods exist. The entries elimination creates can be counted before any arithmetic runs, from the graph alone.

sparsity · fill 110²10⁴10⁶10⁸10¹⁰10¹²10¹⁴10⁻¹⁷10⁻¹⁴10⁻¹¹10⁻⁸10⁻⁵10⁻²10¹condition number κ(A)relative errorforward errorbackward errorpredicted: κ · u8×8, 20 seeds per κ; dashed is the worstthe problem worsens, not the method

A small residual is not a small error

Substituting the answer back and finding that it fits is the most natural check there is, and it verifies the wrong thing. A residual of 10⁻¹⁷ is entirely compatible with an answer whose second digit is wrong.

error · backward error 1234567891010⁻¹⁸10⁻¹⁵10⁻¹²10⁻⁹10⁻⁶10⁻³1indexsingular valuecutoff, σ₁ · 10⁻¹⁰numerical rank 10gap 8.2·10⁶an opiniontrue rank 410×10, built with 4 nonzero valuesrank is a decision

Rank is a decision

A floating-point matrix does not have a rank. It has a spectrum of singular values, and somewhere in that spectrum is a place where the values stop being signal and start being noise. Deciding where is a judgement, and the evidence for it is a gap.

spectra · rank 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.

iterative · krylov 110¹10²10³10⁴10⁵10⁶10⁷00.250.50.751amplification of the input perturbationfraction of directions at or belowκ = 10·10⁵worst found 7.6·10⁵6×6, 200 directionsmedian reaches 0.29 of κ

The condition number is an amplifier

κ is usually introduced as a definition and then quoted. It is a measurement: perturb the input by a known amount, look at how much the output moves, and the largest ratio you can find is the number.

error · conditioning 10⁻⁴10⁻³10⁻²10⁻¹110¹10²10³10⁻¹⁶10⁻¹³10⁻¹⁰10⁻⁷10⁻⁴10⁻¹relative change in the coefficients, along the worst directionrelative increase in the residualcoefficients doubled39% change, fit unmoved in the sixth digit308×: the third digit movesκ(A) = 3.6·10⁶. Exact arithmetic would pick one point on this floor. It would not raise it.24 points, degree 9, monomial basisthe data leaves them free

The valley with no bottom

A degree-nine fit's coefficients can be moved by a third of their own size before the residual changes in the sixth significant figure. The arithmetic did not lose those digits. The data never contained them.

leastsquares · fitting 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.

iterative · preconditioning 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.

iterative · gmres 10¹10²10³10⁴10⁵10⁶10⁷10⁸10⁹10¹⁰110⁴10⁸10¹²10¹⁶condition number κ(A)× short of a double solvebf16fp16fp32bf16fp16fp32the reference is soundreference solve, worst backward error10⁻¹⁶bfloat16 threshold κ256fp32 threshold κ1.7·10⁷eight refinement steps, residual always in doubleflat at 1 means it reached double

Where the hardware went

bfloat16 carries eight mantissa bits, which puts its refinement threshold at a condition number of 256. That is not an exotic matrix. It is an ordinary one, and past it the method still improves the answer by a factor of four hundred while getting nowhere near a usable one.

arithmetic · mixed precision 00.250.50.751-1-0.500.51mode frequency θ / πdamping factor per sweeppredicted±0.333measureda coarse grid seestwo routes to one factorsmoothing factor, scanned0.33smoothing factor, closed form0.33worst mode disagreement4.4·10⁻¹⁶63 interior points, one sweepthe left-hand end is what the coarse grid is for

The error smoothing cannot reach

One weighted Jacobi sweep multiplies every mode of the error by a number, and the number is a sine. Half the modes are cut by three or better, and the other half come back at 0.999 — which is not a failure of the method but the fact the whole of multigrid is built on.

iterative · multigrid 3000000120000-21000000-0.5-1.500001.5-0.5000000-2T = ZᵀAZthe highlighted boxes each hold one conjugate pair, and no real rotation removes themthe form, and that it is one‖A − ZTZᵀ‖/‖A‖1.8·10⁻¹⁵‖ZᵀZ − I‖2.5·10⁻¹⁵worst eigenvalue error2.7·10⁻¹⁵surviving subdiagonal26×6, spectrum chosen before the matrix was builtquasi-triangular is as far as the reals go

The form a real matrix can reach

A real matrix with complex eigenvalues has no real triangular form, and the reason is one line — a real triangular matrix has a real diagonal, and a similarity does not move the spectrum. What it has instead is triangular except for one two-by-two block per conjugate pair, and the count is decided by the matrix rather than by where the iteration stopped.

spectra · real schur 110¹10²10³110¹10²10³off-diagonal entry ccondition number of the eigenvalue√(1 + c²)decoupled: 1measuredthree routes, one number‖A − ZTZᵀ‖/‖A‖1.7·10⁻¹⁵closed form100computed 1/|yᵀx|100worst measured movement46four eigenvalues, two conditioning numbersthe symmetric case has one, and it is 1

A condition number for one eigenvalue

In the symmetric case every eigenvalue has condition number exactly one. In this four-by-four matrix two of them have condition number 100.005 and the other two have exactly 1, and the number belongs to the eigenvalue rather than to the matrix.

spectra · eigen conditioning

Named alongside it

The objects these essays reach for when they reach for this one.

ResidualOrthogonalitySingular valuesConjugate gradientsJacobi iterationKrylov subspacePerturbationSparsityBackward errorCholeskyFill-inForward error

All concepts