Residuals, Forward Error and Backward Error

Definition Let $A \in \mathbb{F}^{m \times n}$ and $b \in \mathbb{F}^m$, and let $\hat{x} \in \mathbb{F}^n$ be a proposed solution of $Ax = b$. The residual of $\hat{x}$ is

It measures the mismatch in the equation, in the output space, and it can be computed from the data without knowing any exact solution. In exact arithmetic $r = 0$ exactly when $\hat{x}$ solves the system. For $b \ne 0$ the relative residual is $\lVert r \rVert / \lVert b \rVert$. Now let $A$ be square and invertible and $b \ne 0$, so that the exact solution $x = A^{-1}b$ is nonzero. The forward error of $\hat{x}$ is $\hat{x} - x$, measured absolutely by $\lVert \hat{x} - x \rVert$ or relatively by

The forward error lives in the input space, and it usually cannot be computed, because $x$ is unknown. A residual and an error are different kinds of quantity, with different units in general, so they should not be compared as bare numbers without a sensitivity relation. A residual computed in floating-point arithmetic also carries its own rounding error, so a tiny computed residual has to be read together with that error.

Definition The backward error of a proposed solution $\hat{x}$ is the size of the smallest change of the data, within a stated class of allowed changes, that makes $\hat{x}$ an exact solution. The class and the way changes are measured must be stated, because different models give different numbers.

In the right-hand-side model, $A$ stays fixed and we require $A\hat{x} = b + \delta b$. The only such change is $\delta b = -r$, so the backward error is $\lVert r \rVert$, and relative to $b \ne 0$ it equals the relative residual $\lVert r \rVert / \lVert b \rVert$.

In the matrix model with the Euclidean norm, $b$ stays fixed, $\hat{x} \ne 0$, and we require $(A + \delta A)\hat{x} = b$, that is $\delta A\,\hat{x} = r$. Every such $\delta A$ has $\lVert \delta A \rVert_2 \ge \lVert r \rVert_2 / \lVert \hat{x} \rVert_2$, because $\lVert r \rVert_2 \le \lVert \delta A \rVert_2 \lVert \hat{x} \rVert_2$. The rank-one matrix

satisfies $\delta A\,\hat{x} = r$ and has $\lVert \delta A \rVert_2 = \lVert r \rVert_2 / \lVert \hat{x} \rVert_2$, so it attains the minimum. Relative to $\lVert A \rVert_2$ the backward error is $\lVert r \rVert_2 / (\lVert A \rVert_2 \lVert \hat{x} \rVert_2)$, which has a different denominator from the first model. If $\hat{x} = 0$ and $r \ne 0$, no change of $A$ alone can help, since every matrix sends $0$ to $0$. Allowing changes of both $A$ and $b$, or bounding each entry separately, defines still other backward errors.

Definition Let $A \in \mathbb{F}^{n \times n}$ be invertible, let $b \ne 0$ and $x = A^{-1}b$, and let $\hat{x}$ have residual $r = b - A\hat{x}$. In any vector norm on $\mathbb{F}^n$ with its induced matrix norm and condition number $\kappa(A) = \lVert A \rVert \, \lVert A^{-1} \rVert$,

Subtracting $Ax = b$ from $A\hat{x} = b - r$ gives $A(\hat{x} - x) = -r$, so $\hat{x} - x = -A^{-1}r$ and $\lVert \hat{x} - x \rVert \le \lVert A^{-1} \rVert \, \lVert r \rVert$. Dividing by $\lVert x \rVert$ and using $\lVert b \rVert \le \lVert A \rVert \, \lVert x \rVert$ gives the upper bound. For the lower bound, $\lVert r \rVert = \lVert A(\hat{x} - x) \rVert \le \lVert A \rVert \, \lVert \hat{x} - x \rVert$ and $\lVert x \rVert = \lVert A^{-1}b \rVert \le \lVert A^{-1} \rVert \, \lVert b \rVert$. Multiplying $\lVert \hat{x} - x \rVert \ge \lVert r \rVert / \lVert A \rVert$ by $1/\lVert x \rVert \ge 1/(\lVert A^{-1} \rVert \, \lVert b \rVert)$ gives it. The relative residual is the relative right-hand-side backward error, so in this model the relative forward error is at most the condition number times the relative backward error. In the Euclidean norm both bounds are attained for suitable $b$ and $\hat{x}$, but for a particular $\hat{x}$ the error can lie anywhere between them. The proof uses $A^{-1}$, so it does not apply to a singular system.

Definition Let $A \in \mathbb{F}^{n \times n}$ be invertible, with a declared norm. The system $Ax = b$, with $A$ fixed and $b$ subject to error, is ill-conditioned when $\kappa(A)$ is large, so that, for some right-hand sides, small relative changes in $b$ can cause relative changes in $x$ up to $\kappa(A)$ times larger. There is no universal threshold. If the data carry relative errors up to $\eta$, solutions consistent with the data can differ from $x$ by relative amounts up to $\kappa(A)\,\eta$, and for some right-hand sides and perturbations that worst case is reached. When $\kappa(A)\,\eta$ is near $1$ or larger, the bound no longer guarantees even one correct digit. Ill-conditioning belongs to the problem, its norm and its error model. No algorithm removes it, although a poor algorithm can add errors of its own.

The size of the determinant is not a test of conditioning. Writing $A = U\Sigma V^{*}$ with unitary $U$ and $V$, whose determinants have absolute value $1$, gives

In absolute value the determinant is the product of the singular values, while the condition number is a ratio of two of them. Scaling $A$ by $c \ne 0$ multiplies the determinant by $c^n$ and leaves the condition number unchanged.

Definition The floating-point rounding model describes arithmetic with a unit roundoff $u$. A real number $y$ in the normal range is stored as $\operatorname{fl}(y) = y(1 + \delta)$ with $\lvert \delta \rvert \le u$, and each basic operation $\circ$, one of $+$, $-$, $\times$ and $\div$, applied to stored numbers $s$ and $t$ returns

Overflow and underflow are excluded from this simple model. Rounding to nearest with $d$ significant decimal digits gives $u = \tfrac{1}{2} \cdot 10^{1 - d}$, and IEEE double precision has $u = 2^{-53} \approx 1.1 \cdot 10^{-16}$.

An algorithm for square systems $Ax = b$ is backward stable, normwise and under this model, if for every input it returns a vector $\hat{x}$ that exactly solves a nearby system

with $\lVert \delta A \rVert \le c_n u \lVert A \rVert$ and $\lVert \delta b \rVert \le c_n u \lVert b \rVert$, where $c_n$ is a modest constant that depends only on the dimension. If $A$ is invertible, $b \ne 0$ and $x = A^{-1}b$ is the exact solution, combining this with conditioning gives

which follows from $A(\hat{x} - x) = \delta b - \delta A\,\hat{x}$, the inequality $\lVert \hat{x} \rVert \le \lVert x \rVert + \lVert \hat{x} - x \rVert$ and $\lVert b \rVert \le \lVert A \rVert \, \lVert x \rVert$. So backward stability guarantees a small forward error only when $\kappa(A)\,c_n u$ is small. It controls the error that the algorithm adds, while the condition number decides what that error does to the answer. Proving that a whole algorithm is backward stable requires bounding every rounded operation it performs. A small residual in one example proves nothing of the kind.

Definition Let $A \in \mathbb{F}^{m \times n}$ have singular values $\sigma_1 \ge \cdots \ge \sigma_p \ge 0$, where $p = \min(m, n)$. Its exact rank is the number of nonzero singular values, counted with no tolerance. For a tolerance $\tau > 0$, its numerical rank is the number of singular values greater than $\tau$. By the Eckart-Young theorem in the spectral norm, this is the smallest rank of a matrix within spectral distance $\tau$ of $A$,

If exactly $k$ singular values exceed $\tau$, the truncated SVD $A_k$ has rank $k$ and $\lVert A - A_k \rVert_2 = \sigma_{k+1} \le \tau$, with $\sigma_{p+1}$ read as $0$, while every $B$ of rank less than $k$ has $\lVert A - B \rVert_2 \ge \sigma_k > \tau$. The tolerance has to come from the problem, for example from the size of the measurement errors in $A$ or from the rounding errors of the arithmetic, not from a wish for tidy output. It is an absolute size in the units of $A$, so rescaling $A$ without rescaling $\tau$ can change the count. Replacing $A$ by a matrix of lower rank also changes the problem. It can make a solution less sensitive to noise, but the reduced model no longer represents every direction of $A$.