Numerical Linear Algebra
where the algebra stops being the arithmetic
Theme

The thread: Whose fault is it — page 13

Essays 289 to 303 of the 303 on this theme, in the same order.
smallest groupsecond groupthird grouplargest groupfirst form, leadinglost0.030.046·10⁻¹⁶first form, trailing7·10⁻¹⁵0.03lostlostreversed, leading7·10⁻¹⁶0.040.04lostsecond form, leading0.030.03lost7·10⁻⁹second form, trailing3·10⁻¹²0.04lostlostscaled at the second root, leadinglostlostlost8·10⁻⁸scaled at the second root, trailing4·10⁻¹⁶10⁻¹³lostlostscaled at the third root, leadinglostlost3·10⁻¹⁴6·10⁻¹⁶scaled at the third root, trailing3·10⁻⁸lostlostlostgreen: under 10⁻¹² · blue: 10⁻⁹ · violet: 10⁻⁵ · red: worse · bold: a seed's count was not fouran interior root needs its own scaling The eigenvalue problem that is not linear

An interior root needs its own scaling

A cubic matrix polynomial with three groups of eigenvalues needed only two reductions: the leading one kept the middle group at every separation. A quartic has two interior groups, and the prediction was that the leading reduction would keep both, so that two reductions would still do. At eight decades between groups it keeps neither — both interior groups are wrong in the second digit — and no unscaled route keeps the second group at all. Two reductions still suffice, but each must be scaled at an interior tropical root: the trailing one at the second root for the lower two groups, the leading one at the third for the upper two. And the second companion pencil, predicted to keep the middle groups as the reversed polynomial does, loses them as the first pencil does.

 inner iterations in allfirst rule462second rule690second, capped675123456789101112130100200300400500600700Newton stepinner iterations so farfirst rulesecond rulesecond, cappeddotted: the step before the residual rosethe gap opens before the reset When the problem arrives again

The step before the reset

On the strongly nonlinear problem, inexact Newton's second forcing rule took two more Newton steps than the first and half as much work again, and the blame went to its safeguard, which let the forcing term jump back to 0.9 when the residual rose. Capping that jump at ten times the last forcing term was predicted to remove the two steps and bring the rule within five per cent of the first. It removes the steps, at three targets of four. It saves 0.7 per cent of the work, and leaves the rule 45 per cent dearer than the first. The reset was the cheap part: four loose steps for 55 inner iterations. The cost was the step before it, solved to a forcing term of 0.024 for 109 iterations, where the first rule spent 5 — and on the ill-conditioned problem, where the residual never rose at all, the cap costs ten per cent.

 spread used, % from the truethe last trace's spread, trace 10-12times the last ratio, trace 10-25guarded at a turn, trace 10-180.70.750.80.850.90.9511.051.11.15trace in the sequencespread used ÷ true spread48101216the last trace's spreadtimes the last ratioguarded at a turndotted: the turn · dashed: the true spreadthe turn's direction is the error's sign Randomised, and the guarantee that changes kind

A turn is seen one trace late

A carried spread multiplied by the ratio of the last two spreads restored a drifting sequence's coverage, and the prediction for a sequence whose drift turns was a dip of four or five points on the trace after the turn, gone a trace later, and a guard that removes it for nothing. Run on the turn the prediction described — the spread rising and then falling — there is no dip: the extrapolation promises more rise, uses a spread eleven per cent too large, and covers 70.5 per cent against a table of 68.3. Run on the mirror turn, the spread falling and then rising, the dip is there and larger than predicted: a spread a quarter too small and 59.5 per cent coverage. Either way it is gone a trace later and the sequence's average stays within a point. The guard costs nothing and does not help, because no rule that reads past spreads can see a turn before the trace that shows it.

 longest ÷ shortesttwo rows, 12²1.5two rows, 24²1.311.11.21.31.41.5longest ÷ shortest critical path12²16²20²24²12²16²20²24²g1g2g312²16²20²24²grid, one rowgrid, level setgeometric graphgrid, two rowsdashed: every shuffle gives the same factorfree unless the separator is thick Sparsity, and what elimination costs

A clique has no order

Ordering a dissection's pieces with their separator in view cut the critical path by a quarter, and the next proposal was to order the separator's own vertices by how much fill each carries in. On every separator the dissections here produce — the grid's rows, its level-set diagonals, and the level sets of three irregular graphs of 300 vertices — twenty-four shuffles of the separator's order give the same factor to the last entry. The reason is one line of graph theory: once the pieces are eliminated, a separator whose vertices all touch one connected piece is a clique, and a clique's order changes nothing. The order matters only on a separator that is not a clique, and the one built here to be that — two grid rows thick — loses to a separator of one row under every order tried. The separator's size is the decision; its order is not one.

 at ε = 10⁻⁴tread at one, separable ÷ mixture7.7-6-5-4-3-210⁻⁷10⁻⁶10⁻⁵10⁻⁴10⁻³10⁻²10⁻¹1fraction of the isotropic kernel, log₁₀relative width of the treadtread at one, separable symboltread at one, mixturetread near 6.5, separabletread near 6.5, mixtureevery line has slope one: width in proportion to εthe width moves to the crowded tread Structure, and the solver that cannot see it

The part a circulant can cancel

A separable kernel with a small isotropic part was preconditioned by the circulant sampled from the whole mixture's symbol. The cheaper choice keeps only the separable part's symbol and leaves the isotropy as error, and the prediction was that it widens the same treads by about the same amount — the same step count to a thousandth of isotropy — and falls behind past a hundredth. The step count is the same to within one on every grid up to ε = 10⁻³, and the cheaper preconditioner is behind on every grid from ε = 0.1. But it does not widen the same treads. It puts its error on the tread at one, where sixty-four of the hundred eigenvalues sit, widening it 7.6 times as fast as the mixture's symbol does, and it leaves the tread the mixture widens most nine times narrower. The counts agree because conjugate gradients pays for the widest tread's logarithm, not for which tread it is.

 1251020shifted symmetricsmallest pivot, plainsmallest pivot, diagonal firstharmonic mean of the spectrumGram, 48 rowssmallest pivot, plainsmallest pivot, diagonal firstharmonic mean of the spectrumGram, 30 rowssmallest pivot, plainsmallest pivot, diagonal firstharmonic mean of the spectrumover the smallest eigenvaluethirty matrices a family, logarithmic axisthe pivot follows the mean, not the bottom Elimination, and the swap

A definite pivot sits at the mean

On a positive-definite matrix every Cholesky pivot is at or above the smallest eigenvalue, so a shift read off the smallest pivot always has to be pushed down, and it was predicted to need a median of four pushes and to be usable unpushed on none. The second half holds on all ninety matrices measured. The first does not: shrinking by a quarter at a time takes a median of six to ten steps, because the smallest pivot sits within a factor of two of the spectrum's harmonic mean, and the harmonic mean sits four to nine times above the smallest eigenvalue. The push was a device for the indefinite case, which has no floor to aim at. A definite matrix does, and halving towards zero gets below the spectrum in a median of two to four.

 0816243240485664reflectors appliedworst ‖Qᵀ(Qx) − x‖02·10⁻¹⁵4·10⁻¹⁵6·10⁻¹⁵8·10⁻¹⁵96 rows192 rows384 rowsmedian of five matrices, 64 columns eachlonger vectors, a higher plateau, reached as soon Orthogonality, measured

The plateau is set by the length

A Householder factor applied through its reflectors, never formed, loses orthogonality over its first twenty or so reflectors and then stops. The saturation was predicted to be a property of the vector length: on 384-long vectors the growth should continue longer before it stops. It does not. On 96-, 192- and 384-long vectors the departure reaches nine-tenths of its plateau within twenty reflectors, and on the longest it gets there soonest. What the length sets is the plateau itself — 3.6, 5.0 and 6.7·10⁻¹⁵ — which is the unit roundoff times the square root of the length, times a constant between 3.1 and 3.4. Twice the reflectors on the same length move it by nine per cent.

 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.

 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.

 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.

 10²10⁻⁴10⁻³10⁻²10⁻¹grid points nworst nodal errorupwindtunedcentralthe same tuning, another problemtuned ÷ central at n = 3149tuned ÷ central at n = 12796central's error at the finest grid10·10⁻⁵exact on the problem it was derived fromand harmful on the one beside it Iterating, instead of factorising

A parameter that is also a price

ξ = coth(Pe) − 1/Pe is the fraction of h/2 that makes a boundary-layer solution exact at every node. On a problem with no layer in it, the error the same scheme commits is ξ times upwinding's — 0.2511 against a ξ of 0.2504, 0.7461 against 0.7448 — so the number that buys the exactness is also the invoice.

 0183654729010812610⁻¹⁶10⁻¹³10⁻¹⁰10⁻⁷10⁻⁴10⁻¹products with Asizeresidual boundtrue errorbounded memorybasis vectors kept8products with A140worst error in the k wanted7.1·10⁻¹⁵the bound is free and the error is notand the basis never grows Eigenvalues, singular values, rank

The same budget, spent five ways

A restarted method has one budget — products with A — and two ways to spend it, in many short cycles or a few long ones. At about a hundred and forty products the answer is the same to a factor of seven whichever split is chosen, and the residual bound the method reports spans ten orders of magnitude across the same five runs.

 -213284300.250.50.751rotation of the anisotropy (degrees)factor / couplingusableconvergence factoraxis couplingdiagonal couplingthe standard answer, and the anglefactor at 0°0.19factor at 45°0.8axis ÷ diagonal coupling at 45°2the hierarchy reads the matrixand the matrix lost the direction Iterating, instead of factorising

How much direction there was to lose

At 45° the nine-point stencil hands smoothed aggregation the same wrong hierarchy at every anisotropy — six strong neighbours per interior point, 121 aggregates, the identical partition from ε = 10⁻⁴ to 0.099. The convergence factor that one hierarchy produces runs from 0.802 to 0.581 over the same range.

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

All themes