Numerical Linear Algebra
where the algebra stops being the arithmetic
Theme

The thread: Claims that can fail — page 2

Essays 25 to 48 of the 213 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.

 forward error, relative to a solution of exactly (1, 1)no pivoting · as given1 0 interchangesno pivoting · rows scaled1 0 interchangespartial · as given0 1 interchangepartial · rows scaled1 0 interchangesscaled partial · as given0 1 interchangescaled partial · rows scaled0 1 interchangecomplete · as given0 1 interchangecomplete · rows scaled0 1 interchangethe same problem twicepartial, as given10⁻¹⁸partial, rows scaled1its relative residual10⁻¹⁷complete, rows scaled10⁻¹⁸the two systems have the same solutionand one pivot rule cannot see it Elimination, and the swap

The pivot that reads the units

Partial pivoting compares the entries of a column and takes the largest. Those entries carry units, so the comparison depends on them — and there is a row scaling, on the standard two-by-two that pivoting exists to fix, which makes partial pivoting perform the identical catastrophic elimination it was introduced to prevent, with no interchange at all.

 02468101210⁻¹⁴10⁻¹²10⁻¹⁰10⁻⁸10⁻⁶10⁻⁴10⁻²1iteration‖r‖ / ‖b‖Laplaciancyclic shiftno progress at allevery eigenvalue of the shift is on the unit circleand it predicts nothing Iterating, instead of factorising

The spectrum that predicts nothing

For a symmetric matrix the eigenvalues govern how fast an iteration converges. Drop symmetry and they stop governing anything — there is a matrix whose eigenvalues are as evenly spread as eigenvalues can be, on which GMRES makes no progress at all until the last possible step.

 00.250.50.75110⁻¹110¹share of the noise placed in the matrixleast-squares error ÷ total least-squares errorequally accuratetotal leastsquares aheadordinary leastsquares aheadthe model, not the methodadvantage, all noise in b0.28advantage, all noise in A2.5seeds at each share40the same total noise at every pointand only where it sits changes Least squares, and the road not to take

When the matrix is wrong too

Every least-squares problem here has assumed A is exact and b is not, and moved b onto the column space of A. Where both were measured, the smallest correction that makes the system consistent moves the matrix as well — and on the problems where that answer is more accurate, it has the larger residual, by construction rather than by luck.

 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⁻¹⁶.

 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.

 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.

 det mod p, as a fraction of p3713a prime that divides the answerdet A3·10⁴primes swept25unlucky5rate0.2det, in bits15singular mod p is not singularand one residue cannot tell them apart Exact arithmetic, and what it costs instead

A prime that divides the answer

A modular elimination reports a singular matrix and is telling the truth — over the field with p elements the matrix is singular. Over the rationals it is not. Nothing in the residue distinguishes the two cases, no quantity is small enough to be suspicious, and the wrong answer is a correct computation of a different question.

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

 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.

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

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

 545545545455554545545545545454555555454545545545545455554545545545rows and columns, in the order the points arrivedark: kept dense · light: two thin factors, rank printedthe strong partitionblocks112kept dense46largest rank5numbers stored2.7·10⁴relative compression error3.4·10⁻¹⁰the picture is decidedbefore a number is read Neither sparse nor dense

Which pairs are allowed to be small

A hierarchical representation is a partition of the matrix into blocks, and the rule that produces it reads four numbers per pair of index clusters and not one entry of the matrix. On a 256-square it yields 112 blocks, 66 of them stored as two thin factors, none of rank above five.

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

 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.

 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.

 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.

 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.

 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.

 567891003691215log₂ nlargest rank in the partitionweak: touching pairs compressedstrong: touching pairs refusedthe test bounds a rank and costs storagestrong, blocks250weak, blocks94strong, numbers6.8·10⁴weak, numbers6.1·10⁴strong ⁄ weak1.1the better partitionis the more expensive one Neither sparse nor dense

The test that costs what it saves

The partition that refuses to compress a touching pair keeps every rank at five while the other lets them climb from nine to thirteen. It also stores more numbers at every size measured — 67,968 against 61,440 at n = 512 — and which of those two facts matters is a question about how large the problem is going to get.

 10¹10²10³03691215what one rebuild is worth, in preconditioned iterationscheapest number of members between rebuildsa dense Cholesky here is worth 6.7 iterationsringed: where the growth rule beats every fixed periodcheapest period, by setup costsetup worth 51setup worth 102setup worth 203setup worth 506setup worth 1008setup worth 20012how long to keep itis a question about what it cost to build When the problem arrives again

What a rebuild is worth

One sequence, one drift, one preconditioner — and six different right answers, because the cheapest rebuild period depends on what a rebuild cost to build and on nothing else. The optimum walks from every member to every twelfth as the setup gets dearer, and the free rule that reads the iteration count beats it in the middle of that range and loses at both ends.

All themes