A nearby problem of the wrong kind

The sentence this whole site is built on is that a good algorithm returns the exact answer to a nearby problem. Every essay here that reports an error reports which of the two factors it came from, and the backward error — the distance to the nearby problem — is the factor that belongs to the algorithm.

A hundred and eighteen essays have measured that distance. Not one of them has asked what the nearby problem is.

For a general matrix the question has no content: any matrix is a matrix, and a perturbation of A is a perturbation of A. For a matrix with structure it has a great deal, and this essay is about a case where the answer is uncomfortable.

What a Toeplitz system is #

A Toeplitz matrix is constant along its diagonals: T[i][j] depends only on i − j. An n×n one is therefore 2n − 1 numbers in an n² object, and it is what a stationary process, a convolution, a time-invariant filter or a discretised translation-invariant kernel produces.

A fast solver for it takes O(n²) operations instead of O(n³) — the same trade a matrix that is one row makes one structure further in — and the reason is exactly the count above: Levinson’s recursion works on the 2n − 1 numbers and never touches an entry off the first row. It builds the solution of the k×k leading system from the (k−1)×(k−1) one using a single reflection coefficient, and it never forms a factor of anything.

That last clause is where this essay starts. There is no factorisation to be the exact factorisation of a nearby matrix, so the standard backward error analysis — the one that says LU with partial pivoting is stable because the computed factors are the exact factors of A + ΔA — has nothing to be about.

The measurement needs a known answer #

The family here is the Kac–Murdock–Szegő matrix ρ^|i−j|, and it is used because this site already has its closed-form inverse — an answer that is known rather than computed. So the true solution of every system below is known in exact terms rather than computed, and every forward error is a measurement.

That matters more than usual, because the whole essay is about distinguishing between two kinds of smallness, and a comparison of two float computations could not.

Levinson on this family agrees with elimination and with the closed-form inverse: relative forward error below 10⁻⁹ at ρ = 0.5, 0.9 and 0.99, on a 12×12 system, and the two solvers land in the same place. Nothing is wrong with the algorithm, which is the point.

The smallest perturbation of any kind #

Rigal and Gaches settled this in 1967 and the answer is a formula rather than a search. The smallest ‖ΔA‖ with (A + ΔA)x̂ = b is exactly ‖r‖/‖x̂‖, and it is attained by

ΔA = r x̂ᵀ ⁄ ‖x̂‖²

which is rank one.

That is the object every backward-stability claim on this site has been about — including the growth factor’s — and it has never been drawn. It is the top matrix in the hero: dense, every entry nonzero, and constant along none of its diagonals — the measured departure from being Toeplitz is 0.885, on a scale where one means the diagonals carry as much variation as the matrix itself.

So the reassurance in its usual form reads: the computed answer is the exact solution of a nearby problem, which is not a problem of the kind that was posed, and which has n² degrees of freedom where the original had 2n − 1.

The smallest perturbation of the right kind #

The honest question is what the smallest perturbation that is Toeplitz costs, and it is a constrained least-squares problem rather than a formula.

Write ΔA as a combination of the basis matrices of the structure — one per constant diagonal — and the requirement ΔA·x̂ = r is linear in the coefficients. Weighting them by how many entries each diagonal has turns the whole thing into an ordinary minimum-norm problem, so the structured backward error is a least-squares answer rather than an estimate.

On the 40×40 system at ρ = 0.999:

• smallest of any kind: 2.0·10⁻¹⁷;
• smallest that is Toeplitz: 1.15·10⁻¹⁴;
• smallest that is symmetric Toeplitz: 5.6·10⁻¹².

Insisting on the structure costs two and a half million at the worst point of the family, and insisting on the symmetry as well costs another factor of five hundred.

Two things about that figure are worth separating.

Neither number is alarming. 5.6·10⁻¹² is a small perturbation of a matrix by any standard, and a small residual is not a small error is about the gap in the other direction — and the computed solution really is the exact solution of a genuine symmetric Toeplitz system a hair away from the one posed. The reassurance survives; it is merely worth orders less than the number that is printed.

And the two solvers are indistinguishable on the reported number. Levinson’s normwise backward error and Gaussian elimination’s are within a factor of two at every ρ, which is what “both are backward stable” means and is what a library would report. On the structured quantity they separate. The number that can tell them apart is the one nobody computes.

Why the constrained problem is a least-squares problem #

The construction is worth a paragraph, because it is the reason the number above is a minimum rather than an example.

A symmetric Toeplitz perturbation is a sum over its n constant diagonals: ΔA = Σ tₖ Eₖ, where Eₖ has ones on the k-th diagonal and its mirror and zeros elsewhere. The requirement is ΔA·x̂ = r, which is n equations linear in the n unknowns t. And the size of the perturbation is

‖ΔA‖²_F = Σₖ mₖ tₖ²

where mₖ is how many entries the k-th diagonal has, which is a weighted sum of squares. Substituting uₖ = √mₖ·tₖ turns it into an ordinary minimum-norm problem — minimise ‖u‖ subject to Gu = r — whose answer is the pseudoinverse applied to r.

Two things fall out of writing it that way and both are used above.

The minimum is attained, so the number reported is the smallest structured perturbation and not one somebody found. A larger one would prove nothing; a smaller one does not exist. That is the difference between a bound and a measurement three errors and one number turns on.

And feasibility is checkable, which is a refusal rather than a number. If Gu = r has no solution — if the residual is not in the range of the map from structured perturbations to their action on x̂ — then no perturbation of that kind explains the answer at all, and the routine reports so rather than returning the least-squares approximation as though it were a perturbation. That is what happens at the last rung of the ladder, and it is why the number there is absent rather than large.

And the gap between the two widens as the family gets harder, which the size sweep alone cannot show.

Only the structured number moves. Across ρ = 0.5, 0.7, 0.8, 0.9 and 0.95 the unconstrained backward error sits at 2.17·10⁻¹⁷, 1.67·10⁻¹⁷, 1.99·10⁻¹⁷, 1.72·10⁻¹⁷ and 1.89·10⁻¹⁷ — flat, and at the rounding level throughout — while the smallest Toeplitz perturbation runs 3.83·10⁻¹⁵, 6.83·10⁻¹⁵, 6.44·10⁻¹⁵, 1.01·10⁻¹³ and 2.65·10⁻¹³, growing by a factor of 69.

Each halving of 1 − ρ does not cost the same amount. Eight-tenths to nine-tenths multiplies the ratio by 18; nine-tenths to 0.95 multiplies it by 2.4. The blow-up is real and it is not a clean power of 1 − ρ, which matters because the bounds usually quoted for this are stated as one.

So the number a solver reports is the one that does not notice. The unstructured backward error is 10⁻¹⁷ on every one of those five matrices and would be quoted as backward stability on all of them; the question the caller actually asked — is there a nearby Toeplitz system this answer solves exactly — has an answer four orders of magnitude worse at ρ = 0.95 than at ρ = 0.5. The reassurance does not merely fail to improve as the problem hardens, it is constant by construction, because the space of unconstrained perturbations is large enough to absorb whatever the family does. A quantity that cannot get worse is not measuring the difficulty.

Which is the reassurance getting weaker with size #

That last observation is the one worth carrying out of this essay, because it runs against how a backward error is usually read.

A backward-stable algorithm’s error bound grows slowly with n — a modest polynomial, usually — and that is taken to mean the reassurance holds up as problems get larger. For the structured question it does the opposite. The constrained perturbation is choosing among 2n − 1 numbers to explain a residual that lives in n dimensions, while the unconstrained one chooses among n².

That is an argument, and it has an exponent. Fifteen right-hand sides at each size on the same family:

   n     unstructured    Toeplitz      symmetric     T ⁄ unstr   sym ⁄ unstr
   8      2.90·10⁻¹⁷    2.72·10⁻¹⁶    1.44·10⁻¹⁴        9.4        5.0·10²
  16      2.00·10⁻¹⁷    3.98·10⁻¹⁶    3.12·10⁻¹⁴       20          1.6·10³
  32      2.07·10⁻¹⁷    1.86·10⁻¹⁵    1.16·10⁻¹³       90          5.6·10³
  64      2.25·10⁻¹⁷    3.09·10⁻¹⁵    1.02·10⁻¹³      140          4.5·10³

The first column does not move. 2.90·10⁻¹⁷ to 2.25·10⁻¹⁷ across a factor of eight in n, which is what a mild polynomial looks like when it is mild — and it is the number a library reports.

The ratio to the structured quantity fits n^1.37 for Toeplitz and n^1.26 for symmetric Toeplitz. So over a range where the reported backward error is flat to within a factor of 1.3, the distance to the nearest problem of the same kind grows more than tenfold.

The exponent is worth one remark, because it is not the one the degree-of-freedom count suggests. The unconstrained perturbation has n² parameters and the constrained one 2n − 1, a ratio growing like n; the measured gap grows like n^1.4, which is faster. The extra comes from the constrained problem having to satisfy n equations with 2n − 1 parameters — mildly underdetermined, so its minimum-norm answer is close to being pinned — while the unconstrained one has n² parameters for the same n equations and is free to spread the perturbation as thinly as it likes. The two are competing on different margins, and only one of them has room to improve as n grows.

So the discount grows with the problem. A reader who checks a reported backward error at n = 8 and finds it reassuring has, at n = 64, the same reassuring number standing for a distance more than ten times larger. The reassurance does not merely say less than it appears to; it says progressively less as the problem gets bigger, which is the opposite of what a bound polynomial in n leads one to expect.

What to report instead #

The measurement makes the recommendation specific rather than a gesture, and it has three parts of increasing cost.

Report the residual and the solution norm — the quantity a method can see. ‖r‖/‖x̂‖ is the unstructured backward error exactly — Rigal and Gaches, no search — and it costs one matrix–vector product. It is what libraries print and it is worth printing, because it is the one quantity that bounds the others from below.

Report the structure the matrix was stored in. A caller who knows the matrix came in as 2n − 1 numbers knows the reported figure is a lower bound for a question about a different object, and the essay’s table says by roughly how much. That costs nothing at all and is a change to a docstring.

And compute the structured number when the verdict matters. It is one least-squares solve against an n × n matrix built from the structure basis, so it costs about a solve — a constant factor on top of a fast method that was chosen to be fast. That is a real cost and it is the right one to pay when the question is which of two solvers is better on this class, because the table above shows the two solvers separating there and nowhere else.

The third is the one this essay is arguing for and the second is the one that would help most people, which is the usual ordering: the expensive fix is the interesting one and the free fix is the one that is missing.

Large structured problems are exactly the ones where the fast algorithm is worth having, which is to say: the reassurance is weakest where the algorithm is most necessary.

And at the tightest description there is none #

There is one more rung and it is the one that changes the shape of the argument.

The Kac–Murdock–Szegő matrix holds one number. Not 2n − 1 and not n — the whole matrix is ρ^|i−j|, and ρ is what the model contains. Ask for the smallest perturbation of the model that makes the computed solution exact, and the answer is that there is none.

No value of ρ whatever has the computed x̂ as its exact solution. The closest one leaves a large share of the residual unexplained, and the structured backward error at that rung is not a large number, it is infinite.

That is not a defect of the solver and it is not a surprise once stated: a one-parameter family is a curve in the space of matrices, and a computed solution generically corresponds to no point on a curve. What is worth noticing is that the sentence — the exact answer to a nearby problem — has quietly stopped being true, at the description of the problem that the person who posed it would recognise.

What a caller could actually do #

The structured backward error costs a least-squares solve with 2n − 1 unknowns and n equations, which is O(n³) if done naively and is therefore more expensive than the O(n²) solve it is auditing. Nobody is going to compute it routinely.

Two cheaper things are available and neither is standard.

Compare the two solvers. Elimination and the fast algorithm on the same system give two answers, and the difference between them is a lower bound on how far either is from the truth. It costs a slow solve, so it is a check for the first run rather than for every run — which is exactly how a reference implementation is normally used, and this collection’s habit of two routes to a number applied to a case where the second route is merely the expensive one.

Or check the residual against the structure. The residual r and the computed x̂ are both to hand, and whether r lies in the range of the structured perturbation map is a rank question about an n×(2n−1) matrix — cheaper than the minimisation and enough to distinguish “there is a nearby problem of the right kind” from “there is not”. On the last rung of the ladder that check is what says the answer is absent.

The general point is the one this site keeps arriving at: the diagnostic quantity exists, is computable, and is not what gets printed. What is printed is the normwise number, because it has a closed form and costs nothing, and the closed form is available precisely because the question it answers is the unconstrained one.

What this does not say #

Three limits, because the finding is easy to over-read.

It is not that Levinson is unstable. It is backward stable in the normwise sense, its forward error tracks the closed-form answer, and on this family it is as accurate as elimination. Every claim about it in this collection stands.

It is not that structured backward errors are always large. They are bounded below by the unconstrained one and above by nothing in general, and on this family they stay at 10⁻¹² — small. There are algorithms and structures where the structured quantity is enormous or infinite, and there are pairs where it is barely larger; measuring is the only way to know which.

And it is not an argument for the slow method. Elimination’s structured backward error is also orders above its normwise one — the shape belongs to the question, not to the fast algorithm. Whatever is being given up here was given up by asking a normwise question about a structured problem, and no choice of solver recovers it.

The pattern across three fields #

This collection has now measured three quantities that a summary number conceals, and they are the same shape in three different places.

A condition estimate is a lower bound that libraries print as though it were the value, and the matrices where it under-reports are the ones where the number mattered.

A residual is small on a computation that solved a different problem, and a small residual is not a small error.

And now a backward error is a distance to a problem of a kind nobody posed, and the distance to the nearest problem of the right kind is orders larger.

In all three the printed number is the one with a closed form, the useful one requires an extra computation, and the extra computation is skipped for a reason that is entirely defensible in isolation. What the three have in common is worth naming as a rule of thumb rather than a discovery: a quantity that is free to compute is answering the question that has a closed form, and that is rarely the question being asked.

The refusal #

The assertion is fed a Levinson recursion on a Toeplitz matrix that is not positive definite.

The recursion divides by α, the residual variance of the order-k prediction problem, and on a positive definite matrix that quantity is positive and decreasing. On an indefinite one it goes negative, the next reflection coefficient is a ratio of two numbers whose signs mean nothing, and the recursion continues to the end and returns a vector.

The vector has the right length and the wrong everything else, and it comes back without a warning, because nothing in the loop tests a quantity that is positive on every matrix the method was designed for. It is the same shape as elimination without a pivot search: an algorithm that runs to completion on an input it was not written for and returns an answer of the correct type.

What is next #

This essay measured what the structure costs when it is imposed on the perturbation. The next one measures what it buys when it is imposed on the problem — the condition number of a matrix described by fewer numbers, which falls as the description shrinks, and where the fall happens turns out not to be where the linear algebra is.