LU Decomposition - Executable Engineering: Interactive Numerical Methods in Python

import numpy as np
import matplotlib.pyplot as plt
from scipy.linalg import lu, solve
from scipy.sparse import diags
from scipy.sparse.linalg import spsolve

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

What is LU Decomposition?¶

Any square matrix can be decomposed into the product of a lower triangular matrix ( $\mathbf{L}$ ) and an upper triangular matrix ( $\mathbf{U}$ ):

\mathbf{A} = \mathbf{L}\mathbf{U}

(1)

Where:

$\mathbf{L}$ has 1 s on the diagonal and the multipliers used during elimination below the diagonal.
$\mathbf{U}$ is exactly the upper-triangular matrix that remains after forward Gaussian elimination!

This means LU Decomposition is literally just Gaussian Elimination, but instead of throwing away the elimination steps, we save them inside the $\mathbf{L}$ matrix.

Assembling the L Matrix¶

Let’s look back at the $3 \times 3$ system we solved in the Gaussian Elimination chapter. Our original coefficient matrix was:

\mathbf{A} = \begin{bmatrix} 4 & 3 & -5 \\ -2 & -4 & 5 \\ 8 & 8 & 0 \end{bmatrix}

(2)

During elimination, we used three multipliers:

$m_{2,1} = -0.5$
$m_{3,1} = 2$
$m_{3,2} = -0.8$

The lower triangular matrix $\mathbf{L}$ is constructed by simply placing these exact multipliers into their respective $(i, j)$ positions below the main diagonal (with 1s on the diagonal):

\mathbf{L} = \begin{bmatrix} 1 & 0 & 0 \\ m_{2,1} & 1 & 0 \\ m_{3,1} & m_{3,2} & 1 \end{bmatrix} = \begin{bmatrix} 1 & 0 & 0 \\ -0.5 & 1 & 0 \\ 2 & -0.8 & 1 \end{bmatrix}

(3)

The upper triangular matrix $\mathbf{U}$ is simply the final matrix produced at the end of the forward elimination phase:

\mathbf{U} = \begin{bmatrix} 4 & 3 & -5 \\ 0 & -2.5 & 2.5 \\ 0 & 0 & 12 \end{bmatrix}

(4)

If we multiply these two matrices together, we perfectly recover our original $\mathbf{A}$ matrix! $$ \mathbf{L}\mathbf{U} = \begin{bmatrix} 1 & 0 & 0 \ -0.5 & 1 & 0 \ 2 & -0.8 & 1 \end{bmatrix}

\begin{bmatrix} 4 & 3 & -5 \\ 0 & -2.5 & 2.5 \\ 0 & 0 & 12 \end{bmatrix}

(5)

= \begin{bmatrix} 4 & 3 & -5 \ -2 & -4 & 5 \ 8 & 8 & 0 \end{bmatrix} = \mathbf{A} $$

Why store the elimination steps?¶

If we decompose $\mathbf{A}$ into $\mathbf{L}\mathbf{U}$ , solving $\mathbf{A}\mathbf{x} = \mathbf{b}$ becomes a two-step process:

\begin{aligned} (\mathbf{L}\mathbf{U})\mathbf{x} &= \mathbf{b} \\ \mathbf{L}(\mathbf{U}\mathbf{x}) &= \mathbf{b} \end{aligned}

(6)

Let $\mathbf{y} = \mathbf{U}\mathbf{x}$ . First, solve $\mathbf{L}\mathbf{y} = \mathbf{b}$ using forward substitution (since $\mathbf{L}$ is lower triangular).
Then, solve $\mathbf{U}\mathbf{x} = \mathbf{y}$ using backward substitution (since $\mathbf{U}$ is upper triangular).

The advantage: If you need to solve the same system for 100 different $\mathbf{b}$ vectors (e.g. testing different loads on a bridge), you only have to perform the expensive $\mathcal{O}(n^3)$ elimination once! The forward and backward substitutions only take $\mathcal{O}(n^2)$ .

PLU Decomposition in Python¶

In reality, standard LU decomposition fails if a zero pivot is encountered. Robust algorithms use pivoting (row swapping), which introduces a Permutation matrix ( $\mathbf{P}$ ):

\mathbf{A} = \mathbf{P}\mathbf{L}\mathbf{U}

(7)

Let’s decompose a matrix using Scipy!

A = np.array([[ 4, -2,  1],
              [-2,  4, -2],
              [ 1, -2,  4]])

# Calculate the PLU decomposition
P, L, U = lu(A)

print("Permutation Matrix P:\n", P)
print("\nLower Triangular L:\n", L)
print("\nUpper Triangular U:\n", U)

# Reconstruct A to verify
print("\nReconstructed A (P @ L @ U):\n", P @ L @ U)

Computational Complexity¶

Although forward and back substitution both occur in $\mathcal{O}(n^2)$ , the actual $\mathbf{L}\mathbf{U}$ decomposition itself still requires $\mathcal{O}(n^3)$ operations.

However, there are a number of major reasons why this is the preferred approach over standard Gaussian elimination:

The $\mathcal{O}(n^3)$ decomposition only needs to occur once, and the system can then be solved for multiple $\mathbf{b}$ vectors efficiently in $\mathcal{O}(n^2)$ .
There are highly parallelized algorithms for distributed, GPU computing which can speed up the decomposition wall-time significantly.
The $\mathbf{L}$ and $\mathbf{U}$ factors preserve sparsity patterns inherently better than calculating exact inverses.
The factors are useful for other operations (e.g., calculating the determinant of $\mathbf{A}$ is simply the product of the diagonals of $\mathbf{U}$ ).

Direct solvers scaling with $\mathcal{O}(n^3)$ is terrible. This implies doubling the number of variables increases the computational time by a factor of ~8!

Here sparsity comes to your rescue. Depending on the type of sparsity, one can typically get to algorithms ~ $\mathcal{O}(n^2)$ or better (depending on the structure of the sparsity)!

Sparse LU and the danger of “Fill-in”¶

If solving multiple $\mathbf{b}$ vectors is the goal, why don’t we just calculate the exact inverse $\mathbf{A}^{-1} = \mathbf{U}^{-1} \mathbf{L}^{-1}$ once and just do $\mathbf{x} = \mathbf{A}^{-1}\mathbf{b}$ ?

The answer is Fill-in.

Fill-in is the phenomenon where, during computation, zeros in a sparse matrix become non-zeros. This destroys the sparsity pattern, requires massive amounts of memory, and drastically slows down the computation.

Let’s look at what happens when you try to take the inverse of a sparse banded matrix.

# Create a sparse tridiagonal matrix
n = 10
A = np.zeros((n, n))
np.fill_diagonal(A, 2)
np.fill_diagonal(A[1:], -1)
np.fill_diagonal(A[:, 1:], -1)

# Calculate its exact inverse
A_inv = np.linalg.inv(A)

# Plot their sparsity patterns side-by-side
fig, (ax1, ax2) = plt.subplots(1, 2, figsize=(10, 4))
ax1.spy(A)
ax1.set_title("Original Banded Matrix A (Sparse!)")
ax2.spy(A_inv)
ax2.set_title("Inverse Matrix A^-1 (Completely Dense!)")
plt.show()

Yikes!¶

The exact inverse of a sparse banded matrix is almost always completely dense. If this matrix was $10,000 \times 10,000$ , storing the inverse would instantly crash your computer’s RAM!

The LU Solution¶

This is why we use LU decomposition!

The LU decomposition of a banded matrix produces $\mathbf{L}$ and $\mathbf{U}$ matrices that are also banded!
LU decomposition completely prevents catastrophic fill-in, preserving our memory and computational speed.

# Compute the PLU decomposition of the tridiagonal matrix A
P, L, U = lu(A)

# Plot the sparsity patterns of L and U
fig, (ax1, ax2, ax3) = plt.subplots(1, 3, figsize=(15, 4))

ax1.spy(A)
ax1.set_title("Original Matrix A (Banded)")

ax2.spy(L)
ax2.set_title("Lower Triangular L (Banded!)")

ax3.spy(U)
ax3.set_title("Upper Triangular U (Banded!)")

plt.show()

Standard Software Solvers¶

In practice, you rarely write your own LU decomposition solver.

Numpy and Scipy both default to highly optimized (P)LU decomposition under the hood when you call their solve functions:

np.linalg.solve(A, b)
scipy.linalg.solve(A, b)

A = np.array([[2, -1, 0],
              [-1, 2, -1],
              [0, -1, 2]])
b = np.array([1, 2, 3])

x_numpy = np.linalg.solve(A, b)
print("Solution using numpy.linalg.solve:\n", x_numpy)

x_scipy = solve(A, b)
print("\nSolution using scipy.linalg.solve:\n", x_scipy)

Sparse Software Solvers¶

If you know your matrix is sparse, you must explicitly use a sparse solver to prevent it from being treated as dense.

The advent of distributed computing motivated incredible algorithms suited for massive sparse systems:

PARADISO (PARallel Direct SOlver)
SuperLU (Supernodal LU)
UMFPACK (Unsymmetric-pattern MultiFrontal method)

Scipy provides scipy.sparse.linalg.spsolve, which routes to SuperLU or UMFPACK automatically! Let’s see the speed difference.

# Generate a 200x200 sparse tridiagonal matrix
n = 200
main_diag = np.full(n, 2)
upper_diag = np.full(n - 1, -1)
lower_diag = np.full(n - 1, -1)
A_sparse = diags([lower_diag, main_diag, upper_diag], offsets=[-1, 0, 1], format='csr')

b = np.random.rand(n)

# Convert to dense format for comparison
A_dense = A_sparse.toarray()

print("Timing the Sparse Solver:")
%timeit spsolve(A_sparse, b)

print("\nTiming the Dense Solver:")
%timeit solve(A_dense, b)