MIT 18.06SC Lecture 4: LU Decomposition

Linear Algebra
MIT 18.06
Matrix Factorization
Author

Chao Ma

Published

October 7, 2025

Context

My lecture notes | Exercises notebook

Gilbert Strang’s fourth lecture introduces one of the most important matrix factorizations: LU decomposition, which factors any invertible matrix \(A\) into the product of a Lower triangular matrix and an Upper triangular matrix. This factorization is the foundation of efficient numerical linear algebra.

What is LU Decomposition?

Goal: Factor any invertible matrix \(A\) as the product of: - \(L\) = Lower triangular matrix (with 1’s on diagonal) - \(U\) = Upper triangular matrix (the result of elimination)

\[ A = LU \]

Why is this useful?

  1. Efficient solving: \(Ax = b\) becomes two simpler triangular solves:

    Step 1 - Forward substitution: Solve \(Lc = b\) for \(c\)

    Step 2 - Back substitution: Solve \(Ux = c\) for \(x\)

    How this works:

    Since \(A = LU\), we have \(Ax = LUx = b\). Let \(Ux = c\), then: \[ LUx = Lc = b \]

    Forward substitution (solving \(Lc = b\)):

    Since \(L\) is lower triangular with 1’s on the diagonal, we can solve for \(c\) step by step: \[ \begin{aligned} c_1 &= b_1 \\ c_2 &= b_2 - m_{21}c_1 \\ c_3 &= b_3 - m_{31}c_1 - m_{32}c_2 \\ &\vdots \end{aligned} \]

    Each \(c_i\) depends only on previously computed values, so we solve forward from \(c_1\) to \(c_n\).

    Back substitution (solving \(Ux = c\)):

    Since \(U\) is upper triangular, we solve backward from \(x_n\) to \(x_1\): \[ \begin{aligned} x_n &= \frac{c_n}{u_{nn}} \\ x_{n-1} &= \frac{c_{n-1} - u_{n-1,n}x_n}{u_{n-1,n-1}} \\ &\vdots \end{aligned} \]

    Result: We’ve solved \(Ax = b\) without ever explicitly computing \(A^{-1}\)!

  2. Reusable factorization: When \(A\) is fixed but \(b\) changes, we can reuse \(L\) and \(U\)

    • Factorization: \(O(n^3)\) operations (done once)
    • Each solve: \(O(n^2)\) operations
    • Huge savings for multiple right-hand sides!

  3. Foundation of numerical computing: Used in MATLAB, NumPy, and all scientific computing libraries

How Elimination Creates U

The Elimination Process

Starting with \(A\), we apply elimination matrices \(E_{21}, E_{31}, E_{32}, \ldots\) to get upper triangular \(U\):

\[ E_{32} E_{31} E_{21} A = U \]

Example (3×3 case):

\[ A = \begin{bmatrix} 2 & 1 & 1 \\ 4 & -6 & 0 \\ -2 & 7 & 2 \end{bmatrix} \]

Step 1: Eliminate below first pivot (rows 2 and 3)

\[ E_{21} = \begin{bmatrix} 1 & 0 & 0 \\ -2 & 1 & 0 \\ 0 & 0 & 1 \end{bmatrix}, \quad E_{31} = \begin{bmatrix} 1 & 0 & 0 \\ 0 & 1 & 0 \\ 1 & 0 & 1 \end{bmatrix} \]

Step 2: Eliminate below second pivot (row 3)

\[ E_{32} = \begin{bmatrix} 1 & 0 & 0 \\ 0 & 1 & 0 \\ 0 & -1 & 1 \end{bmatrix} \]

Structure of Elimination Matrices

An elimination matrix \(E_{ij}\) eliminates the entry at position \((i,j)\) by subtracting a multiple of row \(j\) from row \(i\).

General form: \[ E_{ij} = I - m_{ij} \mathbf{e}_i \mathbf{e}_j^T \]

where: - \(m_{ij}\) = multiplier = \(\frac{A_{ij}}{\text{pivot at } (j,j)}\) - \(\mathbf{e}_i\) = \(i\)-th standard basis vector - The \((i,j)\) entry of \(E_{ij}\) is \(-m_{ij}\)

Key properties: 1. Lower triangular (operates below diagonal) 2. Determinant = 1 (doesn’t change volume) 3. Easy to invert: \(E_{ij}^{-1} = I + m_{ij} \mathbf{e}_i \mathbf{e}_j^T\) (just flip the sign!)

Inverting to Get L: The Key Insight

From elimination, we have:

\[ E_{32} E_{31} E_{21} A = U \]

Multiply both sides by the inverses (in reverse order):

\[ A = E_{21}^{-1} E_{31}^{-1} E_{32}^{-1} U = LU \]

where: \[ L = E_{21}^{-1} E_{31}^{-1} E_{32}^{-1} \]

The Beautiful Result

When elimination matrices are multiplied in the right order, their inverses combine beautifully:

\[ L = \begin{bmatrix} 1 & 0 & 0 & \cdots \\ m_{21} & 1 & 0 & \cdots \\ m_{31} & m_{32} & 1 & \cdots \\ \vdots & \vdots & \ddots & \ddots \end{bmatrix} \]

The multipliers \(m_{ij}\) (used during elimination) directly fill in the entries of \(L\) below the diagonal!

No extra computation needed — just save the multipliers as you eliminate.

Computational Complexity

Operation Counts

For an \(n \times n\) matrix:

Step Operations Order
Elimination (find U) \(\frac{n^3}{3} + O(n^2)\) \(O(n^3)\)
Forward substitution \((Lc = b)\) \(\frac{n^2}{2}\) \(O(n^2)\)
Back substitution \((Ux = c)\) \(\frac{n^2}{2}\) \(O(n^2)\)

Why \(\frac{n^3}{3}\)?

At step \(k\), we update an \((n-k) \times (n-k)\) submatrix:

\[ \text{Total operations} = \sum_{k=1}^{n-1} (n-k)^2 \approx \int_0^n x^2 \, dx = \frac{n^3}{3} \]

When is LU Worth It?

Single solve: \(Ax = b\) costs \(O(n^3)\) either way

Multiple solves: If solving \(Ax = b_1, Ax = b_2, \ldots, Ax = b_m\): - Without LU: \(m \times O(n^3)\) - With LU: \(O(n^3)\) (once) + \(m \times O(n^2)\) ✅

Huge savings when \(A\) is fixed but \(b\) changes!

Hands-On Exercises

Let’s practice LU decomposition with concrete examples.

Show code 
import numpy as np

print("✓ Libraries imported successfully")
✓ Libraries imported successfully

Exercise 1: Manual LU Decomposition (2×2)

Compute the LU decomposition of \(A = \begin{bmatrix} 2 & 3 \\ 4 & 7 \end{bmatrix}\) by hand.

Steps: 1. Perform elimination to get \(U\) 2. Record the multiplier \(m_{21}\) to build \(L\) 3. Verify \(A = LU\)

Show code 
from IPython.display import display, Markdown, Latex

# Original matrix
A = np.array([[2, 3],
              [4, 7]])

display(Markdown("**Original matrix A:**"))
display(Latex(r"$$A = \begin{bmatrix} 2 & 3 \\ 4 & 7 \end{bmatrix}$$"))

# Compute multiplier m21
m21 = 4/2  # row2[0] / row1[0]

display(Markdown("<hr />"))
display(Markdown("**Step 1: Compute multiplier**"))
display(Latex(f"$$m_{{21}} = \\frac{{4}}{{2}} = {m21}$$"))

# Build L matrix
L = np.array([[1, 0],
              [m21, 1]])

display(Markdown("<hr />"))
display(Markdown("**Step 2: Build L matrix**"))
display(Latex(r"$$L = \begin{bmatrix} 1 & 0 \\ m_{21} & 1 \end{bmatrix} = \begin{bmatrix} 1 & 0 \\ 2 & 1 \end{bmatrix}$$"))

# Build U matrix (result after elimination)
# After: row2 = row2 - m21*row1
# [2, 3]        [2, 3]
# [4, 7]  -->   [0, 1]  (because 7 - 2*3 = 1)
U = np.array([[2, 3],
              [0, 1]])

display(Markdown("<hr />"))
display(Markdown("**Step 3: Build U matrix (after elimination)**"))
display(Markdown("Row 2 → Row 2 - 2 × Row 1"))
display(Latex(r"$$U = \begin{bmatrix} 2 & 3 \\ 0 & 1 \end{bmatrix}$$"))

display(Markdown("<hr />"))
display(Markdown("**Verification: $A = LU$**"))
display(Latex(r"$$LU = \begin{bmatrix} 1 & 0 \\ 2 & 1 \end{bmatrix} \begin{bmatrix} 2 & 3 \\ 0 & 1 \end{bmatrix} = \begin{bmatrix} 2 & 3 \\ 4 & 7 \end{bmatrix} = A \quad \checkmark$$"))

Original matrix A:

\[A = \begin{bmatrix} 2 & 3 \\ 4 & 7 \end{bmatrix}\]

Step 1: Compute multiplier

\[m_{21} = \frac{4}{2} = 2.0\]

Step 2: Build L matrix

\[L = \begin{bmatrix} 1 & 0 \\ m_{21} & 1 \end{bmatrix} = \begin{bmatrix} 1 & 0 \\ 2 & 1 \end{bmatrix}\]

Step 3: Build U matrix (after elimination)

Row 2 → Row 2 - 2 × Row 1

\[U = \begin{bmatrix} 2 & 3 \\ 0 & 1 \end{bmatrix}\]

Verification: \(A = LU\)

\[LU = \begin{bmatrix} 1 & 0 \\ 2 & 1 \end{bmatrix} \begin{bmatrix} 2 & 3 \\ 0 & 1 \end{bmatrix} = \begin{bmatrix} 2 & 3 \\ 4 & 7 \end{bmatrix} = A \quad \checkmark\]

Key observation: The multiplier \(m_{21} = 2\) goes directly into position \((2,1)\) of \(L\)!

Exercise 2: LU Decomposition (3×3)

Perform LU decomposition on:

\[ A = \begin{bmatrix} 2 & 1 & 1 \\ 4 & -6 & 0 \\ -2 & 7 & 2 \end{bmatrix} \]

Goal: Find \(L\) and \(U\) such that \(A = LU\)

Show code 
from IPython.display import display, Markdown, Latex

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

display(Markdown("**Original matrix A:**"))
display(Latex(r"$$A = \begin{bmatrix} 2 & 1 & 1 \\ 4 & -6 & 0 \\ -2 & 7 & 2 \end{bmatrix}$$"))

display(Markdown("<hr />"))
display(Markdown("### Step 1: Eliminate column 1"))

# Calculate multipliers for column 1
m21 = A[1, 0] / A[0, 0]  # 4/2 = 2
m31 = A[2, 0] / A[0, 0]  # -2/2 = -1

display(Markdown("**Multipliers:**"))
display(Latex(f"$$m_{{21}} = \\frac{{4}}{{2}} = {m21}, \\quad m_{{31}} = \\frac{{-2}}{{2}} = {m31}$$"))

# Create A1 after first elimination
A1 = A.copy().astype(float)
A1[1] = A1[1] - m21 * A1[0]  # row2 - 2*row1
A1[2] = A1[2] - m31 * A1[0]  # row3 - (-1)*row1

display(Markdown("**After eliminating column 1:**"))
display(Markdown("- Row 2 → Row 2 - 2 × Row 1"))
display(Markdown("- Row 3 → Row 3 - (-1) × Row 1"))
display(Latex(r"$$A^{(1)} = \begin{bmatrix} 2 & 1 & 1 \\ 0 & -8 & -2 \\ 0 & 8 & 3 \end{bmatrix}$$"))

display(Markdown("<hr />"))
display(Markdown("### Step 2: Eliminate column 2"))

# Calculate multiplier for column 2
m32 = A1[2, 1] / A1[1, 1]  # 8/(-8) = -1

display(Markdown("**Multiplier:**"))
display(Latex(f"$$m_{{32}} = \\frac{{8}}{{-8}} = {m32}$$"))

# Create U (final upper triangular)
U = A1.copy()
U[2] = U[2] - m32 * U[1]  # row3 - (-1)*row2

display(Markdown("**After eliminating column 2:**"))
display(Markdown("- Row 3 → Row 3 - (-1) × Row 2"))
display(Latex(r"$$U = \begin{bmatrix} 2 & 1 & 1 \\ 0 & -8 & -2 \\ 0 & 0 & 1 \end{bmatrix}$$"))

display(Markdown("<hr />"))
display(Markdown("### Build L from multipliers"))

# Build L from multipliers
L = np.array([[1, 0, 0],
              [m21, 1, 0],
              [m31, m32, 1]])

display(Markdown("The multipliers directly fill in $L$:"))
display(Latex(r"$$L = \begin{bmatrix} 1 & 0 & 0 \\ m_{21} & 1 & 0 \\ m_{31} & m_{32} & 1 \end{bmatrix} = \begin{bmatrix} 1 & 0 & 0 \\ 2 & 1 & 0 \\ -1 & -1 & 1 \end{bmatrix}$$"))

display(Markdown("<hr />"))
display(Markdown("### Verification: $A = LU$"))

display(Latex(r"$$LU = \begin{bmatrix} 1 & 0 & 0 \\ 2 & 1 & 0 \\ -1 & -1 & 1 \end{bmatrix} \begin{bmatrix} 2 & 1 & 1 \\ 0 & -8 & -2 \\ 0 & 0 & 1 \end{bmatrix} = \begin{bmatrix} 2 & 1 & 1 \\ 4 & -6 & 0 \\ -2 & 7 & 2 \end{bmatrix} = A \quad \checkmark$$"))

Original matrix A:

\[A = \begin{bmatrix} 2 & 1 & 1 \\ 4 & -6 & 0 \\ -2 & 7 & 2 \end{bmatrix}\]

Step 1: Eliminate column 1

Multipliers:

\[m_{21} = \frac{4}{2} = 2.0, \quad m_{31} = \frac{-2}{2} = -1.0\]

After eliminating column 1:

  • Row 2 → Row 2 - 2 × Row 1
  • Row 3 → Row 3 - (-1) × Row 1

\[A^{(1)} = \begin{bmatrix} 2 & 1 & 1 \\ 0 & -8 & -2 \\ 0 & 8 & 3 \end{bmatrix}\]

Step 2: Eliminate column 2

Multiplier:

\[m_{32} = \frac{8}{-8} = -1.0\]

After eliminating column 2:

  • Row 3 → Row 3 - (-1) × Row 2

\[U = \begin{bmatrix} 2 & 1 & 1 \\ 0 & -8 & -2 \\ 0 & 0 & 1 \end{bmatrix}\]

Build L from multipliers

The multipliers directly fill in \(L\):

\[L = \begin{bmatrix} 1 & 0 & 0 \\ m_{21} & 1 & 0 \\ m_{31} & m_{32} & 1 \end{bmatrix} = \begin{bmatrix} 1 & 0 & 0 \\ 2 & 1 & 0 \\ -1 & -1 & 1 \end{bmatrix}\]

Verification: \(A = LU\)

\[LU = \begin{bmatrix} 1 & 0 & 0 \\ 2 & 1 & 0 \\ -1 & -1 & 1 \end{bmatrix} \begin{bmatrix} 2 & 1 & 1 \\ 0 & -8 & -2 \\ 0 & 0 & 1 \end{bmatrix} = \begin{bmatrix} 2 & 1 & 1 \\ 4 & -6 & 0 \\ -2 & 7 & 2 \end{bmatrix} = A \quad \checkmark\]

Key observation: All three multipliers \((m_{21}, m_{31}, m_{32})\) go directly into their corresponding positions in \(L\):

\[ L = \begin{bmatrix} 1 & 0 & 0 \\ m_{21} & 1 & 0 \\ m_{31} & m_{32} & 1 \end{bmatrix} = \begin{bmatrix} 1 & 0 & 0 \\ 2 & 1 & 0 \\ -1 & -1 & 1 \end{bmatrix} \]

Note: In practice, numerical libraries like SciPy provide scipy.linalg.lu() which computes LU decomposition efficiently and includes automatic row permutation (pivoting) for numerical stability.

Show code 
from scipy.linalg import lu
from IPython.display import display, Markdown, Latex

A = np.array([[2, 1, 1],
              [4, -6, 0],
              [-2, 7, 2]], dtype=float)

# SciPy returns P, L, U where PA = LU (P is permutation matrix)
P, L_scipy, U_scipy = lu(A)

display(Markdown("**SciPy's LU decomposition:**"))
display(Markdown("SciPy returns $P$, $L$, $U$ where $PA = LU$ ($P$ is a permutation matrix)"))

display(Markdown("<hr />"))
display(Markdown("**Permutation matrix P:**"))
display(Latex(r"$$P = \begin{bmatrix} 0 & 1 & 0 \\ 1 & 0 & 0 \\ 0 & 0 & 1 \end{bmatrix}$$"))
display(Markdown("(This swaps rows 1 and 2 for numerical stability)"))

display(Markdown("<hr />"))
display(Markdown("**Lower triangular L:**"))
display(Latex(r"$$L = \begin{bmatrix} 1 & 0 & 0 \\ 0.5 & 1 & 0 \\ -0.5 & 1 & 1 \end{bmatrix}$$"))

display(Markdown("<hr />"))
display(Markdown("**Upper triangular U:**"))
display(Latex(r"$$U = \begin{bmatrix} 4 & -6 & 0 \\ 0 & 4 & 1 \\ 0 & 0 & 1 \end{bmatrix}$$"))

display(Markdown("<hr />"))
display(Markdown("**Verification: $PA = LU$**"))

# Note: If P = I (identity), then our manual L and U should match
if np.allclose(P, np.eye(3)):
    display(Markdown("✓ No row swaps needed! Our manual $L$ and $U$ match SciPy."))
else:
    display(Markdown("⚠ **Row swaps were performed** (pivot strategy for numerical stability)."))
    display(Markdown("SciPy chose the largest pivot to minimize rounding errors."))
    display(Markdown("Our manual decomposition is valid but uses a different pivot order."))

SciPy’s LU decomposition:

SciPy returns \(P\), \(L\), \(U\) where \(PA = LU\) (\(P\) is a permutation matrix)

Permutation matrix P:

\[P = \begin{bmatrix} 0 & 1 & 0 \\ 1 & 0 & 0 \\ 0 & 0 & 1 \end{bmatrix}\]

(This swaps rows 1 and 2 for numerical stability)

Lower triangular L:

\[L = \begin{bmatrix} 1 & 0 & 0 \\ 0.5 & 1 & 0 \\ -0.5 & 1 & 1 \end{bmatrix}\]

Upper triangular U:

\[U = \begin{bmatrix} 4 & -6 & 0 \\ 0 & 4 & 1 \\ 0 & 0 & 1 \end{bmatrix}\]

Verification: \(PA = LU\)

⚠ Row swaps were performed (pivot strategy for numerical stability).

SciPy chose the largest pivot to minimize rounding errors.

Our manual decomposition is valid but uses a different pivot order.