Numerical Linear Algebra
where the algebra stops being the arithmetic
Concept

Cancellation — where it appears

The loss of significance when two nearly equal numbers are subtracted, which destroys the digits that differ and keeps the ones that do not. It is not a rounding error — every operand and the operation are exact — and it is the reason two algebraically identical formulas can differ by every digit.

Named by 12 essays across 8 fields — each of them below, with the objects they name alongside it.

κ of the suma component of b − Ax1.01·10¹⁷ad − bc, near-degenerate3.6·10¹⁶qᵢᵀqⱼ, an orthogonality check7.39·10¹⁵pᵀAp, a curvature409zᵀAz, a trace probe41.8rᵀr, a residual norm1measured, not assumedhighest10¹⁷lowest1above 10¹⁰3terms64sums of squares are safeand nobody decides anything from one

The vector that hides it

Every quick demonstration of a parallel sum uses positive numbers, and positive numbers are the one family where the effect is absent. Measured on six inner products this site already computes, the summation condition number runs from exactly 1 to 10¹⁷ — and the safe end is where nobody makes a decision.

machine · Summation 1112131415193111.365129.731148.096166.461probes takenrunning estimate of the tracenormal±1one probe, no errorthe exact trace99±1 variance, this matrix0±1 variance, rotated57normal variance545the same spectrum in a general basiscosts the ±1 probe its whole advantage

Counting what cannot be looked at

The trace is n additions and one of the most expensive quantities in the subject to estimate, because the matrices whose trace is wanted are never stored. Hutchinson's estimator is unbiased with one line of algebra — and its variance depends on which random vector is used, by a factor that is a property of the matrix, and on a diagonal matrix one choice is exact from the first probe and the other is not.

randomised · Trace estimation 10⁻²10⁻¹110¹10²10³10⁴10⁻⁴10⁻³10⁻²10⁻¹1backward error, in units of ushare of systems above ituCramereliminationCramer, control24-bit arithmeticworst Cramer, in u395worst elimination, in u1.2control, worst Cramer4.9κ of the worst system3.4·10⁶one derivation, two computationsand only one of them is stable

A rule that is correct and unusable

Cramer's rule gives every component of the solution in closed form, in terms of determinants, and it is a theorem. On two-by-two systems whose rows are nearly parallel it returns an answer with a backward error of 458 units of roundoff where elimination returns 1.3 — on a matrix whose condition number is 32,000 and which elimination solved perfectly.

elimination · Determinant 10⁻⁹10⁻⁶10⁻³110³10⁶13212937κ · usignificand bitsκu = 1a bound was provedthe method refusedit never returns a wrong boundlargest κu with a proof0.45smallest κu without one0.89cases refused, of the grid13a refusal is not a wide bound — it is no bound at alland it is the only failure mode here

A bound that is proved

Every error statement on this site so far is a measurement of one run. Interval arithmetic makes a different kind of claim — the answer lies in this set, for this input, with no probability attached — and its failure mode is that it returns nothing at all. On a Hilbert system it proves a bound 23 times the error it bounds, and one size later it refuses.

arithmetic · Interval ‖H − Hₙ₋₁‖∞ ÷ σₙpositive residues2.0000mixed residues2.0000clustered poles2.0000four poles2.0000wide spectrum2.0000‖ΔH‖ 3.89·10⁻⁴ σ 1.94·10⁻⁴‖ΔH‖ 3.92·10⁻⁴ σ 1.96·10⁻⁴‖ΔH‖ 0.00502 σ 0.00251‖ΔH‖ 8.41·10⁻⁴ σ 4.21·10⁻⁴‖ΔH‖ 4.6·10⁻⁴ σ 2.3·10⁻⁴2four systems, four ratios, one number: 2.0000one state, one equalityand the sum is that, repeated

The state that is removed is not a mode

Balanced truncation removes one state and pays exactly twice one Hankel singular value. The natural reading is that the state removed is the model's least important mode and that the σ is that mode's own size — and on five systems that reading over-estimates by between 1.21 and 7.46, never once under.

reduction · Balanced truncation 10⁻¹⁷10⁻¹⁵10⁻¹³10⁻¹¹10⁻⁹10⁻⁷10⁻⁵10⁻³10⁻¹10¹10³10⁵σ — how far the periodic wrap is from singularrelative forward error10⁻¹10⁻²10⁻³10⁻⁴10⁻⁵10⁻⁶10⁻⁷10⁻⁸10⁻⁹10⁻¹⁰10⁻¹¹10⁻¹²periodic correctionκ(C)·uelimination on Tn = 64, median of five answersperiodic correction, σ = 10⁻⁸1.2·10⁻⁴elimination on T, σ = 10⁻⁸1.1·10⁻¹⁴κ(C) at σ = 10⁻⁸4·10⁸κ(T) at σ = 10⁻⁸1712the same matrix T in every columnonly the wrap it is solved through changes

The circulant the problem did not contain

A matrix that differs from a circulant in two corner entries can be solved through the circulant, by a transform and a two-by-two correction, and the cost claim is exact. The accuracy claim is not. On tridiag(−1, 2 + σ, −1), whose condition number stops at 1,712, the correction is wrong by 1.2·10⁻⁴ at σ = 10⁻⁸ while elimination is right to 1.1·10⁻¹⁴ — because the periodic neighbour is singular at σ = 0 and the two-by-two system inherits that. Solved by Cramer's rule, as here, the two amplifications multiply; a later measurement found that a pivoted solve of the same two-by-two system removes the second.

structure · Circulant -7-6-5-4-3-210⁻¹⁶10⁻¹⁴10⁻¹²10⁻¹⁰10⁻⁸10⁻⁶10⁻⁴10⁻²110²10⁴log₁₀ of the twist's distance from landing, in turnsrelative forward errorsimple zero, landed on at 0simple zero, landed on at 0.3double zero, d = 2n = 64, twist 10⁻⁵ of a turn from landingpair: error, κ(A) = 4.1·10⁶5.4·10⁻¹²single: error, κ(A) = 3.8·10⁶3·10⁻⁸double: error, κ(A) = 4.2·10¹²7515the same condition number of the wrapa different number of modes near zero

Two near-zeros cost less than one

Solve a well-conditioned tridiagonal matrix through a nearly singular wrap and the correction's accuracy is not set by how singular the wrap is. At κ = 4·10⁷ one wrap returns the answer to 8.8·10⁻¹¹ — better than κ·u — and another, at κ = 3.8·10⁷, returns it to 3.3·10⁻⁶. The difference is how many of its samples sit near the symbol's zeros. A real wrap lands on a conjugate pair, a rank-two correction absorbs the pair exactly, and its two-by-two system has condition number 1.00 — which mattered because that system was solved by Cramer's rule; solved with pivoting, the single landing costs what the pair costs.

structure · Circulant 10⁻¹²10⁻¹⁰10⁻⁸10⁻⁶10⁻⁴10⁻²10⁻¹⁴10⁻¹10¹²10²⁵10³⁸10⁵¹10⁶⁴δ, the gap between consecutive eigenvaluesrelative error in eᴬan answer with no correct digitsV f(Λ) V⁻¹scaling and squaring‖A(δ) − A₀‖exact eigenvalues throughoutκ(V) at the smallest δ3.3·10⁸²eigen route2.9·10⁶⁵scaling and squaring4.1·10⁻¹²distance to the limit7·10⁻¹²the eigenvalues are the diagonaland they are exact at every stop

The series that has to be squared back

The Taylor series for the matrix exponential is not wrong — every term is computed correctly — and on Moler and Van Loan's two-by-two its largest term is 5.4 million times the answer it sums to. The method that replaces it scales the matrix down and squares the result back, and both halves of that sentence cost: too few squarings and the approximant is out of range, too many and each one doubles the rounding.

spectra · Matrix function 10⁻¹⁷10⁻¹⁵10⁻¹³10⁻¹¹10⁻⁹10⁻⁷10⁻⁵10⁻³10⁻¹10⁻¹⁴10⁻¹²10⁻¹⁰10⁻⁸10⁻⁶10⁻⁴10⁻²1εrelative error in J(x)vforwardcentralcancellationtruncationagainst a derivative that is exactforward floor1.3·10⁻¹⁰central floor1.1·10⁻¹²truncation slope, forward1truncation slope, central2no ε reaches the roundoffand the analytic derivative is free of the choice

An operator with no entries

At the sizes where linear algebra is expensive the matrix does not exist. What exists is a subroutine that returns Av. Every Krylov method survives that unchanged; every algorithm that reads an entry disappears. And the derivative such a code computes is accurate to ten digits instead of sixteen, which turns out to cost nothing at all.

iterative · Matrix-free 1234567810⁻¹⁷10⁻¹⁴10⁻¹¹10⁻⁸10⁻⁵10⁻²10¹singular value, largest firstrelative errorone-sided Jacobizero-shift QRshifted QReigenvalues of BᵀBagainst a rational bisectionσₘᵢₙ, exactly2.1·10⁻³⁰worst, one-sided Jacobi4.4·10⁻¹⁶worst, zero-shift QR2.2·10⁻¹⁶worst, eigenvalues of BᵀB1a relative error is a ratioand the denominator is the answer

Small compared to what

This site's own singular value routine has carried a sentence since the month it was written — that one-sided Jacobi computes the small singular values to high relative accuracy and the standard method does not. It has never been measured here, because measuring it needs a σ that is known rather than computed. A bidiagonal matrix and a Sturm count in exact rationals supply one.

spectra · Relative accuracy 01428425670849811210⁻²²10⁻¹⁸10⁻¹⁴10⁻¹⁰10⁻⁶10⁻²conjugate gradient iterationrelative residualthe unit roundoff, 1.11·10⁻¹⁶the answer's residualthe residual reportedtwo residuals, one runreported, at its best6.9·10⁻²¹the answer's, at its best5.1·10⁻¹⁰unit roundoff1.1·10⁻¹⁶largest iterate on the way9.3·10¹³iterations drawn110the recurrence remembers every roundingand the stopping test is written in it

The residual the method reports

Conjugate gradients prints a relative residual of 6.9·10⁻²¹. The unit roundoff is 1.1·10⁻¹⁶, so that is not a small residual and not a large one — it is not a residual. The vector the method is holding at that step has ‖b − Ax‖/‖b‖ = 5.1·10⁻¹⁰, and nothing in the run says so.

iterative · Residual gap 01020304050607010⁻⁶10⁻³110³10⁶10⁹k‖Aᵏ/k!‖‖e^A‖ = 2.6largest term 1.4·10⁷what the series throws awaylargest term1.4·10⁷‖e^A‖2.6digits cancelled away5.4·10⁶error after the sum5.2·10⁻⁹every term is computed correctlyand the sum has lost seven digits

The error the method already knows

Summing the exponential's Taylor series throws away a known number of digits, and the number is on the machine while the sum is being formed. The largest term divided by the answer, times the unit roundoff, tracks the relative error that comes out — to within a factor of nine, across fourteen orders of magnitude of it — and nothing reports it.

spectra · Matrix function

Named alongside it

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

Condition numberUnit roundoffForward errorBackward errorCapacitance matrixCatastrophic cancellationCirculant matrixDiscrete laplacianExact ground truthFlop countKrylov subspaceMatrix exponential

All concepts