The thread: Measured, not assumed — page 18

Least squares, and the road not to take

The degree that is safe to overshoot

The rules that choose a Tikhonov parameter miss by factors of millions on one draw in twenty. Transplanted to the degree of a polynomial fit, in a basis orthonormal on the data, the same rules never cost more than 2.7 times the best degree's error in three hundred draws. The reason is the shape of the valley they search: six degrees too few costs from 44 to 16,000 times the best error, forty degrees too many costs about twice it. The one rule with a tail, the discrepancy principle, has its threshold half a standard deviation above the residual it is waiting for.

Elimination, and the swap

Where the multipliers go

Bunch–Kaufman bounds the growth in D and not the entries of L, and the warning attached to that is that everything which later uses the factors inherits the size of L. On a matrix built to make those entries 1.2 over ε, they reach 1.2·10¹⁰ while the solve's backward error stays at 1.6·10⁻¹⁶, |L||D||Lᵀ| stays at 6.7 times ‖A‖, and a step of refinement changes nothing. The large multipliers are where the rule has put the matrix's ill-conditioning. A direction of negative curvature read from those factors finds 7·10⁻¹⁶ of the curvature that is there; the bounded rule's finds 13%.

Exact arithmetic, and what it costs instead

The field decides it, usually

A matrix whose rank depends on the field it is read over was built, the first time, from its invariant factors outward, because random integer matrices never seemed to show the effect. Random 0/1 matrices show it at almost every size that is not tiny. At twenty rows, 99.8% of them are invertible over the rationals, 29% modulo two, 56% modulo three — and 71% of the ones the rationals call invertible are singular modulo two. Modulo two they obey, corank by corank, the law for uniformly random matrices over that field; modulo three and five, which their entries cannot fill, they converge to that field's law anyway.

Randomised, and the guarantee that changes kind

A fade made of drops

A one-nonzero sketch's median error fell by about forty per cent for every factor of ten in how far a coherent matrix had been mixed toward an incoherent one, and nothing explained the rate. Followed one draw at a time, no draw fades at that rate. Each holds its coherent error — as large as the singular value of the direction its hash lost — and then drops, within one to three decades of mixing, never faster than one decade of error per decade of mixing. The median's steady slope is where the drops happen to fall. Change the spectrum and they fall elsewhere: at a decay of 0.9 there is no slope, only a cliff.

Regularisation, and the answer that is chosen

A third penalty on a flat floor

Two penalties at once were worth three per cent at most, and the explanation offered was that one penalty already does the work — which predicts that a third buys less still and that the best choice sits on a face of the parameter cube. The third buys a median of exactly nothing on all five signals and 0.66% on its best draw. But on twelve draws of forty the best triple does use all three, and on every one of them the nearest point with a penalty switched off is within that same 0.66%. The minimum is not on a face; it is on a floor so flat that where it lands is noise. Choosing the right single penalty is worth a factor of two.

Structure, and the solver that cannot see it

Where one step stops being enough

One step of iterative refinement took the squared Laplacian's circulant-wrap solve to elimination's accuracy at every size, and the account was that each step multiplies the error by the wrap's condition number times the rounding, so one step suffices while that product is small. Raised to higher powers, the band tests the account and half of it holds: each step does contract by about κ(wrap)·u. The other half fails. The fifth power at sixteen points needs three steps with κ(wrap)·u near 10⁻⁶, where the third power at 128 points needs one with thirty times more, because the first solve starts up to a thousand times further from the answer than κ(wrap)·u says.

Where the flop count stopped predicting the time

The circle between two eigenvalues

A contour count's error is not approximately governed by the nearest eigenvalue; it is exactly one closed-form term per eigenvalue, and summing those terms reproduces the quadrature to a millionth at 1,720 radius and point pairs. Three things follow. The rate is the ratio of the two moduli the circle sits between, not a distance, so two circles 0.2 from their nearest eigenvalue converge three times apart. Ten digits cost about 21 points divided by log₁₀ of that ratio, 71 points with twelve eigenvalues inside and 3,476 with six. And the best circle is not halfway: at the geometric mean of two eigenvalues on one ray their two terms are equal and opposite, and ten digits cost 27 points where the midpoint needs 84.

Orthogonality, measured

The factor nobody forms

A blocked Householder factorisation's orthogonal factor, multiplied out, departs from orthogonality half as far in blocks of sixteen as one reflector at a time, and that was read as blocking buying a factor of two. Libraries do not multiply it out. Applied to vectors through its stored blocks — which is how every caller uses it — the same factor departs by 3.1 to 4.0·10⁻¹⁵ at every block size from one to sixty-four, and stops growing after about twenty reflectors instead of adding them up. The factor of two was the price of forming the product, and a factor that is never formed never pays it.

Sparsity, and what elimination costs

Pieces ordered blind

One level of nested dissection followed by minimum degree lengthened the factorisation's critical path, and the blame moved from the separator to the halves: minimum degree's path through a 12 × 24 half was longer than through the whole 24 × 24 grid. A 12 × 24 grid on its own has a path of 12,913, a third of the square's. What made the half expensive is that it was ordered as if the separator were not there. Count the separator's vertices in every degree and number them last, and the depth-one path falls from 53,508 to 41,050 — minimum degree's own — while the pieces' own arithmetic does not change at all. At depth two the same ordering beats full dissection.

Methods that were designed apart

The first step's filter, kept for the whole run

Conjugate gradients' first step applies a ramp — Tikhonov's slope up to a corner, then one — and the essay that found it proposed the ramp as a fixed method, predicting that it would match the iteration below half the answer better than any sharpened Tikhonov and lose the match near the answer. The first half holds by more than predicted: 0.044 from the iteration below half, where the best member of the family was 0.106. The second does not: near the answer the ramp is within 0.016, closer than any filter measured. And two numbers read off the iteration's own Ritz values — its tail's slope and its smallest Ritz value — reproduce it past the top as well, where every fixed filter falls ten to seventeen per cent below it.

Reduction, and what a model is for

Where a ringing mode stops being a pair

A lightly damped mode puts two nearly equal Hankel singular values into balanced truncation, and the reduction error was read as twice the larger of the first pair cut — right to two per cent at damping 0.01, wrong by half at 0.1. Read instead as the pair's two values added, the error is right to within six per cent at every order up to damping 0.05, and the pair's split turns out to be twice the damping. What it cannot do is say when the reading fails. The split predicts nothing: the widest split at damping 0.1 has the best reading. The failures come from somewhere the pair does not look — the other dropped modes' response under the peak, and a cut between two singular values that agree to under one per cent, which left a reduced model with a pole a thousand times closer to the axis than any pole of the system it came from.

Where the flop count stopped predicting the time

The short leaf belongs at the bottom

A recursive LU that cuts every dimension at a multiple of the cache's edge moves 14 to 21 per cent fewer words than the pure recursion, except where one leaf per dimension has to be short — and the one-column leftovers cost the most. Moving the short leaf to the top, where it becomes a thin strip of a large product, was predicted to cure that. It makes it worse on all twenty-nine sizes that have a leftover, by up to three hundredths. The measurement that matters more is the edge itself: one column under it the aligned cut keeps most of its saving, one column over it every size costs two and a half to three times the pure recursion. Halving to the same base case barely notices either mistake.

Neither sparse nor dense

One build tells the leaf

The rule that picks a hierarchical matrix's cheapest leaf — halve it while the rank at half the leaf is below a quarter of the leaf — reads ranks that exist only after the matrix has been built at every candidate leaf. Fed instead the rank a code can guess from the tolerance alone, one rank per two decades, it names the best leaf on twenty-three cells of thirty, and on the rank-one kernel it cannot see it stores up to seventy-two per cent too much. Doubling the guess for an oscillation, as earlier measurements suggested, makes it worse: eighteen of thirty. Building once, at a single leaf, and holding that build's mean rank fixed names the best leaf on twenty-nine — one more than the rule that builds at all five — and never stores three per cent more than the least.

When the index is a tuple

A fit with no answer to find

Half of all random 3 × 3 × 2 tensors, and most larger ones, have no rank-three decomposition, because their pencil has a complex pair. A rank-n fit to one of them does not wander and does not stall. Two starts settle at the same error to five digits, and that error is the distance from the tensor to the surface where its pencil has a double eigenvalue — found with no fitting at all, and matched to within one and a half per cent. Meanwhile the fit's terms grow without limit, like the square root of the sweep count, while the fitted pencil's two closest eigenvalues close on each other at exactly the rate the terms grow. The error has an answer; the decomposition does not.

The matrix a constraint makes

One hyperplane hides nothing

Iterative refinement against a saddle-point matrix gives its verdict when the growing direction overtakes the minimum's slowest decay, and a right-hand side with that direction removed delays it by a fixed number of steps a decade. With two negative curvatures, the prediction was that hiding both is paid for by the faster, and that hiding only the faster leaves the slower one's race to run. The first half holds to half a per cent on three pairs, and to one and a half on two whose curvatures differ by a factor of two and of 1.2. The second does not: hiding either direction alone moves the verdict by a fixed handful of steps, the same at two decades as at sixteen, because the direction left in view is already growing. Only the intersection of the two hyperplanes hides anything — and it hides less than either curvature can alone, 1,127 steps at worst against 2,907, with random right-hand sides' slowest at 79 against 417.

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.

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.

The arithmetic underneath

The far well in a byte

Half precision's subnormals carried a two-well chain's far well ten octaves deeper than flush to zero, and the essay that measured it predicted that on the eight-bit formats the band would be worth almost nothing and starting high almost everything. Measured on nine chains, the band is worth what it always was, two or three octaves; starting at the top of the range is worth ten and fourteen; and at three significant bits the far well is lost on eight chains of nine with no exponent limit at all, because the precision fails before the range does.

Elimination, and the swap

A curvature direction the factors cannot refine

A direction of negative curvature read from Bunch–Kaufman's factors held 7·10⁻¹⁶ of the curvature that was there, and the bounded rule's held 14 per cent. The prediction was that a few steps of inverse iteration with the same factors would recover it from either. One step leaves both under a thousandth, and two put both on positive curvature, at the eigenvalue nearest zero — because a solve amplifies the smallest eigenvalue in magnitude, not the most negative. What recovers the curvature is the matrix, not its factors: Lanczos from either direction reaches ninety-nine per cent in six to nine products at every coupling, and power iteration from Bunch–Kaufman's direction has not reached a tenth after sixty.

The matrix that is a graph

Each spectrum hears the other's pairs

Pairs of graphs with the same Laplacian spectrum share a spanning-tree count, so the critical group could only tell them apart by how that count factors — and did on two pairs in five. Turn the census round, to the 733 pairs of eight-vertex graphs with the same adjacency spectrum, and the tree count is free to differ. It does on 717 of them, by a median of six per cent. The critical group separates exactly those 717 and not one more: on all sixteen pairs whose counts agree, including a non-cyclic one, the two graphs have the same group. And the Laplacian spectrum separates all 733, as the adjacency spectrum separated every Laplacian pair.

Least squares, and the road not to take

The scale that only moved a pivot

Multiply the constraint rows of a saddle-point system until its multipliers are the size of its solution, and the extra error the route was blamed for — 4.6·10⁻⁴ against the null-space route's 1.8·10⁻⁸ — falls to 3.3·10⁻⁸. The prediction holds and its reason does not. A scale of ten does what a scale of 6·10⁵ does; hold the elimination's row order fixed and nine decades of scale move the error by less than a factor of five. What the scale changed was which row partial pivoting took at the second step, and taking the constraint rows first does the same job with no scale at all.

Orthogonality, measured

The worst residual belongs to the route

Every stable least-squares route's error rises with the residual at about a fiftieth of the bound κ²uρ, and the question left open was whether a residual aimed at the weak directions closes the gap. It can be aimed exactly: the map from residual to error is, on these problems, one direction and rounding. Aimed, Householder's constant rises from a median of 0.014 to 0.11 — √(m − n) = 6.9 times a typical direction, at 32, 64 and 128 rows to three figures — and stops a factor of nine short of the bound. The direction is each route's own: aimed at one, another route draws a fifth of its worst. And aimed at the appended block, the one-per-cent rule that made it worth its extra block falls to half a per cent, with the block twelve times behind Householder on the same data.

When the problem arrives again

A target the residual can promise

A caller who wants inexact Newton's answer to six digits usually passes 10⁻⁶ as the residual tolerance, and on five problems at four targets that misses by 7 to 35,000 times. Divide the target by a condition number and it is met on every cell where a residual can certify it — and the condition number need not be known: the Ritz values conjugate gradients produces on the way give the Jacobian's own, settled by the sixth Newton step and within eleven per cent of the truth at the end, and two to nineteen times below the linear part's. Stopped that way, with the usual floor under the forcing term, the run lands within a fifth of its target and is cheaper on 16 of 19 cells than the unfloored run stopped by someone who knows the answer.

Structure, and the solver that cannot see it

The symbol that builds the staircase

A block-circulant preconditioner left the isotropic kernel's spectrum a ramp of thirty-eight clusters where the separable kernel's was a staircase of ten, and the proposal was a circulant sampled from the isotropic kernel's own symbol, to gather the ramp into treads of its own. It gathers the separable kernel instead: six clusters at every grid from 4 to 10 and every correlation, sixty-four of a hundred eigenvalues at exactly one, and eight conjugate-gradient steps at every side from 6 to 16 while Chan's circulant climbs from 18 to 28. On the isotropic kernel it shortens the ramp from 38 clusters to 30 and the step count by about a sixth, and both still grow with the grid. The treads come from a Kronecker product, not from a small difference: on both kernels the matrix and its symbol circulant differ in dozens of directions.