The step that returns the digits

Where the matrix is worst, and where the answer is not swept a ladder filter’s nodal matrix along the frequency axis and found its condition number predicting nothing. The solution vector was right to three units of round-off at every frequency; the response taken out of it lost four decades into the stopband, and the loss was predicted by one quantity that mentions no matrix: the round-off divided by the response. A small component of a well-solved vector, it concluded, is a small number carrying the absolute error of the large ones.

That is a diagnosis, and it invites a repair. The digits the arithmetic did not have applied the obvious one to a different network — more digits, a pair of doubles carrying 31 — and found it recovered the bridge’s lost digits completely. Thirty-one digits is the expensive cure, twenty times slower. There is a cheap one that numerical analysis has used since the first computers: solve, compute how wrong the answer is, solve for the correction, add it. Iterative refinement costs one matrix–vector product and one pair of triangular solves with a factorisation already in hand.

Applied to the ladder filter, one step of it in plain double precision returns every digit. Applied to the bridge, any number of steps returns nothing. The difference is a precise statement about where each network’s digits were lost.

Solve, check, correct #

The step is short to state. Having solved $A x = b$ for an approximate $x$, compute the residual $r = b - A x$ — how far the approximate solution is from satisfying the equations — then solve $A\,d = r$ for a correction, using the elimination already done, and replace $x$ with $x + d$. If the elimination were exact, $d$ would be exactly the error and the corrected $x$ exact. It is not, but the error in $d$ is a small fraction of $d$, so each step multiplies the error by something small.

The classical form computes the residual in more precision than the solve, because the residual is a difference of nearly equal quantities — $b$ and $A x$ agree to most of their digits by construction. The figure uses the opposite, cheaper form: the residual in plain double, the same precision as everything else. The comparison is against the same network eliminated from the start in a pair of doubles carrying 31 digits, which shares nothing with the double solve but the matrix.

Plainly solved, the filter’s response is right to a unit of round-off in the passband and loses accuracy as the response falls: $1.1\times10^{-16}$ at 10 Hz, $9.3\times10^{-15}$ at 3 kHz, $1.1\times10^{-12}$ at 10 kHz, where the output is $6.78\times10^{-7}$ of the input. That is the earlier essay’s finding, on the earlier essay’s line. One step of refinement brings every frequency to within 2.4 units of round-off. No extra digits were used anywhere.

One step, and the residual’s precision does not matter #

The refinement converges in one step. At 10 kHz the error goes from $1.1\times10^{-12}$ to exactly zero against the 31-digit answer after the first step, then sits at a few units of round-off on later steps as the correction itself rounds; at 3 kHz, $9.3\times10^{-15}$ to $1.7\times10^{-16}$. Computing the residual in a pair of doubles changes the floor from a few units of round-off to zero and nothing else — it does not make the first step work any better, because the first step was not limited by the residual’s precision.

That is the opposite of what the textbook account of refinement would lead one to expect. Refinement in working precision is usually described as useless for accuracy, because a residual computed in the same precision as the solve is mostly round-off; what saves it here is that the problem was never the precision. The filter’s matrix is well conditioned enough — its condition number peaks at $1.59\times10^{7}$ — that a correction solved from even a rough residual is accurate relative to itself, and the correction is small. The plain solve’s trouble was a particular kind of error in the elimination, and a single honest residual exposes it.

Stable by the norm, not by the row #

The kind of error is the whole explanation, and it has a name.

A solve’s backward error is the smallest change to the data for which the computed answer would be exact. Measured the normal way — the size of the residual against the size of the whole matrix times the whole solution — the plain solve’s is $1.4\times10^{-18}$ at every frequency: it is exactly the answer to a problem changed by round-off. Gaussian elimination with pivoting is normwise backward stable, and the figure confirms it. But “changed by round-off” in the norm allows a change of round-off relative to the largest entry in every row, and in the stopband rows the entries that matter are small. Measured row by row, against each row’s own entries and solution components, the plain solve is exact only for data changed by $5.3\times10^{-13}$ at 10 kHz: the elimination has, in effect, perturbed the small stopband equations by an amount that is round-off for the input node and four decades more than round-off for the output.

A solve that is componentwise backward stable gives every component to round-off when the problem is well conditioned, small components included. One step of refinement in working precision makes Gaussian elimination componentwise backward stable — a result due to Skeel — and the dashed curve shows it happening: $9.8\times10^{-17}$ at every frequency after one step. The response’s digits were lost in the elimination’s rows, not in the matrix’s condition and not in the arithmetic’s width; the residual is computed against the original equations, row by row, and so it sees exactly the rows the elimination treated carelessly.

Why the condition number allowed it #

The condition number was no use as a predictor of the stopband’s loss, and the earlier essay showed that at length. It is exactly the right quantity for the question refinement asks, which is whether a correction computed with the same flawed elimination is good enough to improve on what it corrects. Each step multiplies the error by roughly the condition number times the round-off, $\kappa\varepsilon$, and refinement converges whenever that is less than one. For this filter $\kappa$ peaks at $1.59\times10^7$ at the band edge, so $\kappa\varepsilon$ is $3.5\times10^{-9}$: each step shrinks the error by eight decades, and one step is more than enough at every frequency.

That gives the condition number its proper job back. Where the matrix is worst, and where the answer is not found it overstating the solve’s error by many decades; here it bounds the rate at which the error can be removed, and the bound is loose in the safe direction. A network whose $\kappa$ approached $10^{16}$ — the feedback amplifier of the matrix that is ill, and the answer that is not, near the gain at which the solver refuses — would be one where refinement in double could not converge, because a correction would be as wrong as the thing it corrected.

The digits the arithmetic did not have reported the paired-double solve as twenty times slower than a plain one. One refinement step costs a matrix–vector product, of the order of the matrix’s number of entries, and two triangular solves with factors already in hand; for the filters here it adds a few per cent to the solve, and it recovers the same digits in the stopband that the twenty-fold slower route did. The expensive route is needed only where the cheap one cannot reach — which, by the next section, is where the digits were never in the solve.

Deeper stopbands, more to return #

The deeper the stopband, the smaller the response and the more the plain solve loses: $3.0\times10^{-16}$ at third order, $1.1\times10^{-12}$ at fifth, $1.4\times10^{-10}$ at ninth, where the response at ten times the cutoff is $4.9\times10^{-12}$ of the input. The even orders lose less than their odd neighbours, because their stopbands are shallower at this frequency, and every order stays under the round-off divided by its response. After one step every order is at round-off. A ninth-order filter’s stopband, which the plain solve reports to four significant figures, is reported to sixteen after a single matrix–vector product and two triangular solves — a cost of a few per cent on top of the solve.

Every derivative, and the one that is zero computed sensitivities with a transposed solve on the same factorisation; refinement is the same economy put to a different use. The factorisation is the expensive part of any solve, and once it is done, both a derivative and a correction are nearly free.

The bridge, which cannot be refined #

The bridge of the matrix that is ill, and the answer that is not loses one digit per decade of imbalance and has none left at a part in $10^{15}$. The earlier essay found that a pair of doubles recovers every digit. Refinement does not.

Plainly solved and refined twice in double, the bridge’s output error is the same at every imbalance — $5.4\times10^{-13}$ at $\delta = 10^{-4}$, $8.9\times10^{-5}$ at $10^{-12}$ — and on the line of round-off divided by $\delta$. Refinement found nothing to correct, because there was nothing wrong with the solve: both node voltages were already right to their last bit. The output is their difference, two numbers near 0.5 that agree to all but $\delta$ of their value, and two doubles near 0.5 cannot carry a difference smaller than their own round-off. The loss is in storing the answer as two node voltages, which happens after the solve and is invisible to its residual.

Carrying the node voltages as pairs of doubles through the refinement removes that loss, and the dots show what is left: an error still on the same line. Measured against the network the matrix actually holds, the paired answer is exact. The network the matrix holds is not the one that was drawn: its conductances are $1/R_3$ rounded to a double and its diagonal entry is $G_3 + G_4$ rounded to a double, and each rounding moves the bridge’s output by about round-off over $\delta$ — the same size as the first loss. The bridge loses its digits twice, both times in storage: once storing its data as conductances, once storing its answer as node voltages. The earlier essay’s pair of doubles recovered them because it built the matrix in pairs too.

Three places a digit can go #

Taken with the earlier essays, the networks now separate cleanly by where their digits were lost.

In the elimination. The filter’s stopband. The matrix and the data are fine, the arithmetic is wide enough, and the elimination perturbed some rows more than round-off relative to their own size. One step of refinement in working precision recovers everything, for a few per cent of the cost of the solve.

In the storage of the answer. The bridge’s output, a small difference of two stored node voltages. No refinement of the voltages can help; only carrying them in more digits, or choosing unknowns whose values are the small quantity itself — the kind of change the node that is not in the circuit measured when it moved the reference node, which here would make the bridge’s output a node voltage that refinement could then refine.

In the storage of the data. The bridge’s conductances, which carry an error of round-off over δ before anything is solved. Nothing downstream recovers this; the answer the arithmetic computes is exactly right for a network that differs from the drawn one. And beyond all three lies the floor that outlives the arithmetic: the filter synthesis, whose loss is in the data in the deeper sense that the problem itself is ill posed, and where more digits made things worse.

The habit that separates them is an old one, and the one that one step, computed twice applied to a transient: compute the same quantity two ways, and attribute the difference. A refinement step is a second route that costs almost nothing and splits the first category from the other two.

Keeping what the arithmetic threw away #

There is a family resemblance worth naming, because the same move appears in digital filtering. The rounding a section keeps found that a fixed-point filter section which remembers the error each rounding makes, and adds it back into the next computation, loses none of the accuracy a plain section loses to its dead band: the error was known exactly at the moment it was made, and discarding it was the only thing wrong. Refinement is the same move in a linear solve. The residual is the error the elimination made, measured against the original equations; the plain solve throws it away; one step keeps it and feeds it back.

In both cases the repair costs almost nothing because it uses information the computation already had, and in both cases it recovers everything the plain computation lost to its own procedure and nothing it lost elsewhere. The filter section’s error feedback cannot correct the quantisation of its input, which arrived already rounded; refinement cannot correct the bridge’s stored conductances, which arrived already rounded. The two repairs draw the same line between a loss made by the algorithm and a loss carried in by the data.

What the refinement does not settle #

The figures use dense Gaussian elimination with partial pivoting, the solver used for every network solved here. A sparse solver that orders its eliminations to reduce fill-in may perturb the stopband rows differently, and the size of the plain solve’s componentwise error is a property of the elimination order as much as of the network; the refinement’s effect — componentwise stability in one step — should not depend on it, but the size of what it recovers would. And the networks are small: a few dozen unknowns. A refinement step on a large circuit simulation costs one more back-substitution per frequency, which is cheap, but the matrix–vector product for the residual must be formed from the original matrix, which a solver that factorises in place no longer has unless it kept a copy.

The numbers #

The fifth-order ladder filter’s response at ten times its cutoff is $1.1\times10^{-12}$ wrong solved once in double, and within 2.4 units of round-off after one refinement step in the same precision; holding the residual in a pair of doubles changes only the floor. The plain solve is normwise backward stable to $1.4\times10^{-18}$ and componentwise only to $5.3\times10^{-13}$ in the stopband rows; one step makes it $9.8\times10^{-17}$. At ninth order the plain error is $1.4\times10^{-10}$ and the refined one round-off. The bridge refines to nothing: its node voltages were already exact, and its digits are lost storing their difference and storing its conductances, each by about round-off over $\delta$.

Still open: a reference that refines, a sweep that reuses, and a residual that warns #

The bridge with its reference at a bridge node. Moving the reference to one side of the bridge makes the output a node voltage. Whether refinement in double then recovers the output to round-off — leaving only the data’s storage error — would turn the second category into the first by a change of unknowns, and measure what that costs in the matrix’s condition.

Refinement across a sweep. A frequency sweep solves a matrix that changes slowly from point to point. A factorisation at one frequency, refined with residuals at the next, might converge in a few steps without a new factorisation, and the number of steps it needs against the frequency step is the cost curve of a cheaper sweep.

The residual as an alarm. A componentwise backward error is computed from the residual and costs no extra solve. Whether its size, row by row, predicts which components of a plain solve have lost digits — and so whether a simulator could refine only when the alarm sounds — is a test this essay’s figures already contain the data for.