Numerical Linear Algebra
where the algebra stops being the arithmetic
Theme

The thread: Claims that can fail — page 4

Essays 73 to 96 of the 247 on this theme, in the same order.
the diagonal of the hat matrix, hᵢ = aᵢᵀ(AᵀA)⁻¹aᵢ · dashed: its average p/m = 0.200p/m10the leave-one-out residual: eᵢ/(1 − hᵢ), and forty refitsbars: closed form · dots: refitted without that pointone number, two fieldsΣ hᵢ, exactly p10largest leverage0.5closed form against refits8.9·10⁻¹³1 − h of the first row0.5y appears in the residualand nowhere in the leverage Least squares, and the road not to take

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.

 110¹10²01020304050noise ÷ thicknesstrials mirrored, %a coin: 50%t = 10⁻²t = 10⁻³t = 10⁻⁴per cent mirroredσ/t = 3, mean of three7.7σ/t = 10, mean of three34σ/t = 100, mean of three4720 points, 400 trials a stop, one seed per thicknessthe ratio decides, not the thinness Orthogonality, measured

A rotation that comes back mirrored

Align twenty noisy points and the nearest orthogonal matrix to the answer is a reflection in 7.7 per cent of trials at noise three times the set's thickness and a third of them at ten — at thicknesses of 10⁻², 10⁻³ and 10⁻⁴ alike. The determinant fix is never a small correction. It moves the answer by exactly 2, it costs exactly 4σ₃ of residual, and it leaves the rotation's error at half the noise however thin the set becomes.

 110¹10²10³10⁴10⁻¹⁴10⁻¹¹10⁻⁸10⁻⁵10⁻²10¹sweeprelative error, and the largest term's sizerising: the swamp's largest rank-one termfalling, slowly: its errorfalling, once: a fit with an answera plateau with a rising floorswamp error0.0014swamp term10term growth2.2benign error9.7·10⁻¹⁵benign term growth1the error alone cannot tellthe size of the terms can When the index is a tuple

An iteration that walks out of the set

Every sweep of alternating least squares is the exact minimiser of its own subproblem, so the objective can only fall. What it cannot do is converge, when the target's nearest rank-r point is not in the rank-r set — and a plateau at a small residual looks identical to slow convergence unless the size of the terms is plotted beside it.

 12345678910¹10³10⁵10⁷number of indicesrandom numbers drawna dense Gaussian sketchdashes: the tensor's own entriesa Khatri–Rao sketcha random matrix nobody can afforddense at d = 81.5·10⁸structured4032the tensor's entries1.7·10⁷dense ⁄ structured3.7·10⁴crossing at d2the sketch outgrows its tensorand the structured one does not Randomised, and the guarantee that changes kind

Sketching what is never unfolded

A range finder multiplies its matrix by a few random vectors. For a mode-k unfolding those vectors have nᵈ⁻¹ entries, so the random object is the size of the tensor divided by n — and by six indices it is larger than the tensor it is sketching.

 0481216202402468101214unknowns on the separatorcolumns above 10⁻⁸the same matrix, renumberedin the separator's own orderthe ordering the geometry hands overseparator 73separator 236renumbered, largest11the block, largest11share of the square stored0.52the fill is totaland it is not independent Sparsity, and what elimination costs

The fill that is not independent

Eliminate both halves of a grid and what is left on the separator is 100 per cent nonzero — the sparsity field's result, unchanged. Its off-diagonal block is 11 by 12 and six columns describe it to eight digits. Renumber the separator and the same block needs all eleven.

 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.

 s10 bitss20 bitss30 bitss40 bitss514 bitsinvariant factors, in bitstwo routes, and a conserved quantitydet A-2.3·10⁴Π invariants2.3·10⁴SNF widest15HNF widest24det of the transform-1an algorithm and a definitionagreeing as integers Exact arithmetic, and what it costs instead

What a determinant does not determine

Two integer matrices can have the same determinant, the same rank and the same size, and define genuinely different maps. What separates them is a list of integers each dividing the next — computed here twice, once by unimodular elimination and once from the gcds of every minor, which share no algorithm at all.

 2468101210⁻¹⁷10⁻¹⁴10⁻¹¹10⁻⁸10⁻⁵10⁻²log₁₀ τ‖Bx − d‖ / ‖d‖sketched with the objectiveweighted, not sketchedkept out of the sketchone power instead of twokept out — feasibility1.7·10⁻¹⁶sketched at τ = 10⁸1.7·10⁻⁸unsketched at τ = 10⁸1.7·10⁻¹⁶objective ÷ optimum1.2a sketch preserves a normand a constraint is not one Randomised, and the guarantee that changes kind

The half of a problem a sketch may touch

A sketch guarantees that a norm is preserved to within a factor. An equality constraint is a statement that a quantity is zero, and no multiplicative guarantee says anything about zero. Sketch a constrained problem written as a weighted one and the constraint is not destroyed — it is demoted, from a violation of 1/τ² to one of ε/τ, exactly half the exponent.

 11.31.61.92.210⁻¹⁶10⁻¹³10⁻¹⁰10⁻⁷10⁻⁴10⁻¹10²log₁₀ quadrature pointsdistance from the true counthalf an eigenvaluean integer, eventuallytrue count2at 4 points2at 128 points2finest error1.4·10⁻¹³the integral is an integerand a rounding hides how far it was Where the flop count stopped predicting the time

The last digit is the cheapest

Every cost curve measured here has the same shape: the first digits are cheap and the last ones are not. One method inverts it. Doubling the work buys twice as many digits as the previous doubling did, so the price of a digit halves every time it is paid.

 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.

 -0.0833-0.167-0.0833-0.1671-0.167-0.0833-0.167-0.0833R · A · P, normalised to a unit centrestored entries per rowlevel 0 · 31×314.87/rowlevel 1 · 15×158.22/rowlevel 2 · 7×77.37/rowlevel 3 · 3×35.44/rowlevel 4 · 1×11.00/rowstill a stencil, still annihilates a constantentries in an interior row9weight outside the 3×30row sum0the isotropic model problemnine, at every level below the first Iterating, instead of factorising

The coarse problem is a different problem

In one dimension the Galerkin coarse operator is the coarse discretisation, entry for entry — this site asserted it. In two dimensions a five-point operator produces a nine-point coarse one, so the recursion solves a different discretisation at every level below the first, and converges at 0.20 a cycle regardless.

 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.

 02244668811013215417610⁻¹⁶10⁻¹⁴10⁻¹²10⁻¹⁰10⁻⁸10⁻⁶10⁻⁴10⁻²iterationchange between iteratestwo rates, one curveα0.85λ₂(P)1predicted rate0.85measured0.84iterations176against the solve1.4·10⁻¹⁶the upper dashed line is αᵏthe curve is on the other one The matrix that is a graph

A chain with no stationary vector

A page with no outgoing links loses forty per cent of the walker's probability in six hundred steps. A directed cycle never converges at all. And on a graph whose links only run one way, the entire rank of half the vertices is exactly one minus the teleportation parameter.

 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.

 10¹10³10⁵10⁷10⁹10¹¹10¹³10⁻¹⁷10⁻¹⁴10⁻¹¹10⁻⁸10⁻⁵10⁻²10¹condition number κ(A)‖QᵀQ − I‖classical once, Householdermodified once, Householderclassical twice, Householderclassical twice, Cholesky QRκ²uκuat κ = 10⁸classical once, Householder0.0042modified once, Householder5.7·10⁻⁹classical twice, Householder3.1·10⁻¹⁵classical twice, Cholesky QR10⁻¹⁵64×16 in blocks of 4, three seedsHouseholder inside does not help between Orthogonality, measured

A stable block is not a stable basis

Block Gram–Schmidt orthogonalises twice over — between blocks, and inside each one. Householder inside the blocks does not stop the classical between-block step losing orthogonality like κ², 4.2·10⁻³ at κ = 4.3·10⁷, and a second pass does not stop Cholesky QR inside the blocks breaking down at κ = 10⁸. Each level fails only on ill-conditioning placed at its own level, and one variant holds 3·10⁻¹⁵ on every placement.

 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 same lattice, twicewhat the reduction may not changedet, before1det, after1defect, before7.1defect, after1LLL steps1the determinant is the invariantand the defect is what is being reduced Exact arithmetic, and what it costs instead

A basis that describes its lattice badly

The same set of points has infinitely many bases, they are all correct, and they are not equally useful. One measurement separates them — the product of the vectors' lengths over the lattice determinant — and the determinant is the invariant the reduction may not change, which is what makes the reduction checkable.

 0816243211.11.21.31.4terms added, each followed by a truncationerror ⁄ best rank-k erroroptimalterms with nothing in commona subspace that driftsthe rounding nobody should have feareddrifting, worst excess1independent, worst excess1a linear bound would say32energy discarded, first1.8·10⁻⁷energy discarded, last0.03thirty-two roundingsand four per cent Neither sparse nor dense

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.

 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.

 θx (frequency across x)θy0π/2π0π/2πthe coarse grid'sunder 0.2under 0.40under 0.60under 0.80under 0.95under 1.01damping per sweeptwo routessmoothing factor, scanned1closed form130×30 frequency cellsthe marker is the mode nothing removes Iterating, instead of factorising

A direction the smoother cannot see

Give the Laplacian a strong direction and multigrid stops working — from 0.2016 a cycle to 0.9565 — with every component unchanged and the condition number identical to twelve digits. The problem did not get harder. The link between the method's two halves broke.

 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.

 11.31.61.92.210⁻¹⁶10⁻¹³10⁻¹⁰10⁻⁷10⁻⁴10⁻¹10²log₁₀ quadrature pointsdistance from the true counthalf an eigenvaluean integer, eventuallytrue count2at 4 points2at 128 points2finest error1.4·10⁻¹³the integral is an integerand a rounding hides how far it was Randomised, and the guarantee that changes kind

Counting what is inside a circle

A trace of a matrix nobody wants to form, integrated around a contour, gives an integer — how many eigenvalues are inside. It converges exponentially, it is estimated with random probes, and the probe block is a ceiling that the answer does not mention.

 10⁻¹¹10⁻⁹10⁻⁷10⁻⁵10⁻³10⁻¹10⁻¹²10⁻¹⁰10⁻⁸10⁻⁶10⁻⁴10⁻²ε the blocks were compressed at‖A − LLᵀ‖ ⁄ ‖A‖the factorisationthe knob, on a factorisationtruncations34residual at 10⁻⁴2.2·10⁻⁵residual at 10⁻¹⁰9.9·10⁻¹²slope1worst ratio to the representation1.8ten decadesand a slope of one Neither sparse nor dense

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.

All themes