Numerical Linear Algebra
where the algebra stops being the arithmetic
Theme

The thread: Measured, not asserted — page 3

Essays 49 to 72 of the 110 on this theme, in the same order.
how far is it from A to an orthogonal matrix?smaller is nearer · the polar factor minimises this in every unitarily invariant normpolar factor U1.8554QR, signs fixed2.1265QR as returned3.8226200 drawn at randomκ = 10polar factor1.9QR, signs fixed2.1QR as returned3.8best of 200 random2.7‖A − QR‖ is the same either wayand ‖A − Q‖ is not Orthogonality, measured

The nearest orthogonal matrix

Every field that has to clean up a drifted rotation reaches for QR, and QR does not answer the question. The nearest orthogonal matrix is the orthogonal factor of the polar decomposition — nearer by about a tenth, and, more to the point, the same matrix whatever order the columns were written in. QR's answer changes completely.

 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.

 -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 Eigenvalues, singular values, rank

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.

 10¹10²10³10⁴10⁵10⁶10⁷10⁸10⁹10¹⁰110⁴10⁸10¹²10¹⁶condition number κ(A)× short of a double solvebf16fp16fp32bf16fp16fp32the reference is soundreference solve, worst backward error10⁻¹⁶bfloat16 threshold κ256fp32 threshold κ1.7·10⁷eight refinement steps, residual always in doubleflat at 1 means it reached double The arithmetic underneath

Where the hardware went

bfloat16 carries eight mantissa bits, which puts its refinement threshold at a condition number of 256. That is not an exotic matrix. It is an ordinary one, and past it the method still improves the answer by a factor of four hundred while getting nowhere near a usable one.

 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 Sparsity, and what elimination costs

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.

 10²020406080100120unknownsiterations2D, no preconditioner2D, preconditioned1D, preconditionedone construction, two dimensions2D steps at 16 unknowns102D steps at 100 unknowns211D steps at 100 unknowns10the same averaging, the same transformand a count that no longer stops growing Structure, and the solver that cannot see it

Two dimensions, and the cluster that thins

The same kernel, the same averaging, the same transform — applied along two axes instead of one. In one dimension the preconditioned step count is 7, 10, 10, 10; on square grids with the same unknown counts it is 10, 18, 20, 21, and the share of the spectrum near one falls from 56% to 17%.

 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.

 10²10³10⁴10⁵10⁶10⁷10⁸10⁻¹⁴10⁻¹¹10⁻⁸10⁻⁵10⁻²condition number‖QᵀQ − I‖one passsweeptwicetreewhat the second pass removesfitted slope, one pass2fitted slope, two passes0.96rounds, two passes6rounds, the sweep24one pass squares the condition numberand two do not Where the flop count stopped predicting the time

Doing it twice

Cholesky QR squares the condition number — a fitted slope of 1.95 in κ against the Householder sweep's 1.00. Run the identical routine a second time on the Q it returned and the slope is 0.93, the orthogonality is at or below the sweep's at every κ, and the price is one more all-reduce.

 10²10⁴10⁶10⁸10¹⁰110²10⁴10⁶10⁸10¹⁰condition number of the matrixcondition number seen by the iterationκ(A), unchangedκ(AR⁻¹)a bound with no κ(A) in itκ(AR⁻¹), every κ(A)2.2κ(SU), the other route2.2κ(A), across the sweep10⁸κ(AR⁻¹), across the sweep1the sketch never sees the spectrumand the spectrum cancels out of the answer Randomised, and the guarantee that changes kind

The sketch that is not the answer

Sketch-and-solve throws away the original problem and keeps the small one's answer, which is why its answer moves with the seed. Use the same sketch as a preconditioner instead and the condition number the iteration sees is the same number at every κ from a hundred to ten billion — identically the same, to nine digits, because the spectrum cancels out of it.

 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¹¹110¹10²10³10⁴condition number of the matrixgrowth factorbound 2^11partial pivotingCholeskyno pivot to gain fromCholesky growth, every κ1Cholesky interchanges0partial pivoting, worst9the bound, 2^112048both eliminations reach the same growthand only one of them had to swap to get there Elimination, and the swap

A factorisation with nothing to pivot for

Cholesky's growth factor is not bounded by one. It is equal to one, at every size and every condition number, and the two-line reason is why the algorithm needs no pivoting at all — not "usually gets away without it". Its only failure is the square root of a non-positive number, which is exactly the test for definiteness, and in floating point that test moves with the precision.

 10¹10³10⁵10⁷10⁹10¹¹10¹³10¹⁵10⁻¹⁸10⁻¹⁵10⁻¹²10⁻⁹10⁻⁶10⁻³1κ₂(A), the matrix that was updatedrelative forward errorSherman–Morrisondirect solve of A + uvᵀone answer, two routesκ of the answer's matrix1κ of the matrix replaced10·10¹³update formula's error2.5·10⁻⁴direct solve's error1.1·10⁻¹⁶the question's condition number is 1at every point on this axis Least squares, and the road not to take

A correction cheaper than the problem

Sherman and Morrison's formula updates a solved system for a rank-one change to the matrix, at 4n² operations instead of (2/3)n³. It is exact algebra. On a problem whose updated matrix is the identity — condition number one, the easiest system there is — it returns a forward error of 2.5·10⁻⁴ where a direct solve returns 10⁻¹⁶.

 159131721252910⁻¹⁶10⁻¹³10⁻¹⁰10⁻⁷10⁻⁴10⁻¹steprelative error in the orthogonal factorNewtonNewton, scaledNewton–Schulzone fixed point, three costsscaled Newton, steps7Newton–Schulz, steps28Newton at step 657scaled Newton at step 64.4·10⁻¹³a Newton step needs an inverseand a Schulz step needs two products Orthogonality, measured

An iteration that only multiplies

Newton's iteration for the polar factor needs an inverse every step. Newton–Schulz needs only matrix products — nothing that reads an entry, nothing that pivots — and it converges if and only if every singular value is below √3. At 1.73205 it converges and at 1.73206 it returns an orthogonal matrix that is not the answer, with a residual of 5·10⁻¹⁶ and nothing to say so.

 00.250.50.751-1-0.500.51mode frequency θ / πdamping factor per sweeppredicted±0.333measureda 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.

 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.

 3000000120000-21000000-0.5-1.500001.5-0.5000000-2T = ZᵀAZthe highlighted boxes each hold one conjugate pair, and no real rotation removes themthe form, and that it is one‖A − ZTZᵀ‖/‖A‖1.8·10⁻¹⁵‖ZᵀZ − I‖2.5·10⁻¹⁵worst eigenvalue error2.7·10⁻¹⁵surviving subdiagonal26×6, spectrum chosen before the matrix was builtquasi-triangular is as far as the reals go Eigenvalues, singular values, rank

The form a real matrix can reach

A real matrix with complex eigenvalues has no real triangular form, and the reason is one line — a real triangular matrix has a real diagonal, and a similarity does not move the spectrum. What it has instead is triangular except for one two-by-two block per conjugate pair, and the count is decided by the matrix rather than by where the iteration stopped.

 3456789101100.250.50.751bits in the exponent fieldeach curve as a fraction of its own maximumbfloat16fp16rangeprecisionto 617 decadesto 3.9 digitswhat the split buysbfloat16: largest number3.4·10³⁸fp16: largest number6.6·10⁴a bit of exponent doubles the rangea bit of significand adds a third of a digit The arithmetic underneath

The other half of a format

fp16 and tf32 have the same eleven significand bits and their largest numbers are 65,504 and 3.4·10³⁸. For two phases this site simulated the significand alone, so it was obliged to report them as the same format — which is a claim, and a false one.

 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%.

 10²11.251.51.7522.25processorstraffic saved, as a factorthe √4 the law promisesbreak-evenmeasureda limit is not a sizesaving at p = 640.88saving at p = 5761.4what the law promises2memory, as a factor4a loss at sixty-four processorsand 72% of the law at five hundred Where the flop count stopped predicting the time

Memory bought with messages

Holding four copies of the data instead of one is supposed to cut a matrix multiplication's communication by √4. Measured on a machine of 64 processors it costs 14% more traffic; at 576 it saves 44%, which is 72% of what the law promises. The memory is exactly four times, and that part is not asymptotic.

 D from PAPᵀ = LDLᵀ — the shaded pairs are 2×2 pivots10⁻⁶0.749······0.749·········1.2·10⁻⁶1.4······1.4·········2.1·10⁻⁶0.549······0.549·········3.6·10⁻⁶0.614······0.614·three rules, one matrix‖PAPᵀ − LDLᵀ‖, blocks5.8·10⁻¹⁷‖PAPᵀ − LDLᵀ‖, diagonal3.1·10⁻¹¹growth, blocks1.3growth, diagonal5·10⁵the zero block is what the problem saysand one rule does not need it to be nonzero Elimination, and the swap

When symmetry is not enough

The matrix [[0, 1], [1, 0]] is symmetric, nonsingular and perfectly conditioned, and there is no diagonal entry to pivot on. Every factorisation restricted to symmetric interchanges and one-by-one pivots fails on it, at any depth of searching, because every entry it could search is zero. The repair is to take two variables at once.

 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.

 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.

 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.

 15731pointsfinest gridone unknown — the recursion bottoms out in a divisionthe coarse operator, two ways‖RA_hP − A_2h‖/‖A_2h‖10⁻¹⁸unknowns / finest grid1.7cycle cost, in fine sweeps12each coarse point reaches three fine ones½, 1, ½ — and the restriction is its transpose Iterating, instead of factorising

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.

All themes