The thread: Two routes to a number — page 9
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.
Elimination, and the swapThe push a pivot needs belongs to the matrix
An indefinite factorisation's most negative pivot looks like an estimate of the most negative eigenvalue, and on one planted family Bunch–Kaufman's, pushed a tenth further, was a usable shift for inverse iteration. On thirty random symmetric matrices neither rule's pivot is: unpushed it is nearest the smallest eigenvalue on 12 and 5 of them, pushed a tenth on 13 and 9, and only a push of two serves 29. The push each matrix needs runs from 0.59 to 2.27 for Bunch–Kaufman — whose pivot is more negative than every eigenvalue on eight matrices and under half the smallest on others — and from 0.85 to 2.01 for the bounded rule. What does work needs no push at all: multiply the shift by a quarter until the factorisation of A − σI has no negative pivot, which Sylvester's law says is exactly when σ is below the spectrum. It lands every time, in a median of two or three factorisations.
When the problem arrives againTwo forcing terms, one floor
Inexact Newton's adaptive forcing rule comes in two versions, and every measurement of the floor under it used the second, which reads the ratio of the last two residuals. The first reads how wrong the linear model was at the last step, and it was natural to hope that a rule watching the model's own error would not solve its last step past what the stopping test needs. It does. Run without the floor on five problems at four forward targets, it over-solves its last step by a median of 34 times, against the second rule's 95; the floor saves 17.5 per cent of its inner work, against 16.0. Floored, the first rule is the cheaper by 4.6 per cent overall, and it buys that with Newton steps: two more on the median cell, a third less inner work on the strongly nonlinear problem, and 4.5 per cent more on the easy one.
Structure, and the solver that cannot see itSix steps were six eigenvalues
A circulant sampled from a separable kernel's symbol held conjugate gradients at eight steps on every grid, because the preconditioned spectrum was six values. Add a fraction ε of the isotropic kernel and the prediction was that the six treads would widen, the count would stay near eight while they stayed under one per cent, and then climb towards the isotropic count. A ten-thousandth of isotropy leaves six clusters at the one-per-cent rule and takes the count from 8 to 13–15; a millionth takes it to 11. The eight steps were finite termination on six distinct eigenvalues, not convergence on six clusters, and once the eigenvalues are distinct what a tread costs is set by its width against the solver's tolerance: at a tolerance of 10⁻⁴ a ten-thousandth of isotropy costs two steps. The widths grow in proportion to ε, and summed tread by tread they account for the climb to within four steps. Past a tenth the mixture is harder than either kernel alone — 37 steps on a 16 × 16 grid against 35 for the isotropic kernel — while its condition number is lower than both.
Where the flop count stopped predicting the timeA centre that chooses which pair cancels
A contour count's error is one closed-form term per eigenvalue, and two real eigenvalues on one ray cancel at their geometric mean — about the origin. Moved along the real axis, the centre was predicted to buy a further factor of two by putting real eigenvalues on two rays, the far one cancelling with a sign flip. The flip is parity: an eigenvalue behind the centre changes sign only at odd point counts. The factor is not two. Counting the delay problem's first two and three real eigenvalues, the best centre costs 22 and 21 points against the origin's 41 and 78; counting all four, 32 against 341, because about the origin the circle has to pass between a real eigenvalue and a complex pair and about 2.3 it does not. The centre does not add a cancellation. It chooses which pair cancels and pushes everything else away.
Methods that were designed apartA power the spectrum does not set
Conjugate gradients' filter past the answer is a ramp raised to a power m between 1.27 and 1.50, read off its own Ritz values, and the power was predicted to be the blur's: computable from the decay of the singular values to within 0.1, so the iteration's filter could be built without running it. It is not. Summed over the singular values above the iteration's own corner, the decay gives anything from 0.76 to 4.3, because the Ritz values are not those singular values — fewer of them early in the run, each standing for a cluster, and more of them late, the leading ones found again and again. The power is the spacing between the smallest Ritz value and its two neighbours. And nothing depends on it: placed by GCV's first minimum, the power ramp at every m from 1 to 3 lands within a tenth of its best on 63 to 66 draws of 72.
Exact arithmetic, and what it costs insteadEvery entry under the determinant
The Hermite form's swell was put down to its transform — the record of how the answer was reached — while the form's own entries stayed bounded. Instrumented separately, the form's intermediates swell exactly as much: 3,371 bits against the transform's 3,366 for a random 48 × 48 matrix whose determinant has 219. Carried out modulo the determinant instead, every entry stays within the determinant's width, the answer is identical on every matrix to n = 80, and it is almost always the identity with the whole determinant in its last pivot and last column. What the bound does not buy is speed: counted in bit operations the modular route does more work until just past n = 72, because its steps multiply numbers the determinant's size where plain elimination multiplies wide numbers by small quotients.
The matrix a constraint makesA tolerance is not an accuracy
A fixed perturbed block predicted what an inexact inner solve would do to a triangularly preconditioned saddle-point system: a two-step residual growing like the inner tolerance, a constant a few times the model's, and a step rule exact to one step on nine runs in ten. A real inner iteration under flexible GMRES delivers something else. Conjugate gradients stopped at a residual tolerance hands back errors ten to a billion times smaller than it was asked for, the outer residual follows that error with a slope near one and a constant ten to a hundred times the model's, the rule counts four runs in five — and the real solve still costs no more outer steps, while loosening it never once saved inner work.
Iterating, instead of factorisingThe 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.
Randomised, and the guarantee that changes kindA drift put back at the table's price
A trace estimator that carries its spread from the trace before loses one step of drift, and the drift is visible in the spreads themselves, so the repair proposed was to multiply the carried spread by the ratio of the last two. Over 400 drifting sequences it works: averaged over the third to sixteenth traces it covers 67.8 to 68.9 per cent at one standard error at every drift, where the plain carry falls to 65.4 and the table says 68.3. It costs what the table says closing that gap costs — 10.5 per cent more products at the fastest drift against a price of 11.9 — so at that drift it spends more than a fresh pilot every trace. Its factor is noisy, scattering the spread it uses two and a half times as widely as the plain carry does, and on a sequence that does not drift that noise costs nothing at all. And the single-trace comparison the proposal named turned out to be smaller than its own draw noise.
The matrix that is a graphNo subtraction, and no waiting
The power iteration keeps a positive vector's small entries correct because it never subtracts, and it pays for that in steps. On a Markov chain made of two clusters joined by a coupling ε it needs about 7/ε of them — 6,900 at ε = 10⁻³, 69,000 at 10⁻⁴ — and then stalls at a fixed point of its own rounded arithmetic about u/(4ε) from the answer; at ε = 10⁻⁶ it is still 45 per cent wrong after 200,000. Gaussian elimination needs no waiting and subtracts: on the same chains its error grows like 1/ε, and on a chain whose entries fall to 10⁻¹⁸⁴ it returns the small ones as noise, eight of them negative. The Grassmann–Taksar–Heyman elimination does neither. Every entry of every chain comes back to a relative 10⁻¹⁵, in a third of n³ operations, because the stationary vector is well conditioned under the perturbations it makes.
Methods that were designed apartThe copies cost steps, not the filter
Past the answer, conjugate gradients holds more Ritz values than singular values above its corner, because rounding lets it find the leading ones again and again. Those copies were a suspect for why the filter's power could not be read off the spectrum. Reorthogonalising every basis vector removes them, and the run reaches the same iterate — the same best error to a per cent — in up to four times fewer steps. The power barely moves, the operator's sum still misses it on 62 draws of 72, and without copies there are still fewer Ritz values than singular values above the corner: the clusters, not the copies, broke the prediction.
Exact arithmetic, and what it costs insteadThe swell was one-sided
The Smith form can be computed modulo the determinant, and it was predicted to return the right invariant factors, to stay within twice the determinant's width, and to cost about what the modular Hermite form costs instead of the far larger count of the plain Smith routine. The first three hold. The comparison does not: the plain Smith routine never needed a far larger count. With column operations as well as row operations, its smallest-entry pivot is a unit after a step or two, every entry it holds is a minor of the matrix, and nothing it writes is wider than the determinant — at every size to sixty-four, at a fiftieth of the modular route's work. The swell the Hermite form suffered belonged to doing elimination with one hand.
Where the flop count stopped predicting the timeCount the upper half, read the lower
Moving a contour's centre off the real axis was predicted never to pay for a count closed under conjugation: one member of each pair comes closer, the other goes away, and the slower one always gets slower. For one circle the prediction holds at every offset measured, from 52 points on the axis to 567 three units above it. But one circle is not the only way to count a set like that. On the delay problem, the four real eigenvalues on a circle of their own and each band's upper cluster on a small circle of its own cost 55 points against 52 for one band, 68 against 97 for two, and 115 against 275 for six. What pays is the mirror: for a problem with real coefficients the count about the conjugate centre is the conjugate of the count already made, equal to 1.6·10⁻¹⁵, so the lower half of the plane costs nothing.
Eigenvalues, singular values, rankSmall 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.
Two errors, and whose fault they areThe second factor was the last row
Wilkinson's polynomial in Chebyshev coefficients loses its roots exponentially from the colleague matrix, and two factors accounted for the loss: the roots' own conditioning, and the size of the matrix's last row against the rest of it. A diagonal similarity by powers of two — balancing, which changes no eigenvalue and rounds nothing — removes the second factor entirely: the balanced error sits within a factor of two and a half of the conditioning at every degree from fourteen to forty, three and a half thousand times better than unbalanced at forty. What is left is the conditioning, still exponential. Balancing the monomial companion matrix takes its norm down by forty-five orders and gains nothing, because there the error was never the matrix's.
Least squares, and the road not to takeThe weights the fit chooses for itself
Total least squares needs each column's noise level, and the tempting shortcut is to read it off the fit: weight, fit, take each column's correction as its noise, reweight, repeat. The prediction was that this converges, to weights that depend on where it starts. It does not converge at all. Each step multiplies every weight by its own column's share of the smallest singular vector, the one self-consistent weighting is unstable — a departure from it doubles every step — and every run leaves for the edge in two to seven steps, landing on the reverse regression or on ordinary least squares by which side it started. The self-consistent weights are nearly the same whatever the noise is: they describe the fit's geometry, and only replicate readings describe the noise.
Randomised, and the guarantee that changes kindPromotion buys rank, not safety
A range finder that grows its basis from cheap sparse probes and certifies it with dense ones was made safe on a coherent matrix by promoting a refusing dense probe into the basis — ten promotions for ten coherent directions, none for none. On matrices mixed between the two, the promotions fade as predicted: ten at the coherent end, nine at a mixture angle of 0.1, none from 0.8. But the danger they guard against does not fade with them. It is gone at an angle of 0.05, where the rule that never promotes already certifies every draw correctly. What promotion buys there is rank — a basis of 30 columns instead of 38 — for four fifths more probe work, and the threshold that decides it is a price, not a safety margin.
Where the flop count stopped predicting the timeThe prefix sits further down
A beam over contraction orders cuts twice at every level — once among the pairings out of each partial order, ranked by the size of their result, and once among the partial orders that survive, ranked by their cost so far — and one tolerance had been widening both. Measured separately, the exhaustive order's pairing is the size ranking's first choice on 34 of 60 seven-tensor networks; its prefix sits a median factor of 2.3 above the cheapest prefix on offer, and needs more room than its pairing on 49. So the width belongs at the order cut, as predicted. What was not predicted is the other half: widening the pairing cut alone is worse than not widening at all, because it hands the cost-so-far ranking more chances to choose a cheap-first prefix.
When the index is a tupleThe penalty walks the same curve
A rank-two fit to a random 2 × 2 × 2 tensor of rank three has no best answer: its error falls towards the tensor's distance to the boundary along a law, d + C/s², while its terms grow without limit. A penalty on the term sizes restores a best answer, and the prediction was that its error at terms three times the tensor's norm would sit within ten per cent of the law's d + C/9 — on the assumption, unproved, that the unbounded fit's path is already the best trade between size and error. On twelve tensors the penalised fit's error at that size is within 1.1 per cent of the law, and every penalised fit with terms past two and a half times the norm lies on the same curve. The assumption holds as far as a penalty can test it. And the penalised fit needs tens of thousands of sweeps to settle where the unbounded fit arrives in at most 144, so the cheap way to a bounded fit is to stop the unbounded one when its error reaches the law's price — not when its terms reach the bound, which on four tensors of twelve happens before the fit has found the boundary at all.
Iterating, instead of factorisingThe 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. The guard against it has been here from the first essay 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.
Methods that were designed apartThe runs part before the copy
Conjugate gradients by its short recurrence and the same iterates computed with every basis vector reorthogonalised were predicted to have the same filter power at the same step, to within 0.01, until the short recurrence found its first copy of a converged Ritz value. Before they part they agree to eleven digits, not two. And they part three to seventeen steps before the copy, at a step the blur decides and the noise does not: one to two steps after the short recurrence's basis loses orthogonality past 10⁻⁴, four to eight after the leading Ritz value converges. The difference starts at the bottom of the spectrum, where the power is made. The copy is the last thing to happen, not the first.
Least squares, and the road not to takeThe rule that divides by one minus the leverage
On a smooth function the rules that choose a polynomial fit's degree were forgiven by a lopsided valley: too few degrees cost thousands of times the best error, too many about twice. On |x| and on a function with a pole just outside the interval, the left wall was predicted to flatten, the best degree to climb, and the information criterion's tight threshold to become a liability. All three happen. The wall falls to a twentieth more, the best degree on |x| runs to the largest the basis has, and the criterion is the costliest rule on |x| at one per cent noise. But the rule that actually fails is leave-one-out: at small noise it stops at degree 26 and costs twelve times the best error, because by degree 26 the leverage of the end points is 0.97 and its formula divides every residual by what is left.
Structure, and the solver that cannot see itThe 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.