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³), 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. 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, 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 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_k E_k, where E_k 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 = Σ_k m_k t_k²

where m_k is how many entries the k-th diagonal has, which is a weighted sum of squares. Substituting u_k = √m_k·t_k 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.

And feasibility is checkable. 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.

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². The gap between what they can achieve widens as n grows, so the distance from the posed problem to the nearest problem of the same kind grows faster than the reported distance does.

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.