Concept

Residual — where it appears

How far a computed answer misses satisfying the equation, which is what a solver can measure and is not the error it wants. It costs one matrix–vector product and it is what every stopping criterion is written on, so the gap between it and the error is where most silent failures live.

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

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
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
024681012141610⁻¹1index kλthe transformthe eigensolverC = F* Λ F is a factorisationworst relative disagreement2.6·10⁻¹⁵‖Cx − b‖/‖b‖ from the transform solve4.7·10⁻¹⁶imaginary part of a real spectrum1.2·10⁻¹⁶n = 16, and the whole matrix is 16 numberseigenvectors known in advance

The matrix that is one row

A circulant of size 16 is sixteen numbers, has no zero entry anywhere, and hands over its entire spectrum in closed form — the discrete Fourier transform of its first column, exactly. An eigensolver spends a sweep of Jacobi rotations over 256 entries arriving at the same answer, and agrees to 1.2·10⁻¹⁵.

structure · Circulant
everything Ax can reachb = (1.1, 0.4, 1.5)Ax, the closest reachable pointr = b − Ax‖Aᵀr‖ / (‖A‖‖r‖)1.7·10⁻¹⁶‖b‖² − ‖Ax‖² − ‖r‖²1.3·10⁻¹⁵‖r‖1.3200 random nearby points of the plane were tried; none is closer.a 3×2 system, Householder QRperpendicularity is checked

The projection and the right angle

The least-squares solution is the one whose residual is perpendicular to everything the columns can reach. That is not a mnemonic — it is an equation, Aᵀr = 0, and the computed answer satisfies it to 10⁻¹⁶.

leastsquares · Least-squares
10²10³10⁴10⁵matrix size ncountoperations, bothwords, unblockedwords, blocked (b = 6)the answer does not move‖PA − LU‖/‖A‖, unblocked2.8·10⁻¹⁶‖PA − LU‖/‖A‖, blocked2.8·10⁻¹⁶difference between them0the dashed curve is both orderings' operation countthe solid pair is what they cost

The same arithmetic at a different price

A blocked and an unblocked elimination perform 72,568 operations each — the same operations, associated differently — choose the same pivots, and return a factorisation identical to the last bit: ‖PA − LU‖/‖A‖ = 4.487946226420872·10⁻¹⁶ in both. One of them moves 41,332 words between fast and slow memory and the other moves 19,476.

cost · Blocking
21-13-3-121-212-443-12as givenrows in the order 1 2 3 443-1201.251.252.502.51.5-30-0.5-0.52after step 1pivot 443-1202.51.5-3000.5400-0.21.4after step 2pivot 2.543-1202.51.5-3000.540003after step 3pivot 0.5‖PA − LU‖/‖A‖0largest multiplier0.75row order 4 3 2 1the pivot is chosen

Elimination is a sequence of choices

Gaussian elimination is taught as a procedure with no decisions in it. There is one decision at every step — which row to use — and every stability property the algorithm has comes from making it well.

elimination · Elimination
0408012016020024010⁻¹⁶10⁻¹³10⁻¹⁰10⁻⁷10⁻⁴10⁻¹iteration‖e‖ ⁄ ‖e₀‖ in the A-normmeasuredκ 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
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
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
21018263442500102030405060708090100110120130140150160170180190200210220230240250260270280290300310320330340significand bitsiterationspreconditioner roundedarithmetic roundedno preconditionersame bits, different casualtyerror, 3-bit preconditioner8.8·10⁻¹³error, 3-bit arithmetic0.16‖A − LLᵀ‖/‖A‖ of the factor0.083a direction may be roundeda measurement may not

The part of a solver that may be rounded

A preconditioner computed and applied with a three-bit significand still returns thirteen correct digits — it costs seventeen extra iterations and nothing else. Round the working arithmetic instead and the step count barely moves while the answer loses exactly the digits the format dropped.

combination · Mixed-precision
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
the matrix43 entriestip eliminated first253 entriestip eliminated last43 entries‖A − LLᵀ‖/‖A‖, tip first1.4·10⁻¹⁶‖A − LLᵀ‖/‖A‖, tip last0dense factor is n(n+1)/2 = 253 · sparse factor is 2n − 1 = 43one row swapped to the endnothing numerical chose between them

Two ends of the same arrow

One matrix, one row moved from the front of the elimination order to the back, and the factor goes from completely dense to no fill at all. Both factorisations are exact to rounding, and nothing numerical chose between them.

sparsity · Fill
0816243240110²10⁴10⁶10⁸10¹⁰10¹²10¹⁴matrix size ngrowth factor max|u| / max|a|the 2ⁿ⁻¹ boundworst of 30 randommedian randomWilkinson's matrix sits on the bound30 Gaussian matrices per sizeat n = 40: bound 5.5·10¹¹, worst 4.8

The bound that is never attained

Partial pivoting's stability guarantee permits the entries to double at every step — a factor of 5.5·10¹¹ at n = 40. The measured growth on random matrices of that size is about three. The gap is eleven orders of magnitude, and the guarantee is still worth having.

elimination · Growth
04812162010⁻¹10⁻⁰.⁵1target rank k‖A − Aₖ‖₂published boundrandomisedσₖ₊₁, optimalhow far apart the three areworst seed spread1.1bound / median at k = 1211median / optimum at k = 12160×60, 6 seeds, oversampling p = 5band is best to worst

Randomisation does not create structure

On a matrix whose singular values are all equal, a rank-ten randomised approximation has error 1.0 — and so does the optimal deterministic one. Neither achieved anything, and only one of them is usually sold with the implication that it might.

randomised · Randomised
H13 x = b, b formed exactly so that x = (1, 2, …, 13)12345678910111213exact1.00002.00003.00133.97865.18864.993010.46720.048921.2657-2.575319.21478.906013.5113computed6.6 correct digits4.8 correct digits3.4 correct digits2.3 correct digits1.4 correct digit0.8 correct digitno correct digitsno correct digitsno correct digitsno correct digitsno correct digits0.6 correct digit1.4 correct digitbackward error2.2·10⁻¹⁷κ = 1.7·10¹⁸The algorithm solved a neighbouring problem perfectly. That problem's answer is this one.right-hand side built in BigInt rationalsthe truth is known

An answer that is known

Almost every demonstration of numerical error estimates the error by computing the same thing more carefully. The Hilbert matrix does not need that: its inverse is a closed form in integers, so the true answer is available exactly and the error is measured rather than approximated.

error · Exact ground truth
the matrix105 entriescorner first — sparsest227 entries, growth 1.9·10¹¹largest first — safe242 entries, growth 1.19the middle factor is the smaller one, and its answer has no correct digitsboth factorisations reproduce the matrix‖PA − LU‖/‖A‖, sparsest3.8·10⁻¹⁷‖PA − LU‖/‖A‖, pivoted5.4·10⁻¹⁷forward error, sparsest3·10⁻⁵forward error, pivoted4.8·10⁻¹⁶red marks are entries elimination createdthe fill argument and the stability argument disagree

Structure and stability stop being separable

The sparsest variable to eliminate on this matrix has a diagonal entry of 10⁻¹². Eliminating it produces the smaller factor, reproduces the matrix to 3.8·10⁻¹⁷ — better than pivoting does — and returns an answer wrong in the fifth digit.

sparsity · Sparse pivoting
κ 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
012345610⁻¹⁶10⁻¹³10⁻¹⁰10⁻⁷10⁻⁴10⁻¹refinement step‖x − x*‖ / ‖x*‖a full double-precision solveresidual in24-bitresidual indoubleone argument apartκ·u of the factorisation6·10⁻⁴double residual, final3.2·10⁻¹³same-precision, final1.3·10⁻⁴30×30, κ = 10⁴, same factors in both runsidentical cost

Buying the accuracy back

Factorise in single precision, then correct the answer using residuals computed in double, and the result is what a full double-precision solve would have given. Compute those residuals in single instead and the identical algorithm, at identical cost, recovers nothing.

arithmetic · Mixed-precision
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
00.250.50.75110⁻¹110¹share of the noise placed in the matrixleast-squares error ÷ total least-squares errorequally accuratetotal leastsquares aheadordinary leastsquares aheadthe model, not the methodadvantage, all noise in b0.28advantage, all noise in A2.5seeds at each share40the same total noise at every pointand only where it sits changes

When the matrix is wrong too

Every least-squares problem here has assumed A is exact and b is not, and moved b onto the column space of A. Where both were measured, the smallest correction that makes the system consistent moves the matrix as well — and on the problems where that answer is more accurate, it has the larger residual, by construction rather than by luck.

leastsquares · Total least-squares
-0.16-0.63-0.493.43-0.025-0.632.3-0.69-0.95-1.70.076-0.49-0.694.3-1.6-1.41.93.4-0.95-1.65.40.0141.63-1.7-1.40.0144.60.055-0.0250.0761.91.60.0553.1A, symmetric→-0.164.600004.65.92.500002.51-1.20000-1.24.90.3800000.385.80.2400000.242H = QᵀAQ, tridiagonal‖A − QHQᵀ‖/‖A‖1.1·10⁻¹⁵below the subdiagonal0worst eigenvalue movement7.1·10⁻¹⁵a similarity, so the spectrum is untouched — and every later step is O(n²) rather than O(n³)one reduction, then every iteration is cheapthe eigenvalues did not move

The form that makes it affordable

One Householder reduction, done once, turns every subsequent iteration of the eigenvalue algorithm from cubic to quadratic cost. It changes no answer at all, which is why it is easy to describe as an optimisation and wrong to.

spectra · The QR algorithm
10⁻³10⁻²10⁻¹1110¹10²10³10⁴pivot threshold τgrowth factor · entries in L+U, ÷ entries in Agrowthfillthe library default‖PA − LU‖/‖A‖ at τ = 0.12.9·10⁻¹⁶growth at τ = 0.138entries in L+U at τ = 0.1372one knob, two measurements, opposite directionsand the default is most of both

A threshold between fill and growth

One number decides how small a pivot an elimination will accept. At 0.001 the factor holds 172 entries and the matrix grows by 1,330; at 1 it holds 260 and grows by 1.2. The libraries ship 0.1, and the measurement says why.

sparsity · Sparse pivoting
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.

iterative · Stationary
00.250.50.751-1-0.500.51mode frequency θ / πdamping factor per sweeppredictedmeasured±0.333a 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

The same problem on a coarser grid

Restriction, the coarse operator and interpolation are three matrices with nine distinct entries between them. Two of the three are each other's transpose, and their product with the fine operator is the coarse discretisation exactly — not approximately, entry for entry, at every level.

iterative · Multigrid

The problem the solver was actually given

A linearisation is exact — it has the polynomial's eigenvalues, with their multiplicities, and the whole loss is arithmetic. A nonlinear eigenvalue problem does not offer that. Every algorithm replaces the function first, and the term that replacement contributes is committed before any number is rounded and appears in no residual.

polynomial · Approximation before linearisation

Influence is decided before the data

The diagonal of the hat matrix sums to the number of columns and the response appears nowhere in it, so a fit has exactly p units of influence to hand out among m observations. The same row at h = 0.5 is a ten-fold outlier on one design and a boundary case on another, and which of those it is was settled before a single measurement was taken.

leastsquares · Leverage

The gap refinement can close

Multiplying by a computed inverse is not backward stable, and refinement at the working precision repairs it. That much is settled. The claim beside it — that the forward error does not move — was read at one conditioning and four corrections too late. Swept over ten, it moves at every one, and it lands on the LU route's own number after a single correction.

elimination · Inversion

The rounding that was not the problem

A rank-k block plus a rank-k block is a rank-2k block, exactly, so every arithmetic in this format truncates after every addition. A Cholesky performed inside it does ninety-eight of those and its residual is 1.14·10⁻⁹ against a representation error of 1.40·10⁻⁹ — the roundings cost nothing measurable.

hierarchy · Recompression

The two numbers a caller has

Choosing between the two least-squares methods is a statement about where the noise is, and the two quantities a caller can compute are both blind to it. The residual separates the answers by 0.14 per cent where their accuracies differ by 14, and κ(A) falls from 3.54 to 2.46 across a sweep in which the error rises by a factor of sixty-two.

leastsquares · Total least-squares

A knob calibrated in residuals

A formatted Cholesky has two numbers in it and only one of them is an accuracy. Across twelve trees — three sizes by four leaf sizes — the leaf moves the truncation count from 0 to 258 and moves the ranks of the blocks not at all, while the residual follows the tolerance at slopes between 1.022 and 1.046 and sits at about a tenth of it throughout.

hierarchy · Recompression

A basis built from the points

A polynomial fit computed in monomials and in an orthogonal basis gives the same curve on exact data, and the valley essay drew the two lying on top of each other. Add 0.1% noise and they separate — by 1.7·10⁻⁵ at degree 40 and 0.004 at degree 48 — because the fitted curve moves with the basis by its condition number times the rounding times the residual. Chebyshev polynomials keep that small only on points spread like their weight; on a sample with a hole in it they reach κ = 1.55·10⁷. A basis orthogonalised against the sample points themselves stays at 1 on every set.

leastsquares · Fitting

The eigenvalues that are answers to nothing

A rational approximant of degree five turns a six-by-six problem into a thirty-six-by-thirty-six one, and thirty-six numbers come back. Six are the answer. The rest are exact eigenvalues of the approximant, lying where the function it approximates is not a real number at all.

polynomial · Spurious spectrum

Five indices are cheaper than two

The same 4,096 unknowns cost 1.049·10⁶ multiplications indexed as a 64 × 64 grid and 1.966·10⁵ indexed as six axes of four. The dense factorisation that ignores the indexing costs 4.581·10¹⁰ at every one of them, and the residual improves in the same direction as the cost.

tensor · Kronecker

The count that is not the budget

A Cholesky performed inside a low-rank format truncates 0, 2, 10, 34 and 98 times as the leaf falls from 128 to 8, and those five integers are the same at every accuracy from 10⁻¹² to 10⁻². Across all ten decades the factorisation's residual stays below the representation's own error at a ratio between 0.81 and 1.00 — with two entries that read 1.83 and 1.78, and neither of them is accumulation.

hierarchy · Recompression

A stopping test is a race

One matrix, one right-hand side, one tolerance, thirteen partition counts — and eleven different iteration counts between 674 and 690. Every run converged, every answer is right to the accuracy asked for, and what differs is the bill.

machine · Stopping test

An accuracy that is a backward error

Every backward error on this site is something an algorithm produced and somebody then measured. This one is a line in the program. Solving with a compressed matrix gives a residual that is the compression's own error, at a slope of 1.000 over ten decades, so the knob that sets the storage sets the backward error directly.

error · Backward error

A stopping rule that follows the run it is given

A preconditioner that reaches the answer four times sooner leaves four steps within 10% of its best instead of sixteen, and a rule that stops by the residual ought to miss so narrow a window more often. Over forty draws of the noise it misses it less: the discrepancy principle stops at 1.030 times the preconditioned run's best against 1.073 times the plain run's. And past the edge it stops within a factor of 1.7 of a run whose own best is 5.5 times Tikhonov's — faithful to a run that has already failed.

combination · Iterative regularisation

One minus a leverage is a subtraction

Every deletion diagnostic divides by 1 − h, and computing it as one minus a computed leverage loses digits in proportion to 1/(1 − h), however accurate the leverage. The complementary block of a QR factor gives the same number as a sum of squares and loses κ(A)·u instead: every digit on a well-conditioned design, and half the digits the subtraction loses on a design whose far point is what made 1 − h small.

leastsquares · Leverage

Bracketing an error nobody can measure

The error of a low-rank Gramian factor is the one quantity a caller cannot compute, because computing it needs the Gramian the factor exists to avoid forming. Two numbers that can be computed sit either side of it — a rational factor known before the run, and a residual known after — and they stay a factor of four apart across a fourfold change of size.

reduction · Gramian decay

Two observations that hide each other

Two observations at the same place, wrong by the same amount, each look harmless when deleted alone, because a fit without one still has the other. Single deletion sees the shared error cut by (1 − 2h)/(1 − h) — measured at 261 times at the far end — and only the pair's two-by-two block of the hat matrix says what the two of them hold.

leastsquares · Leverage

Feasible and wrong

A third constraint that nearly repeats the first takes the best route's answer from 2.96·10⁻¹⁵ to 1.16·10⁻⁴, and the other two routes to no correct digit at all. Every one of those answers satisfies every constraint to 10⁻¹⁵. The quantity a caller checks after a constrained solve is the one quantity here that says nothing.

leastsquares · Constrained least-squares

One line that buys a quarter of the run

The adaptive forcing rule has a floor on it that no published statement of the rule carries: do not solve a step to an accuracy the outer loop will not use. Removing it costs 9 to 27 per cent of the whole inner run. Keeping it costs between 23 and 2,600 times the forward error — accuracy the residual test never asked for and both runs satisfy the test either way. The line is a trade between a residual and an error, and which of the two the caller meant decides whether it is a saving.

sequence · Inexact newton

Three errors and one number

This site's identity has two factors and a division of blame between them. Two fields have now added a third party and a fourth, and only one of the four is a property of anything — the others are decisions, made before the arithmetic, reported by nothing.

error · Backward error

The fifth author

Four authors of a wrong answer have been named on this site and each is a statement about one computation. The fifth is not: it is what separates two computations that are both correct, it is a backward error of measurable size, and no residual, bound or condition number contains it.

error · Backward error

The factor a sparse code keeps anyway

Every deletion diagnostic divides by one minus a leverage, and computing it as a subtraction loses a digit for every decade the leverage is from one. The route that does not subtract needs the orthogonal factor, which a sparse factorisation is supposed not to have. Three repairs that avoid it all fail at exactly a unit of roundoff over the divisor — and the fourth, which reaches the orthogonal factor through the Householder vectors a sparse code keeps in order to solve anything at all, returns the same bits as a stored factor in 900 operations.

leastsquares · Leverage

A corner the penalty can afford

Every smooth reading of the deconvolution's grid needed about forty points and then stopped improving, and the step was the suspect. Give the step one coefficient of its own and forty-eight points reach an error of 0.0070 at 0.1% noise, against 0.118 for the best smooth reading on ninety-six — the step was most of the error. But the same step given two coefficients recovers half as well, and given a doubled node at each edge it recovers worse than no breakpoint at all, while representing the signal to 0.07%. What decides is what the penalty is charged for the corner, and whether the data can say where it is.

regularisation · Regularisation

Nine steps of pessimism

A proved bound is 8 to 26 times the error it bounds, at every precision from 16 to 40 significand bits. A carried interval is (√2)ᵐ times too wide after m re-enclosures. The two cross between eight and nine, so the method everybody warns against is the tighter of the two for a short computation.

arithmetic · Interval

The reading that never moves

Across thirty runs — five grids from 36 to 196 unknowns, six working precisions from 53 significand bits down to 8 — the residual conjugate gradients stops on stays between 1.10·10⁻¹³ and 9.95·10⁻¹³. Over the same thirty runs the error of the answer spans a factor of 2.39·10¹¹, and the step count more than doubles. The one number the run publishes is the only one that responds to neither axis.

machine · Stopping test

A problem with no answer

If two matrices share a null vector then det(A − λB) is identically zero and every λ is an eigenvalue, which means none of them is. Perturb such a pencil by a ten-billionth and a solver returns six numbers with residuals below 10⁻⁹. Change the seed and it returns six different numbers, spread over forty-four, with residuals just as small.

spectra · Pencil

The zero that means it is finished

Every Krylov method ends by dividing by a number the previous step produced, and when that number is zero the recurrence stops. In Arnoldi the stop is the answer — the subspace has closed, the solution is inside it, and the residual is at the unit roundoff. The literature calls it a lucky breakdown, and the adjective is doing real work.

iterative · Breakdown

The same zero, and nothing was found

Change the recurrence by two lines and the divisor stops being a norm. It becomes an inner product of two vectors from two different sequences, and an inner product of two different vectors is zero on a whole hyperplane — with neither vector anywhere near zero, nothing invariant, and nothing converged. The arithmetic event is identical and the meaning is opposite.

iterative · Breakdown

An iterate that must be made smaller

Applying a Kronecker-sum operator to a low-rank iterate multiplies its ranks by d and adding two of them adds their ranks, so a solver in a compressed format cannot keep what it produces. Every step is followed by a truncation — and whether that truncation is a floor on the residual depends on the right-hand side rather than on the truncation.

iterative · Low-rank iteration

A run that is over at step five

A conjugate gradient whose every iterate is cut to a rank budget reaches the floor that budget allows at step 5, 36, 42 or 59, and then does nothing for the rest of the run. Four times the iterations move the floor by a factor of 1.8, and past the answer's own rank they move it the wrong way.

iterative · Low-rank iteration

The same budget, spent five ways

A restarted method has one budget — products with A — and two ways to spend it, in many short cycles or a few long ones. At about a hundred and forty products the answer is the same to a factor of seven whichever split is chosen, and the residual bound the method reports spans ten orders of magnitude across the same five runs.

spectra · Lanczos

Named alongside it

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

Condition numberBackward errorForward errorConjugate gradientsKrylov subspaceStopping criterionExact ground truthLeast-squaresLow-rank approximationUnit roundoffNormal equationsOrthogonal projection

All concepts