Numerical Linear Algebra
where the algebra stops being the arithmetic
Theme

The thread: Whose fault is it — page 4

Essays 73 to 96 of the 267 on this theme, in the same order.
-9-7-5-3-1110⁵10⁶10⁷log₁₀ ε of the preconditionermultiplications for the whole solveno preconditioner: 41 stepsiteration count, rescaledmultiplications spenttwo curves, opposite waysκ21steps, no preconditioner41cheapest ε0.5its rank1tightest ⁄ cheapest2.4the count is what is printedand the cost is what is spent Neither sparse nor dense

The accuracy worth paying for

Used as a preconditioner, a hierarchical representation gets better at every accuracy — the iteration count falls monotonically all the way to the tightest tolerance. The total work does not. Its minimum sits at a rank-one preconditioner on an easy problem and six decades further along on a hard one.

 110¹10²10³10⁻¹⁸10⁻¹⁴10⁻¹⁰10⁻⁶10⁻²10²10⁶numbers that describe the matrixcondition number, and backward errordensetoeplitzsymmetricρ aloneno such problemcondition numberbackward errorone matrix, four descriptionsκ, all n² entries7.9·10⁴as one number, ρ62backward error, unconstrained2·10⁻¹⁷as a symmetric Toeplitz matrix5.6·10⁻¹²fewer numbers, better conditionedand no nearby problem left Structure, and the solver that cannot see it

The condition number of the model

Describe a 40×40 Toeplitz matrix by its 1,600 entries and its condition number is 78,800. Describe it by the one number it actually contains and the condition number is 61.9. The three decades in between are not an approximation or a bound — they are what κ has been over-stating, and the drop is not where the linear algebra is.

 10²10⁴10⁶10⁸10¹⁰10¹²10¹⁴10⁻¹⁸10⁻¹⁵10⁻¹²10⁻⁹10⁻⁶10⁻³1κ₂(A)relative errorforward, A⁻¹bforward, A\bbackward, A⁻¹bbackward, A\bat κ = 10¹⁴η, LU solve2.2·10⁻¹⁷η, via the inverse4.5·10⁻⁵forward, LU solve2.8·10⁻⁴forward, via the inverse0.015one factorisation, two ways to use itand one of them forfeits the backward error Elimination, and the swap

The inverse that is never formed

x = A⁻¹b is how the solution of a linear system is written and it is not how it is computed. The usual reason given is cost — three times the arithmetic. The real reason is that one of the two routes is backward stable and the other is not, and at κ = 10¹⁴ they differ by twelve orders of magnitude in the number that says whose fault a wrong answer is.

 does this matrix look nearly singular?green: the test agrees with the truth · red: it does not · the bar under each number is its magnitude, over sixty-two decades|det A||det A|^(1/n)σₘᵢₙ1/κ = σₘᵢₙ/σₘₐₓ0.1·I at n = 40perfectly conditioned10⁻⁴⁰0.10.11κ = 10¹⁰, |det| = 1nearly singular1110·10⁻⁶10·10⁻¹¹Hilbert at n = 8nearly singular2.7·10⁻³³8.5·10⁻⁵1.1·10⁻¹⁰6.6·10⁻¹¹the two counterexamplesκ of the scaled identity1its determinant10⁻⁴⁰κ of the normalised matrix10¹⁰its determinant1det(cA) = cⁿ det(A)so a determinant carries the units n times over Two errors, and whose fault they are

The number that decides nothing

The determinant is the first scalar anybody attaches to a matrix and the last one worth consulting. A tenth of the identity has a determinant of 10⁻⁶⁰ and a condition number of exactly one. The Hilbert matrix's determinant stops being right at n = 13 and stops being a number at n = 29, and nothing in between reports either.

 -1-0.75-0.5-0.25010⁻¹⁷10⁻¹⁴10⁻¹¹10⁻⁸10⁻⁵10⁻²10¹10⁴log₁₀ σ — the barrier's reduction factorresidual after one reused stepconvergedthe pattern free, the factors notentries moved6off-diagonal0survived at σ = 0.996survived at σ = 0.106 steps0 stepsthe few entries that movedare the ones that dominate When the problem arrives again

What survives one step of the barrier

An interior-point method solves the same system dozens of times with the same pattern and different numbers, and exactly p entries change between one step and the next. The pattern is reusable for ever. The factorisation is reusable for none of them, and the threshold that says so is a reduction factor of about a per cent against schedules that use ten.

 10⁻⁹10⁻⁷10⁻⁵10⁻³10⁻⁸10⁻⁵10⁻²gap between the two eigenvalueshow far it movedthe eigenvectorsthe eigenvaluestheir plane‖E‖ / gapone perturbation, three answerseigenvalue shift, spread over the sweep1plane angle, spread over the sweep1eigenvector angle, spread1.6·10⁵the dashed line is Davis–Kahan's ‖E‖/gaptwo of the three never noticed Eigenvalues, singular values, rank

The gap decides the eigenvector

A symmetric matrix's eigenvalues move by at most the size of the perturbation, whatever the spectrum looks like. Its eigenvectors are governed by a completely different quantity — the distance to the neighbouring eigenvalue — and at a gap of 10⁻⁹ the same perturbation turns them through 27°.

 015304560759010512010⁻²10⁻¹110¹steprelative sizeleast error: 20discrepancy stop: 7errorresidualthe knob is an integerleast error, at step20error there0.14error at step 1206the residual falls at every stepthe error turns and keeps rising Two errors, and whose fault they are

The zero you are allowed to write

A deflation criterion sets a subdiagonal entry to zero because it is small. A drop tolerance discards an entry of a factor because it is small. A truncation discards a singular value because it is small. Three fields, three vocabularies, no shared arithmetic — and plotted as work saved against error accepted, one curve.

 1357910⁻¹⁶10⁻¹³10⁻¹⁰10⁻⁷10⁻⁴log₁₀ b, the couplingrelative errorsmall half, directstructure-preservingsmall half, by 1/λone divisionlarge half, direct5.5·10⁻¹⁵small half, direct10⁻⁷small half, by 1/λ4.9·10⁻¹⁵structure-preserving6.5·10⁻⁶the structure is not decorationit is where half the accuracy is Structure, and the solver that cannot see it

A perturbation that keeps the symmetry

The smallest perturbation that makes a computed answer exact is the backward error. Ask for the smallest one that also keeps the problem's structure and the number can only go up — and measured on a palindromic quadratic it goes up by 1.17, while the structure the computed spectrum has lost is not in either number.

 110¹10⁻⁷10⁻⁵10⁻³distance from the branch point, λ + c|g − r|the eigenvaluesdegree 98 poleschosen, not computedreach0.06left end, from the cut0.2rational, worst8.5·10⁻⁴polynomial, worst0.0068linearisation size, both54committed before the solveand invisible to it The eigenvalue problem that is not linear

An error committed before the arithmetic

Before a nonlinear eigenvalue problem is solved, somebody says where they think the eigenvalues are. That sentence sets the accuracy of everything that follows by five orders, costs nothing to say, and cannot be revised once the approximation built on it is in hand.

 -12-11-10-9-8024681012log₁₀ of how nearly dependent the columns areverdicts that disagreed, of 40724103which one is correctmatrices tested200verdicts disagreed26fused was right9unfused was right17lower part: thefused build was righta sign has no last digitso a verdict has nowhere to hide The answer that depends on the machine

A matrix that is definite on one machine

Two hundred Gram matrices, two conforming builds, and twenty-six of them get different answers to "is this positive definite". The exact verdict, from determinants in BigInt rationals, says the fused build is right nine times and the other one seventeen.

 10⁻¹110¹10²-45-35-25-15-551525s, the smallest of the three interpolation pointsrightmost pole of the reduced modelunstable above this linebalanced truncation, order 3exact, and unusablefull system's pole-3.5placements swept19unstable models4worst pole24their interpolation1.3·10⁻¹⁴balanced truncation-0.84the conditions all holdand the model cannot be run Reduction, and what a model is for

A model that cannot be run

A stable system, reduced by matching its transfer function at three points exactly, comes back with a pole in the right half plane at four of nineteen placements — and matches at all three points to 1.3·10⁻¹⁴ while doing it. The construction did what it promised.

 01234510⁻¹⁸10⁻¹⁵10⁻¹²10⁻⁹10⁻⁶10⁻³1correction steprelative errorLU route: η = 2.2·10⁻¹⁷LU route: forward 2.8·10⁻⁴forward errorbackward errorwhat a correction buysη before refinement4.5·10⁻⁵η after four steps2.8·10⁻¹⁷forward, unchanged7.1·10⁻⁴cost of a step, flops1800the residual is repairableand the accuracy floor is the problem's Elimination, and the swap

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.

 0246810121416182001020304050member of the sequencepreconditioned conjugate gradient iterationskept from the first memberrebuilt every memberone operator, two policieskept, first member18kept, last member52rebuilt, first member18rebuilt, last member13ε at the last member0.033the problem got easierand the kept preconditioner got worse at it When the problem arrives again

The penalty for keeping it is a ratio

A kept incomplete Cholesky costs 40 iterations against a rebuilt one's 10 on 64 unknowns, and 55 against 17 on 256. Across six grids the difference between the two rises by 27 per cent and the ratio between them falls by 19. Neither quantity is free of the problem's size, and the one a policy is paid in is the one that transfers worse.

 conjugate gradientsbest step20best error0.1410% window, last/first6.5Landweberbest step1778best error0.1410% window, last/first901110¹10²10³10⁴10⁵10⁶10⁷10⁻¹110¹10²matrix–vector productsrelative errorCGLS best: step 20Landweber best: step 1,778CGLSLandweberthe same answer at two pricesand a window three orders wide Methods that were designed apart

A step that is not a unit of work

Landweber's iteration reaches conjugate gradients' best answer on the same deconvolution — 0.1414 against 0.1426 — at step 1,778 instead of step 20, and at 0.1% noise at step 56,234 instead of 44. Each step costs the same two products. And within 10% of its best it runs from step 7 to step 6,310, where conjugate gradients runs from 4 to 26: the slow method is the one that forgives a late stop.

 0246810121410⁻¹⁶10⁻¹³10⁻¹⁰10⁻⁷10⁻⁴10⁻¹interior-point iterationrelative error, and μcertified from iterate 1μthe iterate's errorthe crossover's errorone solve, checkediterations15first certified iterate1iterate error there0.22crossover error there4·10⁻¹⁴the active set arrives long before the digitsand a crossover collects them at once The matrix a constraint makes

The active set before the digits

An interior-point method takes fifteen iterations on a quadratic programme with forty constraints, and its iterate has eight correct digits at the eleventh. Take the constraints its diagonal calls active at the first iterate, solve the equality problem they define once, and check the answer against the conditions for optimality. It passes, to thirteen digits. The step's matrix had a condition number of 43 at that iterate, and 7·10¹⁵ at the last.

 the upper pair is distance from the truth; the lower pair is ‖Ax − b‖least squares · error0.08523total least squares · error0.037least squares · ‖Ax − b‖3.965total least squares · ‖Ax − b‖4.078two orderingserror ratio (ls ÷ tls)2.3residual ratio (tls ÷ ls)1seeds40no vector makes the residual smallernot even the one the problem was built from Least squares, and the road not to take

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.

 the invariant planesolid: beforedashed: aftersame perturbation, two questionsthe vectors turned, radians0.029the plane turned, radians7.6·10⁻⁸what left the plane5.6·10⁻⁸drawn in the unperturbed plane's own basisa radius is not determined; the circle is Eigenvalues, singular values, rank

The plane survives what its vectors do not

At a gap of 10⁻⁹ a perturbation of 10⁻⁶ turns the two eigenvectors through half a radian and turns the plane they span through 7.6·10⁻⁸ — a ratio of six million. Ask for the subspace instead of the vectors and a hopeless computation becomes a well-conditioned one, with no change to the arithmetic.

 10⁻¹⁶10⁻¹⁴10⁻¹²10⁻¹⁰10⁻⁸10⁻⁶10⁻⁴10⁻²110⁻⁷10⁻⁵10⁻³10⁻¹10¹relative error in the datarelative error in the answerthe answer's errorthe data's errornothing here is the arithmetic'sdata error10⁻¹²answer error10·10⁻⁴amplification10·10⁸rounding committed0residual0the residual is exactly zeroand the answer is wrong anyway Exact arithmetic, and what it costs instead

An exact answer to a measured problem

The residual is the zero vector, nothing was rounded at any step, and the answer is wrong in its first digit. Data accurate to fourteen places, an exact solve of the system it defines, and an error of 10⁻⁵ — because conditioning was never a statement about arithmetic and removing the arithmetic error removes none of it.

 1110¹10²10³10⁴10⁵the noise level it is told ÷ the true oneerror ÷ the oracle's0.40.50.71.523twice the oraclecliff, ρ = 0.69told the truthworst of 400middle 80%median draweach draw against its own oraclecliff, median ρ0.69told the truth, median1.1told a half, median3201shaded: the middle 80% of drawsbelow the cliff the error has doubled Regularisation, and the answer that is chosen

Thirty-two coefficients instead of a noise level

The discrepancy principle has to be told the noise, and told too little it does not degrade — it falls off a cliff, at 0.80 of the truth when the noise is 10% and at 0.58 when it is 0.001%, exactly where the understatement forces the filter past its best truncation. The missing number is in the data. The root mean square of the last thirty-two coefficients never sends the rule over the cliff at or below 1% noise in four hundred draws, where eight coefficients with the same median do so thirty-five times.

 -1-0.500.511.520eigenvalue of P⁻¹K1 − φ1φS never formedγ1smallest ν0.056furthest from φ, 1 − φ0.57MINRES steps11an exact Schur approximation, bought by changing Hthe golden ratio without S The matrix a constraint makes

Where the augmentation puts the cost

Add γAᵀA to the objective block of a saddle-point system and its Schur complement tends to I/γ, so the cheapest possible approximation becomes the right one and the golden-ratio spectrum arrives — within 7.6·10⁻⁶ at γ = 10⁶. MINRES falls from 21 steps to 6. The inner solve with the augmented block rises from 14 conjugate gradient steps to 43, their product does not fall at all, and the answer loses seven and a half digits on the way.

 0510152025303540012345matrix size ngrowth factorpartial pivotingrook pivotingcomplete pivotingsolid: median of 30 · dashed: worst of themat n = 40partial pivoting, median3.3rook pivoting, median2complete pivoting, median1.6rook, worst of the draw2.8the same matrices under every rulerook 3.3× partial's search Elimination, and the swap

A pivot that searches one row and one column

Rook pivoting looks down a column for its largest entry, along that entry's row for a larger one, and back down that entry's column, until it finds an entry largest in both. On Gaussian matrices of size 64 it keeps the median growth factor at 2.53 against partial pivoting's 4.06 and complete pivoting's 1.88, and it compares 6,987 entries against 2,080 and 89,440. On Wilkinson's matrix it holds the growth at exactly 2 where partial pivoting reaches 9.2·10¹⁸. And on a matrix built to make it walk it compares 113,376 entries — more than complete pivoting.

 08162432404810⁻¹10²10⁵10⁸10¹¹10¹⁴10¹⁷10²⁰degreeκ of the design matrix1/u: past this a double holds nothingmonomialsChebyshevArnoldi on the pointsκ at degree 48monomials1.5·10¹⁷Chebyshev1.6·10⁷Arnoldi on the points1200 two intervals with a gap between themthree bases for one space of polynomials Least squares, and the road not to take

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.

 10²10³10⁻¹⁶10⁻¹³10⁻¹⁰10⁻⁷10⁻⁴10⁻¹steps takenlargest relative error in a coefficientcarried factorrecomputedcarried + one correction+ a second correctionHouseholder QR, same rowsdrift of the factorafter 2875 stepsκ(A) of the window2·10⁴carried factor2.8·10⁻⁸recomputed5.8·10⁻⁹carried + one correction7.1·10⁻¹³+ a second correction3.7·10⁻¹²Householder QR, same rows1.1·10⁻¹²measured against the coefficients in exact rationalsthe refresh repairs the factor, not the answer When the problem arrives again

The repair the drift did not need

A sliding window's carried Cholesky factor drifts 3.9·10⁻¹⁴ from its data, and multiplying by κ(AᵀA) predicts eight lost digits in the coefficients, a stream conditioned at 10¹² losing the answer, and a periodic refresh of the factor as the default repair. Measured against coefficients computed exactly in rationals, all three come out differently. On a stream made ill-conditioned by scaling, the conditioning never reaches the coefficients. On a collinear stream, a freshly recomputed factor is as wrong as the drifted one. And one correction from the window's own rows reaches Householder's accuracy for a fraction of a refresh's cost.

 significand bits10³10⁴10⁵10⁶10⁷10⁸10⁹10¹⁰κ(A)87892634456329625911109435761199127514162024every matrix positive definiteruns producing a false certificate14of runs in total72never above, in significand bits12first κ at eight bits10⁵the comparison was correctand what it proved was not true Two errors, and whose fault they are

Deciding that a zero has arrived

The previous tolerances were offers — accept this much error, save this much work. A detection threshold is not an offer, because both directions are failures. One matrix here has three genuinely near-invariant subspaces, and the constant somebody typed decides which of them the recurrence stops at; at eight significand bits the same kind of constant produces a proof of something false.

All themes