A step that is not a unit of work

The iterative field ended its first account of regularisation on a sentence about integers. A parameter that counts steps put conjugate gradients on a 64-point deconvolution at one per cent noise, found the error least at step 20, and named the step count as the knob: an integer nobody had called a parameter, doing the job that λ does for Tikhonov. Everything measured on the subject since has been measured on that one iteration, including the step at which its textbook description expires.

That leaves a question nobody here has asked. If the step count is the parameter, is a step the same thing in every iteration that regularises by stopping? The cheapest way to find out is to run a second iteration on the identical problem, with the identical noise, and read off where it is best.

The second iteration here is Landweber’s, and its answer is step 1,778.

The oldest regulariser there is, written as one line #

Landweber’s iteration is gradient descent on the least-squares misfit with a fixed step length:

x₍ₖ₊₁₎ = xₖ + ω Aᵀ(b − A xₖ)

Each step multiplies by A once and by Aᵀ once, which is exactly what a step of CGLS costs, so on this problem a count of steps is a count of work for both methods and the comparison needs no conversion factor. The iteration is also a stationary method in the sense of a rate that is known in advance: the same linear map applied to the same vector over and over, with nothing adapted from one step to the next. Richardson’s iteration on the normal equations is the same object under another name.

What makes it worth running beside conjugate gradients is that its filter factors are available in closed form at every step. Starting from zero, after k steps

fₖ(σ) = 1 − (1 − ωσ²)ᵏ

so the iterate is the filtered sum Σ fₖ(σᵢ)(uᵢᵀb/σᵢ)vᵢ with a weight that is a function of σ alone. Conjugate gradients also applies a filter, but its filter is a polynomial whose roots sit at Ritz values the iteration discovers as it runs; Landweber’s has no roots to discover. It is one at the top of the spectrum, near zero at the bottom, and never leaves the interval between.

The step length is not free. Once ωσ₁² exceeds two, the factor (1 − ωσ₁²) has modulus above one, the top component grows geometrically, and the iteration diverges instead of regularising slowly. That is one of the claims the figure refuses, with a step length of 2.2 over σ₁², because a curve drawn from the closed form there would be drawing a formula the iteration does not obey. Everything below uses ω = 1/σ₁² unless it says otherwise.

Eighty-nine times the work for the same answer #

The two best errors in the first figure are not merely close. At every noise level the slider reaches, they agree to better than one per cent, and the steps they take to get there do not:

The error columns could be one column. The step columns share nothing but their direction. At ten per cent noise the slow method is barely slower — seven steps against three, fourteen products against six — and at a tenth of a per cent it takes a hundred and twelve thousand products to reach what conjugate gradients reaches in eighty-eight. The ratio rises at every stop of the table without exception, which is checked as a non-decreasing sequence rather than read off two ends.

The mechanism is in where each filter puts its cutoff. Reading off the index at which each iteration’s filter factor first falls below a half, on this operator at one per cent noise:

and for conjugate gradients, index 11 at step 2, 16 at step 8, 20 at step 16 and 23 at step 20. Conjugate gradients adds a component to what it has recovered every step or two. Landweber adds three to five components for every factor of ten in its step count. The reason is the closed form: for ωσ² small the factor is close to 1 − exp(−ωkσ²), which passes a half where σ² ≈ 0.69/(ωk), so multiplying the step count by ten moves the cutoff down by a factor of ten in σ². This operator’s singular values fall by a factor of ten in σ² every three or four indices near the cutoff, so a decade of steps buys three or four of them.

That is a statement about this spectrum as much as about the method. The blur’s singular values decay exponentially in the index, which is what makes the problem ill-posed rather than ill-conditioned in when the answer is a choice, and exponential decay is exactly what turns a knob that scales σ² into a knob that moves the index logarithmically. A spectrum decaying like a power of the index would give a cutoff moving like a power of the step count instead — slower than conjugate gradients still, but not by the three orders of magnitude seen here. That case is not measured here.

A slow knob is a forgiving knob #

The step count of conjugate gradients has to be found. The essay that introduced it measured a window of 7 to 23 steps inside five per cent of the best error, and found stopping early four times more forgiving than stopping late. A window read in steps is a window read in the knob’s own units, and Landweber’s units are not the same size.

The span of steps within ten per cent of each method’s best, first to last:

From two per cent noise downwards Landweber’s window is between fifteen and a hundred and forty times wider, as a ratio, than conjugate gradients’. At a tenth of a per cent, stopping anywhere in a factor of sixty-seven gives up at most a tenth of the best error; conjugate gradients’ equivalent latitude is a factor of two.

The one-per-cent row is the flattering one and the reason is worth marking, because it is the row the first figure draws. Landweber’s error falls quickly to about 0.155 within a handful of steps and then creeps down a shallow shoulder for three decades, and 0.155 happens to be inside ten per cent of 0.1414. Half a per cent of noise lowers the best error to 0.1225, the shoulder is no longer inside the window, and the window’s first step jumps from 7 to 631. A window statistic depends on where a plateau sits against a threshold somebody chose, which is the same kind of number the zero a reader is allowed to write is careful about. The robust claim is the ratio of ratios, and it is the one the check holds: at least ten at every noise level where Landweber’s best step is past a hundred.

Lopsidedness survives the change of units, measured the way the knob moves. At one per cent, Landweber’s window runs 2.4 decades below its optimum and 0.55 above it; conjugate gradients’ runs sixteen steps below and six above. Both iterations forgive an early stop more than a late one. What differs is that a late stop by a factor of three takes conjugate gradients from step 20 to step 60, far outside its window, and takes Landweber from 1,778 to 5,334, still inside ten per cent.

What a stopping rule does with a slow knob #

None of the optima above is available without the answer, so the practical question is where a rule lands. The discrepancy principle stops at the first step whose residual has fallen to the stated noise level, and it has already been measured on conjugate gradients as a rule that errs early and is right to.

The rule’s stops, with the error at each, across the whole range:

Landweber’s stops are read at the steps its curve is sampled at, forty a decade, so a stop in the thousands is resolved to about six per cent; conjugate gradients’ are exact.

The error a reader receives from the rule is within three per cent between the two methods at six of the seven levels, and within five at the seventh, where ten per cent noise stops Landweber a step past its own best and conjugate gradients on it. That is a property of the rule and the noise more than of which iteration the rule is attached to: both stop where the residual first reaches the stated noise, and both residual curves reach it while the error is still falling.

What the rule costs in work is another matter, and it is not the ratio of the optima. At one per cent the rule’s stop is 2.4 times later on Landweber, not 89 times, because the residual reaches the noise level early on a curve that then spends three decades improving the error in its fourth significant figure. A reader who stops by the rule pays roughly double for Landweber; a reader who insists on the optimum pays ninety times.

The ratio of stops does not stay at two. It reads 1.3, 2.0, 1.8 and 2.4 down to one per cent noise and then 4.5, 33 and 104. Below one per cent Landweber’s residual approaches the noise level along the same logarithmic crawl as its cutoff, so the step at which the residual finally crosses is pushed out with everything else — and the rule starts paying the slow knob’s price in work while still collecting the same error as on the fast one. The two rows at the bottom of the table are where the choice of iteration begins to cost a reader who stops by the rule, and the crossover sits between half a per cent and two tenths of one.

And there is a third way to stop that only the slow knob makes available: a budget. Landweber stopped at a thousand steps, for no reason connected to the data, returns 0.1419 at one per cent noise — within 0.4% of its best. Conjugate gradients run for a hundred and fifty steps has already reached an error of ten. An iteration whose cutoff moves by a few indices a decade is regularised by almost any reasonable budget, and one that moves by an index a step is regularised only by a rule.

The parameter was never the step count #

Landweber’s closed form says what the parameter actually is. For ωσ² small the factor depends on ω and k only through their product, so halving the step length should double the best step and change nothing else.

It does, to one per cent across a factor of 3.8 in the step length, at the same best error to five figures. The step count of Landweber’s iteration is ωk wearing a unit, and the unit is chosen by whoever sets ω. A stopping rule written as “run for N steps” transfers between two Landweber codes only if they share a step length, and between Landweber and anything else not at all.

The longest step length is the one place the product is not the whole story, and the figure shows where. At ω = 1.9/σ₁² the top component’s factor is 1 − 1.9 = −0.9, so its share of the iterate alternates in sign and decays like 0.9ᵏ: the first iterate overshoots the largest singular direction by ninety per cent, the error at step 1 is 0.75 against 0.29 at the unit step length, and it takes about thirty steps for the overshoot to fall below the error the other two curves are carrying by then. After that the dashed curve is the solid one moved left by a factor of 1.9. So a longer step buys its speed back from the bottom of the spectrum and pays for it at the top, and the payment is finished long before the optimum; nothing in the argument depends on the choice between 1 and 1.9.

What does cost something is knowing σ₁ at all. Every step length above is a multiple of 1/σ₁², and σ₁ was read off the SVD the figures already had. A code without that SVD estimates σ₁ by a power iteration, which is more matrix–vector products spent before the first Landweber step, and an underestimate by more than a factor of √2 puts ωσ₁² above two and the iteration diverges. That preamble is not priced in the tables here.

The product has a second reading, and it connects the iteration to the penalty. Where ωkσ² is small the Landweber factor is close to ωkσ², and where σ² is large against λ² the Tikhonov factor σ²/(σ² + λ²) is close to σ²/λ², so the two filters agree in both limits when λ² = 1/(ωk). At the best ωk of 1,778 that λ is 0.0237. The Tikhonov oracle on the same problem, found by sweeping λ directly, is 0.0224, and the four-knob comparison’s coarser grid put it at 0.025. Landweber stopped at its best step is, to within six per cent in λ, Tikhonov at its best λ.

That also accounts for the one number in the first table that is not a tie. Landweber’s best error at one per cent is 0.1414, conjugate gradients’ 0.1426, and the four knobs put Tikhonov’s at 0.1406. The order is the order of how much each filter departs from a monotone weight: Tikhonov’s is monotone by construction, Landweber’s is monotone and slightly sharper, and conjugate gradients’ goes above one and comes back. An overshoot of a fifth at index 17 is a component recovered with its noise amplified, and it costs eight parts in a thousand.

A description that does not expire, and a check that did #

An expiry date the noise does not move measured the step at which conjugate gradients’ filter polynomial stops describing its iterate: step 17 or 18 on this operator, because the Lanczos basis the polynomial is derived from stops being orthogonal. Landweber has no basis. Its closed form is an identity about a matrix power, and a matrix power has nothing to lose orthogonality in.

The lower curve is the claim, and it holds to rounding over three thousand steps. The middle curve is the part worth recording, because it was the first measurement and it said something false.

The check was first written with the filter as the formula is printed: one minus ωσ², raised to the kth power, subtracted from one. It reported a disagreement of 4.6·10⁻¹¹ at step 10, 4.6·10⁻¹⁰ at 100 and 4.6·10⁻⁹ at 1,000 — a drift growing exactly in step with the step count, which is what rounding accumulating one step at a time inside the iteration would look like. It was not in the iteration. Forming 1 − ωσ² keeps the absolute accuracy of ωσ² and loses its relative accuracy wherever ωσ² is small; raising the result to the kth power multiplies that relative error by k; and the filtered sum then divides the lost part by the smallest singular values in it. Rewritten as −expm1(k·log1p(−ωσ²)), which never forms the difference, the same check agrees with the iterate to rounding at every step.

So the two descriptions on the figure fail in opposite ways. Conjugate gradients’ polynomial is a correct identity evaluated on a basis the arithmetic has corrupted, and it expires. Landweber’s is a correct identity evaluated by a formula the arithmetic corrupts, and it drifts — and the iterate was right throughout. Both routes are kept in the check: the stable form is required to agree to rounding, and the printed form is required to drift, so that a future rewrite which reintroduces the subtraction is caught by a check rather than by a reader.

Why an iteration this slow is run at all #

Nothing above makes Landweber the better method on this problem. At its best step it is ninety times the work of conjugate gradients for eight parts in a thousand of error, and at the discrepancy principle’s stop it is two to a hundred times the work for none. The case for it is made of properties the comparison did not test, and they are worth listing because each is a reason the slow knob turns up in practice.

It computes no inner products. A step is a product, a subtraction and a scaled addition; conjugate gradients also needs two inner products a step, and on a machine where every inner product is a collective communication, the message rather than the arithmetic is the cost that decides. How the product counts above translate into messages is not measured here.

Its filter is a fixed function of σ at every step. A constraint on the answer — nonnegativity, a box, a support — can be imposed by projecting each iterate, and the projected iteration is still a descent method on the same misfit. A projected conjugate gradient iteration loses the optimality its polynomial came from. The consequence for the error on a constrained problem is the obvious next measurement and has not been made.

And its window is wide. Where no stopping rule is available — no noise estimate, no residual threshold anybody trusts — a method whose error stays within ten per cent of its best over a factor of nine hundred in steps is regularised by the budget that happened to be spent. That is the property the tables above do establish, and it is the one a reader choosing between the two should weigh against the factor of ninety.

What the comparison has not established #

One operator, and one whose spectrum decays exponentially. The logarithmic knob is a consequence of that decay, argued above from the closed form and measured on this blur. On an operator with power-law decay Landweber’s cutoff would move polynomially in the step count, and the ratio of best steps would be smaller; that is a reasonable expectation and not a measurement.

Landweber’s curve past three thousand steps is the closed form’s. The iteration itself is run to 3,000 and checked there; the ten-million-step end of each figure is the stable formula evaluated at that step. The agreement to rounding at every checked step, and the mechanism of the only drift ever found, are the grounds for trusting the extrapolation, and neither is a proof that nothing new happens at step 10⁶.

Accelerated Landweber is not measured. Semi-iterative variants — the ν-methods and Chebyshev acceleration — reach a given ωk in roughly its square root of steps while keeping a filter that is a fixed function of σ, which sits between the two iterations here on both the cost and the forgiveness axis. Where exactly it sits has not been measured.

And “the same error” is at each method’s oracle. A practical rule lands a few per cent above it on both, as measured above, and the comparison of work at the rule’s stop — a factor of 2.4 at one per cent — is the one a reader actually pays. The factor of 89 is the price of insisting on the optimum.

Where this goes from here #

The step count turned out to be a unit chosen by the iteration. Two measurements follow directly, and both change something the comparison above held fixed.

The preconditioner. On a well-posed system a preconditioner changes how fast an iteration converges and not what it converges to. On this problem the iteration is stopped before it converges, so a preconditioner changes where every step lands — and one good enough, by the usual measure, puts the best iterate before the first step. A preconditioner that arrives past the answer measures that, and measures what a smoothing operator used as a preconditioner does to a signal with an offset.

The smoothness of the answer. Every error in the tables above is scored against the collection’s constructed signal, two bumps and a step, and on that signal Landweber, conjugate gradients and Tikhonov all land within a per cent. Landweber’s filter is an exponential in ωkσ² and Tikhonov’s is a rational function, and on an answer smooth enough to tell them apart that difference decides how fast each method’s error can fall as the data improves. The method that cannot use a smooth answer measures the rates.

And a rule written in the right unit. A discrepancy-principle stop is already free of the step length, because it reads the residual. A budget is not, and a budget written in ωk rather than in steps would transfer between step lengths and would be the natural stopping rule for the slow iteration — its behaviour across problems with different spectra is unmeasured, and it is the practical question the forgiving window raises.