Gaussian elimination

Gaussian elimination#

Given \(A \in \re^{m \times m}\), we want to write it as

\[ A = L U \]

where \(L\) is lower triangular and \(U\) is upper triangular.

import numpy as np

Gaussian elimination without pivoting#

We now implement Algorithm 20.1.

def LU(A):
    m = A.shape[0]
    U = A.copy()
    L = np.eye(m)
    for k in range(m-1):
        for j in range(k+1,m):
            L[j,k] = U[j,k] / U[k,k]
            U[j,k:] = U[j,k:] - L[j,k] * U[k,k:]
    return L, U

Generate a random matrix and compute the decomposition

m = 10
A = 2 * np.random.rand(m,m) - 1
L, U = LU(A)

Test that \(A = L U\) is satisfied, compute \(\max| A - L U |\) over all elements

print('max|A-LU| = ', np.abs(A - L @ U).max())

max|A-LU| =  5.995204332975845e-15

Solution using LU decomposition#

We find the solution of

\[ LUx=b \]

in two steps.

Step 1: Solve

\[ Ly = b \]

using forward substitution

\[ y_i = \frac{1}{L_{ii}} \left[b_i - \sum_{j=0}^{i-1} L_{ij} y_j\right], \qquad i=0,1,\ldots,m-1 \]

Note that \(L_{ii}=1\), so we can drop the division step.

Step 2: Solve

\[ Ux = y \]

using backward substitution

\[ x_i = \frac{1}{U_{ii}}\left[y_i - \sum_{j=i+1}^{m-1} U_{ij} x_j \right], \qquad i=m-1,m-2,\ldots,0 \]

def LUSolve(L,U,b):
    m = L.shape[0]
    # solve Ly = b
    y = np.empty_like(b)
    for i in range(m):
        y[i] = b[i] - np.dot(L[i,0:i], y[0:i])
    # solve Ux = y
    x = np.empty_like(b)
    for i in range(m-1,-1,-1):
        x[i] = (y[i] - np.dot(U[i,i+1:m],x[i+1:m])) / U[i,i]
    return x

Create a random right hand side vector \(b\) and solve \(A x = b\).

b = np.random.rand(m)
x = LUSolve(L,U,b)
print('||Ax-b|| = ',np.linalg.norm(A@x - b, np.inf))

||Ax-b|| =  5.773159728050814e-15

Minimizing memory: We have used two vectors x and y to implement the solution but we can manage with only one. We solve as

Step 1: \(x = L^{-1} b\)
Step 2:

\[ x_i = \frac{1}{U_{ii}}\left[x_i - \sum_{j=i+1}^{m-1} U_{ij} x_j \right], \qquad i=m-1,m-2,\ldots,0 \]

In the second step, we can overwrite into \(x\) during backward substition since \(U\) is upper triangular. The following code does this:

def LUSolve(L,U,b):
    m = L.shape[0]
    x = np.empty_like(b)  # Single vector
    # solve Ly = b
    y = x                 # y is a pointer to x
    for i in range(m):
        y[i] = b[i] - np.dot(L[i,0:i], y[0:i])
    # solve Ux = y
    for i in range(m-1,-1,-1):
        x[i] = (y[i] - np.dot(U[i,i+1:m],x[i+1:m])) / U[i,i]
    return x

Further memory savings can be obtained by modifying the matrix \(A\) in place so that \(L\) is stored in the lower triangular part of \(A\) and \(U\) is stored in upper triangular part of \(A\), rather than creating two new matrices \(L,U\).

Example: failure of Gaussian elimination#

Solve \(A x = b\) with

\[\begin{split} A = \begin{bmatrix} 10^{-20} & 1 \\ 1 & 1 \end{bmatrix}, \qquad b = \begin{bmatrix} 1 \\ 0 \end{bmatrix} \end{split}\]

whose exact solution is

\[\begin{split} x = \frac{1}{1 - 10^{-20}} \begin{bmatrix} -1 \\ 1 \end{bmatrix} \approx \begin{bmatrix} -1 \\ 1 \end{bmatrix} \end{split}\]

A = np.array([[1.0e-20,  1.0],
              [1.0,      1.0]])
b = np.array([1.0, 0.0])
L, U = LU(A)
x = LUSolve(L,U,b)
print('x    = ', x)
print('Ax-b = ', A@x-b)

x    =  [0. 1.]
Ax-b =  [0. 1.]

The solution is totally wrong in the first component. If we interchange the two equations, then

A = np.array([[1.0,      1.0],
              [1.0e-20,  1.0]])
b = np.array([0.0, 1.0])
L, U = LU(A)
x = LUSolve(L,U,b)
print('x    = ', x)
print('Ax-b = ', A@x-b)

x    =  [-1.  1.]
Ax-b =  [0. 0.]

we get the correct solution. This is known as row pivoting.

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. In the code, we store a permutation vector rather than a matrix.

def PLU(A):
    m = A.shape[0]
    U = A.copy()
    L = np.eye(m)
    P = np.arange(m,dtype=int) # Permutation matrix
    for k in range(m-1):
        i = np.argmax(np.abs(U[k:m,k])) + k
        U[[k,i],k:m] = U[[i,k],k:m] # swap row i and k
        L[[k,i],0:k] = L[[i,k],0:k] # swap row i and k
        P[[k,i]] = P[[i,k]]         # swap row i and k
        for j in range(k+1,m):
            L[j,k] = U[j,k] / U[k,k]
            U[j,k:m] = U[j,k:m] - L[j,k] * U[k,k:m]
    return P, L, U

Now we test the above function for LU decomposition.

m = 3
A = np.random.rand(m, m)
P, L, U = PLU(A)
print("A =\n", A)
print("L =\n", L)
print("U =\n", U)
print("P =\n", P)

# Permute rows of A
PA = np.empty_like(A)
for i in range(m):
    PA[P[i],:] = A[i,:]
print("PA-LU =\n", PA - L@U)

A =
 [[0.48340333 0.70438618 0.73757663]
 [0.49333255 0.28432849 0.62347304]
 [0.16637833 0.36625499 0.43186578]]
L =
 [[1.         0.         0.        ]
 [0.97987315 1.         0.        ]
 [0.33725391 0.63498495 1.        ]]
U =
 [[0.49333255 0.28432849 0.62347304]
 [0.         0.42578032 0.12665213]
 [0.         0.         0.14117486]]
P =
 [1 0 2]
PA-LU =
 [[0. 0. 0.]
 [0. 0. 0.]
 [0. 0. 0.]]

From \(L U = P A\) and \(A x = b\) we get

\[ L U x = P A x = P b \]

which can be solved as

\[ L y = P b, \qquad U x = y \]

When we solve, we have to permute the right hand side vector \(b\).

def PLUSolve(L,U,P,b):
    m = L.shape[0]
    # solve Ly = Pb
    y = np.empty_like(b)
    for i in range(m):
        y[i] = b[P[i]] - np.dot(L[i,0:i], y[0:i])
    # solve Ux = y
    x = np.empty_like(b)
    for i in range(m-1,-1,-1):
        x[i] = (y[i] - np.dot(U[i,i+1:m], x[i+1:m])) / U[i,i]
    return x

Example: revisit failure case#

A = np.array([[1.0e-20,  1.0],
              [1.0,      1.0]])
P, L, U = PLU(A)
b = np.array([1.0, 0.0])
x = PLUSolve(L,U,P,b)
print('P    = ',P)
print('x    = ', x)
print('Ax-b = ', A@x-b)

P    =  [1 0]
x    =  [-1.  1.]
Ax-b =  [0. 0.]

Example#

We try a problem with a random matrix and right hand side.

m = 10
A = 2 * np.random.rand(m,m) - 1
P, L, U = PLU(A)
b = np.random.rand(m)
x = PLUSolve(L,U,P,b)
print('||Ax-b|| = ',np.linalg.norm(A@x - b, np.inf))

||Ax-b|| =  6.661338147750939e-16

Here is the permutation vector.

print('P = ',P)

P =  [8 6 5 4 2 9 3 0 7 1]