14 Conditioning and Numerical Stability
In numerical linear algebra, errors in computed solutions arise from two distinct sources: the inherent sensitivity of the mathematical problem, and the accumulation of roundoff errors within the executing algorithm. The condition number quantifies problem sensitivity, whereas numerical stability describes the fidelity of the algorithm. A reliable computational result requires both a well-conditioned problem and a backward stable algorithm.
14.1 Residual, Forward Error, and Backward Error
Let \(\bA\bx = \bb\) be a nonsingular linear system with exact solution \(\bx^* = \bA^{-1}\bb\), and let \(\hat{\bx}\) denote a computed approximation. We assess the quality of \(\hat{\bx}\) through three interrelated quantities.
For a computed approximation \(\hat{\bx}\):
The forward error is the vector difference \(\Delta\bx = \hat{\bx} - \bx^*\), with relative measure \(\frac{\|\hat{\bx} - \bx^*\|}{\|\bx^*\|}\).
The residual is the vector defect \(\br = \bb - \bA\hat{\bx}\), with relative measure \(\frac{\|\br\|}{\|\bA\|\|\hat{\bx}\| + \|\bb\|}\).
The backward error is the smallest perturbation \((\delta\bA, \delta\bb)\) such that \(\hat{\bx}\) is the exact solution to the perturbed system \((\bA + \delta\bA)\hat{\bx} = \bb + \delta\bb\).
Let \(\bA \in \fR^{n \times n}\) be nonsingular, let \(\bb \neq \bzero\), and let \(\bx^*\) satisfy \(\bA\bx^* = \bb\). For any vector \(\hat{\bx}\) with residual \(\br = \bb - \bA\hat{\bx}\), the relative forward error satisfies \[ \begin{align} \frac{\|\hat{\bx} - \bx^*\|}{\|\bx^*\|} \leq \kappa(\bA) \frac{\|\br\|}{\|\bb\|}, \end{align} \] where \(\kappa(\bA) = \|\bA\|\|\bA^{-1}\|\) is the condition number of \(\bA\).
Because \(\bb = \bA\bx^*\), the residual can be written in terms of the forward error as \(\br = \bb - \bA\hat{\bx} = \bA(\bx^* - \hat{\bx})\). Multiplying both sides by \(\bA^{-1}\) yields \(\bx^* - \hat{\bx} = \bA^{-1}\br\). Taking vector norms and applying the submultiplicative property of induced matrix norms gives \[ \begin{align} \|\hat{\bx} - \bx^*\| \leq \|\bA^{-1}\| \|\br\|. \end{align} \] Furthermore, from \(\bb = \bA\bx^*\), we obtain \(\|\bb\| \leq \|\bA\| \|\bx^*\|\), which rearranges to \(\frac{1}{\|\bx^*\|} \leq \frac{\|\bA\|}{\|\bb\|}\). Multiplying these two inequalities produces \[ \begin{align} \frac{\|\hat{\bx} - \bx^*\|}{\|\bx^*\|} \leq (\|\bA\| \|\bA^{-1}\|) \frac{\|\br\|}{\|\bb\|} = \kappa(\bA) \frac{\|\br\|}{\|\bb\|}. \end{align} \]
(A small residual does not guarantee a small error) Consider the diagonal system \[ \begin{align} \bA = \begin{pmatrix} 1 & 0 \\ 0 & 10^{-8} \end{pmatrix}, \qquad \bb = \begin{pmatrix} 1 \\ 10^{-8} \end{pmatrix}. \end{align} \] The exact solution is \(\bx^* = (1, 1)^T\). Suppose an algorithm produces the approximation \(\hat{\bx} = (1, 0)^T\). The residual evaluates to \[ \begin{align} \br = \bb - \bA\hat{\bx} = \begin{pmatrix} 1 \\ 10^{-8} \end{pmatrix} - \begin{pmatrix} 1 \\ 0 \end{pmatrix} = \begin{pmatrix} 0 \\ 10^{-8} \end{pmatrix}, \end{align} \] which satisfies \(\|\br\|_2 / \|\bb\|_2 \approx 10^{-8}\). However, the forward error is \(\hat{\bx} - \bx^* = (0, -1)^T\), yielding a relative forward error of \(\|\hat{\bx} - \bx^*\|_2 / \|\bx^*\|_2 = 1/\sqrt{2} \approx 0.707\) (a seventy percent error). The matrix has condition number \(\kappa_2(\bA) = 10^8\), which amplifies the small residual into an \(O(1)\) error in the solution.
14.2 Conditioning and Backward Stability
The condition number of a nonsingular matrix \(\bA \in \fR^{n \times n}\) with respect to a given norm is \[ \begin{align} \kappa(\bA) = \|\bA\|\|\bA^{-1}\|. \end{align} \] In the Euclidean 2-norm, \(\kappa_2(\bA) = \sigma_{\max}(\bA) / \sigma_{\min}(\bA)\). The condition number measures the worst-case sensitivity of the exact solution \(\bx^*\) to relative perturbations in the input data \(\bA\) and \(\bb\).
For any computed approximation \(\hat{\bx}\) with residual \(\br = \bb - \bA\hat{\bx}\), the smallest relative perturbation size satisfying \((\bA + \delta\bA)\hat{\bx} = \bb + \delta\bb\) with \(\|\delta\bA\| \leq \eta \|\bA\|\) and \(\|\delta\bb\| \leq \eta \|\bb\|\) is given explicitly by \[ \begin{align} \eta(\hat{\bx}) = \frac{\|\br\|}{\|\bA\|\|\hat{\bx}\| + \|\bb\|}. \end{align} \] An algorithm is termed backward stable if for every problem instance, the computed solution satisfies \(\eta(\hat{\bx}) = O(\varepsilon_{\text{mach}})\).
From \((\bA + \delta\bA)\hat{\bx} = \bb + \delta\bb\), rearranging gives \(\br = \bb - \bA\hat{\bx} = \delta\bA\hat{\bx} - \delta\bb\). Taking norms and applying the triangle inequality yields \[ \begin{align} \|\br\| \leq \|\delta\bA\|\|\hat{\bx}\| + \|\delta\bb\| \leq \eta \|\bA\|\|\hat{\bx}\| + \eta \|\bb\| = \eta (\|\bA\|\|\hat{\bx}\| + \|\bb\|), \end{align} \] which establishes that \(\eta \geq \frac{\|\br\|}{\|\bA\|\|\hat{\bx}\| + \|\bb\|}\). To prove that this lower bound is attainable, define the rank-one matrix perturbation \(\delta\bA = \frac{\|\bA\|\|\hat{\bx}\|}{\|\bA\|\|\hat{\bx}\| + \|\bb\|} \frac{\br \hat{\bx}^T}{\|\hat{\bx}\|_2^2}\) and the vector perturbation \(\delta\bb = -\frac{\|\bb\|}{\|\bA\|\|\hat{\bx}\| + \|\bb\|} \br\). Direct substitution confirms that \((\bA + \delta\bA)\hat{\bx} = \bb + \delta\bb\) with perturbations achieving the minimal bound.
(The fundamental theorem of numerical analysis) Combining backward stability with problem conditioning yields the master error estimate of scientific computation: \[ \begin{align} \text{Forward Error} \lesssim \kappa(\bA) \times \text{Backward Error}. \end{align} \] A backward stable algorithm guarantees that the backward error is on the order of machine precision \(\varepsilon_{\text{mach}}\). If the computed solution is inaccurate, the loss of precision is entirely due to the ill-conditioning of the mathematical model \(\kappa(\bA)\), not a defect of the numerical algorithm.
14.3 The Solver Selection Hierarchy
To maximize computational efficiency and numerical stability, direct and iterative solvers should be selected according to the structural hierarchy of the coefficient matrix:
Symmetric Positive Definite (Dense): Cholesky factorization \(\bA = \bL\bL^T\) (\(\frac{1}{3}n^3\) flops, unconditionally backward stable without pivoting).
Symmetric Positive Definite (Large Sparse): Preconditioned Conjugate Gradient (PCG) (\(O(\text{nnz})\) per iteration, \(O(n)\) memory).
General Square (Dense): Gaussian Elimination with Partial Pivoting \(\bP\bA = \bL\bU\) (\(\frac{2}{3}n^3\) flops, backward stable in practice).
General Square (Large Sparse): Preconditioned GMRES or BiCGStab (\(O(\text{nnz})\) per iteration).
Overdetermined Full-Rank (\(m \geq n\)): Householder QR factorization \(\bA = \bQ\bR\) (\(2mn^2 - \frac{2}{3}n^3\) flops, avoids squaring the condition number).
Rank-Deficient or Ill-Conditioned (\(m \geq n\)): Singular Value Decomposition \(\bA = \bU\bsigma\bV^T\) (reveals numerical rank and computes minimum-norm solutions).
14.4 Iterative Refinement
Given a computed solution \(\hat{\bx}^{(0)}\) from an initial LU or Cholesky factorization, iterative refinement improves the solution accuracy:
Compute the residual with high precision: \(\br^{(k)} = \bb - \bA\hat{\bx}^{(k)}\).
Solve the error correction system using the existing factors: \(\bA \bd^{(k)} = \br^{(k)}\).
Update the solution: \(\hat{\bx}^{(k+1)} = \hat{\bx}^{(k)} + \bd^{(k)}\).
If \(\kappa(\bA) \varepsilon_{\text{mach}} < 1\) and the residual is evaluated with extended precision, iterative refinement converges to machine precision accuracy in a few steps.
Let \(\bA = \begin{pmatrix} 1 & 1 \\ 1 & 1 + 10^{-10} \end{pmatrix}\). Compute \(\bA^{-1}\) explicitly and evaluate \(\kappa_\infty(\bA)\).
Suppose \(\bA\bx = \bb\) is solved using a backward stable algorithm on a computer with \(\varepsilon_{\text{mach}} \approx 10^{-16}\). If \(\kappa_2(\bA) = 10^{11}\), estimate the number of correct decimal digits in the computed solution.
Prove that for any two nonsingular matrices \(\bA\) and \(\bB\), \(\kappa(\bA\bB) \leq \kappa(\bA)\kappa(\bB)\).
Show that for an orthogonal matrix \(\bQ\), \(\kappa_2(\bQ\bA) = \kappa_2(\bA)\) for any square matrix \(\bA\).