Stability of Gaussian elimination

\(\newcommand{\re}{\mathbb{R}}\)

Stability of Gaussian elimination#

%config InlineBackend.figure_format = 'svg'
import numpy as np
import matplotlib.pyplot as plt
from scipy.linalg import lu

Gaussian elimination with partial pivoting#

The rows are swapped in the pivoting step and the decomposition is of the form

\[ P A = L U \]

where \(P\) is a permutation matrix.

Now we test on many random matrices whose entries are normally distributed with mean zero and standard deviation \(m^{-1/2}\). We will use scipy LU routine since it will be faster than our own implementation. NOTE: This can take some time to run.

sizes = np.ceil(10**np.linspace(np.log10(5),np.log10(1000),50))
nsamples = 100
mvalues, rhovalues = [], []

for m1 in sizes:
    m = int(m1)
    for r in range(nsamples):
        A = np.random.normal(0.0, 1.0/np.sqrt(m), (m,m))
        P,L,U = lu(A)
        rho = np.abs(U).max() / np.abs(A).max()
        rhovalues.append(rho)
        mvalues.append(m)

Plot growth factor vs matrix size \(m\).

plt.loglog(mvalues, rhovalues, '.')
plt.loglog(mvalues, np.sqrt(mvalues), label='$m^{1/2}$')
plt.xlabel('$m$')
plt.ylabel('Growth factor, $\\rho$')
plt.legend()
plt.grid(True)
../_images/88003ad91d1348736169a97e8f72fa93f92208a062b549fcc28883da29fcef7d.svg

Compare this to Figure 22.1 from Trefethen and Bau.

Contents