When the answer is a choice

Every solve on this site so far has had one answer, and the only question was how accurately it was reached. Even the valley with no bottom, which found that the coefficients of a degree-9 fit can be moved 39% of their own size along the worst singular direction before the residual changes in the sixth digit, was measuring how flat a minimum is — there was still a minimum.

This field is about the problems where that is false. The data is consistent with a range of solutions differing by orders of magnitude, no arithmetic can choose between them, and something outside the data has to.

The choice is the computation. It is not a preliminary to it, it is not a matter of taste, and it cannot be avoided by using better arithmetic — which is what makes this field belong on a site whose premise is that the algorithm’s contribution and the problem’s are separately measurable. Here the algorithm’s contribution is zero and the answer is still wrong by eight orders of magnitude.

The problem #

A one-dimensional deconvolution. A signal is blurred by a Gaussian kernel and the blurred version is measured with 0.1% noise; the task is to recover the signal.

The signal is two smooth bumps and a step, and the step is deliberate: a smooth signal is recovered well by every method here and would hide the differences between them. The blur matrix is symmetric Toeplitz — the same object the structure field is built on, arriving here for a completely different reason — with each row normalised to sum to one, so the operator preserves a constant and the difficulty is entirely about the smoothing.

The exact signal is known because the problem was constructed, so the forward error is known rather than estimated. That is the same move exact.js makes with the Hilbert matrix and it is worth as much here: without it, every parameter-choice rule in the third essay could only be compared with the others.

Ill-posed is not the same as ill-conditioned #

The distinction is the field’s foundation and it is routinely collapsed.

An ill-conditioned problem has a large condition number. It amplifies perturbations, the amplification is computable, and the site has been measuring it since its foundation phase: a 13×13 Hilbert system with κ = 1.7·10¹⁸ returns 0.049 where 8 belongs, and the algorithm is blameless.

An ill-posed problem is worse in a specific way. Its singular values decay to zero smoothly, with no gap anywhere, so there is no rank to find and no threshold to set. Three measurements say so here:

The decay is exponential rather than a power. Fitting log σ against k gives R² = 0.930; fitting it against log k gives 0.610. A matrix whose σ fall like 1/k is ill conditioned and has a numerical rank; one whose σ fall like exp(−k²) has neither.

There is no gap. The largest ratio between consecutive singular values over the usable range is under 4. Rank is a decision establishes that a numerical rank is a decision about a gap and that the gap is worth printing beside it. Here the point is that there is not one — every threshold is as defensible as every other, which is the same statement as the answer is a choice.

And the unregularised answer is worse than returning zero. A Householder least-squares solve — the best-conditioned route this site has, the one the road that squares the problem recommends over the normal equations — returns a relative error of 5.5·10⁸. Not 5.5, not 550. Five hundred and fifty million times the size of the signal it was asked to recover.

κ is 5.7·10¹². The algorithm’s backward error is at the level of rounding. The identity forward error ⪅ condition number × backward error accounts for the whole of it, and the account contains no error on anybody’s part.

One expression, two filters #

Both standard methods are the same sum with a different weight:

x_reg  =  Σₖ fₖ · (uₖᵀb / σₖ) · vₖ

Truncated SVD: fₖ = 1 for the first K terms and 0 after. Tikhonov: fₖ = σₖ²/(σₖ² + λ²), which is nearly one when σ ≫ λ and nearly σ²/λ² when σ ≪ λ. No regularisation at all: fₖ = 1 for every k.

That last line is the diagnosis. The unregularised solution is the same sum with every weight set to one, so its last terms are uₖᵀb divided by a σ of 10⁻¹⁶ — and uₖᵀb down there is not signal, it is the noise. Dividing noise by 10⁻¹⁶ is the whole of the 5.5·10⁸.

Writing the two methods as one expression is not tidying-up. It makes the comparison between them a comparison of weights rather than of methods, and it is enforced by there being exactly one implementation of the sum in lib/regular.js: both routines produce an f and hand it to the same function. Two methods that differ by one line cannot differ in any way the essay has not accounted for.

It is also what makes the figure possible. A plot of two solution vectors would show two curves that both roughly resemble the signal; a plot of two filters shows what the methods actually do, which is decide how much of each singular direction to believe.

The two routes #

The filter-factor sum goes through the SVD. A solver takes a different route entirely: form the stacked least-squares problem

[ A  ]        [ b ]
[ λI ] x  ≈   [ 0 ]

and solve it with a QR. Same answer, no SVD anywhere, and the arithmetic shares nothing with the first route.

They agree to a relative difference of 4.4·10⁻¹⁵ at λ = 0.1, rising to 5.2·10⁻¹² at λ = 10⁻⁴ — which is the two routes to a number habit doing what it is for, and is what makes the filter picture safe to build the field’s whole exposition on.

The stacked route has something the sum does not, and it is worth more than a second opinion. Its condition number has a closed form: the singular values of [A; λI] are √(σₖ² + λ²), so

κ  =  √((σ₁² + λ²)/(σₙ² + λ²))

which is asserted as an equality against the measured κ and agrees to fifteen digits. Regularising does not only change the answer that comes back; it changes the problem handed to the solver. At λ = 0.1 the stacked system’s κ is 10.07 against the original’s 2.9·10¹².

The first version of that assertion used the bound everyone quotes, κ ≤ σ₁/λ, and failed: 10.0659 measured against 10.0161 claimed. The bound is too tight by exactly the +λ² in the numerator. A closed form has no such approximation in it, and this is the reason to prefer one over the memorable version.

What regularisation actually does to the signal #

The honest presentation, and the one a filter plot cannot give.

Even at the best λ available, the step never comes back sharp. That is not a failure of the parameter choice and no better rule fixes it: the blur has destroyed the components that represent an edge, and no filter over the remaining ones can put them back. What less noise buys is a better rounded corner, not a corner.

So regularisation does not remove noise. It removes the ability to represent whatever the filter cut off, and the signal pays too. The essay’s title is meant literally: choosing λ is choosing which features of the answer to give up, and there is no value of it that gives up none.

The two failure directions are not symmetric either, and the figure asserts both. Under-regularised, the answer’s norm is larger than the truth’s — the amplified noise has to go somewhere. Over- regularised, it is smaller. That asymmetry is what the L-curve of the third essay is a picture of.

Why a better solver is not the answer #

The refusal at the end of lib/regular.js is fed the claim that a backward-stable method solves an ill-posed problem accurately, and requires it to fail.

That refusal is doing more work than most on this site, because the claim is the natural reading of everything the earlier fields establish. The road that squares the problem says: do not use the normal equations, use QR. Two Gram–Schmidts says: use Householder, not the classical process. Each of those is a real improvement and each is worth taking.

None of them touches this. The QR route and the normal-equations route on this problem return answers that are both useless, in slightly different ways, and the difference between them is invisible next to the 5.5·10⁸ they share. There is a whole class of problem on which every piece of advice the earlier fields give is correct and irrelevant, and recognising which class a problem is in is the first thing this field is for.

The knob is not the knob #

One more refusal, and it is about the parameter rather than the method. λ = 0 is not the gentle end of the regularisation knob — it is off the end of it, and returns the 5.5·10⁸ catastrophe. The generator refuses it rather than drawing it as a mild case.

That is the same shape as tau = 0 in the sparse pivoting field, where the threshold rule accepts a pivot created by fill and cancellation, ‖PA − LU‖/‖A‖ reaches 2.5·10¹¹ and the solve returns infinities. Both knobs look continuous and both have a point at one end where the method stops being the method. Every real code has a floor for the same reason.

The condition number is not the whole story either #

Worth being precise, because κ = 5.7·10¹² is the number a reader will reach for and it does not by itself say ill-posed.

The 13×13 Hilbert matrix has κ = 1.7·10¹⁸ — six orders of magnitude worse — and its system is perfectly well posed: there is one right answer, this site knows it exactly through exact.js, and the reason the float solve fails to find it is entirely accounted for by the amplification. Give the same system to a solver with enough precision and the answer comes back.

Give this one more precision and it does not. The obstruction is the noise in b, and no amount of arithmetic removes it: the coefficients of the noise in the trailing singular directions are what they are, and they get divided by σ ≈ 10⁻¹⁶ in exact arithmetic just as surely as in floating point. A hundred-digit solve of this problem returns a hundred-digit-accurate answer to the wrong question.

That is the sharpest available statement of the difference. Ill-conditioning is a property of the matrix and precision is its remedy. Ill-posedness is a property of the matrix and the data together, and there is no remedy — only a choice about what to give up.

It is also why this field’s essays report the forward error throughout and the earlier fields report the backward one. A backward error here is uninformative by construction: every method in this field has one at the level of rounding, including the one that returns 5.5·10⁸.

What each filter is good at #

Neither is uniformly better, which is why both are still in use, and the figure’s drag is what shows it rather than an argument.

Truncation is a decision and Tikhonov is a weighting. Truncation keeps K terms whole and discards the rest, so its answer lies exactly in a K-dimensional subspace and can be reported as such — which matters when the components are interpretable. Tikhonov keeps a fraction of everything, so its answer has a component in every direction, most of them tiny.

Truncation has a discrete parameter. There are only n choices and they can all be tried, which is what makes the Picard analysis of the next essay possible at all. Tikhonov’s λ is continuous, and choosing it is the third essay.

And Tikhonov never has to form the SVD. The stacked route above is a QR of a 2n×n matrix, which is affordable at sizes where an SVD is not, and is what any large-scale implementation actually does. The filter picture is the right way to understand both methods and the wrong way to compute either at scale.

The measured errors on this problem, at the matched parameters the figure draws, are within a few per cent of each other. The choice between them is about what the answer is for rather than about accuracy.

What the blur is standing in for #

The problem drawn is a one-dimensional Gaussian deconvolution, and the choice of it is a simplification worth naming.

Every real instance of this is at least two-dimensional. An image restored from a blurred photograph, a computed tomogram, a seismic inversion, an unfolding of a detector response — all of them are integral equations of the first kind discretised, and all of them have singular values that decay to zero with no gap, which is the only property this field uses. What the two-dimensional versions add is size: a 512×512 image is a quarter of a million unknowns, and forming an SVD of the operator is impossible.

That does not change any statement in this field and it changes every implementation of them. The filter picture is still the right way to understand what Tikhonov does; it is not the way anybody computes it at that size. The stacked-QR route scales further, and past that the methods are iterative — conjugate gradients on the normal equations, stopped early, which turns out to be itself a regulariser with its own filter factors, and the iteration count becomes the parameter.

That last observation is the one this field most obviously does not reach, and it is a satisfying one: the iterative field’s whole subject is how fast an error falls, and on an ill-posed problem the error falls and then rises, because the later iterations are reconstructing the noise. Stopping is the choice, and it is the same choice as choosing λ wearing different clothes.

Two thresholds, side by side #

Worth putting next to each other, because this site now has both and they are routinely confused.

Rank is a decision sets a threshold on σ to decide a rank: which directions the matrix genuinely spans. The matrix there has a gap — ten nonzero singular values and a cliff — and the decision is about where in the cliff to cut.

This field sets a threshold on σ to decide a filter: how much of each direction to believe. There is no cliff, so the threshold is not finding a structure in the matrix. It is trading two errors against each other, and the trade depends on the noise in the data rather than on the matrix at all.

The same arithmetic — compare σₖ against a number, keep or discard — answering two different questions, one about an operator and one about a measurement. Which is why the next essay is about the right-hand side rather than about the matrix.

What is left #

Where the answer stops being in the data, which is the next essay: the Picard condition, which is computable from the data alone and says whether there is a solvable problem here at all — and which turns out to identify the wrong index by a consistent margin.

Choosing λ without knowing the answer, which is the third and is where the field’s measurements against an oracle live.

And the general-form problem. Everything here penalises ‖x‖, which is a statement that the answer should be small. The general form penalises ‖Lx‖ for some operator L — usually a derivative, which is a statement that the answer should be smooth — and it is what anybody solving a real deconvolution uses. The reduction to the standard form is a generalised SVD, and none of it is here; the choice of L is a stronger and less examinable assumption than the choice of λ, which makes it a subject rather than an omission.