Matrix Factorizations: LU Decomposition, QR Factorization and SVD

Matrix factorizations write a matrix as a product of simpler ones: the LU decomposition A = LU (lower times upper triangular) records Gaussian elimination, the QR factorization A = QR (orthonormal columns times upper triangular) records Gram–Schmidt orthogonalization, and the singular value decomposition A = UΣV^T exhibits the stretching factors of any matrix. In brief: the matrix [[2, 1, 1], [4, 3, 3], [8, 7, 9]] factors as L = [[1, 0, 0], [2, 1, 0], [4, 3, 1]] times U = [[2, 1, 1], [0, 1, 1], [0, 0, 2]], and its determinant is the product of the diagonal of U: 2 · 1 · 2 = 4.

Factorizations are the working tools of numerical linear algebra: they solve systems, least-squares problems and eigenvalue problems efficiently and stably. This guide explains the three most important ones, shows how to compute each by hand for small matrices, and what they are used for, with worked examples and a table.

What is the LU decomposition?

Gaussian elimination turns A into an upper triangular matrix U by row operations of the type R_i → R_i − m_ijR_j. Collecting the multipliers m_ij below the diagonal of a unit lower triangular matrix L (ones on the diagonal) gives A = LU. If row swaps are needed, one gets PA = LU with a permutation matrix P. Then Ax = b is solved by forward substitution Ly = b followed by back substitution Ux = y, both cheap.

Example 1: computing LU

A = [[2, 1, 1], [4, 3, 3], [8, 7, 9]]. The multipliers: R₂ − 2R₁ = (0, 1, 1), R₃ − 4R₁ = (0, 3, 5), R₃ − 3R₂ = (0, 0, 2). So U = [[2, 1, 1], [0, 1, 1], [0, 0, 2]] and L = [[1, 0, 0], [2, 1, 0], [4, 3, 1]]. Check the third row of LU: 4(2, 1, 1) + 3(0, 1, 1) + (0, 0, 2) = (8, 7, 9) ✓.

Example 2: solving a system

Solve Ax = b with b = (4, 10, 24). Forward: y₁ = 4, y₂ = 10 − 2 · 4 = 2, y₃ = 24 − 4 · 4 − 3 · 2 = 2. Back: 2x₃ = 2 gives x₃ = 1; x₂ + x₃ = 2 gives x₂ = 1; 2x₁ + x₂ + x₃ = 4 gives x₁ = 1. The solution is (1, 1, 1). The factorization is reused for other right-hand sides at a small cost.

Example 3: a matrix that needs a swap

[[0, 1], [1, 1]] has a zero in the pivot position, so no LU exists without a swap; PA = LU with P = [[0, 1], [1, 0]] gives PA = [[1, 1], [0, 1]] = I · U.

What is the QR factorization?

A matrix A with independent columns can be written A = QR, where Q has orthonormal columns (Q^TQ = I) and R is upper triangular with positive diagonal. The columns of Q come from the Gram–Schmidt process applied to the columns of A, and R holds the coefficients r_ij = q_i · a_j. Numerically, Householder reflections or Givens rotations are preferred, because Gram–Schmidt can lose orthogonality.

Example 4: QR of a 2 × 2 matrix

A = [[1, 1], [1, 0]], with the columns a₁ = (1, 1) and a₂ = (1, 0). q₁ = a₁/|a₁| = (1, 1)/√2, r₁₁ = √2. r₁₂ = q₁ · a₂ = 1/√2. The remainder a₂ − r₁₂q₁ = (1, 0) − (1/2, 1/2) = (1/2, −1/2), of length 1/√2 = r₂₂, so q₂ = (1, −1)/√2. Thus Q = (1/√2) [[1, 1], [1, −1]] and R = [[√2, 1/√2], [0, 1/√2]]. Check: the second column of QR is (1/√2)q₁ + (1/√2)q₂ = (1/2)(1, 1) + (1/2)(1, −1) = (1, 0) ✓.

What is the singular value decomposition?

Every m × n matrix A has a decomposition A = UΣV^T, where U (m × m) and V (n × n) are orthogonal and Σ is diagonal (m × n) with the singular values σ₁ ≥ σ₂ ≥ … ≥ 0. The singular values are the square roots of the eigenvalues of A^TA, the columns of V are its orthonormal eigenvectors, and the columns of U follow from u_i = Av_i/σ_i. The number of nonzero singular values is the rank. Geometrically, A maps the unit sphere to an ellipsoid whose semi-axes are the σ_i.

Example 5: a rank-1 matrix

A = [[1, 2], [2, 4]]: A^TA = [[5, 10], [10, 20]] has the eigenvalues 25 and 0, so the singular values are 5 and 0 and the rank is 1. With v₁ = (1, 2)/√5 and u₁ = Av₁/5 = (1, 2)/√5: A = 5 u₁v₁^T = 5 · (1/5) [[1, 2], [2, 4]] ✓. The sum of squares of all entries, 1 + 4 + 4 + 16 = 25 = σ₁² + σ₂².

Example 6: singular values of a shear

A = [[1, 1], [0, 1]]: A^TA = [[1, 1], [1, 2]] has the characteristic polynomial λ² − 3λ + 1, so λ = (3 ± √5)/2 ≈ 2.618 and 0.382. The singular values are σ₁ = (1 + √5)/2 ≈ 1.618 and σ₂ = (√5 − 1)/2 ≈ 0.618 (since 2.618 = 1.618²). Their product 1 equals |det A| ✓, and the ratio σ₁/σ₂ ≈ 2.618 is the condition number.

What are factorizations used for?

• LU: solving many systems with the same matrix, computing determinants (the product of the diagonal of U, up to the sign of the permutation) and inverses.
• QR: least-squares solutions of overdetermined systems, since minimizing |Ax − b| reduces to the triangular system Rx = Q^Tb; also the QR algorithm for eigenvalues.
• SVD: the best rank-k approximation of a matrix (keep the k largest singular values; Eckart–Young theorem), image compression, principal component analysis, the pseudoinverse and measuring numerical rank.

Example 7: a least-squares fit

Fit the line y = c₀ + c₁t to the points (0, 1), (1, 2), (2, 4). The system has the matrix A = [[1, 0], [1, 1], [1, 2]] and b = (1, 2, 4). Solving the normal equations A^TAx = A^Tb: [[3, 3], [3, 5]] x = (7, 10). The determinant is 6, so x = (1/6) [[5, −3], [−3, 3]] (7, 10) = (1/6)(35 − 30, −21 + 30) = (5/6, 3/2). The line is y = 5/6 + 1.5t. QR would solve the same problem without forming A^TA, which is better conditioned numerically.

Table of factorizations

Where do they come from?

The elimination underlying LU is old, but its interpretation as a matrix factorization was formulated by Alan Turing in 1948 in his study of rounding errors. The Gram–Schmidt process is named after Jørgen Gram and Erhard Schmidt (1907); stable QR computations with reflections were introduced by Alston Householder in 1958, and the QR algorithm for eigenvalues by John Francis and Vera Kublanovskaya around 1961. The SVD was discovered for real square matrices by Eugenio Beltrami (1873) and Camille Jordan (1874), and an effective algorithm was given by Gene Golub and William Kahan in 1965.

Which mistakes are common?

Putting the multipliers with the wrong sign into L

If the operation is R_i − mR_j, then L contains +m, since L undoes the elimination.

Ignoring row swaps

A zero pivot requires a permutation. Without it LU may not exist, and small pivots cause numerical instability.

Mixing up Q and R roles

Q has orthonormal columns and is generally not triangular; R is the triangular factor.

Confusing singular values with eigenvalues

Singular values are nonnegative and defined for every matrix; they coincide with the absolute values of the eigenvalues only for symmetric (or normal) matrices.