LU & QR Decomposition: Matrix Factorization

Why Decompose Matrices?

Solving Ax = b directly by computing A^{- 1} is expensive (O(n^{3}) and numerically unstable. Matrix decompositions factor A into simpler matrices that make solving systems faster and more stable.

LU Decomposition

Factor into Lower × Upper triangular. Good for repeated solves with same A, different b. Analogy: Gaussian Elimination.

QR Decomposition

Factor into Orthogonal × Upper triangular. Numerically stable, ideal for least squares. Analogy: Gram-Schmidt.

LU Decomposition

Factor a square matrix into a product of Lower and Upper triangular matrices. It is essentially recording the steps of Gaussian Elimination.

A = LU

L = lower triangular (1s on diagonal), U = upper triangular (pivots on diagonal)

Solving with LU

1. Factor once: A = LU (O(n³))
2. For each b, solve Ly = b (forward substitution, O(n²))
3. Then solve Ux = y (back substitution, O(n²))

Partial Pivoting (PA = LU)

In practice, we permute rows to avoid division by small numbers. This gives PA = LU where P is a permutation matrix.

Interactive Simulator

Step through the LU and QR decomposition process. Watch how matrices are factored at each step.

Matrix Decomposition

Factorizing A into simpler components.

LU (Gaussian)

QR (Gram-Schmidt)

Step 1

Setup

Start with A. Our goal is to transform A into Upper Triangular form (U) using row operations, while recording the multipliers in L.

L (Lower)

1

0

0

1

×

U (Upper)

4

3

6

3

Back

Next Step

QR Decomposition

Factor any matrix (even rectangular) into an Orthogonal matrix times an Upper triangular matrix.

A = QR

Q^{T}Q = I (orthonormal columns), R = upper triangular

Why QR is Stable

Orthogonal matrices preserve lengths: \mid\mid Qx\mid\mid = \mid\mid x\mid\mid This means errors do not amplify during computation (condition number of Q is 1).

Gram-Schmidt Process

The algorithm to construct Q is called Gram-Schmidt. It iteratively subtracts projections to force orthogonality.

1

2

3

4

5

1. Define Vectors

Start with two linearly independent vectors.

Initial Vectors

v₁ Angle15°

v₂ Angle60°

Current Operation

Configure vectors...

Back

Next Step

1. u_{1} = a_{1} then q_{1} = u_{1}/\mid\mid u_{1}\mid\mid
2. u_{2} = a_{2} - (a_{2} \cdot q_{1})q_{1} then q_{2} = u_{2}/\mid\mid u_{2}\mid\mid
3. u_{k} = a_{k} - \sum_{j = 1}^{k - 1}(a_{k} \cdot q_{j})q_{j} then q_{k} = u_{k}/\mid\mid u_{k}\mid\mid

Modified Gram-Schmidt is more numerically stable. Modern libraries use Householder reflections (a sequence of mirrors) instead of projections to compute QR.

Case Study: Bulb Characteristic Fitting

The Problem

You have voltage vs brightness data for a bulb. Fit a polynomial: brightness = a₀ + a₁V + a₂V². This is an overdetermined system (more data points than unknowns).

Using QR for Least Squares

1. Build Vandermonde matrix: A_{ij} = V_{i}^{j}
2. Compute A = QR
3. Solve R\hat{a} = Q^{T}b by back substitution (since R is triangular, this is easy)

Why Not Normal Equations?

For high degree polynomials, A^{T}A becomes ill-conditioned. QR avoids forming this product, preserving numerical stability.

LU vs QR vs SVD

ML Applications

Cholesky (LLᵀ)

For symmetric positive definite matrices (like covariance matrices), Cholesky is 2x faster than LU. Used in Gaussian Processes for sampling.

Eigenvalue Algorithms

The QR Algorithm for eigenvalues repeatedly computes QR decompositions. This is the standard way numpy.linalg.eig works for non-symmetric matrices.

Backpropagation Efficiency

Computing gradients through matrix inverses involves solving linear systems, which typically uses pre-computed LU factors. Frameworks like JAX exploit this for efficiency.

Randomized Linear Algebra

Large scale ML uses randomized QR (sampling columns) to get approximate decompositions. This is much faster than full SVD for massive recommendation system matrices.