Numerical Linear Algebra
where the algebra stops being the arithmetic
Theme

The thread: Two routes to a number — page 3

Essays 49 to 59 of the 59 on this theme, in the same order.
10⁻⁴10⁻³10⁻²10⁻¹100.250.50.751anisotropy εsmoothing factor μpoint + fully-line + fullpoint + semi-y⅓, the one-dimensional answertwo routes, three curvesgap between the repairs2.2·10⁻¹⁶scan against closed form2.2·10⁻¹⁶the dashed curve lies on the solid one beneath itone repair, written two ways Iterating, instead of factorising

Smoothing a whole line at once

Solve every grid line in the strong direction exactly rather than sweeping over it, and the smoothing factor goes from 0.9993 back to 0.3340 — which is the one-dimensional answer, on a problem that is not one-dimensional. The repair replaces one ε in the closed form by a one.

 1357910⁻⁵10⁻⁴10⁻³10⁻²10⁻¹1eigenvalue λamplificationkeptdiscardeda filter on the starting vectormeasured against the polynomial3.9·10⁻¹¹worst kept direction1best discarded direction0.8the roots are the discarded Ritz valuesand the vertical lines are where they sit Eigenvalues, singular values, rank

Restarting is a filter

A restart throws away the Ritz values it does not want and begins again from a new starting vector. Written in the eigenbasis, that vector's components have been multiplied by a polynomial with its roots at the discarded values — measured component by component, and agreeing with the polynomial to rounding.

 point · strong in y0.951point · strong in x0.956y-line · strong in y0.037y-line · strong in x0.967semi-y · strong in y0.107semi-y · strong in x0.9661.00 — no convergenceone problem, seen from two sidesy-line, aimed0.037y-line, turned sideways0.97the ratio2631×31 grid, twelve V-cyclesthe direction is in the code, not in the problem Iterating, instead of factorising

Coarsening in one direction only

Leave the smoother alone and halve only the strong direction, and the smoothing factor is 0.3340 — identical to line relaxation's, at every anisotropy and every weight, to twelve digits. The convergence factors are then a factor of three apart, and at 45° both repairs fail outright.

 the grid the operator came froman edge is a coupling the matrix calls strongkeptinterpolatedwhat the entries decidedcoupling ratio, x against y0.001strong couplings across x0strong couplings along y110rows kept or dropped whole1111×11 grid, θ = 0.25the coarse grid, from the matrix alone Iterating, instead of factorising

The coarse grid the matrix chooses

Given a tridiagonal matrix and no information about a grid, the coarsening keeps every other point and derives the weights ½, 1, ½ — the operators the geometric method was handed. Given the anisotropic operator, it discovers semi-coarsening, in the right direction, without a coordinate.

 1234567891010⁻⁵⁵10⁻⁴⁹10⁻⁴³10⁻³⁷10⁻³¹10⁻²⁵10⁻¹⁹10⁻¹³10⁻⁷10⁻¹cyclesizethe residualrepairedthe reported boundwhat the stopping rule readsreported at the last cycle9.4·10⁻⁴¹the residual there5.7·10⁻⁵with the border recomputed3.8·10⁻¹⁴products, cheap and repaired9a bound with nothing under itand one product a cycle to fix it Eigenvalues, singular values, rank

Keeping the vectors, and losing the bound

Thick restarting keeps the Ritz vectors instead of filtering the starting vector — the same eigenvalues for a third of the products with A. Its residual bound reaches 9.4·10⁻⁴¹ while the residual it bounds sits at 5.7·10⁻⁵, and the eigenvalues are correct to 4.3·10⁻¹⁴ the whole time, so nothing reports it.

 081624324000.250.50.7511.25angle of the strong direction (degrees)ratio, and convergence factorenergy ratio of oneenergy ratioclassicalminimiserwith the constraintwhere the formula is optimalenergy ratio at 0°1energy ratio at 45°1.2classical rate at 45°0.35the minimiser's rate there0.28a ratio of exactly one while the assumption holdsand a diagnostic when it stops Iterating, instead of factorising

The formula that was already optimal

Ask for the interpolation that minimises the energy of its own columns and the answer is the classical AMG formula — to zero at every row of the one-dimensional Laplacian, and to four digits in two dimensions. On the operator rotated to 45° the two part company, and the gap between them is a diagnostic that needs no reference solution.

 darker is smaller: 10⁻¹⁰, 10⁻⁶, 10⁻³, 10⁻¹one spectrum, two matricesspectral radius0.8reach of the 10⁻³ level1.3eigenvalues, all at0.8the circle is |z| = 1and every eigenvalue is well inside it Eigenvalues, singular values, rank

The eigenvalues that are not there

For a normal matrix the resolvent norm is exactly one over the distance to the nearest eigenvalue, so a picture of it carries nothing the spectrum did not. Move one entry above the diagonal and the region a perturbation of 10⁻⁸ can put an eigenvalue into stops being a disc and reaches out past the unit circle, while every eigenvalue stays at 0.8.

 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.

 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.

 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^{At}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⁻⁸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.

All themes