Question

It is important that numerical algorithms be stable and efficient. Efficiency is a property of an algorithm, but stability can be a property of the system itself.

Example

\[\begin{split}\left(\begin{matrix}8&6&4&1\\1&4&5&1\\8&4&1&1\\1&4&3&6\end{matrix}\right)x = \left(\begin{matrix}19\\11\\14\\14\end{matrix}\right)\end{split}\]

A = np.array([[8,6,4,1],[1,4,5,1],[8,4,1,1],[1,4,3,6]])
b = np.array([19,11,14,14])
la.solve(A,b)

array([ 1.,  1.,  1.,  1.])

b = np.array([19.01,11.05,14.07,14.05])
la.solve(A,b)

array([-2.34 ,  9.745, -4.85 , -1.34 ])

Note that the tiny perturbations in the outcome vector \(b\) cause large differences in the solution! When this happens, we say that the matrix \(A\) ill-conditioned. This happens when a matrix is ‘close’ to being singular (i.e. non-invertible).

Condition Number

A measure of this type of behavior is called the condition number. It is defined as:

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

In general, it is difficult to compute.

Fact:

\[cond(A) = \frac{\lambda_1}{\lambda_n}\]

where \(\lambda_1\) is the maximum singular value of \(A\) and \(\lambda_n\) is the smallest. The higher the condition number, the more unstable the system. In general if there is a large discrepancy between minimal and maximal singular values, the condition number is large.

Example

U, s, V = np.linalg.svd(A)
print(s)
print(max(s)/min(s))

[ 15.5457   6.9002   3.8363   0.0049]
3198.6725812

Preconditioning

We can sometimes improve on this behavior by ‘pre-conditioning’. Instead of solving

\[Ax=b\]

we solve

\[D^{-1}Ax=D^{-1}b\]\[where :math:`D^{-1}A` has a lower condition number than :math:`A`\]

itself.

Preconditioning is a very involved topic, quite out of the range of this course. It is mentioned here only to make you aware that such a thing exists, should you ever run into an ill-conditioned problem!