SOCR ≫ DSPA ≫ DSPA3 Topics ≫ 

library(ggplot2)
library(dplyr)
library(tidyr)
library(patchwork)
library(plotly)     # interactive figures and ALL 3-D graphics
library(Matrix)
library(DT)

How this chapter uses graphics

Every two-dimensional figure is drawn with ggplot2 and rendered statically. Immediately after each one, the equivalent plot_ly() code appears in a chunk marked eval=FALSE, echo=TRUE, visible in these notes, ready to paste into a live session for an interactive demonstration.

Every three-dimensional figure is drawn with plot_ly() and evaluated, because rotation carries information a fixed projection destroys. In this chapter that includes the residual-sum-of-squares surface over parameter space, the covariance ellipsoid with its principal axes, the geometry of orthogonal projection, and the fitted regression plane.

1 Learning objectives

After completing this chapter you will be able to:

  1. Construct and manipulate matrices in R, and state the cost of each operation in floating-point operations.
  2. Define rank, null space, and the four fundamental subspaces, and explain why \(\operatorname{rank}(X)=p\) is the condition least squares actually requires.
  3. Derive the orthogonal projection \(\hat{Y}=HY\) and prove that \(H\) is symmetric and idempotent with \(\operatorname{tr}(H)=p\).
  4. Compute and compare the LU, Cholesky, QR, eigen, and singular value decompositions, and choose the right one for a given task.
  5. Quantify numerical conditioning, and explain why forming \(X^\top X\) costs you half your significant digits.
  6. Derive the OLS estimator three ways, calculus, geometry, and maximum likelihood, and state the Gauss–Markov theorem precisely.
  7. Derive \(\operatorname{Var}(\hat\beta)=\sigma^2(X^\top X)^{-1}\) and use it to construct every standard error, \(t\)-statistic, and \(F\)-test in a regression summary.
  8. Diagnose a fitted model using leverage, standardized residuals, and Cook’s distance, and explain each in terms of the hat matrix.
  9. Distinguish multicollinearity from dependence, compute VIFs, and explain why spurious correlation is guaranteed when \(k \gg n\).
  10. Fit and interpret regression trees, model trees, and Bayesian additive regression trees, and state the computational cost of each.

Estimated time: 9–12 hours including exercises. Prerequisites: Chapter 1 (reproducibility conventions, dspa_read()) and Chapter 2 (covariance, correlation, EDA). Readers wanting a refresher on calculus, series, and complex numbers may find the mathematical foundations module useful.

1.1 Chapter 3 Live Webapps

The DSPA3 Labs include a number of interactive webapps, which are build as portable, self-contained, and browser accessible applications. These apps demonstrate various concepts discussed in this chapter, including

2 PART I: LINEAR ALGEBRA AND MATRIX COMPUTING

Linear algebra is the study of linear maps between vector spaces. Its practical importance rests on a single observation: almost every process is locally linear. A smooth function \(f\) near a point \(x_0\) satisfies

\[f(x) = f(x_0) + \underbrace{Df(x_0)}_{\text{a matrix}}(x-x_0) + O\!\left(\|x-x_0\|^2\right),\]

so the first-order behaviour of any differentiable system is a matrix. That is why arbitrarily complicated problems collapse into two canonical forms:

\[\underbrace{A\mathbf{x}=\mathbf{b}}_{\text{solve a system}} \qquad\text{and}\qquad \underbrace{A\mathbf{v}=\lambda\mathbf{v}}_{\text{find invariant directions}}.\]

Everything in this chapter is an elaboration of those two equations.

3 Building and indexing matrices

3.1 Construction

seq1 <- 1:6
m1 <- matrix(seq1, nrow = 2, ncol = 3)   # filled COLUMN-major by default
m1
#>      [,1] [,2] [,3]
#> [1,]    1    3    5
#> [2,]    2    4    6
matrix(seq1, nrow = 2, ncol = 3, byrow = TRUE)   # row-major
#>      [,1] [,2] [,3]
#> [1,]    1    2    3
#> [2,]    4    5    6
set.seed(11)
m3 <- matrix(rnorm(20), nrow = 5)
round(m3, 3)
#>        [,1]   [,2]   [,3]   [,4]
#> [1,] -0.591 -0.934 -0.828  0.012
#> [2,]  0.027  1.324 -0.348 -0.223
#> [3,] -1.517  0.625 -1.538  0.888
#> [4,] -1.363 -0.046 -0.256 -0.592
#> [5,]  1.178 -1.004 -1.150 -0.656

R stores matrices in column-major order, a matrix is a vector with a dim attribute, laid out one column after another. This is the same convention as Fortran and LAPACK, and the reason column operations (colSums, X[, j]) touch contiguous memory and are faster than the corresponding row operations on large matrices.

as.vector(m1)          # reading the underlying storage confirms column-major
#> [1] 1 2 3 4 5 6
attributes(m1)
#> $dim
#> [1] 2 3

diag() is overloaded three ways, keyed on the type of its argument:

diag(c(1, 2, 3))   # vector  -> diagonal matrix
#>      [,1] [,2] [,3]
#> [1,]    1    0    0
#> [2,]    0    2    0
#> [3,]    0    0    3
diag(m1)           # matrix  -> its principal diagonal (works for non-square)
#> [1] 1 4
diag(4)            # scalar  -> 4x4 identity
#>      [,1] [,2] [,3] [,4]
#> [1,]    1    0    0    0
#> [2,]    0    1    0    0
#> [3,]    0    0    1    0
#> [4,]    0    0    0    1
c1 <- 1:5
m4 <- cbind(m3, c1)          # column concatenation; picks up the name "c1"
r1 <- 1:4
m5 <- rbind(m3, r1)
dimnames(m5) <- list(NULL, NULL)   # strip inherited names
round(m5, 3)
#>        [,1]   [,2]   [,3]   [,4]
#> [1,] -0.591 -0.934 -0.828  0.012
#> [2,]  0.027  1.324 -0.348 -0.223
#> [3,] -1.517  0.625 -1.538  0.888
#> [4,] -1.363 -0.046 -0.256 -0.592
#> [5,]  1.178 -1.004 -1.150 -0.656
#> [6,]  1.000  2.000  3.000  4.000

3.2 Indexing

m6 <- matrix(1:12, nrow = 3)
m6
#>      [,1] [,2] [,3] [,4]
#> [1,]    1    4    7   10
#> [2,]    2    5    8   11
#> [3,]    3    6    9   12
m6[1, 2]          # single element
#> [1] 4
m6[1, ]           # row 1 -> DROPS to a vector
#> [1]  1  4  7 10
m6[1, , drop = FALSE]   # row 1 as a 1 x 4 matrix
#>      [,1] [,2] [,3] [,4]
#> [1,]    1    4    7   10
m6[, c(2, 3)]     # two columns
#>      [,1] [,2]
#> [1,]    4    7
#> [2,]    5    8
#> [3,]    6    9

Common misconception: “A[i, ] gives me a row of a matrix.” It gives a vector. R silently drops the dimension when a subscript selects a single row or column, so dim(m6[1, ]) is NULL and any subsequent %*% will treat it as a column vector. In code that must work for both a single row and several, always write drop = FALSE. This is one of the most common sources of “non-conformable arguments” errors.

4 Matrix algebra

4.1 Addition, subtraction, and scalar recycling

Addition and subtraction are elementwise and require identical dimensions.

m7 <- matrix(1:6, nrow = 2)
m8 <- matrix(2:7, nrow = 2)
m7 + m8
#>      [,1] [,2] [,3]
#> [1,]    3    7   11
#> [2,]    5    9   13
m8 - m7
#>      [,1] [,2] [,3]
#> [1,]    1    1    1
#> [2,]    1    1    1
m8 - 1          # scalar recycled over every entry
#>      [,1] [,2] [,3]
#> [1,]    1    3    5
#> [2,]    2    4    6

Cost: \(O(mn)\), one operation per entry, no arithmetic can avoid it.

4.2 Elementwise product versus matrix product

R uses * for the Hadamard (elementwise) product and %*% for the matrix product. They are entirely different operations that happen to share notation in loose writing.

\[(A\circ B)_{ij}=A_{ij}B_{ij} \qquad\text{versus}\qquad (AB)_{ij}=\sum_{\ell=1}^{n}A_{i\ell}B_{\ell j}.\]

The Hadamard product needs \(A\) and \(B\) to have the same shape. The matrix product needs the inner dimensions to match: \(P_{m\times k}=L_{m\times n}R_{n\times k}\).

m8 * m7                      # Hadamard: 2x3 with 2x3
#>      [,1] [,2] [,3]
#> [1,]    2   12   30
#> [2,]    6   20   42
m9 <- matrix(3:8, nrow = 3)  # 3x2
dim(m8); dim(m9)
#> [1] 2 3
#> [1] 3 2
m8 %*% m9                    # (2x3)(3x2) -> 2x2
#>      [,1] [,2]
#> [1,]   52   88
#> [2,]   64  109

Complexity. The definition requires \(mnk\) multiply–add pairs, so the textbook algorithm is \(O(mnk)\), or \(O(n^3)\) for square matrices. This is not optimal: Strassen’s algorithm (1969) achieves \(O(n^{\log_2 7})=O(n^{2.807})\) by trading one multiplication for several additions, and successive refinements have pushed the theoretical exponent below \(2.372\). In practice R calls a tuned BLAS (dgemm), which uses the \(O(n^3)\) algorithm with cache blocking, asymptotically worse but far faster at realistic sizes, because memory traffic dominates flop count.

set.seed(3)
A <- matrix(rnorm(500 * 500), 500)
B <- matrix(rnorm(500 * 500), 500)
system.time(A %*% B)[["elapsed"]]     # BLAS dgemm
#> [1] 0.03

4.3 Inner and outer products

Two vectors combine in two fundamentally different ways.

\[\underbrace{\mathbf{u}^\top\mathbf{v}=\sum_i u_iv_i}_{\text{inner product: a SCALAR}} \qquad\qquad \underbrace{\mathbf{u}\mathbf{v}^\top=\big(u_iv_j\big)}_{\text{outer product: a MATRIX of rank 1}}\]

u <- c(1, 2, 3, 4, 5)
v <- c(4, 5, 6, 7, 8)

u %*% v              # 1x1: R silently makes the left operand a ROW vector
#>      [,1]
#> [1,]  100
crossprod(u, v)      # the explicit, preferred form of the inner product
#>      [,1]
#> [1,]  100
sum(u * v)           # same number
#> [1] 100
u %o% v              # outer product: 5x5
#>      [,1] [,2] [,3] [,4] [,5]
#> [1,]    4    5    6    7    8
#> [2,]    8   10   12   14   16
#> [3,]   12   15   18   21   24
#> [4,]   16   20   24   28   32
#> [5,]   20   25   30   35   40
u %*% t(v)           # identical
#>      [,1] [,2] [,3] [,4] [,5]
#> [1,]    4    5    6    7    8
#> [2,]    8   10   12   14   16
#> [3,]   12   15   18   21   24
#> [4,]   16   20   24   28   32
#> [5,]   20   25   30   35   40
Matrix::rankMatrix(u %*% t(v))[1]   # every outer product has rank 1
#> [1] 1

Common misconception: “multiplying two vectors gives the outer product.” It depends entirely on orientation. \(\mathbf{u}^\top\mathbf{v}\) is a scalar measuring alignment; \(\mathbf{u}\mathbf{v}^\top\) is an \(n\times n\) matrix of rank one. R blurs this because a bare vector has no orientation, in u %*% v it is treated as a row on the left and a column on the right. Write crossprod(u, v) for the inner product and tcrossprod(u, v) for the outer product, and the ambiguity disappears.

Rank-one outer products are the atoms of matrix analysis: the singular value decomposition (§3.6.5) writes any matrix as a weighted sum of them.

4.4 Transpose, crossprod, and why it matters

m8
#>      [,1] [,2] [,3]
#> [1,]    2    4    6
#> [2,]    3    5    7
t(m8)
#>      [,1] [,2]
#> [1,]    2    3
#> [2,]    4    5
#> [3,]    6    7
identical(m8[1, 2], t(m8)[2, 1])
#> [1] TRUE

The combination \(A^\top B\) is so ubiquitous in statistics that R provides a dedicated routine:

set.seed(5)
X <- matrix(rnorm(2000 * 200), 2000)

t1 <- system.time(r1 <- t(X) %*% X)[["elapsed"]]
t2 <- system.time(r2 <- crossprod(X))[["elapsed"]]
c(explicit_transpose = t1, crossprod = t2, identical = all.equal(r1, r2))
#> explicit_transpose          crossprod          identical 
#>               0.02               0.03               1.00

crossprod(X) never materializes \(X^\top\), it calls BLAS dsyrk, which exploits the symmetry of the result and computes only the lower triangle, roughly halving the work. Use crossprod(A, B) for \(A^\top B\) and tcrossprod(A, B) for \(AB^\top\).

5 Vector spaces, rank, and the four subspaces

Everything about the solvability of \(A\mathbf{x}=\mathbf{b}\) and the identifiability of a regression model is a statement about subspaces.

5.1 Span, independence, basis

Given vectors \(\mathbf{a}_1,\dots,\mathbf{a}_n\in\mathbb{R}^m\), their span is the set of all linear combinations,

\[\operatorname{span}\{\mathbf{a}_1,\dots,\mathbf{a}_n\}=\left\{\sum_{j=1}^n c_j\mathbf{a}_j : c_j\in\mathbb{R}\right\}.\]

They are linearly independent if \(\sum_j c_j\mathbf{a}_j=\mathbf{0}\) forces every \(c_j=0\). A basis for a subspace is an independent spanning set; every basis of a given subspace has the same size, called its dimension.

5.2 The four fundamental subspaces

For \(A\in\mathbb{R}^{m\times n}\):

Subspace Definition Lives in Dimension
Column space, aka Image space, \(\operatorname{im}(A)=\operatorname{col}(A)\) \(\{A\mathbf{x}:\mathbf{x}\in\mathbb{R}^n\}\) \(\mathbb{R}^m\) \(r\)
Null space, Kernel space, \(\operatorname{ker(A)}=\mathcal{N}(A)\) \(\{\mathbf{x}:A\mathbf{x}=\mathbf{0}\}\) \(\mathbb{R}^n\) \(n-r\)
Row space \(\operatorname{col}(A^\top)\) \(\{A^\top\mathbf{y}:\mathbf{y}\in\mathbb{R}^m\}\) \(\mathbb{R}^n\) \(r\)
Left null space \(\mathcal{N}(A^\top)\) \(\{\mathbf{y}:A^\top\mathbf{y}=\mathbf{0}\}\) \(\mathbb{R}^m\) \(m-r\)

Here \(r=\operatorname{rank}(A)\), the number of linearly independent columns, which equals the number of linearly independent rows, a fact that is not obvious and is worth pausing on.

Rank–nullity theorem. \[\operatorname{rank}(A)+\dim\mathcal{N}(A)=n.\]

Orthogonality relations. \(\mathcal{N}(A)\perp\operatorname{col}(A^\top)\) and \(\mathcal{N}(A^\top)\perp\operatorname{col}(A)\). The second is the one that generates least squares (§3.9): the residual lives in the left null space, orthogonal to everything the model can produce.

In the example below, .Machine$double.eps represents the machine epsilon for double-precision (floating-point) numbers, which is the smallest positive floating-point number \(x\), such that \(1.0 + x \neq 1.0\) on your computer. In standard IEEE 754 double precision, its value is roughly \(2.220446 \times 10^{-16}\) and essentially measures the limit of precision for decimal numbers in R.

As computers cannot store real numbers with infinite precision, floating-point arithmetic always introduces tiny rounding errors. Thus, a singular value of a matrix that should mathematically be zero might instead show up as a tiny non-zero number, e.g., \(10^{-16}\) or \(10^{-17}\), due to computer noise. The expression max(dim(A)) * .Machine$double.eps * max(sv$d) creates a dynamic tolerance threshold. Any singular value (sv$d) smaller than this threshold is treated as zero (numerical noise). This allows the code to accurately calculate the numerical rank (r) of the matrix and find a reliable basis for the null space.

A <- matrix(c(1, 2, 3,
              2, 4, 6,
              1, 1, 1), nrow = 3, byrow = TRUE)
A
#>      [,1] [,2] [,3]
#> [1,]    1    2    3
#> [2,]    2    4    6
#> [3,]    1    1    1
Matrix::rankMatrix(A)[1]     # row 2 = 2 * row 1, so rank is 2, not 3
#> [1] 2
# A basis for the null space, from the SVD's trailing right singular vectors
sv <- svd(A)
r  <- sum(sv$d > max(dim(A)) * .Machine$double.eps * max(sv$d))
null_basis <- sv$v[, (r + 1):ncol(A), drop = FALSE]
round(null_basis, 4)
#>         [,1]
#> [1,]  0.4082
#> [2,] -0.8165
#> [3,]  0.4082
round(A %*% null_basis, 10)  # A maps the null-space basis to zero
#>      [,1]
#> [1,]    0
#> [2,]    0
#> [3,]    0

Why this matters for regression. A design matrix \(X_{n\times p}\) produces a unique OLS solution iff \(\operatorname{rank}(X)=p\), i.e. \(\mathcal{N}(X)=\{\mathbf{0}\}\). If a column is an exact linear combination of others, a dummy-variable trap, a redundant total, \(n<p\), then infinitely many \(\hat\beta\) give identical fitted values, and no data can distinguish them. lm() responds by reporting NA for the aliased coefficients rather than failing.

set.seed(7)
d <- data.frame(x1 = rnorm(50), x2 = rnorm(50))
d$x3 <- 2 * d$x1 - 3 * d$x2          # exactly collinear by construction
d$y  <- d$x1 + d$x2 + rnorm(50)

coef(lm(y ~ x1 + x2 + x3, data = d)) # x3 is aliased -> NA
#> (Intercept)          x1          x2          x3 
#>    0.173273    1.081103    0.808563          NA

6 Norms, inner products, and orthogonal projection

6.1 Norms

A norm assigns a length. The \(\ell_p\) family on \(\mathbb{R}^n\):

\[\|\mathbf{x}\|_p=\left(\sum_{i=1}^n|x_i|^p\right)^{1/p},\qquad \|\mathbf{x}\|_1=\sum_i|x_i|,\quad \|\mathbf{x}\|_2=\sqrt{\sum_i x_i^2},\quad \|\mathbf{x}\|_\infty=\max_i|x_i|.\]

In practice, the choice of norm is a modelling decision with implications,e.g., \(\ell_2\) gives ridge regression and OLS, while \(\ell_1\) gives LASSO and median regression. The difference between them is the reason \(\ell_1\) produces exact zeros.

The inner product \(\langle\mathbf{u},\mathbf{v}\rangle=\mathbf{u}^\top\mathbf{v}\) induces \(\|\mathbf{x}\|_2=\sqrt{\langle\mathbf{x},\mathbf{x}\rangle}\) and satisfies the Cauchy–Schwarz inequality

\[|\langle\mathbf{u},\mathbf{v}\rangle|\le\|\mathbf{u}\|_2\,\|\mathbf{v}\|_2,\]

with equality iff \(\mathbf{u}\) and \(\mathbf{v}\) are parallel. Applying this to centred data vectors is exactly the statement \(|\rho_{xy}|\le 1\), the correlation coefficient is the cosine of the angle between centred variables:

\[\rho_{xy}=\frac{\langle \mathbf{x}-\bar x\mathbf{1},\; \mathbf{y}-\bar y\mathbf{1}\rangle} {\|\mathbf{x}-\bar x\mathbf{1}\|_2\,\|\mathbf{y}-\bar y\mathbf{1}\|_2}=\cos\theta .\]

set.seed(13)
x <- rnorm(100); y <- 0.6 * x + rnorm(100, sd = 0.8)
xc <- x - mean(x); yc <- y - mean(y)

c(cor          = cor(x, y),
  cosine_angle = sum(xc * yc) / (sqrt(sum(xc^2)) * sqrt(sum(yc^2))),
  angle_degrees = acos(cor(x, y)) * 180 / pi)
#>           cor  cosine_angle angle_degrees 
#>      0.518186      0.518186     58.789351

6.2 Orthogonal projection

Let \(S=\operatorname{col}(X)\subseteq\mathbb{R}^n\) with \(X\) of full column rank \(p\). The orthogonal projection of \(\mathbf{y}\) onto \(S\) is the unique point \(\hat{\mathbf{y}}\in S\) minimizing \(\|\mathbf{y}-\hat{\mathbf{y}}\|_2\). Writing \(\hat{\mathbf{y}}=X\mathbf{b}\) and requiring the residual to be orthogonal to every column of \(X\),

\[X^\top(\mathbf{y}-X\mathbf{b})=\mathbf{0} \;\Longrightarrow\; \mathbf{b}=(X^\top X)^{-1}X^\top\mathbf{y} \;\Longrightarrow\; \hat{\mathbf{y}}=\underbrace{X(X^\top X)^{-1}X^\top}_{H}\mathbf{y}.\]

\(H\) is the hat matrix, it puts the hat on \(\mathbf{y}\).

Proposition. \(H\) is symmetric, idempotent, and \(\operatorname{tr}(H)=p\). Likewise \(M=I-H\) is symmetric idempotent with \(\operatorname{tr}(M)=n-p\), and \(HM=0\).

Proof. Symmetry: \(H^\top=X\big((X^\top X)^{-1}\big)^\top X^\top=H\) since \(X^\top X\) is symmetric. Idempotence: \(H^2=X(X^\top X)^{-1}\underbrace{X^\top X(X^\top X)^{-1}}_{I}X^\top=H\). Trace: using \(\operatorname{tr}(AB)=\operatorname{tr}(BA)\), \(\operatorname{tr}(H)=\operatorname{tr}\!\big((X^\top X)^{-1}X^\top X\big)=\operatorname{tr}(I_p)=p\). \(\blacksquare\)

These three facts carry an enormous amount of downstream weight: the \(p\) in “degrees of freedom”, the \(n-p\) divisor in \(\hat\sigma^2\), the diagonal entries \(h_{ii}\) that define leverage, and the \(F\)-test’s numerator and denominator degrees of freedom are all consequences.

set.seed(19)
n <- 60; p <- 3
Xd <- cbind(1, matrix(rnorm(n * (p - 1)), n))
H  <- Xd %*% solve(crossprod(Xd)) %*% t(Xd)

c(symmetric   = all.equal(H, t(H)),
  idempotent  = all.equal(H %*% H, H),
  trace_H     = round(sum(diag(H)), 10),
  p           = p,
  trace_I_H   = round(sum(diag(diag(n) - H)), 10),
  n_minus_p   = n - p)
#>  symmetric idempotent    trace_H          p  trace_I_H  n_minus_p 
#>          1          1          3          3         57         57

The figure belo shows thaat projection is a geometric action. This one is interactive, the perpendicularity of the residual to the plane is the entire content of least squares, and it is only convincing when you can rotate it. Also, see the Matrices as Motions app.

set.seed(23)
# A plane through the origin spanned by two vectors in R^3
Xg <- cbind(c(1, 0, 0.3), c(0, 1, 0.5))
Hg <- Xg %*% solve(crossprod(Xg)) %*% t(Xg)

yv    <- c(0.9, 0.4, 2.2)     # a point off the plane
yhat  <- as.vector(Hg %*% yv) # its projection
resid <- yv - yhat

gs <- seq(-1.5, 1.5, length.out = 25)
grid <- expand.grid(a = gs, b = gs)
px <- matrix(Xg[1, 1] * grid$a + Xg[1, 2] * grid$b, 25)
py <- matrix(Xg[2, 1] * grid$a + Xg[2, 2] * grid$b, 25)
pz <- matrix(Xg[3, 1] * grid$a + Xg[3, 2] * grid$b, 25)

plot_ly() |>
  add_surface(x = px, y = py, z = pz, opacity = 0.45, showscale = FALSE,
              colorscale = list(c(0, "#BBD5EE"), c(1, "#BBD5EE"))) |>
  add_trace(x = c(0, yv[1]), y = c(0, yv[2]), z = c(0, yv[3]),
            type = "scatter3d", mode = "lines+markers", name = "y",
            line = list(width = 6, color = "black"),
            marker = list(size = 4, color = "black")) |>
  add_trace(x = c(0, yhat[1]), y = c(0, yhat[2]), z = c(0, yhat[3]),
            type = "scatter3d", mode = "lines+markers", name = "Hy (fitted)",
            line = list(width = 6, color = "#1F77B4"),
            marker = list(size = 4, color = "#1F77B4")) |>
  add_trace(x = c(yhat[1], yv[1]), y = c(yhat[2], yv[2]), z = c(yhat[3], yv[3]),
            type = "scatter3d", mode = "lines", name = "residual (I-H)y",
            line = list(width = 6, color = "firebrick", dash = "dash")) |>
  layout(title = "Least squares is orthogonal projection onto col(X)",
         scene = list(xaxis = list(title = "e1"), yaxis = list(title = "e2"),
                      zaxis = list(title = "e3"), aspectmode = "cube"))
c(residual_dot_col1 = round(sum(resid * Xg[, 1]), 12),
  residual_dot_col2 = round(sum(resid * Xg[, 2]), 12),
  pythagoras = all.equal(sum(yv^2), sum(yhat^2) + sum(resid^2)))
#> residual_dot_col1 residual_dot_col2        pythagoras 
#>                 0                 0                 1

The residual is numerically orthogonal to both spanning vectors, and \(\|\mathbf{y}\|^2=\|\hat{\mathbf{y}}\|^2+\|\mathbf{y}-\hat{\mathbf{y}}\|^2\). That Pythagorean identity, applied to centered data, is the ANOVA decomposition \[\underbrace{\mathrm{SST}}_{\text{total dispersion}} =\underbrace{\mathrm{SSR}}_{\text{total dispersion}}+ \underbrace{\mathrm{SSE}}_{\text{total dispersion}}\] (§3.10).

  • SST (Total Sum of Squares): \(\underbrace{\Vert{}\mathbf{Y}-\bar Y\mathbf{1}\Vert{}^2}_{\mathrm{SST}}\) captures the total amount of variability present in the dependent variable (\(\mathbf{Y}\)) around its overall mean (\(\bar{Y}\)), which is the baseline variance the model is trying to explain.
  • SSR (Regression Sum of Squares): \(\underbrace{\Vert{}\hat{\mathbf{Y}}-\bar Y\mathbf{1}\Vert{}^2}_{\mathrm{SSR}}\) is the amount of variability that is successfully explained by the model. Effectively, it calculates the distance between the model’s predicted values (\(\hat{\mathbf{Y}}\)) and the mean (\(\bar{Y}\)). A higher SSR means your model captures a large share of the pattern.
  • SSE (Error Sum of Squares): \(\underbrace{\Vert{}\mathbf{Y}-\hat{\mathbf{Y}}\Vert{}^2}_{\mathrm{SSE}}\) measures the unexplained variability or residuals, which represents the squared distance between the actual observed values (\(\mathbf{Y}\)) and what the model predicted (\(\hat{\mathbf{Y}}\)).

The Coefficient of Determination, \[R^2 = \frac{\mathrm{SSR}}{\mathrm{SST}} = 1 - \frac{\mathrm{SSE}}{\mathrm{SST}},\] represents the proportion of total variance in \(\mathbf{Y}\) that is explained by your predictors, ranging from \(0\) to \(1\) with a higher \(R^2\) indicating a better-fitting model.

The Adjusted \(R^2\) \[R^2_{\text{adj}} = 1 - \frac{\mathrm{SSE}/(n-p)}{\mathrm{SST}/(n-1)},\] which increases when adding more variables (even useless ones), penalizes for adding predictors. It divides SSE and SST by their respective degrees of freedom (\(n-p\) for error, where \(p\) is the number of parameters, and \(n-1\) for total), ensuring variables must genuinely improve the model to boost the score.

7 Inverses, pseudo-inverses, and solving systems

7.1 When does an inverse exist?

For a square \(A\in\mathbb{R}^{n\times n}\), a two-sided inverse \(A^{-1}\) satisfying \(AA^{-1}=A^{-1}A=I_n\) exists iff \(A\) is nonsingular, i.e. \(\operatorname{rank}(A)=n\), equivalently \(\det(A)\ne 0\), equivalently \(\mathcal{N}(A)=\{\mathbf{0}\}\). When it exists it is unique.

For a rectangular \(A\in\mathbb{R}^{m\times n}\) one-sided inverses may exist:

  • a left inverse \(L\) with \(LA=I_n\) exists iff \(\operatorname{rank}(A)=n\) (full column rank), and then \(L=(A^\top A)^{-1}A^\top\) works;
  • a right inverse \(R\) with \(AR=I_m\) exists iff \(\operatorname{rank}(A)=m\) (full row rank), and then \(R=A^\top(AA^\top)^{-1}\) works.

Common misconception: “\(A^{-1}\) means the same thing for rectangular matrices.” One-sided inverses are not unique when \(m\ne n\), there are infinitely many left inverses of a tall full-column-rank matrix. Writing \(A^{-1}\) for one of them invites the reader to cancel it on either side, which is invalid. Use \(A^{+}\) (the pseudo-inverse) when you need a canonical choice.

For \(2\times 2\) the closed form is worth memorizing:

\[A=\begin{pmatrix}a&b\\c&d\end{pmatrix},\qquad A^{-1}=\frac{1}{ad-bc}\begin{pmatrix}d&-b\\-c&a\end{pmatrix},\qquad ad-bc\ne 0 .\]

m10 <- matrix(1:4, nrow = 2)
m10
#>      [,1] [,2]
#> [1,]    1    3
#> [2,]    2    4
det(m10)
#> [1] -2
solve(m10)
#>      [,1] [,2]
#> [1,]   -2  1.5
#> [2,]    1 -0.5
round(m10 %*% solve(m10), 12)
#>      [,1] [,2]
#> [1,]    1    0
#> [2,]    0    1

7.2 The Moore–Penrose pseudo-inverse

For any \(A\in\mathbb{R}^{m\times n}\), singular or rectangular, there is a unique \(A^{+}\in\mathbb{R}^{n\times m}\) satisfying the four Penrose conditions

\[AA^{+}A=A,\qquad A^{+}AA^{+}=A^{+},\qquad (AA^{+})^\top=AA^{+},\qquad (A^{+}A)^\top=A^{+}A .\]

It is computed from the SVD (§3.6.5) as \(A^{+}=V\Sigma^{+}U^\top\), where \(\Sigma^{+}\) inverts the nonzero singular values and transposes. When \(A\) is invertible, \(A^{+}=A^{-1}\); when \(A\) has full column rank, \(A^{+}=(A^\top A)^{-1}A^\top\).

\(A^{+}\mathbf{b}\) is the minimum-norm least-squares solution: among all \(\mathbf{x}\) minimizing \(\|A\mathbf{x}-\mathbf{b}\|_2\), it is the one with smallest \(\|\mathbf{x}\|_2\). That property is what makes it the right default for rank-deficient problems.

Asing <- matrix(c(1, 2, 3,
                  2, 4, 6,
                  1, 1, 1), nrow = 3, byrow = TRUE)   # rank 2
Ap <- MASS::ginv(Asing)

c(penrose_1 = all.equal(Asing %*% Ap %*% Asing, Asing),
  penrose_2 = all.equal(Ap %*% Asing %*% Ap, Ap),
  penrose_3 = all.equal(Asing %*% Ap, t(Asing %*% Ap)),
  penrose_4 = all.equal(Ap %*% Asing, t(Ap %*% Asing)))
#> penrose_1 penrose_2 penrose_3 penrose_4 
#>      TRUE      TRUE      TRUE      TRUE
try(solve(Asing), silent = TRUE) |> class()   # solve() refuses; ginv() does not
#> [1] "try-error"

7.3 Solving \(A\mathbf{x}=\mathbf{b}\)

Consider the system

\[\begin{aligned} a + b + 2c &= 6\\ 3a - 2b + c &= 2\\ 2a + b - c &= 3 \end{aligned} \qquad\Longleftrightarrow\qquad \underbrace{\begin{pmatrix}1&1&2\\3&-2&1\\2&1&-1\end{pmatrix}}_{A} \underbrace{\begin{pmatrix}a\\b\\c\end{pmatrix}}_{\mathbf{x}} =\underbrace{\begin{pmatrix}6\\2\\3\end{pmatrix}}_{\mathbf{b}} .\]

This generalizes the scalar case exactly. To solve \(2x-3=5\), move the constant to the right (\(2x=8\)) and multiply by the reciprocal (\(x=4\)). The matrix version replaces “reciprocal” with “inverse”:

A <- matrix(c(1, 1, 2,
              3, -2, 1,
              2, 1, -1), nrow = 3, byrow = TRUE)
b <- c(6, 2, 3)

x_solve <- solve(A, b)      # PREFERRED: LU factorization, never forms A^{-1}
x_inv   <- solve(A) %*% b   # forms the explicit inverse first

rbind(`solve(A, b)` = x_solve, `solve(A) %*% b` = as.vector(x_inv))
#>                [,1] [,2] [,3]
#> solve(A, b)    1.35 1.75 1.45
#> solve(A) %*% b 1.35 1.75 1.45
round(A %*% x_solve - b, 12)
#>      [,1]
#> [1,]    0
#> [2,]    0
#> [3,]    0

Common misconception: “to solve \(A\mathbf{x}=\mathbf{b}\), compute \(A^{-1}\) and multiply.” Mathematically fine; computationally wrong. solve(A, b) performs an LU factorization and two triangular solves, about \(\tfrac{2}{3}n^3\) flops. solve(A) %*% b computes the full inverse (\(\approx 2n^3\) flops, three times the work), then multiplies. It is also less accurate: the explicit inverse accumulates rounding error in every one of its \(n^2\) entries, and the subsequent multiplication propagates all of them. The rule in numerical linear algebra is: never compute an inverse you only intend to multiply by something.

Lower-Upper (LU) factorization, aka LU decomposition, takes a matrix \(A\) and breaks it down into two special, easier-to-work-with matrices

  • \(L\) (Lower Triangular Matrix), where all entries above the main diagonal are zero.
  • \(U\) (Upper Triangular Matrix), where all entries below the main diagonal are zero.

Then, the original matrix equation \(A \cdot x = b\) is transformed into \[(L \cdot U) \cdot x = b\]

Instead of solving the complex system all at once, LU decomposition splits \(A\) into \(L\) and \(U\) supporting computational solving the problem in two simpler, highly efficient steps

  1. Forward Substitution (\(Ly = b\)): Because \(L\) is lower triangular, the computer can easily solve for a temporary vector \(y\) from the top down.
  2. Back Substitution (\(Ux = y\)): Because \(U\) is upper triangular, the computer can easily solve for your final answer \(x\) from the bottom up.

The main benefit is gain in speed, as doing full Gaussian elimination on a large matrix is computationally expensive. By factoring \(A\) into \(L\) and \(U\) once, we reuse that factorization to solve \(Ax = b\) very quickly for many different vectors \(b\).

set.seed(31)
n <- 700
Ab <- matrix(rnorm(n * n), n); bb <- rnorm(n)

t_solve <- system.time(x1 <- solve(Ab, bb))[["elapsed"]]
t_inv   <- system.time(x2 <- solve(Ab) %*% bb)[["elapsed"]]

c(solve_A_b = t_solve, invert_then_multiply = t_inv,
  speedup = round(t_inv / t_solve, 2),
  residual_solve = norm(Ab %*% x1 - bb, "2"),
  residual_inv   = norm(Ab %*% x2 - bb, "2"))
#>            solve_A_b invert_then_multiply              speedup 
#>          5.00000e-02          1.70000e-01          3.40000e+00 
#>       residual_solve         residual_inv 
#>          3.69212e-12          7.59754e-12

And a word on Cramer’s rule. It expresses \(x_j=\det(A_j)/\det(A)\), where \(A_j\) replaces column \(j\) of \(A\) by \(\mathbf{b}\). It is a beautiful theoretical result and a terrible algorithm: cofactor expansion of a determinant is \(O(n!)\), and even computing the \(n+1\) determinants by LU costs \(O(n^4)\) against \(O(n^3)\) for solving directly. It is also numerically unstable. Use it to prove things, never to compute them.

7.4 Reference: core matrix operations

Expression Meaning Cost
t(X) transpose \(O(mn)\) (often free, a view)
diag(x) diagonal matrix / extract diagonal / identity \(O(n)\) / \(O(n)\) / \(O(n^2)\)
A %*% B matrix product \(O(mnk)\)
crossprod(A, B) \(A^\top B\) without forming \(A^\top\) \(O(mnk)\), halved when \(A=B\)
tcrossprod(A, B) \(AB^\top\) \(O(mnk)\)
solve(A, b) solve \(A\mathbf{x}=\mathbf{b}\) via LU \(\tfrac{2}{3}n^3\)
solve(A) explicit inverse \(2n^3\)
qr.solve(A, b) least squares via QR \(2mn^2-\tfrac{2}{3}n^3\)
chol(A) Cholesky, \(A\) symmetric positive definite \(\tfrac{1}{3}n^3\)
eigen(A) eigendecomposition \(O(n^3)\), \(\approx 10n^3\) symmetric
svd(A) singular value decomposition \(O(mn^2)\), \(\approx 14mn^2\) thin
rowSums / colSums / rowMeans / colMeans margins \(O(mn)\)

8 Matrix decompositions

A decomposition rewrites \(A\) as a product of matrices with special structure, triangular, orthogonal, diagonal. The point is always the same: the factored form makes a hard problem easy.

8.1 LU

Gaussian elimination with partial pivoting factors \(PA=LU\) with \(P\) a permutation, \(L\) unit lower-triangular, \(U\) upper-triangular. Then \(A\mathbf{x}=\mathbf{b}\) becomes two triangular solves, each \(O(n^2)\):

\[L\mathbf{z}=P\mathbf{b}\ \ (\text{forward}),\qquad U\mathbf{x}=\mathbf{z}\ \ (\text{back}).\]

Factorization costs \(\tfrac{2}{3}n^3\); each subsequent right-hand side costs only \(2n^2\). This is why you factor once and solve many times.

lu <- Matrix::lu(Matrix::Matrix(A))
lu_exp <- Matrix::expand(lu)
round(as.matrix(lu_exp$L), 4)
#>        [,1]   [,2] [,3]
#> [1,] 1.0000 0.0000    0
#> [2,] 0.6667 1.0000    0
#> [3,] 0.3333 0.7143    1
round(as.matrix(lu_exp$U), 4)
#>      [,1]    [,2]    [,3]
#> [1,]    3 -2.0000  1.0000
#> [2,]    0  2.3333 -1.6667
#> [3,]    0  0.0000  2.8571
all.equal(as.matrix(lu_exp$P %*% lu_exp$L %*% lu_exp$U), A, check.attributes = FALSE)
#> [1] TRUE

8.2 Cholesky decomposition

If \(A\) is symmetric positive definite (\(\mathbf{x}^\top A\mathbf{x}>0\) for all \(\mathbf{x}\ne 0\)), then \(A=R^\top R\) with \(R\) upper-triangular and positive diagonal. Cholesky decomposition exploits symmetry to halve LU’s work: \(\tfrac{1}{3}n^3\) flops, no pivoting needed, and it is unconditionally stable.

Cholesky doubles as a positive-definiteness test: it fails exactly when the matrix is not positive definite. Every covariance matrix in this book is symmetric positive semi-definite, which is why chol() appears constantly in simulation and Gaussian likelihood code.

set.seed(37)
Z  <- matrix(rnorm(200 * 4), 200)
S  <- crossprod(Z) / 199        # a sample covariance: symmetric PSD
R  <- chol(S)
round(R, 4)
#>        [,1]   [,2]   [,3]    [,4]
#> [1,] 1.0258 0.0004 0.0135 -0.0586
#> [2,] 0.0000 1.0786 0.0152 -0.1092
#> [3,] 0.0000 0.0000 0.9738 -0.0830
#> [4,] 0.0000 0.0000 0.0000  0.9674
all.equal(crossprod(R), S)      # R^T R = S
#> [1] TRUE
# Simulating correlated Gaussians is one Cholesky away
Sigma <- matrix(c(4, 2, 2, 3), 2)
L <- t(chol(Sigma))
set.seed(41)
sim <- t(L %*% matrix(rnorm(2 * 10000), 2))
round(cov(sim), 3); Sigma
#>       [,1]  [,2]
#> [1,] 3.879 1.883
#> [2,] 1.883 2.887
#>      [,1] [,2]
#> [1,]    4    2
#> [2,]    2    3

8.3 QR decomposition

QR factorization (or QR decomposition) is a fundamental mathematical technique in linear algebra that breaks down any matrix \(A\) into two specific matrices:

\[A = QR,\] where

  • \(Q\) (Orthogonal Matrix) is a matrix whose columns are orthonormal (meaning they are perpendicular to each other and each have a length of \(1\)). Geometrically, multiplying by \(Q\) represents a rotation or reflection, which preserves lengths and angles.
  • \(R\) (Upper Triangular Matrix) is a matrix where all the entries below the main diagonal are zero. Multiplying by \(R\) represents scaling and shearing.

While LU factorization works for square matrices, QR factorization is especially powerful for non-square, rectangular matrices, \(n \ge p\), where there are more rows than columns.

Solving Least Squares Problems (Linear Regression): To find the line of best fit for a dataset with more data points (\(n\)) than variables (\(p\)), the standard normal equation approach (\(X^\top X b = X^\top y\)) can be numerically unstable. Using QR factorization (\(X = QR\)) makes solving the least squares problem numerically stable and robust.

Unlike simple methods, like standard Gram-Schmidt orthogonalization, proper QR algorithms, using Householder reflections or Givens rotations, avoid catastrophic round-off computational errors, ensuring the columns of \(Q\) remain strictly orthogonal. The iterative QR algorithm is one of the most famous methods used by software, e.g., NumPy, Julia, to find all the eigenvalues of a matrix simultaneously.

Any \(X\in\mathbb{R}^{n\times p}\) with \(n\ge p\) factors as \(X=QR\) with \(Q\in\mathbb{R}^{n\times p}\) having orthonormal columns (\(Q^\top Q=I_p\)) and \(R\in\mathbb{R}^{p\times p}\) upper-triangular. R uses Householder reflections, which cost \(2np^2-\tfrac{2}{3}p^3\) and are backward stable. (Classical Gram–Schmidt produces the same factors in exact arithmetic but loses orthogonality catastrophically in floating point; modified Gram–Schmidt is better but Householder is better still.)

QR is the workhorse of least squares. Substituting \(X=QR\):

\[\|\mathbf{y}-X\boldsymbol\beta\|_2^2=\|\mathbf{y}-QR\boldsymbol\beta\|_2^2 =\|Q^\top\mathbf{y}-R\boldsymbol\beta\|_2^2+\|(I-QQ^\top)\mathbf{y}\|_2^2 ,\]

because multiplying by an orthogonal matrix preserves \(\ell_2\) length. The second term does not involve \(\boldsymbol\beta\), so the minimizer solves the triangular system

\[R\hat{\boldsymbol\beta}=Q^\top\mathbf{y},\]

by back-substitution in \(O(p^2)\). Crucially, \(X^\top X\) is never formed.

set.seed(43)
Xq <- cbind(1, matrix(rnorm(100 * 3), 100))
yq <- Xq %*% c(2, -1, 0.5, 3) + rnorm(100)

qrX <- qr(Xq)
Q <- qr.Q(qrX); Rm <- qr.R(qrX)

c(orthonormal = all.equal(crossprod(Q), diag(ncol(Xq))),
  reconstructs = all.equal(Q %*% Rm, Xq, check.attributes = FALSE),
  R_upper_triangular = all(abs(Rm[lower.tri(Rm)]) < 1e-12))
#>        orthonormal       reconstructs R_upper_triangular 
#>               TRUE               TRUE               TRUE
beta_qr <- backsolve(Rm, crossprod(Q, yq))
as.vector(beta_qr)
#> [1]  1.921448 -0.949995  0.585666  3.017408
as.vector(coef(lm.fit(Xq, yq)))    # lm.fit does exactly this internally
#> [1]  1.921448 -0.949995  0.585666  3.017408

8.4 Eigen-spectra, Eigen-decomposition

A nonzero \(\mathbf{v}\) with \(A\mathbf{v}=\lambda\mathbf{v}\) is an eigenvector; \(\lambda\) is the corresponding eigenvalue. The eigenvectors are the directions the transformation leaves invariant up to scaling, the axes along which \(A\) merely stretches, compresses, or flips.

The historical thread runs from Euler’s work on rotational motion through Lagrange’s inertia matrices (principal axes) to Hilbert, who fixed the prefix eigen, meaning “own” or “characteristic”.

Collecting eigenpairs into \(V=[\mathbf{v}_1\cdots\mathbf{v}_n]\) and \(\Lambda=\operatorname{diag}(\lambda_1,\dots,\lambda_n)\) gives \(AV=V\Lambda\); if \(V\) is invertible, \(A=V\Lambda V^{-1}\).

Spectral theorem. If \(A\) is real symmetric, its eigenvalues are real, its eigenvectors can be chosen orthonormal, and \(A=Q\Lambda Q^\top\) with \(Q^\top Q=I\). Moreover \(A\) is positive definite iff every \(\lambda_j>0\).

Ae <- matrix(c(4, 1,
               2, 3), nrow = 2, byrow = TRUE)   # NOT symmetric
e <- eigen(Ae)
e$values
#> [1] 5 2
round(e$vectors, 4)
#>        [,1]    [,2]
#> [1,] 0.7071 -0.4472
#> [2,] 0.7071  0.8944
# The correct verification is A V = V Lambda  (equivalently A v_j = lambda_j v_j)
V <- e$vectors; L <- diag(e$values)
round(Ae %*% V - V %*% L, 12)
#>      [,1] [,2]
#> [1,]    0    0
#> [2,]    0    0
round(V %*% L %*% solve(V), 10)      # reconstruction A = V Lambda V^{-1}
#>      [,1] [,2]
#> [1,]    4    1
#> [2,]    2    3

Common misconception: “check eigenvectors with (values * diag(n) - A) %*% vectors.” That expression evaluates \((\Lambda-A)V\), which is not zero in general, the eigen-relation is \(AV=V\Lambda\), and \(\Lambda\) must multiply from the right. The two agree only when \(A\) is a scalar multiple of the identity, where every vector is an eigenvector and the check is vacuous. Verify one pair at a time, or use \(AV-V\Lambda\).

For a symmetric matrix the picture is cleaner and geometrically meaningful:

Sy <- matrix(c(4, 1.5, 1.5, 2), 2)
es <- eigen(Sy, symmetric = TRUE)
c(eigenvalues = es$values,
  orthonormal = all.equal(crossprod(es$vectors), diag(2)),
  reconstructs = all.equal(es$vectors %*% diag(es$values) %*% t(es$vectors), Sy),
  positive_definite = all(es$values > 0))
#>      eigenvalues1      eigenvalues2       orthonormal      reconstructs 
#>           4.80278           1.19722           1.00000           1.00000 
#> positive_definite 
#>           1.00000

A sample covariance matrix is symmetric positive semi-definite, so its eigendecomposition has a direct interpretation: the eigenvectors are the principal axes of the data ellipsoid and the eigenvalues are the variances along those axes. That is the whole content of principal component analysis (Chapter 4), and it is best seen in three dimensions, rotating:

mlb <- dspa_read(
  "https://umich.instructure.com/files/330381/download?download_frd=1",
  "01a_data.txt", reader = read.table, header = TRUE)

M <- as.matrix(mlb[, c("Height", "Weight", "Age")])
ctr <- colMeans(M); Sg <- cov(M)
eg <- eigen(Sg, symmetric = TRUE)

# 1-sd ellipsoid  {x : (x-mu)' S^{-1} (x-mu) = 1}  = mu + Q L^{1/2} u,  |u| = 1
th <- seq(0, 2 * pi, length.out = 60)
ph <- seq(0, pi, length.out = 30)
U <- rbind(as.vector(outer(cos(th), sin(ph))),
           as.vector(outer(sin(th), sin(ph))),
           as.vector(outer(rep(1, 60), cos(ph))))
E <- eg$vectors %*% diag(sqrt(eg$values)) %*% U + ctr

ex <- matrix(E[1, ], 60); ey <- matrix(E[2, ], 60); ez <- matrix(E[3, ], 60)

p3 <- plot_ly() |>
  add_surface(x = ex, y = ey, z = ez, opacity = 0.30, showscale = FALSE,
              colorscale = list(c(0, "#9EC5E8"), c(1, "#9EC5E8"))) |>
  add_trace(x = M[, 1], y = M[, 2], z = M[, 3], type = "scatter3d",
            mode = "markers", name = "MLB players",
            marker = list(size = 2, opacity = 0.35, color = "grey30"))

for (j in 1:3) {
  ax <- ctr + eg$vectors[, j] * sqrt(eg$values[j])
  p3 <- add_trace(p3, x = c(ctr[1], ax[1]), y = c(ctr[2], ax[2]),
                  z = c(ctr[3], ax[3]), type = "scatter3d", mode = "lines",
                  name = sprintf("PC%d (var %.1f)", j, eg$values[j]),
                  line = list(width = 8))
}

p3 |> layout(title = "Covariance ellipsoid and its principal axes (eigenvectors)",
             scene = list(xaxis = list(title = "Height (in)"),
                          yaxis = list(title = "Weight (lb)"),
                          zaxis = list(title = "Age (yr)")))
data.frame(component = paste0("PC", 1:3),
           eigenvalue = round(eg$values, 3),
           prop_variance = round(eg$values / sum(eg$values), 4),
           cumulative = round(cumsum(eg$values) / sum(eg$values), 4))

8.5 Singular value decomposition (SVD)

The SVD is the most general and most useful decomposition. Every matrix \(A\in\mathbb{R}^{m\times n}\) factors as

\[A=U\Sigma V^\top,\qquad U^\top U=I,\quad V^\top V=I,\quad \Sigma=\operatorname{diag}(\sigma_1\ge\sigma_2\ge\cdots\ge\sigma_r>0).\]

No squareness, no invertibility, no symmetry required. Equivalently,

\[A=\sum_{j=1}^{r}\sigma_j\,\mathbf{u}_j\mathbf{v}_j^\top\]

a sum of \(r\) rank-one outer products (§3.2.3), each weighted by its singular value.

Relationships worth knowing: \(\sigma_j=\sqrt{\lambda_j(A^\top A)}\); the columns of \(V\) are eigenvectors of \(A^\top A\) and of \(U\) eigenvectors of \(AA^\top\); \(\operatorname{rank}(A)=\#\{\sigma_j>0\}\); \(\|A\|_2=\sigma_1\); and the first \(r\) columns of \(U\) span \(\operatorname{col}(A)\) while the last \(n-r\) columns of \(V\) span \(\mathcal{N}(A)\), the SVD hands you all four fundamental subspaces at once.

Eckart–Young–Mirsky theorem. Among all matrices of rank at most \(k\), the truncation \(A_k=\sum_{j\le k}\sigma_j\mathbf{u}_j\mathbf{v}_j^\top\) minimizes both \(\|A-A_k\|_2=\sigma_{k+1}\) and \(\|A-A_k\|_F=\big(\sum_{j>k}\sigma_j^2\big)^{1/2}\).

This is the theoretical basis of PCA, latent semantic analysis, image compression, recommender systems, and denoising. Cost: \(\approx 14mn^2\) for the thin SVD when \(m\ge n\).

set.seed(47)
Am <- matrix(rnorm(6 * 4), 6)
s <- svd(Am)

c(reconstructs = all.equal(s$u %*% diag(s$d) %*% t(s$v), Am),
  U_orthonormal = all.equal(crossprod(s$u), diag(4)),
  V_orthonormal = all.equal(crossprod(s$v), diag(4)),
  spectral_norm = all.equal(max(s$d), norm(Am, "2")),
  sigma_vs_eigen = all.equal(s$d^2, eigen(crossprod(Am))$values))
#>   reconstructs  U_orthonormal  V_orthonormal  spectral_norm sigma_vs_eigen 
#>           TRUE           TRUE           TRUE           TRUE           TRUE

Low-rank approximation, demonstrated on real high-dimensional data. The GSE5859 microarray study measured expression for 8,793 genes across 208 lymphoblastoid samples from three HapMap populations (CEU, CHB, JPT).

gene <- dspa_read("https://umich.instructure.com/files/2001417/download?download_frd=1",
                  "exprs_GSE5859.csv", header = TRUE)
info <- dspa_read("https://umich.instructure.com/files/2001418/download?download_frd=1",
                  "exprs_MetaData_GSE5859.csv", header = TRUE)
dim(gene); dim(info)
#> [1] 8793  209
#> [1] 208   4
G <- as.matrix(gene[1:200, 2:101])       # 200 genes x 100 samples
G <- sweep(G, 1, rowMeans(G))            # centre each gene
sg <- svd(G)

energy <- cumsum(sg$d^2) / sum(sg$d^2)
ranks  <- c(1, 2, 5, 10, 25, 50)

approx_err <- sapply(ranks, function(k) {
  Gk <- sg$u[, 1:k, drop = FALSE] %*% diag(sg$d[1:k], k, k) %*% t(sg$v[, 1:k, drop = FALSE])
  c(frobenius_rel = norm(G - Gk, "F") / norm(G, "F"),
    predicted     = sqrt(sum(sg$d[-(1:k)]^2)) / sqrt(sum(sg$d^2)))
})
data.frame(rank = ranks, t(round(approx_err, 5)),
           variance_captured = round(energy[ranks], 4))

Observed and Eckart–Young-predicted errors agree to five decimals, the theorem is not an approximation, it is exact.

scree <- data.frame(k = 1:40, sigma = sg$d[1:40], cum = energy[1:40])

p_a <- ggplot(scree, aes(k, sigma)) +
  geom_col(fill = "steelblue") +
  labs(title = "Singular value spectrum", x = "index j",
       y = expression(sigma[j])) + theme_dspa(10)

p_b <- ggplot(scree, aes(k, cum)) +
  geom_line(linewidth = 0.9, colour = "firebrick") + geom_point(size = 1.3) +
  geom_hline(yintercept = 0.9, linetype = "dashed", colour = "grey40") +
  scale_y_continuous(labels = scales::percent) +
  labs(title = "Cumulative variance captured", x = "rank k", y = NULL) +
  theme_dspa(10)

p_a | p_b

# --- Interactive equivalent ------------------------------------------------
plot_ly(x = 1:40, y = sg$d[1:40], type = "bar", name = "singular values") |>
  add_lines(x = 1:40, y = energy[1:40] * max(sg$d), yaxis = "y2",
            name = "cumulative variance") |>
  layout(title = "Singular value spectrum (GSE5859 submatrix)",
         xaxis = list(title = "index j"),
         yaxis = list(title = "sigma_j"),
         yaxis2 = list(overlaying = "y", side = "right", tickformat = ".0%",
                       range = c(0, 1), title = "cumulative"),
         legend = list(orientation = "h"))

The reconstruction quality is a surface over (gene, sample, expression), so this one is interactive:

mk <- function(k) {
  sg$u[, 1:k, drop = FALSE] %*% diag(sg$d[1:k], k, k) %*% t(sg$v[, 1:k, drop = FALSE])
}
plot_ly() |>
  add_surface(z = G, showscale = FALSE, opacity = 0.95,
              colorscale = "Viridis", name = "full rank") |>
  add_surface(z = mk(3) + 6, showscale = FALSE, opacity = 0.95,
              colorscale = "Viridis", name = "rank 3 (offset +6)") |>
  layout(title = "Gene-expression matrix: full rank (bottom) vs. rank-3 approximation (top)",
         scene = list(xaxis = list(title = "Sample"),
                      yaxis = list(title = "Gene"),
                      zaxis = list(title = "Centred expression")))