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`? This happens all the time in the real world. 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? 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=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 Core Idea
The `L` matrix is the perfect record of the elimination steps we performed to get `A` into its upper triangular form `U`.
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`.
First, solve `Lc = b` for `c`. (Forward substitution - very fast)
Then, solve `Ux = c` for `x`. (Back substitution - very fast)
This is better because both systems involve triangular matrices, which are solved with lightning-fast substitution.
Step 1: Solve `Lc = b` for `c`
Given `L` and `b`, find `c = [c₁, c₂]`. This is forward substitution.
[1201][c1c2]=[39]⟹c=[33]
Step 2: Solve `Ux = c` for `x`
Now use `U` and our new `c` to find `x = [x₁, x₂]`. This is back substitution.
[2013][x1x2]=[33]⟹x=[11]
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: 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.