[,1] [,2]
[1,] -0.4472 -0.6325
[2,] -0.4472 -0.3162
[3,] -0.4472 0.0000
[4,] -0.4472 0.3162
[5,] -0.4472 0.6325 [,1] [,2]
[1,] -2.2361 -6.7082
[2,] 0.0000 3.16239 Matrix Decompositions
This chapter is the payoff for the last four. A decomposition factors a matrix into a product of simpler ones — triangular, orthogonal, diagonal — each of which is easy to compute with.
It is the same move as factoring an integer. \(\mathbf{A}\) tells you little at a glance; \(\mathbf{Q}\mathbf{R}\) or \(\mathbf{U}\boldsymbol{\Sigma}\mathbf{V}^\top\) tells you its rank, its conditioning, what it stretches, what it destroys, and how to solve with it. Essentially every numerical linear algebra routine you will ever call is a decomposition followed by something cheap.
9.1 Why decompose a matrix
Three reasons, in increasing order of importance.
Speed. Solving \(\mathbf{A}\mathbf{x} = \mathbf{b}\) costs \(O(n^3)\). If you factor once, each additional right-hand side costs \(O(n^2)\) (Section 6.6).
Stability. Orthogonal factors have condition number 1 (Section 7.11), so algorithms built from them do not amplify error. This is why QR beats the normal equations.
Insight. The factors mean something. The singular values measure how far a matrix is from singular; the eigenvalues govern long-run dynamics; the leading singular vectors are the directions your data actually varies in. A decomposition is a diagnosis, not just a computation.
The main ones, and what each requires:
| Decomposition | Form | Requires |
|---|---|---|
| LU | \(\mathbf{P}\mathbf{A} = \mathbf{L}\mathbf{U}\) | square |
| QR | \(\mathbf{A} = \mathbf{Q}\mathbf{R}\) | any shape |
| Cholesky | \(\mathbf{A} = \mathbf{R}^\top\mathbf{R}\) | symmetric positive definite |
| Eigen | \(\mathbf{A} = \mathbf{V}\boldsymbol{\Lambda}\mathbf{V}^{-1}\) | square, non-defective |
| SVD | \(\mathbf{A} = \mathbf{U}\boldsymbol{\Sigma}\mathbf{V}^\top\) | any matrix at all |
9.2 LU
Covered in Section 6.6: Gaussian elimination, recorded. \(\mathbf{U}\) is what elimination produces, \(\mathbf{L}\) holds the multipliers used to get there, and in practice a permutation \(\mathbf{P}\) handles row swaps for stability.
It is the workhorse behind solve() for general square systems. Use it when the matrix is square and has no special structure. If it does have structure — symmetry, positive definiteness — something better is available.
9.3 QR
QR factors any \(m \times n\) matrix with \(m \geq n\) as
\[ \mathbf{A} = \mathbf{Q}\mathbf{R} \]
with \(\mathbf{Q}\) having orthonormal columns and \(\mathbf{R}\) upper triangular. The columns of \(\mathbf{Q}\) are an orthonormal basis for the column space of \(\mathbf{A}\) — exactly what Gram–Schmidt produces (Section 7.12), though computed stably with Householder reflections instead.
QR’s headline use is least squares. Substituting \(\mathbf{A} = \mathbf{Q}\mathbf{R}\) into the normal equations and using \(\mathbf{Q}^\top\mathbf{Q} = \mathbf{I}\):
\[ \mathbf{R}^\top\mathbf{Q}^\top\mathbf{Q}\mathbf{R}\hat{\boldsymbol{\beta}} = \mathbf{R}^\top\mathbf{Q}^\top\mathbf{y} \;\Longrightarrow\; \mathbf{R}\hat{\boldsymbol{\beta}} = \mathbf{Q}^\top\mathbf{y} \tag{9.1}\]
A triangular system, solved by back-substitution. Crucially, \(\mathbf{X}^\top\mathbf{X}\) never appears, so the condition number is never squared (Section 6.9).
[1] 0.05 1.99[1] 0.05 1.99This is what lm() does. When someone says least squares is solved “by QR”, this is the equation they mean.
9.4 Cholesky
For a symmetric positive definite matrix (Section 8.9) there is a specialized factorization:
\[ \mathbf{A} = \mathbf{R}^\top\mathbf{R} \]
with \(\mathbf{R}\) upper triangular — a matrix square root. It exists precisely when \(\mathbf{A}\) is positive definite, which makes chol() a test as well as a factorization: if it fails, the matrix was not positive definite.
Worked example. For \(\mathbf{A} = \begin{bmatrix}2&1\\1&2\end{bmatrix}\), we need upper triangular \(\mathbf{R}\) with \(\mathbf{R}^\top\mathbf{R} = \mathbf{A}\). Matching entries: \(r_{11}^2 = 2\) so \(r_{11} = \sqrt{2}\); then \(r_{11}r_{12} = 1\) gives \(r_{12} = 1/\sqrt{2}\); finally \(r_{12}^2 + r_{22}^2 = 2\) gives \(r_{22} = \sqrt{3/2}\).
[,1] [,2]
[1,] 1.4142 0.7071
[2,] 0.0000 1.2247 [,1] [,2]
[1,] 2 1
[2,] 1 2Cholesky is about twice as fast as LU because it exploits symmetry to do half the work, and it needs no pivoting. Whenever you have a covariance matrix, a kernel matrix, or \(\mathbf{X}^\top\mathbf{X}\), this is the right tool.
n <- 600
Z <- matrix(rnorm(n * n), n)
S <- crossprod(Z) + diag(n) # symmetric positive definite
bb <- rnorm(n)
chol_solve <- function(S, b) {
R <- chol(S)
backsolve(R, backsolve(R, b, transpose = TRUE))
}
c(
lu = unname(system.time(solve(S, bb))["elapsed"]),
chol = unname(system.time(chol_solve(S, bb))["elapsed"])
) lu chol
0.106 0.048
NoteIn machine learning
Cholesky is how you sample from a multivariate normal. To draw \(\mathbf{x} \sim \mathcal{N}(\boldsymbol{\mu}, \boldsymbol{\Sigma})\), factor \(\boldsymbol{\Sigma} = \mathbf{R}^\top\mathbf{R}\), draw standard normal \(\mathbf{z}\), and set \(\mathbf{x} = \boldsymbol{\mu} + \mathbf{R}^\top\mathbf{z}\). It is also the inner loop of Gaussian process regression, where the same kernel matrix is factored once and reused.
9.5 Eigendecomposition
Covered in Section 8.5: \(\mathbf{A} = \mathbf{V}\boldsymbol{\Lambda} \mathbf{V}^{-1}\), simplifying to \(\mathbf{Q}\boldsymbol{\Lambda}\mathbf{Q}^\top\) when \(\mathbf{A}\) is symmetric.
Its limitation is worth restating, because it motivates everything that follows: the eigendecomposition requires a square matrix, and even then may not exist. A \(100 \times 3\) design matrix has no eigenvalues at all.
9.6 The singular value decomposition
The SVD removes every restriction. Every matrix — square or rectangular, full rank or not, real or complex — factors as
\[ \mathbf{A} = \mathbf{U}\boldsymbol{\Sigma}\mathbf{V}^\top \tag{9.2}\]
where, for \(\mathbf{A}\) of shape \(m \times n\):
- \(\mathbf{U}\) is \(m \times m\) orthogonal; its columns are the left singular vectors
- \(\boldsymbol{\Sigma}\) is \(m \times n\) diagonal with non-negative entries \(\sigma_1 \geq \sigma_2 \geq \cdots \geq 0\), the singular values
- \(\mathbf{V}\) is \(n \times n\) orthogonal; its columns are the right singular vectors
The singular values are ordered, always real, and always non-negative. No caveats, no exceptions.
The connection to eigenvalues is direct. Since \(\mathbf{A}^\top\mathbf{A} = \mathbf{V}\boldsymbol{\Sigma}^\top\boldsymbol{\Sigma} \mathbf{V}^\top\), the right singular vectors are the eigenvectors of \(\mathbf{A}^\top\mathbf{A}\) — symmetric and positive semi-definite, so the spectral theorem applies — and
\[ \sigma_i = \sqrt{\lambda_i(\mathbf{A}^\top\mathbf{A})} \]
Worked example. Take \(\mathbf{M} = \begin{bmatrix}0&2\\1&0\end{bmatrix}\), which swaps the axes and stretches one of them.
\[ \mathbf{M}^\top\mathbf{M} = \begin{bmatrix}1&0\\0&4\end{bmatrix} \]
Eigenvalues 4 and 1, so \(\sigma_1 = 2\) and \(\sigma_2 = 1\), with right singular vectors \(\mathbf{v}_1 = \mathbf{e}_2\) and \(\mathbf{v}_2 = \mathbf{e}_1\). Then \(\mathbf{u}_i = \mathbf{M}\mathbf{v}_i / \sigma_i\) gives \(\mathbf{u}_1 = \mathbf{e}_1\) and \(\mathbf{u}_2 = \mathbf{e}_2\).
[1] 2 1round(sv$u, 4) [,1] [,2]
[1,] -1 0
[2,] 0 -1round(sv$v, 4) [,1] [,2]
[1,] 0 -1
[2,] -1 0 [,1] [,2]
[1,] 0 2
[2,] 1 0Signs may differ from the hand computation — \(\mathbf{u}_i\) and \(\mathbf{v}_i\) can both be negated without changing the product.
The singular values answer questions nothing else does as cleanly:
| Question | Answer from the SVD |
|---|---|
| What is the rank? | number of non-zero \(\sigma_i\) |
| How close to singular? | how small is \(\sigma_{\min}\) |
| Condition number | \(\sigma_1 / \sigma_{\min}\) |
| Best rank-\(k\) approximation | keep the top \(k\) |
9.7 Geometry of the SVD
Equation 9.2 says every linear map, however complicated, is really three simple ones:
Rotate (\(\mathbf{V}^\top\)), stretch along the axes (\(\boldsymbol{\Sigma}\)), rotate again (\(\mathbf{U}\)).
That is the whole content. There is no fourth kind of thing a matrix can do.
Concretely: \(\mathbf{A}\) maps the unit circle to an ellipse. The right singular vectors \(\mathbf{v}_i\) are the directions in the input that end up as the ellipse’s axes; the left singular vectors \(\mathbf{u}_i\) are those axis directions in the output; and the singular values are the axis lengths.
G <- matrix(c(2, 1, 0, 1), nrow = 2, byrow = TRUE)
sg <- svd(G)
th <- seq(0, 2 * pi, length.out = 400)
circ <- cbind(cos(th), sin(th))
ell <- t(G %*% t(circ))
draw_plane(
curves = list(circ, ell),
curve_color = c(amds_gray, amds_colors[1]),
vectors = list(
v1 = sg$v[, 1], v2 = sg$v[, 2],
s1u1 = sg$d[1] * sg$u[, 1],
s2u2 = sg$d[2] * sg$u[, 2]
),
color = c(
amds_gray, amds_gray,
amds_colors[2], amds_colors[3]
),
line_type = c("dashed", "dashed", "solid", "solid"),
label_at = "mid"
)Compare Figure 8.2, where the matrix was symmetric and the eigenvectors were already the ellipse axes. For a non-symmetric matrix the eigenvectors are not perpendicular and are not the axes — the singular vectors are. This is why the SVD is the more reliable description.
9.8 Low-rank approximation
Write Equation 9.2 as a sum of rank-one pieces:
\[ \mathbf{A} = \sum_{i=1}^{r} \sigma_i\,\mathbf{u}_i\mathbf{v}_i^\top \tag{9.3}\]
Each term is an outer product — a rank-one matrix — weighted by its singular value. Since the \(\sigma_i\) decrease, the terms matter less and less. Truncating after \(k\) terms gives \(\mathbf{A}_k\), and the Eckart–Young theorem says something remarkable:
\(\mathbf{A}_k\) is the best possible rank-\(k\) approximation to \(\mathbf{A}\), in both the Frobenius and spectral norms.
Not a good heuristic — provably optimal. No other rank-\(k\) matrix is closer.
Here is a \(40 \times 40\) matrix built from two smooth patterns plus noise:
[1] 27.432 4.477 0.590 0.575 0.525 0.516Two large singular values, then a floor of small ones — the noise. The structure is essentially two-dimensional.
draw_line(
x = 1:20,
y = sz$d[1:20],
xlab = "index",
ylab = "singular value"
)[1] 0.17535021 0.07018404 0.06690933 0.06364460draw_heatmap(Z, show_colorbar = FALSE)draw_heatmap(lowrank(sz, 2), show_colorbar = FALSE)The rank-2 version stores \(2(40 + 40 + 1) = 162\) numbers instead of 1600, a tenfold saving, and it is cleaner than the original — discarding small singular values discards noise.
NoteIn machine learning
Low-rank approximation is one idea wearing many hats. PCA is the SVD of a centered data matrix. Latent semantic analysis is the SVD of a term–document matrix. Recommender systems factor a sparse ratings matrix into user and item factors. LoRA fine-tunes a large model by learning a low-rank update to its weight matrices. All of them are Equation 9.3, truncated.
9.9 Principal component analysis from the SVD
PCA finds the directions in which data varies most. It is the SVD of the centered data matrix, and nothing more.
Center \(\mathbf{X}\) so each column has mean zero, then take \(\mathbf{X} = \mathbf{U}\boldsymbol{\Sigma}\mathbf{V}^\top\):
- The columns of \(\mathbf{V}\) are the principal component directions — orthogonal by construction, which is why component scores are uncorrelated.
- \(\mathbf{U}\boldsymbol{\Sigma}\) holds the scores, the data in the new coordinates.
- The variance along component \(i\) is \(\sigma_i^2/(n-1)\).
That last identity follows from the covariance matrix being \(\mathbf{X}^\top\mathbf{X}/(n-1) = \mathbf{V}\boldsymbol{\Sigma}^2\mathbf{V}^\top/(n-1)\) — the eigendecomposition of the covariance, obtained without ever forming it.
from_svd from_cov
1.80217 1.80217 The figure is drawn with equal axis scaling, which matters: PCA directions are only perpendicular in a picture that does not distort angles.
WarningWatch out
Two traps.
Center your data. Without centering, the first component points at the mean rather than the direction of greatest variance. prcomp() centers by default; svd() does not.
Think before scaling. If your variables are in different units, whichever has the largest numbers will dominate. Scaling to unit variance fixes that but throws away real information when the units are comparable. prcomp(scale. = TRUE) is a modeling choice, not a default.
9.10 The pseudoinverse
What should \(\mathbf{A}^{-1}\) mean when \(\mathbf{A}\) is not square, or is singular? The Moore–Penrose pseudoinverse \(\mathbf{A}^{+}\) answers this by inverting through the SVD:
\[ \mathbf{A} = \mathbf{U}\boldsymbol{\Sigma}\mathbf{V}^\top \quad\Longrightarrow\quad \mathbf{A}^{+} = \mathbf{V}\boldsymbol{\Sigma}^{+}\mathbf{U}^\top \tag{9.4}\]
where \(\boldsymbol{\Sigma}^{+}\) inverts each non-zero singular value and leaves the zeros alone. Inverting what is invertible, ignoring what is not.
It does exactly what you would want in both awkward cases:
- Overdetermined (tall, full column rank): \(\mathbf{A}^{+}\mathbf{b}\) is the least squares solution, identical to \((\mathbf{A}^\top\mathbf{A})^{-1}\mathbf{A}^\top \mathbf{b}\) but computed stably.
- Underdetermined or rank deficient: among the infinitely many solutions, \(\mathbf{A}^{+}\mathbf{b}\) returns the one with smallest norm.
Worked example. The rank-2 matrix from Section 7.8, with a right-hand side chosen to lie in its column space:
pinv <- function(A, tol = 1e-10) {
s <- svd(A)
keep <- s$d > tol * max(s$d)
v <- s$v[, keep, drop = FALSE]
u <- s$u[, keep, drop = FALSE]
v %*% (t(u) / s$d[keep])
}
Ad <- matrix(
c(1, 2, 3, 2, 4, 6, 1, 1, 2),
nrow = 3, byrow = TRUE
)
bd <- drop(Ad %*% c(1, 0, 0))
xp <- drop(pinv(Ad) %*% bd)
round(xp, 6)[1] 0.666667 -0.333333 0.333333Check that it is a solution, and that it is the smallest one. The obvious solution \((1,0,0)\) has norm 1; any \((1,0,0) + t(1,1,-1)\) is also a solution, since \((1,1,-1)\) spans the null space:
[1] 0 0 0 pinv obvious
0.8164966 1.0000000 [1] 0Smaller norm, and orthogonal to the null space — those two facts are the same fact. The minimum-norm solution is the one with no wasted component in the directions the matrix ignores.
WarningWatch out
The tol argument is doing real work. Deciding which singular values count as zero is a judgment call, and on noisy data the answer is not obvious. This is the practical face of the warning in Section 7.8: rank is discrete, floating point is not.
9.11 Choosing a decomposition
| Situation | Use | R |
|---|---|---|
| Square system, no structure | LU | solve(A, b) |
| Symmetric positive definite | Cholesky | chol(A) |
| Least squares | QR |
lm(), qr.solve()
|
| Symmetric, want eigenvalues | Eigen | eigen(A, symmetric = TRUE) |
| Rank, conditioning, PCA, anything rectangular | SVD |
svd(A), prcomp()
|
| Rank deficient or underdetermined | Pseudoinverse | via svd()
|
Two rules of thumb worth internalizing.
Exploit structure when you have it. Cholesky beats LU on positive definite matrices; LU beats SVD on well-conditioned square ones. Using a more general tool than necessary costs time.
When in doubt, use the SVD. It always exists, never lies about rank or conditioning, and degrades gracefully. It is the most expensive of these — roughly an order of magnitude more than LU — but on anything that fits in memory that is usually a price worth paying for an answer you can trust.
9.12 Summary
| Decomposition | Form | Best for |
|---|---|---|
| LU | \(\mathbf{P}\mathbf{A} = \mathbf{L}\mathbf{U}\) | general square systems |
| QR | \(\mathbf{A} = \mathbf{Q}\mathbf{R}\) | least squares |
| Cholesky | \(\mathbf{A} = \mathbf{R}^\top\mathbf{R}\) | symmetric positive definite |
| Eigen | \(\mathbf{A} = \mathbf{Q}\boldsymbol{\Lambda}\mathbf{Q}^\top\) | symmetric matrices, dynamics |
| SVD | \(\mathbf{A} = \mathbf{U}\boldsymbol{\Sigma}\mathbf{V}^\top\) | everything else |
9.13 Exercises
1. Compute the QR decomposition of cbind(1, 1:4) and verify \(\mathbf{Q}^\top
\mathbf{Q} = \mathbf{I}\) and \(\mathbf{Q}\mathbf{R} = \mathbf{A}\).
TipSolution
[,1] [,2]
[1,] 1 0
[2,] 0 1 [,1] [,2]
[1,] 1 1
[2,] 1 2
[3,] 1 3
[4,] 1 4\(\mathbf{Q}\) is \(4 \times 2\) here — orthonormal columns, but not square, so \(\mathbf{Q}\mathbf{Q}^\top \neq \mathbf{I}\). Only \(\mathbf{Q}^\top\mathbf{Q}\) is the identity. This is the “thin” QR, which is all least squares needs.
2. Try chol() on \(\begin{bmatrix}2&3\\3&2\end{bmatrix}\). Explain what happens.
TipSolution
[1] 5 -1[1] "the leading minor of order 2 is not positive"It fails. The matrix is symmetric but has a negative eigenvalue, so it is indefinite (Section 8.9) and has no Cholesky factor. The failure is useful: chol() doubles as a positive definiteness test, and the error is the answer rather than a bug.
3. For the \(40 \times 40\) matrix Z, plot relative approximation error against \(k\). How many components do you need for 1% error, and why does the curve flatten?
TipSolution
[1] 0.1754 0.0702 0.0669 0.0636 0.0608 0.0579 0.0550 0.0522 0.0494 0.0467
[11] 0.0442 0.0416 0.0392 0.0370 0.0347draw_line(
x = ks, y = errs,
xlab = "k", ylab = "relative error"
)The error falls sharply to \(k=2\) and then flattens. Only two components carry signal; beyond that each additional one removes a sliver of noise, and since the noise is spread evenly across the remaining ~38 dimensions, no small number of them helps much. You never reach 1% error, because about 7% of the matrix is noise.
The shape of this curve is how you choose \(k\) in practice — look for the elbow, not a fixed threshold.
4. Confirm that prcomp() and svd() agree on the data Xd from Section 9.9.
TipSolution
[1] 1.342449 0.335768[1] 1.342449 0.335768 PC1 PC2
a1 0.684971 0.728571
a2 0.728571 0.684971 [,1] [,2]
[1,] 0.684971 0.728571
[2,] 0.728571 0.684971Identical, up to sign. prcomp() is svd() on centered data, with the scaling absorbed into sdev. Taking absolute values sidesteps the arbitrary sign of each singular vector.
5. Verify that the condition number equals \(\sigma_1/\sigma_n\) for the ill-conditioned matrix from Section 6.10.
TipSolution
ratio kappa
82.03781 82.03781 The same number. kappa(exact = TRUE) computes the SVD and takes the ratio — which is why the condition number is the honest measure of near-singularity and the determinant is not.
6. Use the pseudoinverse to solve an underdetermined system: one equation, \(x_1 + x_2 + x_3 = 3\), in three unknowns. Which of the infinitely many solutions do you get?
TipSolution
[1] 1 1 1pinv_norm alt_norm
1.732051 3.000000 You get \((1, 1, 1)\) — the minimum-norm solution. Both \((3,0,0)\) and \((0,3,0)\) satisfy the equation, but they have norm 3 against \(\sqrt{3} \approx 1.73\).
Geometrically, the solutions form a plane in \(\mathbb{R}^3\), and \((1,1,1)\) is the point on that plane closest to the origin. That is the same minimum-norm principle behind ridge regression (Section 6.11): when the data cannot choose, prefer the smallest answer.