A backward error the answer does not feel

An equation whose unknown is a matrix solved $AX + XB = C$ two ways, by Bartels and Stewart’s reduction of $A$ and $B$ to Schur form and by the $n^2 \times n^2$ Kronecker system that writes the equation as one linear system, and found the two agreeing. The elimination the matrix does not need found that the Kronecker form, after the same reduction, has nothing left to eliminate. Both essays measured the routes against each other and against sep, the equation’s conditioning. Neither asked what the residual both routes leave at the rounding level actually guarantees.

For a linear system it guarantees a great deal. A small relative residual means the computed answer is the exact answer of a nearby system, which is backward stability, and the error is then at most the condition number times the residual. For a matrix equation “nearby” means something narrower: nearby $A$, $B$ and $C$, three matrices perturbed separately, and an $n^2 \times n^2$ perturbation of the Kronecker matrix is mostly not of that form. Higham worked out what the gap can be. The smallest perturbation of $A$, $B$ and $C$ that makes a computed $\hat X$ exact can be larger than the relative residual by a factor

$$\mu = \frac{(\alpha + \beta)\lVert \hat X\rVert_2 + \gamma}{\sqrt{(\alpha^2 + \beta^2)\,\sigma_{\min}(\hat X)^2 + \gamma^2}} ,$$

with $\alpha, \beta, \gamma$ the norms of $A, B, C$, which is large exactly when $\hat X$ is ill conditioned and $C$ is small beside $(\alpha + \beta)\lVert \hat X\rVert$. This essay builds equations where both happen and measures three things: the residual, the backward error, and the error in the answer.

A family that cancels on purpose #

The equations are of order 8. $A$ is symmetric with eigenvalues 1 to 8 in a random orthogonal basis $Q$; $B = -A + \delta I$; and the exact solution is $X = Q\,\mathrm{diag}(x)\,Q^{\mathsf T}$, which commutes with $A$, with $x_i$ spaced geometrically so that $\kappa(X)$ is $10^t$ for $t$ from 0 to 8. Then $AX + XB = \delta X$ exactly, so $C = \delta X$ is as small beside $AX$ as $\delta$ makes it: half of it, or $10^{-2}$, $10^{-4}$, $10^{-6}$. The equation’s sep is $\delta$ itself, since $\lambda_i(A) + \mu_i(B) = \delta$ for every $i$, so $\delta$ also sets how sensitive the answer is. Each equation is solved by both routes, on three random bases, and the medians are reported.

The structured backward error is computed exactly. The perturbations $\Delta A, \Delta B, \Delta C$ satisfy one linear constraint, $\Delta A\,\hat X + \hat X\,\Delta B - \Delta C = C - A\hat X - \hat X B$, and the smallest of them, measured as $\lVert[\Delta A/\alpha,\ \Delta B/\beta,\ \Delta C/\gamma]\rVert_F$, is the minimum-norm solution of that constraint in vectorised form — a 64-row system at this order, which is small enough to solve outright.

The residual says stable #

Both routes leave relative residuals at the rounding level on all twenty equations: Bartels and Stewart’s between $1.7 \cdot 10^{-16}$ and $5.2 \cdot 10^{-16}$, the Kronecker system’s between $4.1 \cdot 10^{-17}$ and $6.3 \cdot 10^{-17}$, about four times smaller, as a backward-stable solve of one linear system should be. By that measure both routes are stable everywhere and the Kronecker one is the better.

The backward error disagrees. With $X$ well conditioned it is four times the residual on both routes, a constant that comes from the norms in the definitions. As $\kappa(X)$ grows with $C = 10^{-6}X$, Bartels and Stewart’s backward error goes from $1.5 \cdot 10^{-15}$ to $4.8 \cdot 10^{-14}$, $2.3 \cdot 10^{-12}$ and $3.9 \cdot 10^{-11}$ at $\kappa(X) = 10^4$, $10^6$ and $10^8$, while its residual stays at $2 \cdot 10^{-16}$. The Kronecker system’s goes to $1.0 \cdot 10^{-10}$. At the extreme its backward error is 2.6 times Bartels and Stewart’s, where its residual was four times smaller: the two measures rank the two routes in opposite orders. On the dial, larger $C$ flattens all of it; at $C = 10^{-2}X$ the backward errors stop at a few times $10^{-14}$.

Bounded, and not attained #

The figure at the top of the page shows the ratio of backward error to residual for Bartels and Stewart’s route, against $\kappa(X)$ for the four sizes of $C$, with $2\sqrt 2\,\mu$ dashed above each. The factor $2\sqrt 2$ is what the Frobenius norms used here contribute to Higham’s bound, and with it the bound holds on every one of the 120 solves: no ratio exceeds $2\sqrt 2\,\mu$, and on the well-conditioned equations, where $\mu$ is about 1.4, the ratio sits exactly at the bound, at 4.0.

Away from there the bound is loose by a varying amount. For Bartels and Stewart’s route the ratio reaches 185, 9,300 and 212,000 at $\kappa(X) = 10^8$ for $C = 10^{-2}X$, $10^{-4}X$ and $10^{-6}X$, while $\mu$ is 2,850, 285,000 and 28 million. The Kronecker route comes closer to its bound, within a factor of fourteen at the extreme. So $\mu$ is the right shape — it grows with $\kappa(X)$ and with $1/\delta$, and the measured ratios grow with them — but on a given equation it says how bad the backward error could be, not how bad it is.

Where the backward perturbation has to go #

Why the backward error grows is visible in where the smallest perturbation sits. With $X$ well conditioned it is split evenly between $A$ and $B$, 50 per cent each, and none on $C$: perturbing $A$ or $B$ by a little changes $AX + XB$ by a little in every direction, which is the cheapest way to absorb a residual. As $\kappa(X)$ grows, a perturbation of $A$ reaches the residual only through $\Delta A\,\hat X$, and along the directions where $\hat X$ is small that product is small whatever $\Delta A$ is. At $\kappa(X) = 10^8$ the residual has components along those directions that $\Delta A$ and $\Delta B$ cannot reach at any reasonable size, so the perturbation moves to $C$: 83 per cent of it, against 4 per cent on $A$ and 12 on $B$. And $C$ is a millionth of the size of $AX$, so a perturbation of it measured relative to $\lVert C \rVert$ is a million times larger than the same absolute perturbation measured against $A$. That is the backward error’s growth, in two factors: the directions only $C$ can reach, and the smallness of $C$ that makes reaching them expensive.

Why the smaller residual has the larger backward error #

The same picture explains why the two routes rank in opposite orders. What the backward error charges for is not the size of the residual but the part of it that lies where $\hat X$ is small on both sides — in the block of the residual between $\hat X$’s small left and right singular directions, which $\Delta A\,\hat X$ and $\hat X\,\Delta B$ cannot reach and only $\Delta C$ can.

Measured on one basis at $\kappa(X) = 10^8$ and $C = 10^{-6}X$, with five of $\hat X$’s singular values below a thousandth of the largest, Bartels and Stewart’s residual has norm $6.9 \cdot 10^{-15}$ and only $8.7 \cdot 10^{-17}$ of it in that block — a share of its squared size under a thousandth. The Kronecker system’s residual is four and a half times smaller, $1.5 \cdot 10^{-15}$, and $4.4 \cdot 10^{-16}$ of it lies in the block, 8.6 per cent of its squared size and five times Bartels and Stewart’s in absolute terms. At $C = 10^{-2}X$ the Kronecker share is 20 per cent and Bartels and Stewart’s still a thousandth.

That basis is the kindest to Bartels and Stewart’s route of the three, and the other two do not change the ranking. Over all three bases and all four sizes of $C$, Bartels and Stewart’s residual puts between 0.02 and 2.2 per cent of its squared size in the expensive block, and the Kronecker system’s between 6.9 and 41 per cent — at least fourteen times as large a share on every one of the twelve equations. The Kronecker residual is the smaller of the two in total on all twelve, by factors of 2.3 to 9.7, and the larger in the block on all twelve, by factors of 1.5 to 5. Neither share depends on the size of $C$ in any orderly way; what depends on the size of $C$ is how much each unit of residual in that block costs, which is the second factor of the section above.

The reason is in how each route rounds. Bartels and Stewart’s reduces $A$ and $B$ to Schur form, which for this symmetric family is their eigenbasis, which is also $X$’s: its rounding errors are committed in the coordinates where $X$ is diagonal, and they land mostly along $X$’s large directions, where the residual is cheap to explain. Gaussian elimination on the Kronecker matrix rounds in coordinates that know nothing of $X$, and spreads its residual over every direction, the expensive block included. A smaller residual spread evenly costs more than a larger one kept away from the expensive block, and two condition numbers of one matrix is the same lesson for a single linear system: what a perturbation costs depends on its shape as well as its size.

The answer does not feel it #

The forward error is the last measurement and it settles what the backward error means here. Against the exact $X$, Bartels and Stewart’s solution has relative error $6.4 \cdot 10^{-15}$ to $1.7 \cdot 10^{-14}$ with $C = 0.5X$, $1.6 \cdot 10^{-13}$ to $1.3 \cdot 10^{-12}$ with $C = 10^{-2}X$, $4.7 \cdot 10^{-11}$ to $1.0 \cdot 10^{-10}$ with $C = 10^{-4}X$, and $1.2 \cdot 10^{-9}$ to $9.0 \cdot 10^{-9}$ with $C = 10^{-6}X$. Within each size of $C$ the error does not rise with $\kappa(X)$ at all; it falls a little. Across sizes of $C$ it rises by about a hundred for every factor of a hundred in $\delta$, which is the unit roundoff over sep, the conditioning an equation whose unknown is a matrix found governing the equation. The Kronecker route’s forward errors are smaller again, between $1.3 \cdot 10^{-15}$ and $4.6 \cdot 10^{-10}$, and they too ignore $\kappa(X)$.

So the structured backward error rose by a factor of twenty-six thousand across the $\kappa(X)$ sweep at $C = 10^{-6}X$ and the forward error did not move. The two are consistent, because 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: the perturbations the backward error is made of lie 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. The answer is insensitive to exactly the perturbation the backward error is forced to use.

What the measurement says about each measure #

The residual is a statement about the Kronecker system, and both routes are backward stable for it: residuals of $10^{-16}$ to $10^{-17}$ on every equation. That is a true statement about a problem nobody posed, as a backward-stable answer to a problem nobody asked put it for linearised eigenproblems: the perturbations that make it true are of the $n^2 \times n^2$ matrix, not of $A$, $B$ and $C$.

The structured backward error is a statement about the equation, and it can be large for a solver that is, by every other measure, behaving well. A perturbation that moves every coefficient found the same thing one step further: restricting the perturbations to the coefficients anybody believes can make a backward error large for reasons that are about the restriction and not the computation.

The forward error is the one a user wants, and here neither backward measure predicts it on its own. The residual says the routes are equally good and the Kronecker one slightly better; the backward error says they degrade by five and six orders of magnitude; the forward error says nothing degrades with $\kappa(X)$, and the Kronecker route is better by a factor of 2.7 to 36. A code that reported the structured backward error as an accuracy warning would warn on exactly the equations where it need not. A small residual is not a small error is this collection’s oldest warning, and here its converse holds as well: a large backward error is not a large error either, when the perturbation it is made of is one the answer cannot see.

What a code can usefully report is the forward error’s own estimate: the unit roundoff times $\lVert A\rVert_2 + \lVert B\rVert_2$ over sep, which here predicts the forward error within a factor of ten either way on every equation and both routes — measured, between 0.11 and 7.4 times it. An estimate of sep costs a few solves with the triangular Schur factors Bartels and Stewart’s route has already computed, which is what library codes for this equation do. The residual and the structured backward error are both correct numbers, and neither is the one a user wants.

Why sep is the gap here, and what that leaves out #

In this family sep and the smallest eigenvalue sum are the same number, $\delta$, because $A$ and $B$ are symmetric: for normal coefficients the Kronecker matrix is normal too and its smallest singular value is its smallest eigenvalue in modulus. An equation whose unknown is a matrix built non-normal pairs to pull the two apart, and on those the forward error follows sep and not the gap. Here there is nothing to pull apart, which is deliberate: it isolates the backward error’s growth from the conditioning’s, so that everything the backward error does across the $\kappa(X)$ sweep happens at a fixed sep, and a fixed sep is why the forward error holds still.

When the structured backward error is the right measure #

None of this makes the structured backward error a wrong number. It is the right measure when the question is whether the computed $\hat X$ is the exact answer to data within its own uncertainty — when $A$, $B$ and $C$ came from measurements with relative error bars of their own. Then a backward error of $3.9 \cdot 10^{-11}$ relative to $C$ says the computed answer is exact for a $C$ perturbed by four parts in a hundred billion, and if $C$ was measured to six digits that is far inside its error bar and the computation added nothing. The measure answers “could this answer have come from my data”, and the forward error answers “how far is this answer from the one my data defines”; on this family the first can be two hundred thousand times the residual while the second stays where sep puts it. They are different questions, and a report that gives one in place of the other answers the wrong one with a correct number.

What this family does not show #

One family, built so that $X$ commutes with $A$ and $C$ is an exact multiple of $X$; that makes everything exact and also special, since a perturbation of $C$ along $X$’s own small directions is the whole of what the backward error has to do. With $X$ not commuting with $A$, the residual’s components along $X$’s small directions could be reachable through $A$ and $B$ in part, and the backward error would grow less. One order, 8, small enough to form the backward error’s normal equations outright; at larger orders the backward error can only be estimated. Symmetric $A$ with well-separated real eigenvalues, so that the Schur form is diagonal up to rounding and Bartels and Stewart’s route is at its best; non-normal $A$ and $B$, the case an equation whose unknown is a matrix used to separate sep from the eigenvalue gap, would add the Schur reduction’s own error to everything here.

Still open: a family that does not commute, and the condition number that matches #

Not commuting. The prediction with a sign is that with $X$ replaced by a random matrix of the same singular values, so that $AX + XB$ no longer cancels exactly and $C$ is formed from it, the backward error over the residual at $\kappa(X) = 10^8$ and $\delta = 10^{-6}$ falls below a thousand on both routes, and the share of the smallest perturbation on $C$ below a half.

The structured condition number. If the forward error is the backward error times a structured condition number, that number must fall by the factor the backward error rises. 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.