The residual knew the error

A backward error the answer does not feel built Sylvester equations $AX + XB = C$ on purpose to separate two numbers that a linear system keeps together. With $A$ symmetric with eigenvalues one to eight, $B = -A + \delta I$ and $X$ commuting with $A$ — so that $C = \delta X$ exactly — it swept the conditioning of $X$ from 1 to $10^8$ and the size of $C$ from half of $X$ down to $10^{-6}X$. Both routes to the answer, Bartels and Stewart’s Schur reduction and the $n^2 \times n^2$ Kronecker system, left residuals at the rounding level everywhere. The structured backward error — the smallest perturbation of $A$, $B$ and $C$, each measured against its own size, that makes the computed $X$ exact — did not: it climbed by a factor of 26,000 across the sweep at the smallest $C$. And the forward error against the exact $X$ did not move at all.

The essay’s explanation was that “a forward error is bounded by a condition number times a backward error, and the condition number that goes with these structured perturbations falls as fast as the backward error rises”. It asked for that condition number. “The prediction is that the structured condition number for perturbations of $A$, $B$ and $C$ measured relative to their own norms, computed exactly at this order, falls by at least ten thousand across the $\kappa(X)$ sweep at $C = 10^{-6}X$, and that its product with the backward error is within a factor of ten of the measured forward error on every equation.”

It is computed here, and it does the opposite. And the number that does track the error turns out to be the one the earlier essay set aside in the first paragraph: the residual.

The condition number of an equation in three matrices #

Perturb $A$, $B$ and $C$ by $\Delta A$, $\Delta B$ and $\Delta C$ and the answer moves by $\Delta X$ with, to first order, $A\,\Delta X + \Delta X\,B = \Delta C - \Delta A\,X - X\,\Delta B$. Written with the Kronecker matrix $P = I \otimes A + B^{\mathsf T} \otimes I$ that an equation whose unknown is a matrix introduced, that is $\operatorname{vec}\Delta X = -P^{-1}H e$, where $e$ stacks the three perturbations divided by the norms $\alpha$, $\beta$ and $\gamma$ of $A$, $B$ and $C$, and $H = [\alpha (X^{\mathsf T} \otimes I),\ \beta (I \otimes X),\ -\gamma I]$. That $H$ is the same matrix whose minimum-norm solve defined the backward error, which is why the two belong together: the backward error is the smallest $\|e\|$ with $He$ equal to the residual, and the structured condition number is the largest $\|P^{-1}He\|/\|X\|$ per unit of $\|e\|$, the norm of $P^{-1}H$ over the size of $X$.

At order eight $P^{-1}H$ is a 64 × 192 matrix and its norm can be computed directly: 192 solves with $P$, then the largest singular value. It is computed here from the exact $X$ on every equation of the family, three bases each, and set beside the backward errors and forward errors the earlier essay measured, for both routes.

That is affordable only because the order is small. The elimination the matrix does not need priced the Kronecker matrix of order $n^2$ and found the structure that lets Bartels and Stewart avoid ever forming it; the structured condition number has no such shortcut, since it is the norm of $P^{-1}$ applied to three blocks, one of them $n^2$ wide, and at order eight that is 192 solves with a 64 × 64 matrix. At order a hundred it would be thirty thousand solves with a matrix of order ten thousand. A number that expensive has to earn its place, and the question here is whether it does.

It rises #

At $C = 10^{-6}X$ the structured condition number is $7.1 \cdot 10^{6}$ when $X$ is perfectly conditioned and $2.0 \cdot 10^{7}$ when $\kappa(X) = 10^8$. It does not fall; it rises, by a factor of 2.8. The same factor appears at every size of $C$: 13.7 to 38.6 at $C = 0.5X$, 714 to 2,013 at $10^{-2}$, $7.1 \cdot 10^{4}$ to $2.0 \cdot 10^{5}$ at $10^{-4}$. Almost all of the change happens between $\kappa(X) = 1$ and $10^2$, and from there it is flat. Over the same sweep the backward error goes from $1.5 \cdot 10^{-15}$ to $3.9 \cdot 10^{-11}$, so the two curves in the figure run in opposite directions where the prediction needed them to run in parallel and opposite.

The size of the number is no mystery once it is computed. Each level of $C$ sets it at about $(\alpha + \beta)/\operatorname{sep}(A, B)$: with $\alpha$ and $\beta$ near 14.3 each and $\operatorname{sep}$ equal to $\delta$ for this family — the smallest $|\lambda_i(A) + \lambda_j(B)|$ is $\delta$, at $i = j$ — that is $2.9 \cdot 10^{7}$ at $\delta = 10^{-6}$, and the structured number sits a factor of 1.4 below it at the ill-conditioned end. Perturbing $A$ by a relative $\eta$ changes $AX$ by about $\eta\alpha\|X\|$, and the inverse of the Kronecker matrix multiplies that by up to $1/\operatorname{sep}$. That direction is always available. The structured condition number is a supremum over all perturbations of a given size, and the supremum is reached by perturbing $A$ or $B$ along the direction $P^{-1}$ amplifies most, which has nothing to do with how $X$ is conditioned.

So the product is a bound, and a loose one #

The product of the backward error and the structured condition number is a valid first-order bound on every solve: never below the error, its smallest ratio 1.13. But it is a bound whose slack grows with exactly the quantity the family turns up. At $C = 0.5X$ it is within a factor of six for Bartels and Stewart and of 57 for the Kronecker system. At $C = 10^{-6}X$ its median over the three bases is 1.2 times the error at $\kappa(X) = 1$, 10.5 at $10^2$, 512 at $10^4$, $3.2 \cdot 10^4$ at $10^6$ and $6.3 \cdot 10^5$ at $10^8$ for Bartels and Stewart; for the Kronecker system $4.6 \cdot 10^6$ at the far end. The worst single solves overstate by $3.2 \cdot 10^6$ and $9.3 \cdot 10^6$. The prediction asked for a factor of ten.

The two numbers in the product describe different directions. The backward error is the size of the most economical perturbation, and the earlier essay found where that perturbation lands: almost entirely on $C$, along the directions where $X$ is small, because $\Delta A\,X$ and $X\,\Delta B$ both pass through $X$ and cannot reach there cheaply. The condition number is the response to the least economical perturbation of the same size, a change in $A$ along the direction $P^{-1}$ amplifies most. Multiplying them supposes that the cheapest perturbation is also the most damaging one. Here it is close to the least damaging, and the product is wrong by the ratio between the two directions’ effects — six decades at the corner of the family.

The residual’s estimate #

The other estimate needs nothing structured at all. The forward error is $P^{-1}$ applied to the residual: $\operatorname{vec}(\hat X - X) = P^{-1}\operatorname{vec}(R)$ exactly, with $R = C - A\hat X - \hat X B$, by linearity, whatever perturbation of the data anybody imagines produced $\hat X$. So $\|\hat X - X\|$ is at most $\|R\|/\operatorname{sep}(A, B)$, and dividing by $\|X\|$ and writing the residual relative to its own scale $(\alpha + \beta)\|X\| + \gamma$ gives the relative residual times $((\alpha + \beta)\|X\| + \gamma)/(\operatorname{sep}\,\|X\|)$ — the unstructured condition number, which for this family is 55.6, 2,860, $2.86 \cdot 10^5$ and $2.86 \cdot 10^7$ at the four sizes of $C$ and does not depend on $X$’s conditioning at all.

That estimate is within a factor of 0.96 to 9.65 of the forward error on all 120 solves. At $C = 10^{-6}X$ with Bartels and Stewart it is 1.2, 2.9, 3.4, 3.9 and 4.2 times the error across the $\kappa(X)$ sweep, and 3.0 to 4.4 for the Kronecker system. Drag the dial up through the sizes of $C$ and the three curves of the figure behave the same way at every level: the forward error and the residual’s estimate run together, flat across $\kappa(X)$, and the backward error’s estimate leaves them as $X$ grows ill conditioned, by more the smaller $C$ is. At $C = 0.5X$ all three of Bartels and Stewart’s lie within a factor of six.

Plotted directly, error against estimate, the 120 solves make one band seven decades long. The errors run from $7.5 \cdot 10^{-16}$ at $C = 0.5X$ to $1.1 \cdot 10^{-8}$ at $C = 10^{-6}X$, each size of $C$ a cluster two decades further along, and every cluster sits between equality and ten times it. The conditioning of $X$, which moved the backward error’s estimate by six decades, moves points only within their cluster. That is what an estimate is for: a reader handed only the residual and $\operatorname{sep}(A, B)$, knowing nothing about how $X$ was built or which route computed it, could place the error within a decade on every one of these equations.

Bracketing an error nobody can measure put two computable numbers either side of an error that could not be computed directly, and asked how far apart they had to be. Here the error can be computed — the family has an exact $X$ — and the residual’s estimate sits above it by a factor between one and ten on every solve, which is as tight a single computable number as a reader could ask for. The structured product sits above it too, always, and is useless for the same purpose, because a bound that can be a million times too large says nothing about the error it bounds.

It is not quite a bound, since the residual’s own size is rounded and the inverse’s action is not always at its worst, and the lowest of the 120 sits at 0.96. But it is an estimate in the sense the prediction wanted: it says what the error is, to within a factor of ten, on every equation of a family built to make the backward error say something else.

The smaller residual, the worse estimate #

The two routes rank differently, and the earlier essay found the reason. The Kronecker system’s residual is about four times smaller than Bartels and Stewart’s, and its backward error is larger wherever $\kappa(X) \ge 10^6$ and $C$ is small, because more of its residual lies where only $C$ can reach. The backward error’s estimate therefore does worse on the route with the smaller residual — $4.6 \cdot 10^6$ times the error against $6.3 \cdot 10^5$ at the far corner — while the residual’s estimate ranks them correctly: the Kronecker system’s forward error is smaller, and so is its residual. A reader choosing between the two routes by their backward errors would choose the less accurate one at every $\kappa(X)$ from $10^6$ up; choosing by their residuals they would choose the more accurate one everywhere.

What the earlier explanation got right #

The earlier essay’s account of why the forward error did not move was half right. It said that the backward error’s perturbation “lies on $C$ along the directions where $X$ is small, and moving $C$ in those directions moves the answer only in those directions, where it is small to begin with”. That is a statement about a direction, and it is true. What does not follow is that a condition number captures it, because a normwise condition number is a supremum over directions and cannot know which direction a particular backward perturbation took. A condition number that did know would be a function of the residual — it would have to be told where the perturbation lies — and at that point it is simpler to apply $P^{-1}$ to the residual and skip the backward error altogether.

A condition number sent to infinity met the same shape in an interior-point method: a matrix driven deliberately toward singularity, a normwise condition number heading off without limit, and an answer that kept its digits because the errors the factorisation made lay in directions the answer was not sensitive to. The normwise number there was a supremum over directions the computation never visited. The structured condition number here is the same kind of object, and the only cure in both places is to look at the direction the error actually took — which, for a linear equation, is what applying the inverse to the residual does.

A small residual is not a small error is the site’s warning against reading a residual as an answer’s accuracy, and it stands: the residual here is $10^{-16}$ relative and the error $10^{-9}$, seven decades apart, because the residual has to be multiplied by a condition number of $2.9 \cdot 10^{7}$. What this family adds is that the multiplication is honest. The residual, times the one condition number that does not depend on anything but $A$ and $B$, is the error to within a factor of ten. It is the backward error, the quantity that is supposed to be the more meaningful one, that needs a second number to translate it, and that second number is wrong by up to six decades.

When the backward error is still the right question #

None of this makes the structured backward error a bad measure. It answers a question about the data, not the answer: whether the computed $X$ is the exact solution of an equation whose coefficients are within their own uncertainty of the given ones. If $A$, $B$ and $C$ come from measurement, that is the question that decides whether the computation added anything to what the data already allowed, and the exact answer to a nearby problem is the long argument for asking it. On this family it says that Bartels and Stewart’s $X$ needs a perturbation of $3.9 \cdot 10^{-11}$ — far more than rounding — to be exact, and that is true and worth knowing about the algorithm.

Cost points the same way. The residual’s estimate needs the residual, which costs two matrix products, and $\operatorname{sep}(A, B)$, the smallest singular value of the Kronecker matrix, which is itself expensive to compute exactly and is in practice estimated from the Schur forms both routes already have. The structured product needs the backward error — a minimum-norm solve with an $n^2 \times 3n^2$ matrix — and the structured condition number on top. The cheaper number is the better one.

What it does not do is predict the error, and the measurement shows the gap is not a constant to be calibrated but a function of the problem: a factor of one on a well-conditioned $X$ and of millions on a badly conditioned one with a small right-hand side. A condition number scaling cannot move found Skeel’s componentwise number smaller than the normwise one by an amount that depends on the matrix in the same way; there the smaller number was the one that told the truth, and the same is true here of the residual’s estimate against the structured product.

What eight-by-eight equations do not show #

One family, built so that $X$ commutes with $A$ and $C$ is exactly $\delta X$. In a family where $X$ does not commute — the earlier essay’s other open prediction, still untested — the backward error’s perturbation need not lie so cleanly on $C$, and the product could be much closer to the error; the residual’s estimate does not depend on that. Order eight only, which makes $P^{-1}H$ small enough to form; at realistic orders the structured condition number is estimated rather than computed, and an estimate of a supremum can only be looser. And the forward errors are measured against an $X$ that is exact by construction, which a real application does not have.

Still open: a family that does not commute, and a residual for the Schur route #

Not commuting. The earlier essay’s prediction for a non-commuting $X$ still stands: with $X$ replaced by a random matrix of the same singular values, the backward error over the residual at $\kappa(X) = 10^8$ and $C = 10^{-6}X$ falls below a thousand on both routes. This essay adds a prediction with a sign about the estimates: there the structured product overstates the forward error by less than a hundred on every equation, while the residual’s estimate stays within ten — so that the gap measured here is a property of the commuting family’s geometry, not of structured estimates in general.

A residual for the Schur route that costs nothing. Bartels and Stewart solve a triangular equation $SY + YT = F$ in the Schur bases, and its residual there is available without forming $A\hat X$. The prediction is that that residual, rotated back, is within a factor of two of the true residual on every equation here, so that the residual’s estimate of the error can be had for the price of the solve itself.