Conditioning

try:
    import executable_engineering as exe 
except ImportError:
    %pip install -q executable_engineering
    import executable_engineering as exe

import numpy as np

Stability and Uncertainty¶

Previously, we discussed how roundoff and truncation errors introduce uncertainty into numbers stored on a computer. The conditioning of a function relates how uncertainty in the inputs propagates to uncertainty in the outputs.

For a general function $f(\mathbf{x}) = \mathbf{b}$ :

Well-conditioned: Small changes in inputs ( $\mathbf{x}$ ) lead to proportionally small changes in outputs ( $\mathbf{b}$ ).
Ill-conditioned: Small changes in inputs could result in massive, explosive changes in the outputs, and vice-versa.

The condition of a function is a fundamental property of the mathematical transformation itself, not just the specific input/output values.

Recall however that numerical methods have inherent uncertainty in numbers due to finite precision representation. Ill-conditioned functions can catastrophically disrupt finding solutions of linear systems!

Approaching Singularity¶

In the solvability section, we saw two strict cases based on the determinant $|\mathbf{A}|$ :

$|\mathbf{A}| = -30$ and the system has a unique exact solution.
$|\mathbf{A}| = 0$ and the system is unsolvable (parallel lines).

But what happens in the “grey area” when $|\mathbf{A}| \rightarrow 0$ ?

Let’s adjust the price of tables in our fundraising problem to be $20-\epsilon$ .

\begin{aligned} 20 c + (20-\epsilon) t &= 700 \\ c + t &= 20 \end{aligned}

(1)

The determinant is $|\mathbf{A}| = 20(1) - (20-\epsilon)(1) = \epsilon$ . Let’s plot this when $\epsilon$ is small!

# Change epsilon to see how wildly the solution point swings!
epsilon = 10

A = [[20, 20 - epsilon], 
     [1, 1]]
b = [700, 20]

fig = exe.visual_solve_2d(A, b, method='exact')
fig.show()

The Consequence of Small $\epsilon$ ¶

If we solve this algebraically by substituting $c = 20-t$ into the first equation:

\begin{aligned} 20(20 - t) + (20-\epsilon) t &= 700 \\ -\epsilon t &= 300 \\ t &= -\frac{300}{\epsilon} \end{aligned}

(2)

Notice what happens as $\epsilon$ gets smaller:

$\epsilon$	$t$	$c$
-10	30	-10
1	-300	320
0.1	-3000	3020
1e-15	-3e+17	3e+17

This should sound alarm bells! Technically, for any $\epsilon \neq 0$ , the lines intersect and the system is invertible. However, the solution point swings wildly toward infinity as the lines become nearly parallel.

This is simple example of an ill-conditioned system (matrix). If $\mathbf{A}$ is ill-conditioned, the transformation depends highly on precision. When numerical methods are applied to such systems, they inherantly add round-off error which makes the solution unstable (and also interferes with algorithms that check the interations for convergence).

Distortion, not Amplification¶

A common misconception is that ill-conditioning is just “amplification” of error. But if a function simply amplified error, its inverse would reduce error.

This is not the case! The conditioning of a forward function and its inverse are quantifiably the same. Conditioning represents the maximum distortion or stretching caused by the transformation in any direction.

Well-conditioned matrix distortion¶

For a linear transformation $\mathbf{A}\mathbf{x} = \mathbf{b}$ , we can visualize this distortion by applying the matrix $\mathbf{A}$ (and its inverse $\mathbf{A}^{-1}$ ) to a perfect unit circle.

Let’s look at a well-conditioned matrix:

\mathbf{A} = \begin{pmatrix} 2 & 0 \\ 0 & 2 \end{pmatrix}

(3)

Notice how the circle is scaled but remains perfectly proportional. The inverse simply shrinks it back. The mapping is smooth and nowhere would two values of $x$ overlap.

A_well = [[2, 0], 
          [0, 2]]

print('Condition number: ', np.linalg.cond(A_well))
fig1 = exe.visualize_conditioning(A_well)
fig1.show()

Ill-conditioned matrix distortion¶

Now, let’s look at an ill-conditioned matrix where the rows are nearly parallel:

\mathbf{A} = \begin{pmatrix} 1 & 1 \\ 1 & 1.1 \end{pmatrix}

(4)

A_ill = [[1, 1], 
         [1, 1.1]]

print('Condition number: ', np.linalg.cond(A_ill))
fig2 = exe.visualize_conditioning(A_ill)
fig2.show()

The Condition Number ( $\kappa$ )¶

The severe stretching seen above is quantified by the Condition Number ( $\kappa$ ). It is the maximum ratio of the relative error in $\mathbf{x}$ to the relative error in $\mathbf{b}$ .

Let’s perturb our inputs and outputs by a small error $\delta \mathbf{x}$ and $\delta \mathbf{b}$ :

\begin{aligned} \mathbf{x} & \rightarrow \mathbf{x} + \delta \mathbf{x}\\ \mathbf{b} & \rightarrow \mathbf{b} + \delta \mathbf{b} \end{aligned}

(5)

Derivation¶

Substitute the perturbed vectors into our linear system $\mathbf{A} \mathbf{x} = \mathbf{b}$ :

\begin{aligned} \mathbf{A} (\mathbf{x} + \delta \mathbf{x}) &= \mathbf{b} + \delta \mathbf{b} \\ \mathbf{A} \mathbf{x} + \mathbf{A} \delta \mathbf{x} &= \mathbf{b} + \delta \mathbf{b} \end{aligned}

(6)

Since we know $\mathbf{A} \mathbf{x} = \mathbf{b}$ , these terms cancel out, leaving us with a direct relationship between the uncertainties:

\begin{aligned} \mathbf{A} \delta \mathbf{x} &= \delta \mathbf{b} \\ \delta \mathbf{x} &= \mathbf{A}^{-1} \delta \mathbf{b} \end{aligned}

(7)

Ratio of Relative Errors¶

We want to find the ratio of the relative error in the input to the relative error in the output using the norm (magnitude) of the vectors:

\frac{\|\delta \mathbf{x}\|}{\|\mathbf{x}\|} \Bigg/ \frac{\|\delta \mathbf{b}\|}{\|\mathbf{b}\|} = \frac{\|\mathbf{A}^{-1}\delta \mathbf{b}\|}{\|\delta \mathbf{b}\|} \cdot \frac{\|\mathbf{A} \mathbf{x}\|}{\|\mathbf{x}\|}

(8)

The maximum possible value of this ratio across all possible vectors is defined as the condition number $\kappa$ :

\kappa = \max \left( \frac{\|\mathbf{A}^{-1}\delta \mathbf{b}\|}{\|\delta \mathbf{b}\|} \right) \cdot \max \left( \frac{\|\mathbf{A} \mathbf{x}\|}{\|\mathbf{x}\|} \right)

(9)

\kappa = \|\mathbf{A}^{-1}\| \|\mathbf{A}\|

(10)

NB: This result introduces the concept of the norm (magnitude) of a matrix (defined below) and an identity $\|\mathbf{A}^{-1}\delta \mathbf{b}\| \leq \|\mathbf{A}^{-1} \| \cdot \|\delta \mathbf{b}\|$ .

Application¶

From this definition, it is obvious that the condition number is identical for both the forward ( $\mathbf{A}$ ) and backward ( $\mathbf{A}^{-1}$ ) transformations!

It guarantees bounds on the relative error in both directions:

\frac{\|\delta \mathbf{b}\|}{\|\mathbf{b}\|} \leq \kappa \frac{\|\delta \mathbf{x}\|}{\|\mathbf{x}\|} \quad \text{and} \quad \frac{\|\delta \mathbf{x}\|}{\|\mathbf{x}\|} \leq \kappa \frac{\|\delta \mathbf{b}\|}{\|\mathbf{b}\|}

(11)

If $\kappa$ is large (ill-conditioned), small errors in your data $\mathbf{b}$ can result in massive, unpredictable errors in your computed solution $\mathbf{x}$ and due to finite precision, small errors will always be present!!

Calculation via Eigenvalues¶

The formal definition requires calculating the inverse $\mathbf{A}^{-1}$ , which is exactly what we are trying to avoid doing blindly! We want to know the condition number before attempting to solve the system.

Geometrically, the maximal stretching of a matrix is defined by its largest eigenvalue ( $\lambda_{\max}$ ). The maximal stretching of the inverse is the smallest eigenvalue of the original matrix ( $\lambda_{\min}$ ).

Therefore, for symmetric matrices, we can calculate $\kappa$ using the eigenvalues:

\kappa = \frac{|\lambda_{\max}|}{|\lambda_{\min}|}

(12)

In numpy this is calculated using the $np.linalg.cond()$ function.