Stability of the LU factorization

Numerical analysis precept: consider some algorithm $A(X)$. Suppose the algorithm breaks down when $X=X_0$ (e.g., division by 0). Then the numerical error is typically large when $X\sim X_0$.

Consider this case:

$$A=\left(\begin{array}{cc}0& 1\\ 1& \pi \end{array}\right)$$

LU breaks down immediately because $a_{11}=0$.

What happens when we have a small pivot? Consider now:

$$\begin{array}{c}A=\left(\begin{array}{cc}ϵ& 1\\ 1& \pi \end{array}\right),\\ L=\left(\begin{array}{cc}1& 0\\ {ϵ}^{-1}& 1\end{array}\right),U=\left(\begin{array}{cc}ϵ& 1\\ 0& \pi -{ϵ}^{-1}\end{array}\right)\end{array}$$

• Floating point numbers on a computer can only be represented by a finite number of digits.
• Consider $\pi -{ϵ}^{-1}$.
• When $ϵ\ll 1$, ${ϵ}^{-1}\gg \pi$.
• For $ϵ$ sufficiently small, $\pi -{ϵ}^{-1}\equiv -{ϵ}^{-1}$ on a computer.
• Because of numerical roundoff errors, on a computer, we get

$$U=\left(\begin{array}{cc}ϵ& 1\\ 0& \cancel{\pi }-{ϵ}^{-1}\end{array}\right)$$

for $ϵ$ sufficiently small. $LU$ is no longer equal to $A$!

$$\left(\begin{array}{cc}1& 0\\ {ϵ}^{-1}& 1\end{array}\right)\left(\begin{array}{cc}ϵ& 1\\ 0& -{ϵ}^{-1}\end{array}\right)=\left(\begin{array}{cc}ϵ& 1\\ 1& 0\end{array}\right)\ne A$$

If we solve the linear system using numerical approximations to $L$ and $U$, we will get the wrong result.

A computer program cannot store all the digits of ${ϵ}^{-1}$. If ${ϵ}^{-1}\approx 10^{20}$, $\pi$ appears around digit 20. If you store only 16 digits of ${ϵ}^{-1}$, there is no room to store $\pi$.

Consequence: ${ϵ}^{-1}\equiv {ϵ}^{-1}-\pi$!

Solving linear systems using LU, Triangular factorization, LU algorithm, Existence of LU