Numerical Linear Algebra
where the algebra stops being the arithmetic
Theme

The thread: Two routes to a number — page 6

Essays 121 to 144 of the 151 on this theme, in the same order.
027548110813510⁻⁶10⁻⁴10⁻²110²10⁴10⁶power‖Aᵏ‖Kreiss constant 6760e · n · K‖Aᵏ‖ρᵏtwo routes to one peakspectral radius0.8peak of ‖Aᵏ‖2·10⁴Kreiss constant6757e · n · K1.1·10⁵everything here decays in the endand one of these curves says how much first Eigenvalues, singular values, rank

A spectral radius that grows first

ρ(A) below one guarantees that the powers of A go to zero and says nothing about what they do on the way. Here they rise by a factor of twenty thousand before turning over, and the peak is bracketed above and below by a constant computed from the resolvent norms outside the unit circle — two routes to one number, one through the plane and one through the powers.

 110¹10²10³110¹10²10³the accumulating quantity, relative to its first valuethe error, relative to its first valuethe bounds: slope 1what all three do: slope ½three mechanisms, one exponenta left-to-right sum0.49a chain of rotations0.55a residual recurrence0.51every bound's slope1spread of the three0.067a bound is a sum of the roundingsand the roundings have signs The arithmetic underneath

Three walks and one bound

A left-to-right sum, a chain of three thousand rotations and a conjugate gradient residual recurrence share no arithmetic and no vocabulary. Each has a standard bound that is linear in whatever it accumulates against. All three come out at a half — 0.486, 0.554 and 0.507 — and nothing is rescaled.

 00.10.20.30.40.50.60.70.80.9110⁻¹⁶10⁻¹³10⁻¹⁰10⁻⁷10⁻⁴10⁻¹decay rate of the spectrumrelative spread of the rank-k errorthe seedthe machinetwo kinds of variationseeds drawn8partitionings7seed spread0.31machine spread9.8·10⁻¹⁶ratio3.1·10¹⁴one of these is recordedand it is the large one The answer that depends on the machine

The variation that comes with a seed

A randomised low-rank approximation's error moves by 31% between draws and by 10⁻¹⁵ between partitionings of one draw. In the one field on this site whose answer already comes as a band, the machine is inside the width of the line — and it is still there.

 a triangle and a square, joinedK₂,₃ with a pendant edge1234560123456indexeigenvaluethe same, and not the samevertices each6edges each7spanning trees12spectra differ by5.3·10⁻¹⁵highest degree, left4highest degree, right3one has a trianglethe other is bipartite The matrix that is a graph

The spectrum is not the graph

Two graphs on six vertices with the same Laplacian characteristic polynomial — as integer polynomials, not to fourteen digits. One contains a triangle; the other is bipartite. Every method in this field that reads only the spectrum is answering about the class.

 10⁻¹⁷10⁻¹⁴10⁻¹¹10⁻⁸10⁻⁵10⁻²-1-0.500.51relative backward error acceptedfraction of the work not doneQR deflation criterionincomplete Cholesky droplow-rank truncationthree fields, no shared arithmeticslope, deflation criterion0.039slope, drop tolerance0.051slope, rank truncation0.23widest apart, as a ratio5.7a tolerance is an offerand the three offers are one offer Two errors, and whose fault they are

A tolerance is priced by the problem

Three tolerances from three fields sit on one pair of axes and agree to within a factor of 5.74. That factor is the ratio of the two curves that cannot move. Change the only problem in the comparison and the third curve's fitted slope swings from 0.188 to 0.040 while the printed spread does not shift by a digit.

 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.

 κ of the suma component of b − Ax8.09·10¹⁷ad − bc, near-degenerate3.6·10¹⁶qᵢᵀqⱼ, an orthogonality check3.26·10¹⁶zᵀAz, a trace probe95.7pᵀAp, a curvature73.7rᵀr, a residual norm1measured, not assumedhighest8.1·10¹⁷lowest1above 10¹⁰3terms128sums of squares are safeand nobody decides anything from one The answer that depends on the machine

Two machines, one certificate

Nothing a solver returns says which of its answers you got. Four things could be reported instead — the summation condition number, the partition count, an exactly accumulated residual and a directed-rounding interval — and each costs about one pass over data the routine already has in hand.

 Re λIm λwhat survives the arrowsvertices18arcs18worst row sum0worst column sum0largest |Im λ|0.98asymmetry1the null vector is still exactand nothing else about the spectrum is real The matrix that is a graph

A Laplacian that is not symmetric

Point the edges and the matrix stops being symmetric. Its row sums are still exactly zero, so the null vector survives; everything built on the quadratic form does not, and the eigenvalues need a plane rather than a line. Asymmetry permits that and does not force it, which the smallest case here demonstrates by being asymmetric and real.

 finite eigenvalues (degree of det Q)13at infinity (2n − degree)1at infinity, by the rank of M12n, if M were nonsingular14a degree, not a decisiondegree of det Q13at infinity1by the rank of M1singular-value gap∞the count is a degreeand the other route is a judgement The eigenvalue problem that is not linear

One mass removed, and one eigenvalue gone

A coordinate with no inertia reads like a coordinate that has been deleted, and a chain of eight masses with one of them removed would then be a chain of seven, with fourteen eigenvalues. It has fifteen. The massless coordinate is still there, still carrying a damper, and it contributes a first-order equation rather than none.

 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.

 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.

 10¹10⁴10⁷10¹⁰10¹³10¹⁶10⁻¹⁷10⁻¹⁴10⁻¹¹10⁻⁸10⁻⁵10⁻²10¹κ(B)relative error, and asymmetryvia B⁻¹Avia Choleskyasymmetry of B⁻¹Au · κ(B)against a spectrum known exactlyslope, via B⁻¹A0.92slope, via Cholesky0.98worst ratio between them2.3asymmetry of B⁻¹A1.1the symmetry claim is trueand it is not about the accuracy Eigenvalues, singular values, rank

Two matrices and one problem

Ax = λBx is what a finite element model, a structural vibration and a constrained optimisation actually produce, and it is not the one-matrix problem with a change of variables. Everybody is told not to form B⁻¹A because it is not symmetric. That is true, the departure from symmetry is about one, and it is not what decides the accuracy.

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

 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.

 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.

 051015202530354010⁻³10⁻¹10¹GMRES steprelative residualthe two curves are the same curvea number re-derived, not carriedworst reported/actual factor1.8at step40‖VᵀV − I‖ of the basis1.4reported at the last step0.041actual at the last step0.074the same family of methodsand only one of them lies Iterating, instead of factorising

The number that is re-derived

GMRES prints a residual it never computes from its answer either. On the matrix that sends a conjugate gradient recurrence 7.3·10¹⁰ wrong, and on two others chosen to be worse, its number is never more than a factor of 2.86 out — while the basis it is computed from has lost orthogonality entirely. The disease is not iterative methods, and it is not floating point.

 -4-2.35371-0.7074290.9388572.585144.231435.8777102468101214shift σν(σ)an answer that is an integershifts2000disagreements0eigenvalues14steps14the marks are a Jacobi decompositionand the staircase never saw one Eigenvalues, singular values, rank

An eigenvalue count that cannot be slightly wrong

Every spectral computation on this site returns floats with errors in them. Counting eigenvalues below a shift by the signs of an unpivoted elimination returns an integer, and an integer cannot be 6.9999999997 — so the answer is exactly right, or wrong by a whole eigenvalue, and where the second happens is a band of measurable width.

 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.

 0102030405010⁻¹⁷10⁻¹⁴10⁻¹¹10⁻⁸10⁻⁵10⁻²10¹10⁴10⁷decades of gradingworst relative errorthe answer is gonevia BᵀBone-sided Jacobizero-shift QRone axis, four routesBᵀB at the narrowest grading3.6·10⁻¹⁴and at the widest1.9·10⁷Jacobi, worst over the sweep1.5·10⁻¹⁵zero shift, worst1.1·10⁻¹⁵the definition is not a methodand squaring buries what it squares Eigenvalues, singular values, rank

A threshold the matrix does not set

Two numbers come out of a relative-accuracy comparison and they belong to different things. The size of the matrix moves the constant of the routes that never fail, by a factor of 2.7 between n = 4 and n = 10; it does not move the point where the route through BᵀB stops returning an answer, which sits between ten and eleven decades of grading at every size drawn.

 00.250.50.75100.250.50.751xuexact at ε = 0.005exact at ε(1 + Pe)the upwind answeran exact answer to a different question‖upwind(ε) − central(ε(1+Pe))‖/‖·‖0distance to the problem it solves0.026distance to the problem posed0.36the added diffusion is h/2 = 0.001953, whatever ε isso refining removes it Iterating, instead of factorising

A different equation on every grid

Upwinding is the exact discretisation of a convection–diffusion problem with diffusion ε + h/2, entry for entry, at a relative difference of between 0 and 1.26·10⁻¹⁶ on every mesh from 15 points to 511. The equation it is exact for is chosen by the mesh and not by ε — the added diffusion is 0.01563 on a 31-point grid whether ε is 0.2 or 0.001.

All themes