LU Decomposition

The Efficiency Expert

So far, we have treated Gaussian Elimination as a process—a series of steps we perform to solve a single system `Ax = b`.

But what if you have to solve multiple systems with the same coefficient matrix `A`, but different result vectors `b`?

  • Ax=b1Ax = b_1
  • Ax=b2Ax = b_2
  • Ax=b3Ax = b_3

This happens all the time in the real world, for example, when simulating a structure under different loads. Running the entire, slow Gaussian Elimination process for each `b` would be incredibly inefficient. The elimination steps only depend on `A`, so we'd be wastefully re-doing the exact same row operations every single time.

The central question is: Can we save the elimination steps themselves? Can we capture the entire forward elimination process in a matrix?

The answer is yes. This is the motivation for the LU Decomposition.

What is LU Decomposition?
LU Decomposition is a matrix factorization. It is a way of breaking down a single, complex matrix `A` into the product of two simpler, more useful matrices:
A=LUA = LU
  • `L` is a Lower triangular matrix. It has 1s on the diagonal and all zeros above the diagonal.
  • `U` is an Upper triangular matrix. It is the exact same row echelon form that we get from running Gaussian Elimination on `A`.

The magic is this: The `L` matrix is the perfect record of the elimination steps we performed to get `A` into its upper triangular form `U`.

The `L` Matrix: Storing the Multipliers
Remember our elimination operations, like `Row 2 → Row 2 - 2 * Row 1`? The number `2` here is the multiplier we used. The `L` matrix is constructed by simply placing these multipliers (with their sign flipped) into the corresponding positions below the diagonal.

Example:

Let's take a simple matrix `A` and find its LU decomposition.

A=[2145]A = \begin{bmatrix} 2 & 1 \\ 4 & 5 \end{bmatrix}
Step 1: Perform forward elimination to get `U`.
  • We need to eliminate the `4` in the second row.
  • The operation is `R₂ → R₂ - 2*R₁`. The multiplier is 2.
  • [4,5]−2⋅[2,1]=[0,3][4, 5] - 2 \cdot [2, 1] = [0, 3]

So, our Upper triangular matrix `U` is:

U=[2103]U = \begin{bmatrix} 2 & 1 \\ 0 & 3 \end{bmatrix}
Step 2: Construct `L` from the multipliers.
  • `L` will be a lower triangular matrix with 1s on the diagonal.
  • The multiplier we used was 2. It was used to operate on Row 2 using Row 1. So, we place this `2` in the Row 2, Column 1 position of `L`.
L=[1021]L = \begin{bmatrix} 1 & 0 \\ 2 & 1 \end{bmatrix}

You can multiply them to verify: `A = LU`

The Payoff: Solving `Ax = b` with Super Speed
We substitute `LU` for `A`: `LUx = b`. Now, we cleverly split this one hard problem into two easy problems. Let's define an intermediate vector `c` such that `Ux = c`.
  1. First, solve `Lc = b` for `c`.
  2. Then, solve `Ux = c` for `x`.

This is better because both systems involve triangular matrices, which are solved with lightning-fast substitution.

Let's see it in action.

Let `A` be our matrix from before, and let `b = [3, 9]`. We want to solve `Ax = b`.

Step 1: Solve `Lc = b` for `c = [c₁, c₂]`.
[1021][c1c2]=[39]\begin{bmatrix} 1 & 0 \\ 2 & 1 \end{bmatrix} \begin{bmatrix} c_1 \\ c_2 \end{bmatrix} = \begin{bmatrix} 3 \\ 9 \end{bmatrix}

From the first row: `1*c₁ + 0*c₂ = 3` → `c₁ = 3`.
From the second row: `2*c₁ + 1*c₂ = 9` → `2(3) + c₂ = 9` → `c₂ = 3`.

So, `c = [3, 3]`. That was fast.

Step 2: Solve `Ux = c` for `x = [x₁, x₂]`.
[2103][x1x2]=[33]\begin{bmatrix} 2 & 1 \\ 0 & 3 \end{bmatrix} \begin{bmatrix} x_1 \\ x_2 \end{bmatrix} = \begin{bmatrix} 3 \\ 3 \end{bmatrix}

From the second row: `3*x₂ = 3` → `x₂ = 1`.
From the first row: `2*x₁ + 1*x₂ = 3` → `2*x₁ + 1 = 3` → `x₁ = 1`.

The final solution is `x = [1, 1]`.

The upfront work of finding `L` and `U` allows us to solve for any new `b` vector with blistering speed. This is the method computers actually use to solve large systems of equations.

Pitfalls and Properties: The Fine Print

Pitfall #1: The Need for Row Swaps

What happens if we have a zero in a pivot position? The standard `A = LU` decomposition will fail. The solution is a slightly modified form called the **PA = LU Decomposition**. `P` is a **Permutation Matrix**—a special type of orthogonal matrix that simply records the row swaps.

Properties and When It Can Be Used

  • **Existence:** An LU decomposition for a square matrix `A` exists if and only if the matrix can be brought to row echelon form without any row swaps.
  • **Uniqueness:** If it exists, and `A` is invertible, the `L` and `U` matrices are unique.
  • **Determinant:** `det(A) = det(L)det(U) = 1 * det(U)` (product of diagonals of `U`). This is often the fastest way to compute a determinant.
  • **Primary Use Case:** The ideal scenario for using LU decomposition is when you need to solve `Ax=b` for a fixed `A` and many different `b` vectors.

Up Next: We have mastered the art of solving `Ax=b`. It's time to step back and analyze the deeper structure of the matrix `A` itself. We will formally introduce the Four Fundamental Subspaces.