Matrix Norms, Condition Numbers and Numerical Stability: Examples and Counterexamples

For $A_1$, the column sums and row sums are $\{1,4\}$, so [L2] gives $\parallel A_1{\parallel }_1=\parallel A_1{\parallel }_{\mathrm{\infty }}=4$; its inverse is $\mathrm{diag}(1,1\mathrm{/}4)$, with column and row sums $\{1,1\mathrm{/}4\}$, so $\parallel A_1^{-1}{\parallel }_1=\parallel A_1^{-1}{\parallel }_{\mathrm{\infty }}=1$.

For $A_2$, the column sums are $\{2,2+\epsilon \}$ and the row sums $\{2,2+\epsilon \}$, so [L2] gives $\parallel A_2{\parallel }_1=\parallel A_2{\parallel }_{\mathrm{\infty }}=2+\epsilon$.

By [L3], $A_2^{-1}={\epsilon }^{-1}\left(\begin{array}{cc}1+\epsilon & -1\\ -1& 1\end{array}\right)$, whose column sums and row sums are both $\{(2+\epsilon )\mathrm{/}\epsilon ,2\mathrm{/}\epsilon \}$, the first being larger; so [L2] gives $\parallel A_2^{-1}{\parallel }_1=\parallel A_2^{-1}{\parallel }_{\mathrm{\infty }}=(2+\epsilon )\mathrm{/}\epsilon$.

By [L1] and step 1.1, ${\kappa }_1(A_1)={\kappa }_{\mathrm{\infty }}(A_1)=4\cdot 1=4$: for every right-hand side $b\mathrm{\ne }0$ a relative perturbation of $b$ changes the solution by at most four times as much.

By [L1], ${\kappa }_1(A_2)={\kappa }_{\mathrm{\infty }}(A_2)=(2+\epsilon {)}^2\mathrm{/}\epsilon$, and at $\epsilon =1\mathrm{/}100$ this is $(201\mathrm{/}100{)}^2\cdot 100={201}^2\mathrm{/}100=40401\mathrm{/}100$.

The two systems of the same shape therefore have condition numbers $4$ and $40401\mathrm{/}100\approx 404$: the nearly singular one amplifies relative data perturbations about a hundred times more, and the contrast comes from the small determinant $\epsilon$, not from large entries.

The inverse of $H_2$ is $H_2^{-1}=\left(\begin{array}{cc}4& -6\\ -6& 12\end{array}\right)$: multiplying, $(H_2{H}_2^{-1}{)}_{00}=4-3=1$, $(H_2{H}_2^{-1}{)}_{01}=-6+6=0$, $(H_2{H}_2^{-1}{)}_{10}=2-2=0$ and $(H_2{H}_2^{-1}{)}_{11}=-3+4=1$.

By [L2], $\parallel H_2{\parallel }_{\mathrm{\infty }}=1+1\mathrm{/}2=3\mathrm{/}2$ and $\parallel H_2^{-1}{\parallel }_{\mathrm{\infty }}=\max \{4+6,6+12\}=18$, so [L1] gives ${\kappa }_{\mathrm{\infty }}(H_2)=(3\mathrm{/}2)\cdot 18=27$.

The inverse of $H_3$ is $H_3^{-1}=\left(\begin{array}{ccc}9& -36& 30\\ -36& 192& -180\\ 30& -180& 180\end{array}\right),$ as a direct multiplication verifies: the first column against the rows of $H_3$ gives $9-18+10=1$, $9\mathrm{/}2-12+15\mathrm{/}2=0$, $3-9+6=0$, and the other six entries similarly reduce to the identity matrix.

By [L2], $\parallel H_3{\parallel }_{\mathrm{\infty }}=1+1\mathrm{/}2+1\mathrm{/}3=11\mathrm{/}6$, and the row sums of $H_3^{-1}$ are $9+36+30=75$, $36+192+180=408$ and $30+180+180=390$, so $\parallel H_3^{-1}{\parallel }_{\mathrm{\infty }}=408$; [L1] gives ${\kappa }_{\mathrm{\infty }}(H_3)=(11\mathrm{/}6)\cdot 408=748$.

Steps 1.2 and 1.4 exhibit condition numbers $27$ and $748$ at orders $2$ and $3$: the Hilbert family is badly conditioned from the start.

For the node pair $(t_0,t_1)=(100,101)$, the determinant is $101-100=1$, and [L3] gives $V(t_0,t_1{)}^{-1}=\left(\begin{array}{cc}101& -100\\ -1& 1\end{array}\right)$.

By [L2], $\parallel V(100,101){\parallel }_{\mathrm{\infty }}=\max \{101,102\}=102$ and $\parallel V(100,101{)}^{-1}{\parallel }_{\mathrm{\infty }}=\max \{201,2\}=201$, so [L1] gives ${\kappa }_{\mathrm{\infty }}=102\cdot 201=20502$.

For the centered and scaled pair $(s_0,s_1)=(-1\mathrm{/}2,1\mathrm{/}2)$, the determinant is $1\mathrm{/}2-(-1\mathrm{/}2)=1$, and [L3] gives $V(s_0,s_1{)}^{-1}=\left(\begin{array}{cc}1\mathrm{/}2& 1\mathrm{/}2\\ -1& 1\end{array}\right)$.

By [L2], $\parallel V(s_0,s_1){\parallel }_{\mathrm{\infty }}=\max \{3\mathrm{/}2,3\mathrm{/}2\}=3\mathrm{/}2$ and $\parallel V(s_0,s_1{)}^{-1}{\parallel }_{\mathrm{\infty }}=\max \{1,2\}=2$, so [L1] gives ${\kappa }_{\mathrm{\infty }}=(3\mathrm{/}2)\cdot 2=3$.

The change of variable $t=s+201\mathrm{/}2$ maps $s_0$ to $100$ and $s_1$ to $101$, and interpolating the same data $(t_i,y_i)$ in the coordinate $t$ is the system $V(t_0,t_1)c=y$ while interpolating it in the coordinate $s$ is $V(s_0,s_1)c^{\mathrm{\prime }}=y$ with $c_0^{\mathrm{\prime }}=c_0+(201\mathrm{/}2)c_1$ and $c_1^{\mathrm{\prime }}=c_1$; the two systems represent the same polynomial task, yet their condition numbers are $20502$ and $3$.

Steps 1.2 and 1.4, read through step 2.1, show that centering and scaling the nodes improves the conditioning of the same interpolation problem by four orders of magnitude.

The residual is $r=b-A\widehat{x}=(1,\epsilon )-(1,0)=(0,\epsilon )$, so $\parallel r{\parallel }_2=\epsilon$, while $\parallel b{\parallel }_2=\sqrt{1+{\epsilon }^2}$; the relative residual is $\epsilon \mathrm{/}\sqrt{1+{\epsilon }^2}$, which tends to $0$ as $\epsilon \to 0^{+}$.

The error is $e=x-\widehat{x}=(0,1)$, so $\parallel e{\parallel }_2=1$ and the relative forward error is $1\mathrm{/}\sqrt{2}$, independent of $\epsilon$.

Steps 1.1 and 1.2 contradict the claimed implication: the relative residual can be made arbitrarily small while the relative forward error stays at $1\mathrm{/}\sqrt{2}$.

The bounds of [L1] are consistent with the witness and explain it: $\parallel e{\parallel }_2=1=(1\mathrm{/}\epsilon )\cdot \epsilon =\parallel A^{-1}{\parallel }_2\parallel r{\parallel }_2$ attains the first bound exactly, and the condition number $1\mathrm{/}\epsilon$ of [L2] is what converts the tiny residual $\epsilon$ into order-one forward error.

Steps 2.1 and 2.2 refute the false statement: a tiny residual guarantees only a tiny forward error up to the condition number, not by itself.

The residual is $r=b-A\widehat{x}=(0,\epsilon )$ and $\parallel r{\parallel }_2=\epsilon$; also $\parallel A{\parallel }_2=1$, $\parallel \widehat{x}{\parallel }_2=1$ and $\parallel b{\parallel }_2=\sqrt{1+{\epsilon }^2}$.

The forward error is $\parallel x-\widehat{x}{\parallel }_2=\parallel (0,1){\parallel }_2=1$ and the relative forward error $1\mathrm{/}\sqrt{2}$, order one even as $\epsilon \to 0^{+}$.

The conditioning conversion. For the problem map $f(b)=A^{-1}b$ the absolute condition number is ${\kappa }_{\mathrm{a}\mathrm{b}\mathrm{s}}=\parallel A^{-1}{\parallel }_2=1\mathrm{/}\epsilon$, attained in the second coordinate, and the computed answer is exact for the perturbed input $\tilde{b}:=A\widehat{x}=(1,0)=b+h$ with $h=(0,-\epsilon )$. By [A1], $f(b+h)-f(b)=A^{-1}h=(0,-1)$, and direct substitution gives $1\mathrm{/}\sqrt{2}={\kappa }_{\mathrm{r}\mathrm{e}\mathrm{l}}\parallel h{\parallel }_2\mathrm{/}\parallel b{\parallel }_2$.

By [L1], the normwise backward error is ${\eta }_2(\widehat{x})=\epsilon \mathrm{/}(1+\sqrt{1+{\epsilon }^2})\le \epsilon \mathrm{/}2=u\mathrm{/}2$. Thus the one constant $C=1\mathrm{/}2$ works for every input in the declared singleton class $\mathcal{I}$, so the algorithm is normwise backward stable there in the specialized sense of Forward and backward stability for a problem family under an arithmetic model.

Steps 2.1, 1.2 and 1.3 show that this backward-stable algorithm on an ill-conditioned problem class has order-one forward error even though its backward error is order $u$; [L2] records the general first-order conditioning principle illustrated exactly here. Backward stability controls the data perturbation, and only conditioning controls the answer.

The two formulas agree exactly: $z\mathrm{/}(1+w)=(1-w)(1+w)\mathrm{/}(1+w)=1-w$, because $z=1-w^2$ and $1+w>0$.

For $z={10}^{-8}$, the exact value is $f(z)=1-w\approx 5\cdot {10}^{-9}$: indeed $(1-5\cdot {10}^{-9}{)}^2=1-{10}^{-8}+25\cdot {10}^{-18}>1-z$, so $w<1-5\cdot {10}^{-9}$, while $(1-6\cdot {10}^{-9}{)}^2=1-12\cdot {10}^{-9}+36\cdot {10}^{-18}<1-z$, so $w>1-6\cdot {10}^{-9}$; hence $1-w\in (5\cdot {10}^{-9},6\cdot {10}^{-9})$.

The direct formula. By [L1] the computed values are $\mathrm{fl}(w)=w(1+{\delta }_1)$ and $\mathrm{fl}(1-\mathrm{fl}(w))=(1-w(1+{\delta }_1))(1+{\delta }_2)$ with $\mathrm{\mid }{\delta }_1\mathrm{\mid },\mathrm{\mid }{\delta }_2\mathrm{\mid }\le u$; subtracting the exact value $1-w$ leaves an error of absolute value $\mathrm{\mid }{\delta }_2(1-w)-{\delta }_{1}w(1+{\delta }_2)\mathrm{\mid }\le u(1-w)+uw(1+u)$.

The reformulated formula. By [L1], the computed values are $\mathrm{fl}(1+w)=(1+w)(1+{\delta }_3)$ for the denominator and $\mathrm{fl}(z\mathrm{/}\mathrm{fl}(1+w))=\frac{z}{(1+w)(1+{\delta }_3)}(1+{\delta }_4)$ for the quotient, with $\mathrm{\mid }{\delta }_3\mathrm{\mid },\mathrm{\mid }{\delta }_4\mathrm{\mid }\le u$.

By [L2], the relative error of the direct formula is at most $u+uw(1+u)\mathrm{/}(1-w)$; with $w\approx 1$ and $1-w\approx 5\cdot {10}^{-9}$ this is about $u(1+2\cdot {10}^8)$, roughly $2\cdot {10}^{8}u$: with $u=2^{-53}$ the relative forward error can reach about $2\cdot {10}^{-8}$, losing about eight of the roughly sixteen decimal digits of relative precision available in binary64.

By [L2], the relative error of the reformulated value is $\mathrm{\mid }(1+{\delta }_4)\mathrm{/}(1+{\delta }_3)-1\mathrm{\mid }=\mathrm{\mid }{\delta }_4-{\delta }_3\mathrm{\mid }\mathrm{/}\mathrm{\mid }1+{\delta }_3\mathrm{\mid }\le 2u\mathrm{/}(1-u)$, which is of order $u$, not of order ${10}^{8}u$: no subtraction of nearby quantities occurs, and the reformulation evaluates $f$ to full relative accuracy.

Steps 2.1 and 2.2 compare the same function evaluated two ways: the direct formula amplifies roundoff by the cancellation factor $1\mathrm{/}(1-w)\approx 2\cdot {10}^8$, while the rationalized formula keeps the error at order $u$.

With $n=2$ and $2u<1$, [L1] applies and yields $s_2=x_0{y}_0(1+{\delta }_0)+x_1{y}_1(1+{\delta }_1)$ with $\mathrm{\mid }{\delta }_0\mathrm{\mid },\mathrm{\mid }{\delta }_1\mathrm{\mid }\le {\gamma }_2$, which for $x=(1,2)$ and $y=(3,4)$ reads $s_2=3(1+{\delta }_0)+8(1+{\delta }_1)$.

The exact dot product is $x\cdot y=1\cdot 3+2\cdot 4=11$, and the absolute summand total is $\mathrm{\mid }{x}_0{y}_0\mathrm{\mid }+\mathrm{\mid }{x}_1{y}_1\mathrm{\mid }=3+8=11$, while $\parallel x{\parallel }_1\parallel y{\parallel }_1=(1+2)(3+4)=21$.

By the normwise bound of [L1] applied to the numbers of step 1.2, $\mathrm{\mid }{s}_2-11\mathrm{\mid }\le {\gamma }_2\cdot 11\le {\gamma }_2\cdot 21$; substituting $u=2^{-53}$ gives ${\gamma }_2=2^{-52}\mathrm{/}(1-2^{-52})$, so $\mathrm{\mid }{s}_2-11\mathrm{\mid }\le 11\cdot 2^{-52}\mathrm{/}(1-2^{-52})$, a bound of order ${10}^{-15}$.

Step 1.1 verifies the componentwise perturbation form and step 2.1 the normwise bound for this short vector pair: the computed dot product agrees with the exact dot product of slightly perturbed summands and lies within $11{\gamma }_2$ of $11$.

The columns of $A$ are orthogonal and have norms $1$ and $\epsilon$. Writing $A=U\mathrm{\Sigma }{V}^{\mathsf{T}}$ with $U=[e_1{e}_2]$, $V=I_2$, and the $2\times 2$ matrix $\mathrm{\Sigma }=\mathrm{diag}(1,\epsilon )$ shows that the singular values are the diagonal entries $1$ and $\epsilon$ of $\mathrm{\Sigma }$.

By [L2], Gram-Schmidt applied to the columns $a_1=(1,0,0{)}^{\mathsf{T}}$ and $a_2=(0,\epsilon ,0{)}^{\mathsf{T}}$ returns the already orthogonal unit vectors $q_1=(1,0,0{)}^{\mathsf{T}}$ and $q_2=(0,1,0{)}^{\mathsf{T}}$ after scaling by their norms. Writing $Q=[q_1{q}_2]$, one has $A=QR$ with $R=Q^{\mathsf{T}}A=\mathrm{diag}(1,\epsilon )$.

Hence ${\kappa }_2(A)=1\mathrm{/}\epsilon ={10}^8$, and $A^{\mathsf{T}}A$ is the $2\times 2$ diagonal matrix $\mathrm{diag}(1,{\epsilon }^2)$, whose singular values are $1$ and ${\epsilon }^2$.

By step 1.2, $R=\mathrm{diag}(1,\epsilon )$, so its singular values are $1$ and $\epsilon$ and therefore ${\kappa }_2(R)=1\mathrm{/}\epsilon ={10}^8$.

By [L1], ${\kappa }_2(A^{\mathsf{T}}A)={\kappa }_2(A{)}^2={10}^{16}$: forming the normal equations squares the conditioning of the least-squares problem.

Steps 3.1 and 2.2 compare the two routes on the same badly scaled problem: the normal-equation matrix carries condition number ${10}^{16}$, the square of ${\kappa }_2(A)$, while the Gram–Schmidt triangular factor $R$ carries ${10}^8$, the original conditioning of the problem.

The inverse of the upper triangular matrix $A$ is $A^{-1}=\left(\begin{array}{ccc}1& -2\mathrm{/}3& 1\mathrm{/}6\\ 0& 1\mathrm{/}3& -1\mathrm{/}12\\ 0& 0& 1\mathrm{/}4\end{array}\right),$ obtained by back substitution from $AX=I_3$.

By [L2], the column sums of $A$ are $1,5,5$ so $\parallel A{\parallel }_1=5$, and the column sums of $A^{-1}$ are $1,1,1\mathrm{/}2$ so $\parallel A^{-1}{\parallel }_1=1$; [L1] gives ${\kappa }_1(A)=5\cdot 1=5$.

By [L2], the row sums of $A$ are $3,4,4$ so $\parallel A{\parallel }_{\mathrm{\infty }}=4$, and the row sums of $A^{-1}$ are $11\mathrm{/}6,5\mathrm{/}12,1\mathrm{/}4$ so $\parallel A^{-1}{\parallel }_{\mathrm{\infty }}=11\mathrm{/}6$; [L1] gives ${\kappa }_{\mathrm{\infty }}(A)=4\cdot 11\mathrm{/}6=22\mathrm{/}3$.

The rescaled matrix is $DA=\left(\begin{array}{ccc}10& 20& 0\\ 0& 3& 1\\ 0& 0& 4\end{array}\right)$, with inverse $(DA{)}^{-1}=A^{-1}{D}^{-1}=\left(\begin{array}{ccc}1\mathrm{/}10& -2\mathrm{/}3& 1\mathrm{/}6\\ 0& 1\mathrm{/}3& -1\mathrm{/}12\\ 0& 0& 1\mathrm{/}4\end{array}\right)$.

By [L2], the column sums of $DA$ are $10,23,5$ so $\parallel DA{\parallel }_1=23$, and the column sums of $(DA{)}^{-1}$ are $1\mathrm{/}10,1,1\mathrm{/}2$ so $\parallel (DA{)}^{-1}{\parallel }_1=1$; [L1] gives ${\kappa }_1(DA)=23\cdot 1=23$.

Steps 1.2 and 1.3 give ${\kappa }_1(A)=5\mathrm{\ne }22\mathrm{/}3={\kappa }_{\mathrm{\infty }}(A)$: the same matrix has different condition numbers in different induced norms.

Since $D$ is invertible, the systems $Ax=b$ and $(DA)x=Db$ have exactly the same solution set, yet step 1.5 shows ${\kappa }_1(DA)=23\mathrm{\ne }5={\kappa }_1(A)$: diagonal rescaling changes the condition number even though it changes nothing about the set of solutions.

Steps 2.1 and 2.2 refute both halves of the false statement: condition number depends on the declared norm and on the scaling of the equations.

At the tolerance ${\tau }_1=3\mathrm{/}4$, the only singular value strictly above the threshold is ${\sigma }_1=1$, so [L1] gives ${\mathrm{rank}}_{{\tau }_1}(A)=\mathrm{\#}\{j:{\sigma }_j>3\mathrm{/}4\}=1$.

At the tolerance ${\tau }_2=1\mathrm{/}4$, the singular values strictly above the threshold are ${\sigma }_1=1$ and ${\sigma }_2=1\mathrm{/}2$, so [L1] gives ${\mathrm{rank}}_{{\tau }_2}(A)=\mathrm{\#}\{j:{\sigma }_j>1\mathrm{/}4\}=2$.

The singular value ${\sigma }_2=1\mathrm{/}2$ lies strictly between the two tolerances, ${\tau }_2=1\mathrm{/}4<{\sigma }_2<{\tau }_1=3\mathrm{/}4$, which is exactly why it is counted at ${\tau }_2$ but not at ${\tau }_1$.

Steps 1.1, 1.2 and 1.3 exhibit one matrix and two declared tolerances with numerical ranks $1$ and $2$: numerical rank is a function of the matrix together with the tolerance, not of the matrix alone, which refutes the false statement.

Overflow. The representable numbers $x=(3\mathrm{/}2)\cdot 2^{1023}$ and $y=2$ have product $x\cdot y=3\cdot 2^{1023}$, which exceeds $N_{\max }=(2-2^{-52})2^{1023}$ because $3>2-2^{-52}$; the computed result is $+\mathrm{\infty }$, which is not a real number of the form $(x\cdot y)(1+\delta )$ for any real $\delta$.

Underflow. For the smallest positive subnormal $x=2^{-1074}$, the exact quotient $x\mathrm{/}2=2^{-1075}$ lies strictly between $0$ and the smallest positive subnormal, and round-to-nearest returns $0$; but $(x\mathrm{/}2)(1+\delta )=2^{-1075}(1+\delta )$ is positive for every real $\delta$ with $\mathrm{\mid }\delta \mathrm{\mid }\le u$, since $1+\delta \ge 1-u>0$, so the computed value $0$ is not of the required form.

Invalid operation. The exact value of $0\mathrm{/}0$ is not a real number, so no real $\delta$ can satisfy $\mathrm{fl}(0\mathrm{/}0)=(0\mathrm{/}0)(1+\delta )$; the computed result is $\mathrm{NaN}$, which is not a real number at all.

Each of the three exhibits produces a value outside the model of [L1]: infinity, the lost subnormal result, and not-a-number respectively.

Step 2.1 refutes the false statement: the exclusions in [L1] are not decorative, and any statement proved under the relative model must carry the hypothesis that its operations stay in the normal range.