Numerical Linear Algebra
where the algebra stops being the arithmetic
Theme

The thread: Two routes to a number — page 3

Essays 49 to 72 of the 205 on this theme, in the same order.
024681010⁻¹¹10⁻⁹10⁻⁷10⁻⁵10⁻³10⁻¹rank kept in every moderelative errordashes above: √(Σ tail²), the upper bounddashes below: max tail, a floor under the bestsolid: what the projection returnssmooth: pinned to the upper boundrank 10 error1.1·10⁻¹¹its upper bound1.1·10⁻¹¹the lower bound6.3·10⁻¹²error ⁄ bound1error ⁄ lower1.7inside the boundand sitting on it When the index is a tuple

A decomposition made only of SVDs

Everything the definition of tensor rank loses comes back if the SVD's algorithm is carried across instead of its definition — take the leading left singular subspace of every unfolding and project onto all of them. It exists, it costs d matrix decompositions, and its error is within √d of the best there is.

 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 Iterating, instead of factorising

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.

 192123252729313335373941434500.10.20.30.40.50.60.7position jweight on the true signal at jthe blur's row, 5.89 widethe kernel, 2.82 widecomputed from A and λ alonekernel width at half height2.8components kept, Σfₖ28width × Σfₖ / n1.2deepest negative lobe-0.075no data and no truth went into this curvethe answer is the truth seen through it Regularisation, and the answer that is chosen

A second blur, narrower than the first

A regularised answer is not the truth with the noise taken out. It is the truth seen through a second blur, V F Vᵀ, which depends on the operator and λ and on nothing that was measured. At the best λ for 0.1% noise its rows are 2.82 points wide against the instrument's 5.89, they dip to −0.075 on either side, and their width times the number of components kept stays between 1.10n and 1.27n across seven decades of λ. Two spikes four points apart come back as two; three apart, as one.

 -33-28.2505-23.5009-18.7514-14.0018-9.25229-4.502750eigenvalueQ(μ) ≺ 0one Cholesky, two answersabove the certificate8below it8the gap0.67critical β for this n5.8the spectrum is real by classnot by outcome The eigenvalue problem that is not linear

Every eigenvalue real, and a test that says so

A quadratic eigenvalue problem has no reason to have real eigenvalues. One class does, as a property rather than an outcome, and the proof is a Cholesky that completes. The boundary of the class has a closed form, and at the boundary the arithmetic loses half its digits with nothing ill conditioned anywhere.

 1591317212529333710⁻¹110¹10²bidiagonalisation stepsrelative errorleast without: 20no penaltypenalty insidewhat stopping is worthbest without a penalty0.14and at step 40161best with one0.14and at step 400.14the same floor, reached twiceand only one run stays on it Methods that were designed apart

The step that stops mattering

Regularise the problem the iteration has built rather than the problem it was given, and the error curve stops turning. The unregularised run ends 1,127 times above its own best; the same run with a penalty inside it ends 1.000000000003 times above.

 0246810110²10⁴10⁶10⁸10¹⁰10¹²spread of the row units (decades)condition numberκ∞(DA)cond(DA)Hilbert κ∞Hilbert condone system, two numbersκ∞ at no spread9.8κ∞ at 10 decades1.9·10¹⁰cond, either end7Hilbert, equilibrated1.3·10¹⁰the solution is the same at every spreadand one of these curves knows it Two errors, and whose fault they are

The units the matrix is measured in

One linear system, written twice. The rows of the second are the rows of the first in different units, the solution is identical to the last bit, and the condition number has moved by eight orders of magnitude. One of those two numbers is a fact about the problem and the other is a fact about the notation.

 10⁻⁶10⁻⁵10⁻⁴10⁻³10⁻²10⁻¹10112233445566relative drift ‖E‖ / ‖A₀‖preconditioned conjugate gradient iterationson the small eigenvalueson the large onessame drift, two placessmall end at 10⁻²38large end at 10⁻²9κ(M⁻¹A), small end36κ(M⁻¹A), large end1.3‖E‖/‖A₀‖, both0.01how much the matrix changedis not what the preconditioner cares about When the problem arrives again

Where the drift lands

The standing rule for when a preconditioner has gone stale is to rebuild it once the matrix has changed by more than some fraction of itself. Two drifts of exactly the same relative size cost 19 iterations and 5 on the same matrix, and the quantity that separates them is not in the rule at all — the perturbation is divided by the eigenvalue it lands on.

 01234567810⁻⁵10⁻⁴10⁻³10⁻²10⁻¹110¹μ, the entry above the diagonalsep, and the eigenvalue gapmin |λᵢ + μⱼ|sep(A, B)the spectra never moveeigenvalue gap, throughout2sep at μ = 02sep at μ = 89.5·10⁻⁴amplification there1056solvability is the eigenvaluesand conditioning is not Structure, and the solver that cannot see it

An equation whose unknown is a matrix

AX + XB = C is linear in X, so it has a coefficient matrix, and writing it down is the obvious thing to do. At n = 100 that matrix has a hundred million entries for a problem with ten thousand unknowns, and the algorithm everybody uses instead never forms it. Its conditioning is not the eigenvalue gap either, which is the number a reader is invited to consult.

 27121722273237424710⁻²10⁻¹110¹subspace size kvalueno answer below 8GCV's traceλ from GCVλ from the residuala denominator that is not the problem'strace at k = 41trace at k = 4827λ range across the run2.2subspace before an answer8the divisor moves by twenty-sevenand the answer does not move Methods that were designed apart

A parameter chosen on a smaller problem

Inside a hybrid method the regularisation parameter is chosen on a 25×24 problem rather than a 64×64 one. The rule that reads a residual transfers exactly; the rule that reads a trace is biased by exactly two grid steps at twenty-four steps and one at forty, at every noise level from 10% to 0.1%.

 110¹10²10³10⁴10⁵10⁶10⁷10⁻¹⁷10⁻¹⁶10⁻¹⁵10⁻¹⁴10⁻¹³10⁻¹²10⁻¹¹10⁻¹⁰10⁻⁹10⁻⁸1/(1 − h), the leverage of the removed row‖R̄ᵀR̄ − (G − aaᵀ)‖ / ‖G − aaᵀ‖downdatedrefactorisedat h = 1 − 10⁻⁷κ of the downdated matrix4.3κ of the matrix downdated9.3·10⁶rotation's amplification344downdate residual3.5·10⁻¹⁰a hyperbolic rotation is not orthogonaland that is exactly what it is for Least squares, and the road not to take

The observation that cannot be removed

Removing a rank-one term from a Cholesky factor needs a rotation that is not orthogonal, and the number under its square root is 1 − h, where h is the leverage of the row being removed. The algorithm's breakdown condition and the statistician's warning are the same quantity, arrived at from opposite ends, and neither field states it in the other's language.

 012345610⁻¹⁶10⁻¹³10⁻¹⁰10⁻⁷10⁻⁴10⁻¹round‖AΩ‖ ⁄ ‖A‖ at that roundredrawn each roundhalf redrawnkept: the sketch is now the zero matrixthe randomness is spentkept, round 12.6kept, round 29·10⁻¹⁵redrawn, round 23.1kept, final error0.28redrawn, final error0.11a random matrix used twiceis not random the second time Randomised, and the guarantee that changes kind

The sketch that is spent

Every other object a sequence carries has a shelf life. A random sketch has one use. Deflate what its first round found and apply it again, and it returns the zero matrix — 9.0·10⁻¹⁵ where the first round saw 2.62 — because the input has been made orthogonal to the very draw the guarantee is over.

 the bound211no pivoting127τ = 0.00198τ = 0.003100τ = 0.01106τ = 0.03111τ = 0.1116τ = 0.3117τ = 1138entries in Ua bound, and its slackthe bound, from the graph alone211the worst that occurs138loose by1.5no arithmetic was done to compute the bound — only the pattern of AᵀA and its elimination graphGeorge and Ng: U fits inside chol(AᵀA), whatever the swapsprovable, cheap, and loose Sparsity, and what elimination costs

What the symbolic phase can only bound

Without pivoting, the fill can be computed from the graph and the count is exact — 233 predicted, 233 measured. With pivoting it is 233 predicted and 242 measured, and what survives is a bound that is right at every threshold and loose by 1.7 times at the largest grid drawn.

 01234567891011121310⁻¹110¹10²10³iterationinterpolation point σa fixed point that is a conditionorder4iterations13σ against −λ(Aᵣ)2.8·10⁻¹²H₂ error4.3·10⁻⁵it stopped movingand the condition holds there Reduction, and what a model is for

Interpolating at the model’s own poles

One choice of interpolation points is not arbitrary — the mirrored poles of the model about to be built. It is a fixed point rather than a guess, and when it is reached it beats a method costing O(n³) — by 0.4 per cent, which is the honest size of the whole contest.

 02468101210⁻¹⁷10⁻¹⁴10⁻¹¹10⁻⁸binary64 operations per elementrelative errorruntime orderfixed partitionpre-rounded, one passpre-rounded, two passescompensated (Kahan)exact accumulationwhat determinism costsreproducible policies3cheapest, ops/element3its error2.3·10⁻⁷two passes10⁻¹⁷ordinary, 8 pieces2.7·10⁻¹¹nothing is cheap and exactthe middle of the figure is the answer The answer that depends on the machine

What determinism costs

Six ways to add up a vector, priced in operations per element and in accuracy. Nothing sits in the bottom left of the figure — an answer that is the same on every machine costs between three and twelve operations where an answer that is not costs one.

 2^18modulus, as a power of twoa lattice with one short vectornumerator355denominator1132·max(n, d)², bits18first recovered at18moduli tried42below the bound there are two answersand the algorithm cannot prefer one Exact arithmetic, and what it costs instead

A fraction recovered from one remainder

A solution over the rationals can be computed modulo a prime power and then recovered — the residue determines the fraction uniquely, but only once the modulus is twice the square of the fraction's longer part. Below that there is no partial credit: the algorithm returns a different fraction with the same residue, and it is a perfectly good one.

 10⁻¹⁴10⁻¹³10⁻¹²10⁻¹¹10⁻¹⁰10⁻⁹10⁻¹⁵10⁻¹²10⁻⁹10⁻⁶10⁻³1relative size of the entrywise perturbationrelative forward errorκ∞ · εcond(A,x) · εmeasuredboth bounds holdκ∞(A)1.9·10⁸cond(A, x)4.8ratio of the bounds4·10⁷both curves above the data are boundsand only one of them is a measurement Two errors, and whose fault they are

A condition number scaling cannot move

Skeel's componentwise condition number is invariant under any row scaling — exactly, before any norm is taken, because two diagonal factors cancel entry by entry. It is never larger than the normwise one and can be arbitrarily smaller, and the ratio between them is a diagnostic for which kind of ill-conditioning a matrix has.

 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 Iterating, instead of factorising

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.

 10²10³10⁴10⁻¹⁷10⁻¹⁶10⁻¹⁵10⁻¹⁴10⁻¹³10⁻¹²steps taken‖RᵀR − AᵀA‖ ⁄ ‖AᵀA‖the bound, linear in the steps√k · uevery step safe, the chain notdrift after the run3.9·10⁻¹⁴the bound there3.3·10⁻¹³√k · u there6.1·10⁻¹⁵worst single amplification2.7worst leverage met0.69refreshes1backward stable onceand three thousand times is a different claim When the problem arrives again

Stable once, and three thousand times

A sliding window adds a row and removes one at every step and never looks at the data again. No single step of it amplifies by more than 2.72, no downdate fails, and after three thousand steps the triangular factor in memory is 3.9·10⁻¹⁴ from the matrix it is supposed to be a factor of — six hundred times growth from a per-step bound that says nothing about chains.

 0.46-1.11.5-0.470.2-0.181.90.690.42-0.160.028-0.0160-1.50.160.880.099-0.1100.680.762.30.930.0480-0.250.110.952.11.500001.52.3after reflector 2the highlighted entries are the bulge — the only thing that is not Hessenbergentries below the subdiagonal3reflectors used2reflectors in a whole step5the shifts are never formed — only their sum and productand both of those are real Eigenvalues, singular values, rank

Two shifts that are never formed

The double shift is defined as a factorisation of (A − μI)(A − μ̄I), which nobody computes. What is computed is the first column of that product — three numbers — and the bulge those three numbers create, pushed down the subdiagonal by n − 2 reflectors until it falls off the bottom.

 -2-101234-5-3-1135real partimaginary part4 insidea countable spectruminside the contour4drawn12existing∞worst branch residual1.6·10⁻¹⁵there is no last eigenvalueso the question has to change The eigenvalue problem that is not linear

A problem with infinitely many eigenvalues

Let the matrix depend on λ through something that is not a polynomial and three things stop being true at once. There is no linearisation, there is no characteristic polynomial, and "compute the spectrum" is not a request that can be granted — the only finite question is how many eigenvalues are inside this circle.

 01428425670849810⁻¹⁶10⁻¹⁴10⁻¹²10⁻¹⁰10⁻⁸10⁻⁶10⁻⁴10⁻²iterationchange between iteratestwo rates, one curveα0.85λ₂(P)0.9predicted rate0.77measured0.77iterations108against the solve10⁻¹⁶the upper dashed line is αᵏthe curve is on the other one The matrix that is a graph

A ranking that is an eigenvector

PageRank is the stationary vector of a walk that follows links with probability α and jumps at random otherwise. The iteration and the elimination agree to 4·10⁻¹⁷. What α is set to changes which pages come third, fourth and fifth.

 polar: U₁P − U₂QR: Q₁P − Q₂-3.3·10⁻¹⁶-3.3·10⁻¹⁶2.2·10⁻¹⁶-1.9·10⁻¹⁶5.6·10⁻¹⁷-1.7·10⁻¹⁶6.7·10⁻¹⁶-8.9·10⁻¹⁶2.8·10⁻¹⁶-5.6·10⁻¹⁶-3.6·10⁻¹⁶6.1·10⁻¹⁶-8.3·10⁻¹⁶8.3·10⁻¹⁶-7.2·10⁻¹⁶-10·10⁻¹⁶2.2·10⁻¹⁶4.4·10⁻¹⁶1.7·10⁻¹⁶-4.4·10⁻¹⁶-5.6·10⁻¹⁷5.6·10⁻¹⁷0-7.2·10⁻¹⁶-4.4·10⁻¹⁶3.9·10⁻¹⁶5.6·10⁻¹⁷2.8·10⁻¹⁶-1.1·10⁻¹⁶-3.3·10⁻¹⁶2.2·10⁻¹⁶1.1·10⁻¹⁶-2.8·10⁻¹⁶5.6·10⁻¹⁶-2.8·10⁻¹⁶5.6·10⁻¹⁶0.67-0.35-0.23-0.620.19-0.72-0.770.170.20.290.190.93-0.180.55-0.0290.06-0.460.11-0.17-0.61-0.230.38-0.55-0.380.750.870.661.1-0.36-0.530.42-0.670.180.180.52-0.3Frobenius norms, columns reordered‖U₁P − U₂‖2.8·10⁻¹⁵‖Q₁P − Q₂‖3‖Q‖, for scale2.4κ of the matrix100Frobenius distance between the two answers, on one scale0 to 4polar2.8·10⁻¹⁵QR3.048‖Q‖ = 2.449the column space did not moveand one of the two answers did Orthogonality, measured

A test with no answer in it

A caller with no reference answer can still ask whether a routine answered the right question: reverse the columns, run it again, compare. The polar factor's two answers agree to 10⁻¹⁵ at every conditioning drawn; a QR's differ by 2.353 on matrices whose own norm is 2.449. The test has a floor, and the floor is measurable too.

 26101418222630343810⁻¹⁷10⁻¹⁴10⁻¹¹10⁻⁸10⁻⁵10⁻²10¹stepsizeleast error: 20‖QᵀQ − I‖filter disagreementan identity with an expiry datedisagreement at step 1210⁻¹³disagreement at step 400.014least error at step20exact while the basis is orthogonaland false where the method is best Methods that were designed apart

An expiry date the noise does not move

The polynomial description of conjugate gradients leaves the level of rounding at step 17 or 18 on this operator, at every noise level from 10% to 0.1%. The step worth stopping at moves from 3 to 44 across the same range. They coincide at about 1% noise, which is where the coincidence was first read, and it is a fact about the noise rather than about the method.

 -5-3-11350eigenvalueHZᵀHZK4 negative — Cholesky of H stops at row 30 negative — a minimum on the constraint(10, 4, 0) = In(ZᵀHZ) + (4, 4, 0)one factorisation, no Zpositive, LDLᵀ of K10negative, LDLᵀ of K4negative in H4negative in ZᵀHZ0the count follows the reduced Hessiannot the Hessian The matrix a constraint makes

A minimum the Hessian cannot see

A Hessian with four negative eigenvalues can sit at a constrained minimum, and a Cholesky of it stops at the third row. One symmetric indefinite factorisation of the saddle-point matrix settles the question anyway — ten positive pivots and four negative — without a basis for the null space ever being formed. The count is exact in the algebra and blind in floating point, in a band that grows like κ(A)²; the route through the null space is blind in one that grows like κ(A).

All themes