Condition Number

import numpy as np

The condition number of a matrix is a measure of how close to singular we are.

The definition of the condition number for a matrix \({\bf A}\) is:

\[\mathrm{cond}({\bf A}) = \| {\bf A} \| \| {\bf A}^{-1} \|\]

which requires defining a suitable norm.

For a matrix with a large condition number, then when solving the system \({\bf A}{\bf x} = {\bf b}\), a small change in \({\bf b}\) leads to a large change in \({\bf x}\). This makes algorithms like Gaussian elimination inaccurate for matrices with large condition numbers.

Condition number comes into play with roundoff error, since for every order of magnitude of the condition number, we lose about 1 digit of accuracy in our solution, or comparing an exact solution, \({\bf x}^\star\) with the computed solution, \({\bf x}\):

\[\frac{\| {\bf x}^\star - {\bf x} \|}{\| {\bf x}^\star \|} \approx \mathrm{cond}({\bf A}) \cdot \epsilon_\mathrm{machine}\]

where \(\epsilon_\mathrm{machine}\) is the machine epsilon of the floating point format.

Example

A nice example is given in the paper A simple example of an ill-conditioned matrix by G. J. Tee. Let’s work through that example here.

Let’s start with the matrix:

A = np.array([[11, 10, 14],
              [12, 11, -13],
              [14, 13, -66]], dtype=np.float64)
A

array([[ 11.,  10.,  14.],
       [ 12.,  11., -13.],
       [ 14.,  13., -66.]])

We can compute the determinant:

np.linalg.det(A)

1.0000000000000033

We see that it is non-singular. But it is quite close to being singular. To see this, we can perturb the middle element slightly:

B = A.copy()
B[1, 1] += 1./922.
B

array([[ 11.       ,  10.       ,  14.       ],
       [ 12.       ,  11.0010846, -13.       ],
       [ 14.       ,  13.       , -66.       ]])

np.linalg.det(B)

2.1316282072802905e-14

We see that that small perturbation essentially makes it singular (it would be algebraically singular, but our small number near machine epsilon is due to roundoff).

What is the condition number?

np.linalg.cond(A)

111039.4485341377

That’s quite large.

Let’s consider solving:

\[{\bf A}{\bf x} = {\bf b}\]

with

\[\begin{split}{\bf b} = \left ( \begin{array}{c} 1.001 \\ 0.999 \\ 1.001\end{array} \right )\end{split}\]

b = np.array([1.001, 0.999, 1.001])
x = np.linalg.solve(A, b)
x

array([-0.683,  0.843,  0.006])

Now consider a small perturbation:

\[\begin{split}{\bf b} \rightarrow \left ( \begin{array}{c} 1 \\ 1 \\ 1 \end{array} \right )\end{split}\]

We can think of this as:

\[{\bf A}({\bf x} + \delta{\bf x}) = {\bf b} + \delta {\bf b}\]

And since we are doing a small perturbation to \({\bf b}\), \(\delta {\bf b}\), we would expect the correction to \({\bf x}\), \(\delta {\bf x}\), to also be small.

Let’s solve the new system:

b2 = np.array([1.0, 1.0, 1.0])
x2 = np.linalg.solve(A, b2)
x2

array([ 1.00000000e+00, -1.00000000e+00, -1.66533454e-16])

so a \(1\%\) change in the righthand side, \({\bf b}\) leads to a \(\sim 240\%\) change in the solution, \({\bf x}\). This is not a small difference, and illustrates the behavior of a large condition number.