5.5 Error bounds and iterative refinement
Every method, direct or iterative, produces an approximation \(\tilde {\mathbf x}\). Its residual \(\mathbf r = \mathbf b-A\tilde {\mathbf x}\) can be computed, while its error \(\mathbf x-\tilde {\mathbf x}\) cannot. A small residual seems reassuring; it may not be.
Example 5.5.1. The system \[x_1+2x_2 = 3,\qquad 1.0001x_1+2x_2 = 3.0001\] has the solution \((1,1)^T\). Find the residual of the poor approximation \(\tilde {\mathbf x} = (3,-0.0001)^T\).
Solution. \(\mathbf r = \mathbf b-A\tilde {\mathbf x} = (3-2.9998,\ 3.0001-3.0001)^T = (0.0002, 0)^T\), so \(\|\mathbf r\|_\infty = 0.0002\), while \(\|\mathbf x-\tilde {\mathbf x}\|_\infty = 2\). The two equations describe nearly parallel lines, and points far apart along them almost satisfy both.
Definition 5.5.2. The condition number of a nonsingular matrix \(A\) relative to a natural norm is \(K(A) = \|A\|\,\|A^{-1}\|\). Since \(1 = \|I\| = \|AA^{-1}\|\leq \|A\|\,\|A^{-1}\|\), always \(K(A)\geq 1\). A matrix is well-conditioned if \(K(A)\) is close to \(1\) and ill-conditioned if \(K(A)\) is large.
Theorem 5.5.3. If \(A\mathbf x = \mathbf b\) with \(\mathbf b\neq \mathbf 0\), and \(\tilde {\mathbf x}\) has residual \(\mathbf r = \mathbf b-A\tilde {\mathbf x}\), then \[\|\mathbf x-\tilde {\mathbf x}\|\leq \|A^{-1}\|\,\|\mathbf r\|\qquad \text {and}\qquad \frac {\|\mathbf x-\tilde {\mathbf x}\|}{\|\mathbf x\|}\leq K(A)\frac {\|\mathbf r\|}{\|\mathbf b\|} .\]
Proof. \(A(\mathbf x-\tilde {\mathbf x}) = \mathbf b-A\tilde {\mathbf x} = \mathbf r\), so \(\mathbf x-\tilde {\mathbf x} = A^{-1}\mathbf r\) and \(\|\mathbf x-\tilde {\mathbf x}\|\leq \|A^{-1}\|\,\|\mathbf r\|\). Also \(\|\mathbf b\| = \|A\mathbf x\|\leq \|A\|\,\|\mathbf x\|\), so \(\frac {1}{\|\mathbf x\|}\leq \frac {\|A\|}{\|\mathbf b\|}\); multiplying the two inequalities gives the second bound. □
So a small relative residual guarantees a small relative error only when \(K(A)\) is moderate. For the matrix of Example 5.5.1, \[A^{-1} = \begin {pmatrix}-10000&10000\\5000.5&-5000\end {pmatrix},\qquad \|A\|_\infty = 3.0001,\quad \|A^{-1}\|_\infty = 20000,\quad K_\infty (A) = 60002 .\] The bound gives \(\frac {\|\mathbf x-\tilde {\mathbf x}\|_\infty }{\|\mathbf x\|_\infty }\leq 60002\cdot \frac {0.0002}{3.0001} = 4\), consistent with the actual relative error \(2\).
The condition number also measures the sensitivity of the solution to errors in the data. It can be shown that if \(\tilde {\mathbf x}\) solves the perturbed system \((A+\delta A)\tilde {\mathbf x} = \mathbf b+\delta \mathbf b\) and \(\|A^{-1}\|\,\|\delta A\|<1\), then \[\frac {\|\mathbf x-\tilde {\mathbf x}\|}{\|\mathbf x\|}\leq \frac {K(A)}{1-K(A)\frac {\|\delta A\|}{\|A\|}}\left (\frac {\|\delta A\|}{\|A\|}+\frac {\|\delta \mathbf b\|}{\|\mathbf b\|}\right ).\] In particular, rounding the data to \(t\) significant digits can produce relative errors of about \(K(A)\times 10^{-t}\) in the solution: a system with \(K(A)\approx 10^m\) loses about \(m\) of the digits carried.
Example 5.5.4 (Hilbert matrices). The Hilbert matrix \(H_n\) has entries \(h_{ij} = \frac {1}{i+j-1}\). Its condition numbers grow very rapidly:
| \(n\) | \(3\) | \(5\) | \(8\) | \(10\) | \(12\) |
| \(K_\infty (H_n)\) | \(748\) | \(9.4\times 10^5\) | \(3.4\times 10^{10}\) | \(3.5\times 10^{13}\) | \(3.8\times 10^{16}\) |
Solving \(H_{10}\mathbf x = H_{10}\mathbf 1\) in double precision (about \(16\) digits) by Gaussian elimination with partial pivoting gives a solution whose components differ from \(1\) by up to \(5.9\times 10^{-4}\): about \(13\) of the \(16\) digits are lost, as \(K\approx 3.5\times 10^{13}\) predicts. For \(n = 12\) essentially no correct digits can be expected. Hilbert matrices arise naturally in least-squares fitting by polynomials, which is why that problem is usually solved by other means.
Iterative refinement
If \(\tilde {\mathbf x}\) was computed from an LU factorization, its accuracy can often be improved cheaply: compute the residual \(\mathbf r = \mathbf b-A\tilde {\mathbf x}\) (preferably in higher precision), solve \(A\mathbf y = \mathbf r\) using the same factors (cost \(O(n^2)\)), and set \(\tilde {\mathbf x}\leftarrow \tilde {\mathbf x}+\mathbf y\). Since \(A(\tilde {\mathbf x}+\mathbf y) = A\tilde {\mathbf x}+\mathbf r = \mathbf b\) exactly in exact arithmetic, each step removes most of the remaining error, as long as \(K(A)\) is not so large that the correction itself is meaningless.
Questions on this section
Stuck on something here? Ask below and it stays attached to this topic.