Numerical Linear Algebra
where the algebra stops being the arithmetic
Theme

The thread: Assertions that reject — page 6

Essays 121 to 142 of the 142 on this theme, in the same order.
10⁻¹²10⁻¹⁰10⁻⁸10⁻⁶10⁻⁴10⁻²10⁻¹⁴10⁻¹10¹²10²⁵10³⁸10⁵¹10⁶⁴δ, the gap between consecutive eigenvaluesrelative error in e^Aan answer with no correct digitsV f(Λ) V⁻¹scaling and squaring‖A_δ − A₀‖exact eigenvalues throughoutκ(V) at the smallest δ3.3·10⁸²eigen route2.9·10⁶⁵scaling and squaring4.1·10⁻¹²distance to the limit7·10⁻¹²the eigenvalues are the diagonaland they are exact at every stop Eigenvalues, singular values, rank

The series that has to be squared back

The Taylor series for the matrix exponential is not wrong — every term is computed correctly — and on Moler and Van Loan's two-by-two its largest term is 5.4 million times the answer it sums to. The method that replaces it scales the matrix down and squares the result back, and both halves of that sentence cost: too few squarings and the approximant is out of range, too many and each one doubles the rounding.

 Cheeger's band, for the circulationcirculation · arcsvertices on the smaller sidewhich of the two the theorem is aboutλ₂ of 𝓛0.066λ₂/20.033√(2λ₂)0.36circulation Φ0.12arc conductance0.04their ratio3.1the inequality holdsabout a quantity nobody counts The matrix that is a graph

A conductance the arcs do not measure

Symmetrising a directed Laplacian with respect to its walk recovers everything the arrows took — a real spectrum, a sweep cut, a Cheeger inequality. What it does not recover is the quantity: the inequality bounds the probability that a step of the walk crosses the cut, which on one graph here is three times the weight of the arcs that do.

 012345610⁻¹10²10⁵10⁸10¹¹log₁₀ γ, the change of unitscondition numberthe linearisationthe quadraticone problem, two amplifiersκ(quadratic), first6.6κ(quadratic), last6.6κ(linearisation), last10¹²how far the first moved1the problem is as well conditioned as everand the method is not The eigenvalue problem that is not linear

The number that moves when the problem does

Two quantities are offered as the condition number of one eigenvalue. One is unmoved to eight digits by a change of variable that is exact in both directions, and grows like the square root of the chain's length. The other is inflated by ten orders by that change of variable, and is ten times too large before anything has been done at all.

 ‖QᵀQ − I‖ of the implied Qκ = 10⁸, Cholesky0.37κ = 10⁸, sweep8.5·10⁻⁹κ = 10¹⁰, Cholesky1.3κ = 10¹⁰, sweep3·10⁻⁷κ = 3·10¹⁰, Choleskyrefusedκ = 3·10¹⁰, sweep6.7·10⁻⁶κ = 10¹², Cholesky1.5κ = 10¹², sweep1.4·10⁻⁴the safe run is the one that failsrefusals in the range1wholly non-orthogonal returns2the pivot it refused on-6.5·10⁻¹⁷one of these outcomes is safeand it is the refusal The answer that depends on the machine

The licence is not the boundary

Cholesky QR is licensed by κ²u ≪ 1, which reaches equality at κ = 9.5·10⁷ in double precision. At 10⁸ the factor it returns is already 0.37 away from orthogonal, and it goes on returning factors as far as 10¹³ — refusing at scattered condition numbers in between, at different ones for eight columns and for six.

 036912151810¹10²10³rotationswidththe enclosurethe seta rotation is an isometrymeasured growth a step1.4√2, from the geometry1.4enclosure ÷ set after 201024no rounding error is responsible for any of thisa higher precision does not touch it The arithmetic underneath

Nine steps of pessimism

A proved bound is 8 to 26 times the error it bounds, at every precision from 16 to 40 significand bits. A carried interval is (√2)ᵐ times too wide after m re-enclosures. The two cross between eight and nine, so the method everybody warns against is the tighter of the two for a short computation.

 0481216202410⁻¹⁷10⁻¹⁴10⁻¹¹10⁻⁸10⁻⁵10⁻²10¹matrix–vector products, mrelative error in e^Abforming e^A: 3·10⁻¹⁶crosses at m = 18the vector, not the matrixsteps to the dense answer18dimension100Krylov megaflops0.36dense megaflops2the exponential that is computedis 18×18 Eigenvalues, singular values, rank

The vector was what was wanted

Nobody who computes a matrix exponential wants the matrix. They want eᴬᵗb — one vector, the state of a system at a later time. Twenty matrix–vector products get it to sixteen digits on a hundred-by-hundred problem, without ever forming a hundred-by-hundred exponential, and the exponential that does get computed is twenty by twenty.

 -6-4-20246-6-4-20246real partimaginary partan infinite spectrum, finitely askedinside the contour12probes4moments used3values returned24worst against the closed form6.2·10⁻¹⁵infinitely many eigenvaluesand a question with an answer The eigenvalue problem that is not linear

Where a contour's budget should go

A contour method's ceiling is the number of probes times the number of moments, and the moments are free in solves while the probes are not. Four ways of reaching one ceiling come out four orders apart, the ordering is not monotone, and what separates the best two is not the usual draw but the unlucky one.

 worst nodal error, and nodes outside the interval the equation guarantees0°, streamline only2.39·10⁻¹⁷ · 0 outside0°, with crosswind5.25·10⁻¹⁷ · 0 outside15°, streamline only0.032 · 18 outside15°, with crosswind0.0145 · 9 outside30°, streamline only0.0571 · 50 outside30°, with crosswind0.0247 · 11 outside45°, streamline only0.0661 · 48 outside45°, with crosswind0.0281 · 0 outsidewhat the crosswind term buyserror ratio at 0°0.45error ratio at 15°2.2error ratio at 30°2.3error ratio at 45°2.3free where the scheme was exactand half the error everywhere else Iterating, instead of factorising

The direction the diffusion does not go

Streamline diffusion adds τbbᵀ, a rank-one tensor that annihilates every direction across the flow. That is the design. The price is 18, 50 and 48 nodes where the computed solution leaves the interval the equation guarantees — and half a coefficient of crosswind diffusion halves the error at every angle while costing exactly nothing where the scheme was exact.

 10⁻¹⁷10⁻¹⁵10⁻¹³10⁻¹¹10⁻⁹10⁻⁷10⁻⁵10⁻³10⁻¹10⁻¹⁴10⁻¹²10⁻¹⁰10⁻⁸10⁻⁶10⁻⁴10⁻²1εrelative error in J(x)vforwardcentralcancellationtruncationagainst a derivative that is exactforward floor1.3·10⁻¹⁰central floor1.1·10⁻¹²truncation slope, forward1truncation slope, central2no ε reaches the roundoffand the analytic derivative is free of the choice Iterating, instead of factorising

An operator with no entries

At the sizes where linear algebra is expensive the matrix does not exist. What exists is a subroutine that returns Av. Every Krylov method survives that unchanged; every algorithm that reads an entry disappears. And the derivative such a code computes is accurate to ten digits instead of sixteen, which turns out to cost nothing at all.

 -40-27-14-11225380123456789computed eigenvalueseedevery mark has a residual below 10⁻⁸a small residual, and no answercoefficients of det(A − λB)0worst residual over all seeds1.8·10⁻⁹spread of the answers66seeds drawn8the residual is small at every markand none of the marks means anything Eigenvalues, singular values, rank

A problem with no answer

If two matrices share a null vector then det(A − λB) is identically zero and every λ is an eigenvalue, which means none of them is. Perturb such a pencil by a ten-billionth and a solver returns six numbers with residuals below 10⁻⁹. Change the seed and it returns six different numbers, spread over forty-four, with residuals just as small.

 024681012141618024681012141618distinct eigenvalues in the spectrumstep the recurrence stops atthe step is m, not nn = 30 throughoutspectra drawn8every one breaking at m8worst residual at the breakdown5.6·10⁻¹⁶smallest gain over the step before3.5·10¹⁰an invariant subspace contains the answerand its dimension is what the method costs Iterating, instead of factorising

The zero that means it is finished

Every Krylov method ends by dividing by a number the previous step produced, and when that number is zero the recurrence stops. In Arnoldi the stop is the answer — the subspace has closed, the solution is inside it, and the residual is at the unit roundoff. The literature calls it a lucky breakdown, and the adjective is doing real work.

 ran to the end3173lucky — a subspace closed295serious, cured by a block of two495serious, cured by a longer block29serious, incurable at any length8counted, not estimatedserious, as a fraction0.13of those, cured at two0.93incurable8matrices tried4000measure zero on the realsand an eighth of the integers Iterating, instead of factorising

The same zero, and nothing was found

Change the recurrence by two lines and the divisor stops being a norm. It becomes an inner product of two vectors from two different sequences, and an inner product of two different vectors is zero on a whole hyperplane — with neither vector anywhere near zero, nothing invariant, and nothing converged. The arithmetic event is identical and the meaning is opposite.

 01428425670849811210⁻²²10⁻¹⁸10⁻¹⁴10⁻¹⁰10⁻⁶10⁻²conjugate gradient iterationrelative residualthe unit roundoff, 1.11·10⁻¹⁶the answer's residualthe residual reportedtwo residuals, one runreported, at its best6.9·10⁻²¹the answer's, at its best5.1·10⁻¹⁰unit roundoff1.1·10⁻¹⁶largest iterate on the way9.3·10¹³iterations drawn110the recurrence remembers every roundingand the stopping test is written in it Iterating, instead of factorising

The residual the method reports

Conjugate gradients prints a relative residual of 6.9·10⁻²¹. The unit roundoff is 1.1·10⁻¹⁶, so that is not a small residual and not a large one — it is not a residual. The vector the method is holding at that step has ‖b − Ax‖/‖b‖ = 5.1·10⁻¹⁰, and nothing in the run says so.

 λ = 105 timesλ = 9.55 timesλ = 95 timesλ = 8.55 timesλ = 2.952 timesλ = 2.92 timeseigenvalues that arrived more than once — the matrix has 40 distinct onesa spectrum with the wrong multiplicitiesextra copies, no reorthogonalisation25extra copies, full reorthogonalisation0worst relative error among the copies1.9·10⁻⁸steps taken of 80 asked for, full40no arithmetic error was madeevery one of these is right to eight digits Eigenvalues, singular values, rank

An eigenvalue that arrives twice

A matrix with forty distinct eigenvalues, handed to Lanczos for eighty steps, returns twenty-five extra copies of thirteen of them — the largest arriving five times. Every copy is accurate to 1.9·10⁻⁸ relative. No arithmetic error was made, nothing overflowed, and a caller counting eigenvalues gets the wrong multiplicity from a computation in which no individual number is wrong.

 081624324010⁻¹⁴10⁻¹¹10⁻⁸10⁻⁵10⁻²conjugate gradient steprelative residualabove, on their own scale: rank asked 5, rank kept 3solid: the budgeted residual · light: no budgetbudget 3rank asked for5rank kept3residual, budgeted3.8·10⁻⁴residual, unbudgeted1.6·10⁻¹⁴numbers stored150the step asks for moreat every step Iterating, instead of factorising

An iterate that must be made smaller

Applying a Kronecker-sum operator to a low-rank iterate multiplies its ranks by d and adding two of them adds their ranks, so a solver in a compressed format cannot keep what it produces. Every step is followed by a truncation — and whether that truncation is a floor on the residual depends on the right-hand side rather than on the truncation.

 1611162126313610⁻¹³10⁻¹¹10⁻⁹10⁻⁷10⁻⁵10⁻³10⁻¹10¹indexmagnitude|r_kk|σ_kone factorisation, two verdicts‖AP − QR‖/‖A‖10⁻¹⁵|r_nn|1.1·10⁻¹²σ_min10·10⁻¹³column interchanges33|r_nn| is never below σ_minso the cheap verdict errs one way only Eigenvalues, singular values, rank

A good curve and a bad verdict

The diagonal of a column-pivoted R is famous for the one matrix it is wrong about. On that matrix it is right about thirty-nine of its forty entries — every |r_kk| within a factor of six of the σ_k it stands for — and wrong by 4·10⁶ at the fortieth, which is the only one a rank verdict ever reads.

 01020304050607010⁻⁶10⁻³110³10⁶10⁹k‖Aᵏ/k!‖‖e^A‖ = 2.6largest term 1.4·10⁷what the series throws awaylargest term1.4·10⁷‖e^A‖2.6digits cancelled away5.4·10⁶error after the sum5.2·10⁻⁹every term is computed correctlyand the sum has lost seven digits Eigenvalues, singular values, rank

The error the method already knows

Summing the exponential's Taylor series throws away a known number of digits, and the number is on the machine while the sum is being formed. The largest term divided by the answer, times the unit roundoff, tracks the relative error that comes out — to within a factor of nine, across fourteen orders of magnitude of it — and nothing reports it.

 012345678910⁻¹³10⁻¹⁰10⁻⁷10⁻⁴10⁻¹rank budgetresidual the run stalls atb with no structure at allb constant: the answer is a traintwo ladders, one truncationstructured, rank 10.12structured, rank 45.5·10⁻⁶structured, rank 85.9·10⁻¹⁴unstructured, rank 10.97unstructured, rank 80.47the floor is not the truncation'sit is the answer's Iterating, instead of factorising

A run that is over at step five

A conjugate gradient whose every iterate is cut to a rank budget reaches the floor that budget allows at step 5, 36, 42 or 59, and then does nothing for the rest of the run. Four times the iterations move the floor by a factor of 1.8, and past the answer's own rank they move it the wrong way.

 10⁻¹100.20.40.60.8strength threshold θresidual reduction per cycleθ = εsemi-coarseningkept whole rowsfull coarseningone parameter, two methodsbest factor above ε0.053best factor below ε0.1the ratio across the switch231×31 anisotropic operatora switch, not a dial Iterating, instead of factorising

The switch does not know which side is better

The strength threshold moves the coarsening from full to semi at θ = ε exactly, at every anisotropy. Which of the two converges faster is a separate question with a separate answer, and it changes sign between ε = 0.33 and ε = 0.34 — where nothing whatever happens to the switch.

 10⁻³10⁻²10⁻¹1024681012size of the negative eigenvalue, −λproducts before the test firesharder to find, and milder8 spectra, n = 50products at the largest λ3products at the smallest10smallest share of λ recovered0.14largest0.34the one that hidesis the one that matters least Iterating, instead of factorising

A proof that does not ask how large the matrix is

Proving a Hessian indefinite costs three matrix–vector products when the negative eigenvalue is 3 and nine to eleven when it is a thousandth, and that pair of numbers barely moves across a fourfold range in n. The factorisation that settles the same question costs a third of n³, which grows by a factor of sixty-four over the same range.

 10³10⁶10⁹10¹²10¹⁵10⁻¹⁵10⁻¹²10⁻⁹10⁻⁶10⁻³1largest iterate on the way, ‖x‖difference between the two residualsthe bound, linear in ‖x‖what it does, slope 0.46a bound of one, a walk of a halffitted slope0.46the bound's slope1smallest gap measured3.9·10⁻¹⁵largest gap measured1.5·10⁻¹⁰share of the bound, near end9.1·10⁻⁴share of the bound, far end3.4·10⁻⁹the bound is rightand loose by a square root Iterating, instead of factorising

A walk needs a length

The gap between the two residuals grows as the square root of something, and a square root needs a length. Two quantities are candidates — how far the iterates travelled and how many steps were taken — and only a second sweep separates them. Across a fourfold change in size the iteration count goes from 39 to 96 and the gap goes from 5.04·10⁻¹⁵ to 5.33·10⁻¹⁵.

 110¹00.30.60.91.2trust-region radius Δshare of the exact model decreasethe exact subproblem7 radii, n = 60share at the smallest radius0.95share at the largest0.3products, at most8radii stopped by the curvature4a few products against an eigendecompositionand most of the decrease Iterating, instead of factorising

The certificate that arrives soonest is worth least

The more negative a Hessian's smallest eigenvalue, the sooner conjugate gradients meets a direction of negative curvature — and the less of the exact trust-region decrease that direction turns out to be worth. At λ_min = −10 the step arrives after two products and gets 39.6 per cent; at −10⁻³ the same two products get 89.8, and the whole sweep costs eight.

All themes