Matrix Factorizations and Low-Rank Approximation

Introduction and Context

Matrix factorizations rewrite a matrix as a product of simpler matrices. This is more than an algebraic convenience: it exposes structure, improves numerical stability, reduces computational cost, and supports major AI tasks such as least-squares regression, principal component analysis (PCA), recommender systems, image compression, and latent-factor modelling.

For a matrix A, different factorizations answer different questions. LU factorization is useful for solving square linear systems. QR factorization is especially valuable for least-squares problems. Cholesky factorization exploits symmetry and positive definiteness. Eigendecomposition describes invariant directions of a square matrix. Singular value decomposition (SVD) applies to every real or complex matrix and provides the most informative general-purpose description of its rank, geometry, and approximation quality.

This chapter develops these methods and explains how they should be selected in practical AI systems.

Learning Objectives

By the end of this chapter, you should be able to:

1. Explain the purposes and assumptions of LU, QR, Cholesky, eigendecomposition, and SVD.

2. Use these factorizations to solve systems and least-squares problems.

3. Construct and interpret Moore–Penrose pseudoinverses and compact SVDs.

4. Derive and apply low-rank approximations using the Eckart–Young theorem.

5. Evaluate numerical, computational, and modelling trade-offs in AI applications.

1. Why Factorize a Matrix?

Suppose A is an m × n matrix. Directly computing an inverse is usually unnecessary and often numerically inferior. If a system Ax = b must be solved repeatedly for different vectors b, factorizing A once is generally more efficient. For example, if A = LU, then solving Ax = b becomes two triangular solves:

Ly = b, Ux = y.

Triangular systems can be solved efficiently by forward and backward substitution.

Factorizations also reveal structure:

• Rank indicates the number of independent directions in the data.
• Eigenvectors describe directions preserved by a square transformation.
• Singular vectors describe input and output directions of any linear map.
• Small singular values identify weak or noisy directions.

A factorization is therefore both a computational tool and an interpretive model.

2. LU Factorization

For a square matrix A, an LU factorization has the form

A = LU,

where L is lower triangular, usually with ones on its diagonal, and U is upper triangular. In practical computation, row exchanges are often required:

PA = LU,

where P is a permutation matrix. Pivoting is important because it improves numerical stability and prevents division by very small pivots.

Worked example: LU without pivoting

Let

A = [[2, 1], [4, 3]].

Eliminating the entry below the first pivot uses multiplier 4/2 = 2. Thus

L = [[1, 0], [2, 1]], U = [[2, 1], [0, 1]],

and LU = A.

To solve Ax = [5, 11]^T, first solve Ly = b:

y1 = 5,

2y1 + y2 = 11, so y2 = 1.

Then solve Ux = y:

x2 = 1,

2x1 + x2 = 5, so x1 = 2.

Therefore x = [2, 1]^T.

LU is appropriate for general square systems, but standard LU is not directly designed for rectangular least-squares problems. Pivoting is normally used in software even when a hand calculation does not require it.

3. QR Factorization

A QR factorization expresses a matrix A as

A = QR,

where Q has orthonormal columns and R is upper triangular. For a full-column-rank m × n matrix with m ≥ n, Q is m × n with Q^TQ = I_n, and R is n × n. A full QR may instead use a square orthogonal Q.

Orthogonality is valuable because it preserves Euclidean lengths:

||Qz||_2 = ||z||_2.

It also avoids the numerical problems associated with forming A^TA. For the least-squares problem

min_x ||Ax - b||_2,

with A having full column rank, write A = QR. Then

||Ax-b||_2 = ||QRx-b||_2.

Multiplying by Q^T and separating the component of b in the column space of Q gives

Rx = Q^Tb.

This triangular system yields the least-squares solution. The residual is orthogonal to the columns of A.

Computing QR

Gram–Schmidt constructs orthogonal vectors but classical Gram–Schmidt can lose accuracy when columns are nearly dependent. Modified Gram–Schmidt is better, while Householder reflections are generally preferred in reliable numerical libraries. Givens rotations are useful when updating a factorization or exploiting sparse structure.

Example: least squares through QR

Suppose A has full column rank and A = QR. The least-squares estimate is not obtained by simply solving Ax = b when the system is overdetermined. Instead solve

Rx = Q^Tb.

This approach is usually more stable than solving the normal equations

A^TAx = A^Tb,

because forming A^TA squares the condition number approximately: κ(A^TA) = κ(A)^2 for the 2-norm when A has full column rank.

4. Cholesky Factorization

If A is real, symmetric, and positive definite, then it has a unique factorization

A = LL^T,

where L is lower triangular with positive diagonal entries. Positive definiteness means x^TAx > 0 for every nonzero vector x.

Cholesky is approximately twice as efficient as general LU for a comparable dense symmetric matrix and requires less storage. To solve Ax = b:

1. Solve Ly = b.

2. Solve L^Tx = y.

Cholesky applies to covariance matrices, kernel matrices, and Hessian approximations when the required assumptions hold. It cannot be safely applied to an arbitrary symmetric matrix: symmetry alone does not imply positive definiteness. If a covariance estimate is singular or numerically indefinite, a regularized matrix such as A + λI, with λ > 0, may be used, provided this modelling choice is justified.

5. Eigendecomposition

For a square matrix A, an eigendecomposition seeks

A = VΛV^{-1},

where Λ is diagonal and its entries λ_i are eigenvalues satisfying

Av_i = λ_i v_i.

This form requires enough linearly independent eigenvectors. Some matrices are defective and cannot be diagonalized. For a real symmetric matrix,

A = QΛQ^T,

where Q is orthogonal and all eigenvalues are real. This is the spectral theorem.

Eigendecomposition is central to PCA. If X is a centered data matrix with observations in rows, the sample covariance matrix is

C = (1/(m-1))X^TX.

The eigenvectors of C give principal directions, and the corresponding eigenvalues give variances along those directions. The eigenvectors are ordered by decreasing eigenvalue.

Eigendecomposition is also used for graph-based learning, dynamical systems, and analysis of quadratic forms. However, it is restricted to square matrices and may be unstable or difficult to interpret when eigenvalues are repeated or eigenvectors are poorly conditioned.

6. Singular Value Decomposition

Every real m × n matrix A has an SVD:

A = UΣV^T,

where U is m × m orthogonal, V is n × n orthogonal, and Σ is m × n diagonal in the rectangular sense. The diagonal entries satisfy

σ1 ≥ σ2 ≥ ... ≥ σr > 0,

where r = rank(A), and are called singular values.

The defining relationships are

Av_i = σ_i u_i,

A^Tu_i = σ_i v_i.

The columns v_i are right singular vectors, and u_i are left singular vectors. The singular values satisfy

σ_i^2 = eigenvalues of A^TA,

and also, for nonzero values, of AA^T.

Geometrically, a unit sphere is transformed by A into an ellipsoid. The right singular vectors identify input directions, the singular values give the corresponding stretching factors, and the left singular vectors identify output directions.

Important norms are

||A||_2 = σ1,

and

||A||_F = (Σ_i σ_i^2)^(1/2),

where the Frobenius norm sums the squares of all entries. The rank is the number of nonzero singular values.

Compact SVD

If A has rank r, its compact SVD is

A = U_rΣ_rV_r^T,

where U_r contains only the first r left singular vectors, V_r contains only the first r right singular vectors, and Σ_r is r × r diagonal with positive entries. This representation removes zero singular-value directions and is often substantially smaller than the full SVD.

7. Pseudoinverses

A matrix inverse exists only for a square, nonsingular matrix. The Moore–Penrose pseudoinverse, written A^+, generalizes inversion to rectangular and rank-deficient matrices. If

A = UΣV^T,

then

A^+ = VΣ^+U^T,

where Σ^+ is formed by transposing Σ and replacing each nonzero σ_i by 1/σ_i. Singular values that are zero remain zero in the pseudoinverse.

For a consistent full-column-rank system, A^+b gives the unique least-squares solution. In general, x = A^+b is the minimum-Euclidean-norm least-squares solution. If Ax = b is consistent, it is the minimum-norm exact solution among all solutions.

For full column rank, the equivalent expression is

A^+ = (A^TA)^{-1}A^T,

provided A^TA is invertible. For full row rank,

A^+ = A^T(AA^T)^{-1}.

These formulas are mathematically valid under their stated rank assumptions, but SVD-based computation is usually more reliable for nearly rank-deficient matrices.

8. Low-Rank Approximation

Many datasets are approximately low rank. An image may contain repeated visual patterns, a user–item matrix may be generated by a small number of preferences, and word-document data may be represented by a limited number of semantic factors.

Using the SVD,

A = Σ_{i=1}^r σ_i u_i v_i^T.

The rank-k truncated SVD is

A_k = Σ_{i=1}^k σ_i u_i v_i^T,

where k < r. It stores the strongest k patterns and discards weaker directions.

Eckart–Young theorem

For any matrix A and any rank-k matrix B, where 0 ≤ k < rank(A), the truncated SVD A_k solves

min_{rank(B)≤k} ||A-B||_2

and

min_{rank(B)≤k} ||A-B||_F.

The minimum errors are, respectively,

||A-A_k||_2 = σ_{k+1},

and

||A-A_k||_F = (Σ_{i=k+1}^r σ_i^2)^(1/2).

Thus, among all rank-k approximations, truncated SVD is optimal for both spectral and Frobenius error. If σ_k > σ_{k+1}, the best rank-k approximation is unique; otherwise, the approximation need not be unique, although the minimum error remains fixed.

Choosing k

There is no universally correct rank. Choices may be based on:

• a desired reconstruction error;
• a cumulative explained-variance threshold;
• validation performance in a downstream model;
• storage and latency constraints;
• interpretability and stability.

For PCA, the explained-variance ratio of the first k components is

(Σ_{i=1}^k σ_i^2)/(Σ_{i=1}^r σ_i^2),

when X is centered and singular values are computed from X. This ratio should not be confused with predictive accuracy.

9. Worked Low-Rank Example

Consider a matrix whose singular values are 10, 3, and 1. A rank-1 approximation retains only σ1 = 10.

Its spectral-norm error is σ2 = 3. Its Frobenius-norm error is

(3^2 + 1^2)^(1/2) = √10 ≈ 3.16.

A rank-2 approximation has spectral error σ3 = 1 and Frobenius error 1. The rank-2 model is therefore much more accurate, but it requires an additional singular direction. This illustrates the accuracy–compression trade-off.

For an m × n matrix, a rank-k SVD representation requires approximately k(m+n+1) values, compared with mn values for the original dense matrix. Compression is beneficial when k(m+n+1) is substantially smaller than mn, although metadata, numerical precision, and computational overhead must also be considered.

10. Application Case: Recommender-System Latent Factors

Background

An online retailer records a user–item rating matrix R. Most entries are missing because each user rates only a small fraction of products. Directly treating missing entries as zero confuses “not observed” with “disliked.”

Problem

The retailer wants to predict unseen preferences while reducing storage and computation. It also wants factors that capture broad patterns such as price sensitivity, genre preference, or technical sophistication.

Analysis

A low-rank assumption proposes that ratings are generated by a limited number of latent factors. A basic model is

R_ij ≈ p_i^Tq_j,

where p_i ∈ R^k represents user i and q_j ∈ R^k represents item j. An SVD of a suitably completed or centred matrix provides an initial latent representation. In production, matrix factorization is usually fitted directly to observed entries by minimizing a regularized objective such as

Σ_(i,j observed) (R_ij - p_i^Tq_j)^2 + λ(Σ_i||p_i||_2^2 + Σ_j||q_j||_2^2),

where λ > 0 controls regularization.

Approach and solution

The team centres ratings by user, selects k using validation data, and fits factors only on observed ratings. Regularization limits overfitting. Predictions are evaluated using a held-out test set and ranking measures, not reconstruction error alone.

Outcome

The latent-factor model stores two factor matrices instead of the full rating table, improves ranking quality over a popularity-only baseline, and supports recommendations. However, very small singular values or factors should not automatically be interpreted as meaningful preferences; they may represent noise, sampling effects, or demographic bias.

Learning points

1. Low rank is a modelling assumption, not a guarantee.

2. Missing-data mechanisms must be treated explicitly.

3. Validation should reflect the deployment objective.

4. Regularization and fairness checks are essential.

11. Numerical and Practical Considerations

Conditioning

The condition number in the 2-norm for a full-rank matrix is

κ_2(A) = σ1/σr,

where r = min(m,n) for a full-rank rectangular matrix. A large condition number means that small perturbations in data can cause large changes in computed solutions. Pseudoinverse calculations should use a tolerance to treat sufficiently small singular values as zero.

Randomized and truncated methods

For very large matrices, computing a full SVD may be too expensive. Randomized SVD uses random projections to approximate a dominant singular subspace and is effective when only a modest number of leading components is required. Lanczos and Krylov-subspace methods are also useful for sparse matrices. These methods trade exactness for speed and should be assessed against reconstruction and downstream metrics.

Scaling and preprocessing

SVD and PCA are sensitive to scale. Features measured in large units can dominate. Standardization may be appropriate, but it changes the geometry and must be justified by the application. PCA generally requires centering; applying it to uncentered data answers a different question.

Governance and responsible use

Compression can remove minority patterns and amplify dominant groups. In recommender systems, latent factors can reproduce historical discrimination. Singular vectors are not automatically causal explanations. Systems should document preprocessing, rank selection, missing-data treatment, validation populations, and the consequences of discarded information. Sensitive data should be protected during factorization and storage.

12. Selecting a Factorization

Use LU with pivoting for general square linear systems. Use QR for stable least-squares computation, especially when the matrix is rectangular or poorly conditioned. Use Cholesky when symmetry and positive definiteness are established. Use eigendecomposition for suitable square operators, particularly symmetric covariance matrices. Use SVD when the matrix is rectangular, rank-deficient, ill-conditioned, or when low-rank structure and pseudoinverses are required.

A strong professional workflow avoids explicit matrix inversion, checks dimensions and rank assumptions, uses stable library implementations, monitors condition numbers, and validates approximations on the actual downstream task.

Chapter Summary

Matrix factorizations convert difficult matrix operations into structured computations. LU factorization solves general square systems through triangular solves, normally with pivoting. QR factorization uses orthogonality and is a stable foundation for least-squares problems. Cholesky is efficient for symmetric positive-definite matrices. Eigendecomposition represents suitable square matrices through invariant directions and is central to covariance-based PCA.

SVD is the most general factorization: it applies to every matrix, reveals rank and conditioning, and separates input directions, stretching strengths, and output directions. Its compact form removes zero directions, while the Moore–Penrose pseudoinverse provides minimum-norm solutions for rectangular or rank-deficient systems. Truncated SVD gives the optimal rank-k approximation under both spectral and Frobenius norms by the Eckart–Young theorem. These ideas support compression, PCA, recommender systems, and latent-factor models, but practical use requires attention to scale, missing data, numerical stability, validation, fairness, and interpretability.

Key Terms

Factorization: Representation of a matrix as a product of simpler matrices.

LU factorization: Decomposition of a matrix into lower- and upper-triangular factors.

Pivoting: Row or column exchanges used to improve numerical stability.

QR factorization: Decomposition into an orthonormal factor Q and upper-triangular factor R.

Orthogonal matrix: A square matrix Q satisfying Q^TQ = I.

Least squares: Minimization of the Euclidean residual norm ||Ax-b||_2.

Cholesky factorization: A = LL^T for a symmetric positive-definite matrix.

Positive definite: A matrix satisfying x^TAx > 0 for every nonzero x.

Eigenvalue: A scalar λ for which Av = λv for a nonzero vector v.

Eigendecomposition: Representation of a diagonalizable square matrix using eigenvectors and eigenvalues.

Singular value: A nonnegative stretching factor in the SVD.

SVD: Decomposition A = UΣV^T.

Compact SVD: SVD retaining only nonzero singular directions.

Pseudoinverse: The Moore–Penrose generalized inverse A^+.

Rank: Number of linearly independent directions or nonzero singular values.

Condition number: A measure of sensitivity to perturbations; in the 2-norm, σ1/σr.

Low-rank approximation: Representation of a matrix using a specified maximum rank.

Eckart–Young theorem: The theorem establishing truncated SVD as an optimal rank-k approximation.

Latent factor: An unobserved variable used to explain observed data patterns.

Regularization: A penalty or constraint used to reduce overfitting or instability.

Review Questions

Short Answer Questions

1. What are the roles of L and U in LU factorization?

2. Why is QR generally preferred to normal equations for numerically stable least squares?

3. State the assumptions required for Cholesky factorization.

4. What information do singular values provide about a matrix?

5. Define the Moore–Penrose pseudoinverse.

Long Answer Questions

1. Explain LU, QR, Cholesky, eigendecomposition, and SVD, comparing their assumptions and typical applications.

2. Derive the QR-based least-squares solution for a full-column-rank matrix.

3. Explain the geometric interpretation of the SVD and its relationship to the eigendecomposition of A^TA.

4. State and interpret the Eckart–Young theorem under both spectral and Frobenius norms.

5. Discuss how low-rank factorization can be used in recommender systems, including missing data and regularization.

Analytical / Application Questions

1. A matrix has singular values 12, 4, 2, and 0.5. Calculate the spectral- and Frobenius-norm errors of its rank-2 truncated approximation.

2. You must solve a large least-squares problem with a rectangular, nearly dependent design matrix. Which factorization would you choose, and why would forming A^TA be risky?

3. A PCA system retains 95% of variance but performs poorly on a classification task. Explain at least three reasons this can occur and propose appropriate validation steps.

Multiple Choice Questions with Answers

Question 1

Which factorization is most directly suited to stable least-squares solution for a rectangular matrix?

A. QR factorization

B. Cholesky factorization of A

C. Eigendecomposition of a nonsquare matrix

D. LU without pivoting

Correct Answer: A. QR factorization

Question 2

What condition is required for the standard Cholesky factorization A = LL^T?

A. A must be diagonal

B. A must be symmetric positive definite

C. A must be rectangular

D. A must have rank one

Correct Answer: B. A must be symmetric positive definite

Question 3

In the SVD A = UΣV^T, what do the singular values represent geometrically?

A. Eigenvector angles

B. Input dimensions only

C. Stretching factors

D. The number of rows

Correct Answer: C. Stretching factors

Question 4

For a matrix with singular values 10, 3, and 1, what is the spectral-norm error of its best rank-1 approximation?

A. 1

B. 3

C. 10

D. √10

Correct Answer: B. 3

Question 5

Which expression gives the SVD-based Moore–Penrose pseudoinverse?

A. UΣ^+V^T

B. VΣ^+U^T

C. VΣU^T

D. UΣV^T

Correct Answer: B. VΣ^+U^T

Question 6

For full-column-rank A, the QR least-squares solution is obtained by solving which system?

A. Qx = R^Tb

B. Rx = Q^Tb

C. Ax = Q^Tb

D. R^Tx = Qb

Correct Answer: B. Rx = Q^Tb

Question 7

What does the rank of a matrix equal in terms of its SVD?

A. The largest singular value

B. The sum of singular values

C. The number of nonzero singular values

D. The number of rows minus columns

Correct Answer: C. The number of nonzero singular values

Question 8

What does the Eckart–Young theorem establish?

A. LU is always unique

B. Every matrix is positive definite

C. Truncated SVD gives an optimal rank-k approximation

D. Pseudoinverses require square matrices

Correct Answer: C. Truncated SVD gives an optimal rank-k approximation

Question 9

If a matrix has singular values 5 and 2, what is its 2-norm condition number when it has full rank?

A. 2.5

B. 7

C. 10

D. 0.4

Correct Answer: A. 2.5

Question 10

In a latent-factor recommender model, what do user and item factor vectors represent?

A. Exact observed ratings only

B. Latent patterns used to approximate preferences

C. Permutation matrices

D. Numerical pivoting steps

Correct Answer: B. Latent patterns used to approximate preferences