The thread: Exact ground truth — page 7
The honest residual is not a stopping test
A conjugate-gradient run in reduced working precision reports a residual pinned near 10⁻¹² while its error spans eleven orders, and the repair was to recompute the residual in double. As a report it works: on thirty runs it lands within a factor of a quarter to five of the error. As the thing to stop on it fails outright: below 53 bits it is passed on none of twenty-five runs, which go on to the 300-step cap or until their own recurrence collapses — 7,694 products in all against 605. One clause fixes it. Stop also when two checks in a row fail to halve the double residual, and every run ends within thirty-six products, at the same error, reporting a number that means something — and at 8 and 11 bits on the larger grids, sooner than the recurrence's own test.
The eigenvalue problem that is not linearTwo groups need two reductions
When a quadratic eigenproblem's tropical roots are far apart its eigenvalues fall into two groups, and the remedy offered is a second scaling — one per tropical root, each run keeping the group its root predicts. Raise the damping on an eight-mass chain until the roots are 10²³ apart and measure every eigenvalue: no scaling rescues the small group. Its backward error grows as the unit roundoff times the separation under Fan–Lin–Van Dooren's scaling, under both tropical scalings and under none, to 10⁻⁵ at a separation of 4·10¹¹. What rescues it is the other reduction of the same pencil — inverting the constant term instead of the leading one — which keeps the small group at rounding at every separation and loses the large one instead. Both reductions, each keeping its own group, give all sixteen eigenvalues to 10⁻¹¹ at worst. On the way, the closed form that served as the exact answer turned out to lose the small roots to cancellation, from the same separation.
Exact arithmetic, and what it costs insteadAn order found once knows what the rows know
Searching for a fraction-free elimination's least-cost row order costs forty to two hundred and sixty times the elimination, so it can only pay if one order serves a whole family. The prediction was that it cannot: the best order depends on which rows happen to have short entries, so a reused order should do no better than the natural one. On independent entries and on a shared zero pattern that holds — another member's best order is an ordinary order, at the middle of the 40,320. On a family that shares its row scales it fails completely: another member's order costs a median 1.13 of the least where the natural order costs 1.88. But on every family, including that one, the rule that just takes the smallest pivot entry is cheaper than any reused or trained order — 1.24, 1.08 and 1.01. What an order carries from one matrix to the next is what the family shares — and only when that is a particular matrix, not a visible structure, does an order trained on the family beat the rule.
The matrix a constraint makesA square root the residual does not pay
With the exact Schur complement, the triangular saddle-point preconditioner puts every eigenvalue at one in Jordan blocks of size two, and GMRES finishes in two steps. Solve the first block only to a relative accuracy τ and the blocks split: their eigenvalues leave one like √τ — 0.42 to 0.62 times √τ on eighteen systems across six decades, each pair exactly where a five-by-five eigenproblem puts it — while the rest move like τ. The question left open was whether that square root reaches the iteration. It does not. GMRES's residual after two steps grows like τ, with a slope of 0.86 to 1.27 on every system, because its polynomial has a double root at one and cancels each split pair to second order; every further pair of steps divides the residual by about the same factor; and a rule written in τ alone predicts the step count exactly on 179 of 234 runs and within one on 228.
Eigenvalues, singular values, rankAn 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.
The matrix that is a graphThe group hears a switch only at two
On a regular graph the Laplacian is the degree times the identity less the adjacency matrix, so the two spectra are one piece of information and the tree count comes with them; between cospectral regular graphs the critical group is the only exact invariant left. Made by Godsil–McKay switching on four vertices, 705 distinct cospectral pairs of 3- and 4-regular graphs on 12 to 20 vertices give its answer. It separates 181 of 587 quartic pairs and none of 118 cubic ones. And on every one of the 705 the two groups agree everywhere except at the prime two: a switch is a conjugation by a matrix of halves, the group can see it only through its 2-part, and a group whose tree count has fewer than three factors of two cannot see it at all.
The answer that depends on the machineFive steps is a rate
A stagnation clause that checks the true residual every five steps, and stops when two checks in a row fail to halve it, stopped every reduced-precision conjugate-gradient run at the error it could reach. The worry was that five steps is a property of one fast problem, and the prediction was that on the unpreconditioned Laplacian, several times slower, it would stop some runs early. It does not: unpreconditioned conjugate gradients still halves its residual every five steps on these grids, and the clause stops all twenty runs within a tenth of the least error. Scale the same Laplacian's rows by a decade and a half and it stops nineteen of twenty at step fifteen, up to 775,000 times the attainable error. A clause that sets each interval from the reduction it last saw stops every run of all three kinds within two fifths of it — and on the fast runs pays for that with eight to twelve times as many steps past the best.
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.
Reduction, and what a model is forAn error estimate made of samples
A Loewner model built from samples of a transfer function comes with no error bound, and the number every code reads to choose its order — the pencil's next singular value — is the obvious stand-in. On four systems sampled at 48 points it is the error of the order-k model for the first two orders, within a factor of 2.5, and then falls away from it by about a decade an order: by the fifth it is optimistic by 109, 181 and 457, and on one system by 6,700 at the sixth. An AAA fit of the same degree to the same samples has an error within a factor of six of the Loewner model's at every order, so the shortfall is not the projection's. The model's misfit at its own samples is optimistic the same way. What works is more samples: the error at the 47 midpoints between the samples is within a factor of 2.2 of the true error at every order on every system, and the four midpoints nearest the poles alone within a factor of three.
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.
Eigenvalues, singular values, rankA problem with no answer
If two matrices share a null vector then det(A − λB) is identically zero and every λ is an eigenvalue, which means none of them is. Perturb such a pencil by a ten-billionth and a solver returns six numbers with residuals below 10⁻⁹. Change the seed and it returns six different numbers, spread over forty-four, with residuals just as small.
Structure, and the solver that cannot see itA backward error the answer does not feel
Both ways of solving AX + XB = C leave a residual at the rounding level, and for a linear system that would be backward stability. For the equation it is not: the smallest perturbation of A, B and C that makes the computed X exact can be far larger. On a family built to make it so — X increasingly ill conditioned, C increasingly small beside AX — the structured backward error of Bartels and Stewart's solution climbs to 2.1·10⁵ times its residual, and the Kronecker system's, whose residual is four times smaller, to 2.0·10⁶, within the bound Higham's μ sets. And the error of the answer does not move. Against the exact X it stays where the size of C puts it, 6·10⁻¹⁵ to 1·10⁻⁹, across eight decades of κ(X). The perturbation that makes the backward error large has to sit on C, along directions X is insensitive to, and so the backward error rises and the answer never feels it.
Two errors, and whose fault they areA basis that has to know where the roots are
Wilkinson's polynomial loses its roots in the third figure at degree twenty when they are computed from its monomial coefficients, and the repair offered was to change the basis: in Chebyshev coefficients, with the colleague matrix, the conditioning 'does not degrade exponentially with degree'. Expanded exactly in that basis and rounded once, the integers' roots come back ten orders better at degree twenty — 4.4·10⁻¹³ against 7.6·10⁻³ — and still degrade exponentially, at a third of a decade a degree, to 2.8·10⁻⁴ at forty. Place the same number of roots at the interval's Chebyshev points instead and they stay at 10⁻¹⁴ at every degree, while the monomial route loses them as fast as the integers. The basis does not repair the problem; it has to match where the roots are.
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.
When the problem arrives againA correction that reads every row
A sliding least-squares window moves its answer by the update its entering and leaving rows imply and then corrects it once against all 24 rows it holds. The proposal was to read only the rows of highest leverage, a quarter of them, predicted to recover most of the full correction on a stationary stream. On every noisy stream it does the opposite: the answer ends with a relative error of order one — 2.8 at a noise level of one, where reading every row leaves 3.3·10⁻⁹ and not correcting at all leaves 9.4·10⁻⁴. Choosing the rows at random, by age or in rotation changes nothing. Leaving out one row of 24, the one of least leverage, still leaves 0.29. A correction is a Newton step whose fixed point is the answer only because the residual is orthogonal to all the columns summed over every row; read from some of the rows, its fixed point moves by the omitted rows' share of the residual, and the window settles there. Only on a stream with no noise, where every row's residual is zero at the answer, does a partial correction work.
Elimination, and the swapA doubling that arrives two steps late
A pivot rule that looks one elimination ahead chooses Wilkinson's growth row for row, because the doubling it sets up arrives one step later; looking two ahead escapes it. The argument said a matrix whose damage arrives two steps late would defeat the two-step rule the same way. It exists, and it is two Wilkinson matrices interleaved: W ⊗ I₂. Every look-ahead of depth one or two reproduces the greedy rule's growth on it exactly, 4 to 64 from n = 6 to 14, and depth three gives 2. Three copies defeat depth three and yield to four. The growth each depth is held to is 2^(n/d − 1): a fixed depth stretches the exponential and never removes it, at a cost that rises by a factor of about n for every step it looks.
Least squares, and the road not to takeA unit is a statement about the noise
Total least squares minimises one Frobenius norm over every column of [A b], so the unit a column is written in is a claim about how noisy it is. On a 60-row fit with two measured regressors, rewriting one column in units from 10⁻⁴ to 10⁴ times the recorded ones leaves ordinary least squares exactly where it was and moves total least squares by half again, continuously, between two estimators with their own names: the reverse regression of that column on the others, and the fit that treats it as exact. Its correction is split among the columns as the squares of the coefficients, which the units set and the noise never enters. Dividing each column by its noise level makes every unit give one answer, the best of five estimators on all three placements of the noise tried, and two replicate readings a column are enough to get most of the way there when the recorded units were badly wrong.
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.
Reduction, and what a model is forA pole between the samples
A Loewner model needs a check on its error, and the one predicted to be cheap is four points chosen where the model disagrees most with the model one order higher. With exact samples it is better than predicted: on four systems the true error at those four points is the model's error within two per cent at every order, and the disagreement alone, with no sample at all, within three per cent at every order but the last. With noisy samples the checks keep working on every model that has no pole between the samples, and seventeen of thirty-two models do. Past the order the noise allows, the projection puts a real pole inside the sampled interval, so the error is infinite at one point and small a few per cent of a gap away, and every finite set of checks reports a number between 130 and 49,000 times too small. The test that catches every such model costs nothing: the model's own poles.
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.
The matrix a constraint makesThe multipliers were a set
A crossover that solves the equality problem an interior-point iterate's diagonal suggests, then certifies the answer, refuses a guess with more active rows than unknowns — which a degenerate vertex produces as a matter of course. The two remedies a simplex code and a regularised solver would reach for, choosing a basis or regularising the extra rows, both return the right point from the first or second iterate, and both fail the certificate on most programmes, at almost every iterate. The point is never what fails. The certificate checks the one multiplier vector a solve happens to return, and at a degenerate vertex the multipliers are a set; asking whether the set holds a nonnegative vector certifies the same point at once.
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.
The matrix that is a graphThe stall was in the block
On two dense clusters joined by a coupling ε, the power iteration needs about 7/ε steps and then stalls at about u/(4ε). Aggregation–disaggregation was predicted to need a number of steps independent of ε and to stall at rounding, and with a power step as its smoother it does: fourteen steps at every coupling from a tenth to 10⁻¹², every entry to a relative 10⁻¹⁵. The textbook smoother, block Gauss–Seidel, is faster and stalls at u/ε again — because each block solve forms one minus a self-probability, and a block whose rows sum to 1 − ε is singular to within ε. Solved without that subtraction, the same blocks finish in one to four steps at rounding, and get faster as the coupling weakens.