Numerical Linear Algebra
where the algebra stops being the arithmetic
Theme

The thread: Where the algebra stops being the arithmetic — page 7

Essays 145 to 159 of the 159 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

A function of a matrix is not a function of its entries

Everybody learns that f(A) means diagonalise, apply f to the eigenvalues, undiagonalise. That is a definition, not a method. On a matrix seven picometres from a defective one — with exact eigenvalues and eigenvectors from a closed form — the definition returns an answer wrong by sixty-five orders of magnitude, and a method that never mentions an eigenvalue returns the right one.

 75 aggregates over 225 unknownsthe matrix chose thismean extent along y3mean extent along x1points adopted by pass two15no coordinate enters the methodand the shape follows the coupling Iterating, instead of factorising

Aggregating what the matrix calls strong

The depth phase measured every method it had on the 45°-rotated anisotropic operator — 0.784, 0.883, 0.844 — and diagnosed the failure as being in the discretisation rather than in the hierarchy. Smoothed aggregation is the standard answer to anisotropy. It returns 0.789.

 02468101201change of units, by exponent10⁰10¹10²10³10⁴10⁶10¹⁰10¹⁶10¹⁹10²⁰10⁴⁰10¹⁵⁰as writtenafter scalinga range questionlargest finite value3.4·10³⁸predicted boundary γ1.8·10¹⁹last γ that forms10¹⁹stops where scaling fails0not a poor answerno answer at all The arithmetic underneath

The units that overflow before the answer does

A change of variable that is exact in the algebra requires γ² times a matrix to be a number the format can hold. In binary64 that is a bound nobody meets by accident. In binary32 it arrives at 10¹⁹ and in fp16 at 256, and past it there is no answer rather than a poor one.

 0246810121410⁻¹⁷10⁻¹⁶10⁻¹⁵10⁻¹⁴10⁻¹³10⁻¹²10⁻¹¹10⁻¹⁰10⁻⁹10⁻⁸squarings, srelative error against the closed formbest: s = 2the approximant is out of rangethe squarings amplifyboth sides costerror at s = 06·10⁻¹²best, at s = 26.9·10⁻¹⁶error at s = 153.2·10⁻¹²‖A‖ / 2^s at the best1.1the squarings are exact in the algebraand each one doubles the rounding 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.

 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^Atb — 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.

 λ = 0λ = ∞λ = ∞-0.37690.5654.3836.4282 eigenvalues here, and it is one placecounted exactly, in rationalsfinite eigenvalues4at infinity2degree of det(A − λB)4worst residual, either kind5.5·10⁻¹⁴an eigenvalue is a ratioand a ratio has a direction, not a size Eigenvalues, singular values, rank

An eigenvalue with no value

If the second matrix of a pencil is singular then some of the eigenvalues are infinite, and that is not a degeneracy — it is the algebraic constraints of the model, one per constraint. What survives is a pair of numbers rather than one, and on the line those pairs live on, infinity is an ordinary point with an ordinary residual.

 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.

 1234567810⁻¹⁷10⁻¹⁴10⁻¹¹10⁻⁸10⁻⁵10⁻²10¹singular value, largest firstrelative errorone-sided Jacobizero-shift QRshifted QReigenvalues of BᵀBagainst a rational bisectionσ_min, exactly2.1·10⁻³⁰worst, one-sided Jacobi4.4·10⁻¹⁶worst, zero-shift QR2.2·10⁻¹⁶worst, eigenvalues of BᵀB1a relative error is a ratioand the denominator is the answer Eigenvalues, singular values, rank

Small compared to what

This site's own singular value routine has carried a sentence since the month it was written — that one-sided Jacobi computes the small singular values to high relative accuracy and the standard method does not. It has never been measured here, because measuring it needs a σ that is known rather than computed. A bidiagonal matrix and a Sturm count in exact rationals supply one.

 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.

 012345678-7-5-3-11357conjugate gradient steppᵀAp ⁄ pᵀpλ_min = -0.1positive: a step existsnegative: a certificate existsone matrix, two questionsstep it turns at6quotient there-0.027share of λ_min recovered0.27λ_min, by construction-0.1MINRES steps on the same system37the division that cannot be doneis the answer to a different question Iterating, instead of factorising

The division that cannot be done

Conjugate gradients divides by pᵀAp at every step, and on a matrix that is not positive definite that number can be zero or negative. This site has guarded against it since its first commit and described it as a failure. In the method that made conjugate gradients famous it is the single most valuable object the iteration can produce, and it costs six matrix–vector products.

 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.

 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.

 33.310⁻⁴10⁻³10⁻²10⁻¹1log₁₀ numbers storeddistance to the dominant eigenvalueArnoldi, linearisedprojected quadraticper number heldstorage, linearised2080storage, second-order1040Ritz values, linearised26Ritz values, second-order52half the storageand twice the approximations Iterating, instead of factorising

A Krylov space for a problem that is not linear

A quadratic eigenvalue problem has no matrix to build a Krylov space out of. The recurrence that builds one anyway stores half as many numbers, returns twice as many Ritz values — and stops being a basis at twenty vectors while the answer it gives keeps improving.

All themes