Numerical Linear Algebra
where the algebra stops being the arithmetic
Theme

The thread: Identical algebra, different arithmetic — page 2

Essays 25 to 48 of the 99 on this theme, in the same order.
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.

 -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¹⁰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.

 2345678110²10⁴10⁶10⁸10¹⁰vectors in the basisκ₂ of the basisas derivedsolves spread outone subspace, two spanning setseight moments at one point7.7·10⁹eight points, spread1growth per vector661/u4.5·10¹⁵the same subspaceand only one of them usable Reduction, and what a model is for

A basis that is the same subspace and not the same thing

The interpolation conditions are conditions on a subspace, so any basis of it will do. The one a derivation writes down reaches a condition number of 7.7·10⁹ by its eighth vector, and the rate at which it gets there is set by a number the user chose with no information.

 distinct answersone accumulator303eight pieces72compensated119pre-rounded1exact1worst error 1.17·10⁻⁸worst error 4.61·10⁻⁹worst error 4.96·10⁻¹⁰worst error 4.91·10⁻⁶worst error 0bitwise, or not at allpermutations400one accumulator303pre-rounded1its error4.9·10⁻⁶compensated error5·10⁻¹⁰accuracy and agreement are different propertiesand the accurate one is not the agreed one The answer that depends on the machine

The sum that cannot be wrong

Snap every addend to a common multiple before adding, and every partial sum is exact — so the order stops mattering, by construction rather than by luck. Four hundred permutations return one value where an ordinary reduction returns three hundred and three.

 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.

 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.

 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.

 -5-3.20611-1.412210.3816782.175573.969460log₁₀ |λ||λ| = 1λλ′ = 1, in the algebrapairs6decades of spectrum9.6pairing error, general3.3·10⁻⁹pairing error, structured2.2·10⁻¹⁶a symmetry the solver never knew aboutand the half of the answer it decides The eigenvalue problem that is not linear

A spectrum that comes in reciprocal pairs

A palindromic quadratic reads the same backwards, so λ is an eigenvalue exactly when 1/λ is. A general solver discards that, computes the large half of the spectrum perfectly and the small half to seven digits — and the small half is a division away from being perfect too.

 smooth97.7%hilbert94.5%wave26.6%noise0.4%share of the core's energy on its 6 superdiagonal entrieskeeping only them: 0.152 against 7.94·10⁻⁶keeping only them: 0.235 against 1.77·10⁻⁶keeping only them: 0.857 against 1.34·10⁻¹⁵keeping only them: 0.998 against 0.842orthogonal, and not diagonalsmooth on-diagonal0.98hilbert on-diagonal0.94wave on-diagonal0.27noise on-diagonal0.0044worst slice pair3.5·10⁻¹⁶the slices are orthogonalthe core is not diagonal When the index is a tuple

The orthogonality that cannot be diagonal

A matrix decomposition hands over orthonormal factors and a diagonal middle at once. For three indices the two come apart, and there is no arrangement that has both — so the question stops being which decomposition to use and becomes which of the two properties the computation needs.

 κ(A) = 10⁶ throughout · κ(H) = 100 · the answer is the same answer for every basisorthonormal — κ(Z)1κ(ZᵀHZ)25.6relative error1.07·10⁻¹⁵first m basic — κ(Z)1.99·10⁸κ(ZᵀHZ)3.8·10¹⁶relative error0.0518pivoted basic — κ(Z)2.06κ(ZᵀHZ)31.9relative error6.71·10⁻¹⁶what the choice costsdensity, orthonormal1density, fundamental0.5κ(ZᵀHZ) ÷ κ(Z)², naive0.96error, pivoted choice6.7·10⁻¹⁶every one of them is a basisand one of them loses fourteen digits Orthogonality, measured

The basis nobody chose on purpose

A method that eliminates a constraint has to pick a basis for its null space, and every basis is correct. Their condition numbers are eight orders apart, the reduced problem inherits the square, and the choice is usually made by a one-line rule nobody thought of as a numerical decision.

 -47-37-27-17-7313233301234log₁₀ of the vector's normfp1611 bitsbfloat168 bitstf3211 bitsbinary3224 bitspale: the format's range · blue: √(Σ(xᵢ/m)²)·m · red: √(Σxᵢ²)fp16 and tf32 have the same eleven significand bitsand their bars do not overlap The arithmetic underneath

A norm that overflows before it is a norm

The vector of sixteen thousands has a Euclidean norm of 4,000, which fp16 represents exactly. Written as the square root of the sum of squares it returns infinity, because squaring doubles the exponent — and the expression costs half the format's range on the one computation every iterative method performs at every step.

 10⁻²10⁻¹110¹10²10³10⁴10⁻⁴10⁻³10⁻²10⁻¹1backward error, in units of ushare of systems above ituCramereliminationCramer, control24-bit arithmeticworst Cramer, in u395worst elimination, in u1.2control, worst Cramer4.9κ of the worst system3.4·10⁶one derivation, two computationsand only one of them is stable Elimination, and the swap

A rule that is correct and unusable

Cramer's rule gives every component of the solution in closed form, in terms of determinants, and it is a theorem. On two-by-two systems whose rows are nearly parallel it returns an answer with a backward error of 458 units of roundoff where elimination returns 1.3 — on a matrix whose condition number is 32,000 and which elimination solved perfectly.

 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.

 0246810121410⁻¹⁷10⁻¹⁴10⁻¹¹10⁻⁸10⁻⁵10⁻²10¹10⁴10⁷10¹⁰log₁₀ τ — the weight on the constraintrelative error against the exact answerτ = 1/√uGram–Schmidtnormal equationsHouseholder QRthe ceiling is the method'sHouseholder, τ = 10¹⁴4.8·10⁻¹⁵normal equations14Gram–Schmidt4.9·10¹⁰1/√u6.7·10⁷a constraint is a weight at infinityand the solver decides how far infinity is Least squares, and the road not to take

A constraint is a weight at infinity

Stack an equality constraint on top of a least-squares problem with a large weight and the answer approaches the constrained one like 1/τ². The limit is takeable to any accuracy — and how far it can be taken is a property of the solver, not of the problem. One of them stops at the square root of the precision, and one of them does not stop.

 ℚ6𝔽24𝔽35𝔽56𝔽76𝔽116𝔽136𝔽1016𝔽655376rank, by the ring the entries are read inthe rank of one matrixover ℚ6over 𝔽24over 𝔽35over 𝔽56over 𝔽76nothing is rounded hereand the answer still is not a property of the matrix Exact arithmetic, and what it costs instead

The rank depends on the ring

A floating-point rank is a decision about a threshold. Remove the arithmetic error entirely and the threshold goes away — and the answer still is not a property of the array of numbers, because one integer matrix has rank six over the rationals, five modulo three and four modulo two, with nothing rounded and nothing decided.

 10¹10²10³10⁻⁹10⁻⁷10⁻⁵10⁻³10⁻¹λ, over the spectrum of −A|r(λ)|geometric maxequally spacedgeometricone line of codeshifts6κ of the spectrum389geometric max0.12equally spaced max0.89the factor between7.5the same k solvesand one choice of where Reduction, and what a model is for

Where to put the poles of a rational function

Three times in one field the same question arrives from different directions — ADI shifts, rational approximation of a square root, the decay of a Gramian — and it has one answer. Cluster them geometrically towards wherever the function is difficult, and the alternative that looks reasonable costs orders.

 2022242628303210⁻¹⁷10⁻¹⁴10⁻¹¹10⁻⁸10⁻⁵10⁻²k, where the entries are near 2ᵏrelative error of the determinantas writtenfused: exactproducts need 54 bitsone rounding, the whole answertrue determinant1naive, k = 300fused, k = 301first wrong at k27sizes returning 06both forms conformand the source does not say which The answer that depends on the machine

One multiply the compiler removed

A determinant whose value is exactly 1, computed as exactly 0 by the expression that is written down, and exactly 1 by the same expression with the multiply and the add fused. Both forms conform to IEEE-754, both are legal compilations of the same source, and nothing in the program says which one you have.

 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.

 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.

 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.

 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.

 10²10³10⁴10⁻⁸10⁻⁷10⁻⁶10⁻⁵10⁻⁴terms summedrelative error+∞, −∞ 1.01zero 1.00stochastic 0.50nearest 0.47√n against n, fittednearest, fitted exponent0.47stochastic, fitted exponent0.5toward +∞, fitted exponent1twelve seeds averaged at each sizethe slope is the bias, not the precision The arithmetic underneath

The direction the error leans

The size of one rounding error is set by the precision. How ten thousand of them combine is set by something else entirely — the rounding mode — and the fitted exponents are 0.47 for round-to-nearest and 1.01 for round-toward-infinity, on identical data at identical precision.

 -101-101(x − 1) ÷ 0.001p(x) ÷ 2.68·10⁻¹⁵a square, drawn negativepoints sampled401negative, unfused179negative, fused196signs disagree at98widest gap1.4·10⁻¹⁵unfusedfusedthe true valuethe true curve is the flat oneeverything else is rounding The answer that depends on the machine

A square that evaluates negative

(x − 1)⁶ evaluated near x = 1 comes out negative at 179 of 401 points on one build and 196 on another, and the two disagree about the sign at 98 of them. Neither is nearer the truth: both traces are made entirely of rounding.

All themes