| SOCR ≫ | DSPA ≫ | DSPA3 Topics ≫ |
DSPA3: Data Science and Predictive Analytics (UMich HS650)
Chapter 4: Linear and Nonlinear Dimensionality Reduction
Chapter 4: Linear and Nonlinear Dimensionality Reduction
SOCR/MCAIM (Ivo Dinov)
SOCR/MCAIM (Ivo Dinov)
September 04, 2026
library(ggplot2)
library(dplyr)
library(tidyr)
library(patchwork)
library(plotly) # interactive figures and ALL 3-D graphics
library(MASS)How this chapter uses graphics
Every two-dimensional figure is drawn with
ggplot2and rendered statically. Immediately after each one, the equivalentplot_ly()code appears in a chunk markedeval=FALSE, echo=TRUE, printed 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. This chapter is unusually rich in them, and for good reason: the whole subject is about geometry. A Swiss roll that cannot be rotated is not a Swiss roll, and a 3-D embedding viewed from one fixed angle discards the dimension it was computed to reveal.
1 Learning objectives
After completing this chapter you will be able to:
- Explain why high-dimensional data require reduction, using distance concentration and the notion of intrinsic dimension.
- State and verify the Johnson–Lindenstrauss lemma, and use it to bound how far a dataset can be compressed while preserving pairwise distances.
- Derive PCA two ways, variance maximization and reconstruction-error minimization, and prove they give the same answer.
- Distinguish the first principal component from the least-squares regression line, and explain when they coincide.
- Relate PCA, SVD, and classical MDS exactly, and choose among them on computational grounds.
- Choose the number of components by parallel analysis, broken-stick, or cross-validated reconstruction error rather than by eyeballing a scree plot.
- State the ICA model, its independence assumption, and the two ambiguities plus the one-Gaussian restriction that limit identifiability.
- Contrast factor analysis with PCA as a generative model with a testable fit, and interpret communalities, uniquenesses, and rotations.
- Apply kernel PCA, t-SNE, and UMAP to data on a curved manifold, and state the computational cost of each.
- Quantify embedding quality with trustworthiness and continuity, and list what must not be read off a t-SNE or UMAP plot.
Estimated time: 9–12 hours including exercises. Prerequisites: Chapter 3, especially the SVD (§3.6.5), the spectral theorem (§3.6.4), and Eckart–Young. This chapter is their direct application.
1.1 Chapter 4 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: WHY REDUCE DIMENSIONS
3 Motivating example: two dimensions into one
Suppose we observe the standardized heights of \(n=1{,}000\) pairs of identical twins, each pair a point in \(\mathbb{R}^2\). Twin heights are strongly correlated, so the cloud is a thin ellipse tilted along the diagonal. Such data appear in pediatric neuroimaging-genetics twin studies; we simulate them from a bivariate normal.
\[\begin{pmatrix}Y_1\\Y_2\end{pmatrix}\sim \mathrm{BVN}\!\left(\begin{pmatrix}0\\0\end{pmatrix},\ \begin{pmatrix}1&0.95\\0.95&1\end{pmatrix}\right).\]
set.seed(1234)
n <- 1000
Sigma <- matrix(c(1, 0.95, 0.95, 1), 2, 2)
Y <- MASS::mvrnorm(n, mu = c(0, 0), Sigma = Sigma) # n x 2: rows are twin pairs
colnames(Y) <- c("twin1", "twin2")
head(round(Y, 3), 3)#> twin1 twin2
#> [1,] -1.001 -1.382
#> [2,] 0.226 0.322
#> [3,] 1.314 0.827On orientation. Throughout this chapter, rows
are cases and columns are features, the standard convention in
statistics, and the one prcomp(), dist(), and
cov() all assume. Genomics often transposes this (genes in
rows, samples in columns), so when moving between fields, check which
margin your covariance is computed over: cov(X) gives a
feature × feature matrix from a case × feature
input.
twins <- as.data.frame(Y)
pair_ids <- c(1, 2)
ggplot(twins, aes(twin1, twin2)) +
geom_point(alpha = 0.22, size = 1.1, colour = "grey30") +
geom_point(data = twins[pair_ids, ], aes(colour = factor(pair_ids)),
size = 4.5) +
geom_segment(x = twins$twin1[1], y = twins$twin2[1],
xend = twins$twin1[2], yend = twins$twin2[2],
linetype = "dashed", colour = "grey20") +
scale_colour_manual(values = c("1" = "firebrick", "2" = "darkgreen"),
labels = c("Twin-pair 1", "Twin-pair 2")) +
coord_fixed(xlim = c(-3, 3), ylim = c(-3, 3)) +
labs(title = "Simulated twin heights",
subtitle = sprintf("Correlation %.2f; the cloud is a thin ellipse along the diagonal",
cor(Y[, 1], Y[, 2])),
x = "Twin 1 (standardized height)", y = "Twin 2 (standardized height)",
colour = NULL) +
theme_dspa()
# --- Interactive equivalent ------------------------------------------------
plot_ly() |>
add_markers(x = ~Y[, 1], y = ~Y[, 2], name = "Data scatter",
marker = list(opacity = 0.4)) |>
add_markers(x = Y[1, 1], y = Y[1, 2], name = "Twin-pair 1",
marker = list(color = "red", size = 20,
line = list(color = "yellow", width = 2))) |>
add_markers(x = Y[2, 1], y = Y[2, 2], name = "Twin-pair 2",
marker = list(color = "green", size = 20,
line = list(color = "orange", width = 2))) |>
layout(title = "Scatter plot of simulated twin data",
xaxis = list(title = "Twin 1 (standardized height)"),
yaxis = list(title = "Twin 2 (standardized height)"),
legend = list(orientation = "h"))The quantity we most want to preserve when reducing dimensions is
pairwise distance. stats::dist() computes
distances between the rows of a matrix:
#> [1] 2.100193.1 A transformation that is not a rotation
A natural summary of a twin pair is its average and its difference:
\[z_1=\frac{y_1+y_2}{2},\qquad z_2=y_1-y_2 \qquad\Longleftrightarrow\qquad \mathbf{z}=M\mathbf{y},\quad M=\begin{pmatrix}1/2&1/2\\1&-1\end{pmatrix}.\]
This is a linear map (no translation, so not merely affine), and it is invertible:
\[M^{-1}=\begin{pmatrix}1&1/2\\1&-1/2\end{pmatrix}.\]
#> [,1] [,2]
#> [1,] 0.5 0.5
#> [2,] 1.0 -1.0#> [,1] [,2]
#> [1,] 1 0.5
#> [2,] 1 -0.5#> [,1] [,2]
#> [1,] 1 0
#> [2,] 0 1It does not preserve distances:
D_M <- dist(Z_M)
set.seed(7); idx <- sample(length(D), 4000)
cmp <- data.frame(original = as.numeric(D)[idx], transformed = as.numeric(D_M)[idx])
ggplot(cmp, aes(original, transformed)) +
geom_point(alpha = 0.12, size = 0.7, colour = "steelblue") +
geom_abline(slope = 1, intercept = 0, colour = "firebrick", linewidth = 1) +
coord_fixed(xlim = c(0, 8), ylim = c(0, 8)) +
labs(title = "The average/difference map is not an isometry",
subtitle = "Points off the red line are pairs whose separation changed",
x = "Original distance", y = "Transformed distance") +
theme_dspa()
# --- Interactive equivalent ------------------------------------------------
plot_ly() |>
add_markers(x = ~as.numeric(D)[1:5000], y = ~as.numeric(D_M)[1:5000],
name = "Transformed twin distances", marker = list(opacity = 0.3)) |>
add_trace(x = c(0, 8), y = c(0, 8), mode = "lines", type = "scatter",
line = list(color = "red", width = 4), name = "Preserved distances") |>
layout(title = "Preservation of distances between twins (M transform)",
xaxis = list(title = "Original twin distances", range = c(0, 8)),
yaxis = list(title = "Transformed twin distances", range = c(0, 8)),
legend = list(orientation = "h"))A concrete counterexample makes the failure unambiguous. The points \(v_1=(0,1)^\top\) and \(v_2=(1,0)^\top\) are \(\sqrt2\approx1.414\) apart; after \(M\) they are 2 apart.
v1 <- c(0, 1); v2 <- c(1, 0)
euc <- function(a, b) sqrt(sum((a - b)^2))
c(before = euc(v1, v2),
after = euc(M %*% v1, M %*% v2),
M_orthogonal = isTRUE(all.equal(t(M) %*% M, diag(2))))#> before after M_orthogonal
#> 1.41421 2.00000 0.000004 Isometries: what a rotation preserves
A map between metric spaces that preserves distances is an isometry. For a linear map \(\mathbf{z}=A\mathbf{y}\) the condition is easy to characterize. Writing \(\mathbf{p}=\mathbf{y}_i-\mathbf{y}_j\),
\[d^2(A\mathbf{y}_i,A\mathbf{y}_j)=\|A\mathbf{p}\|^2=(A\mathbf{p})^\top(A\mathbf{p})=\mathbf{p}^\top A^\top A\,\mathbf{p},\]
which equals \(\mathbf{p}^\top\mathbf{p}=d^2(\mathbf{y}_i,\mathbf{y}_j)\) for every \(\mathbf{p}\) if and only if
\[\boxed{\;A^\top A=I\;}\]
i.e. \(A\) is orthogonal. Orthogonal matrices with \(\det A=+1\) are rotations; those with \(\det A=-1\) include a reflection. Either way, lengths and angles are preserved, which is exactly why PCA, whose loadings matrix is orthogonal, can rotate a dataset without distorting it.
Take \[A=\frac{1}{\sqrt2}\begin{pmatrix}1&1\\1&-1\end{pmatrix}.\]
A <- (1 / sqrt(2)) * matrix(c(1, 1, 1, -1), 2, 2)
c(is_orthogonal = isTRUE(all.equal(t(A) %*% A, diag(2))),
determinant = det(A))#> is_orthogonal determinant
#> 1 -1Z <- Y %*% t(A)
D_A <- dist(Z)
c(max_abs_distance_change = max(abs(as.numeric(D) - as.numeric(D_A))))#> max_abs_distance_change
#> 3.55271e-15cmpA <- data.frame(original = as.numeric(D)[idx], rotated = as.numeric(D_A)[idx])
p_left <- ggplot(twins, aes(twin1, twin2)) +
geom_point(alpha = 0.2, size = 0.9, colour = "grey35") +
geom_point(data = twins[pair_ids, ], colour = c("firebrick", "darkgreen"), size = 3.5) +
coord_fixed(xlim = c(-4, 4), ylim = c(-4, 4)) +
labs(title = "Original", x = "Twin 1", y = "Twin 2") + theme_dspa(10)
p_right <- ggplot(as.data.frame(Z), aes(V1, V2)) +
geom_point(alpha = 0.2, size = 0.9, colour = "grey35") +
geom_point(data = as.data.frame(Z)[pair_ids, ],
colour = c("firebrick", "darkgreen"), size = 3.5) +
coord_fixed(xlim = c(-4, 4), ylim = c(-4, 4)) +
labs(title = "After the rotation A",
subtitle = "Same shape, new axes", x = "Rotated axis 1", y = "Rotated axis 2") +
theme_dspa(10)
p_left | p_right
# --- Interactive equivalent ------------------------------------------------
plot_ly() |>
add_markers(x = ~as.numeric(D)[1:5000], y = ~as.numeric(D_A)[1:5000],
name = "Rotated twin distances", marker = list(opacity = 0.3)) |>
add_trace(x = c(0, 8), y = c(0, 8), mode = "lines", type = "scatter",
line = list(color = "red", width = 4), name = "Preserved distances") |>
layout(title = "Preservation of distances between twins (rotation A)",
xaxis = list(title = "Original twin distances", range = c(0, 8)),
yaxis = list(title = "Rotated twin distances", range = c(0, 8)),
legend = list(orientation = "h"))4.1 What the rotation does to the covariance
For \(\mathbf{Y}\sim\mathrm{BVN}(\boldsymbol\mu,\Sigma)\) and \(\mathbf{Z}=A\mathbf{Y}+\boldsymbol\eta\),
\[\mathbf{Z}\sim\mathrm{BVN}\big(\boldsymbol\eta+A\boldsymbol\mu,\ A\Sigma A^\top\big),\]
by the change-of-variables theorem. Affine maps preserve normality; only orthogonal ones preserve distance.
#> [,1] [,2]
#> [1,] 1.95 0.00
#> [2,] 0.00 0.05#> [1] 1.95 0.05\(A\Sigma A^\top=\operatorname{diag}(1.95,\ 0.05)\), and those two numbers are exactly the eigenvalues of \(\Sigma\), since \(\Sigma\)’s eigenvectors are \((1,1)/\sqrt2\) and \((1,-1)/\sqrt2\), which are the rows of \(A\). This rotation diagonalizes the covariance: after it, the two coordinates are uncorrelated, and 97.5% of the total variance (\(1.95/2\)) sits on the first axis alone.
That is PCA, discovered by hand. Everything in §4.4 generalizes this to \(p\) dimensions and finds the diagonalizing rotation automatically.
Uncorrelated and independent are not synonyms. For a jointly Gaussian vector, zero correlation does imply independence, which is why the rotated twin coordinates can be simulated separately below. But the implication fails for variables that are merely marginally normal: \(X\sim N(0,1)\) and \(Y=X\cdot S\) with \(S=\pm1\) at random are each normal and uncorrelated, yet \(|Y|=|X|\) always. This gap is precisely what independent component analysis (§4.9) exploits.
set.seed(2017)
zz <- cbind(rnorm(n, 0, sd = sqrt(1.95)), rnorm(n, 0, sd = sqrt(0.05)))
c(sd_axis1 = sd(zz[, 1]), sd_axis2 = sd(zz[, 2]),
correlation = cor(zz[, 1], zz[, 2]))#> sd_axis1 sd_axis2 correlation
#> 1.372669122 0.223184804 0.0003735314.2 Dropping the second coordinate
Because the rotation concentrated 97.5% of the variance on the first axis, we can discard the second and keep almost all the geometry:
D_1d <- dist(Z[, 1, drop = FALSE])
drop_cmp <- data.frame(full = as.numeric(D)[idx], one_d = as.numeric(D_1d)[idx])
ggplot(drop_cmp, aes(full, one_d)) +
geom_point(alpha = 0.10, size = 0.7, colour = "steelblue") +
geom_abline(slope = 1, intercept = 0, colour = "firebrick", linewidth = 1) +
coord_fixed(xlim = c(0, 8), ylim = c(0, 8)) +
labs(title = "One dimension approximates two",
subtitle = sprintf("Correlation between 1-D and 2-D distances: %.4f; mean relative error %.1f%%",
cor(as.numeric(D), as.numeric(D_1d)),
100 * mean(abs(as.numeric(D_1d) - as.numeric(D)) / as.numeric(D))),
x = "Distance in the original 2-D space",
y = "Distance using only the first rotated axis") +
theme_dspa()
# --- Interactive equivalent ------------------------------------------------
plot_ly() |>
add_markers(x = ~as.numeric(D)[idx], y = ~as.numeric(D_1d)[idx],
name = "Transformed distances", marker = list(opacity = 0.25)) |>
add_trace(x = c(0, 8), y = c(0, 8), mode = "lines", type = "scatter",
line = list(color = "red", width = 4), name = "Preserved distances") |>
layout(title = "Approximate distance preservation in 1-D",
xaxis = list(title = "Original distances", range = c(0, 8)),
yaxis = list(title = "1-D distances", range = c(0, 8)),
legend = list(orientation = "h"))The distances collapse onto the identity line with a small systematic downward bias, reduction can only shorten distances, never lengthen them, since we are projecting. The first rotated coordinate is the first principal component.
5 Why dimension is a problem
Reduction is not merely convenient. In high dimensions, geometry itself misbehaves.
5.1 Distance concentration
Let \(X_1,\dots,X_n\) be IID uniform on the unit cube \([0,1]^d\). As \(d\) grows, the ratio of the spread of pairwise distances to their mean tends to zero:
\[\frac{\operatorname{sd}\big(\|X_i-X_j\|\big)}{E\big(\|X_i-X_j\|\big)}\;\longrightarrow\;0 .\]
Every point becomes almost equidistant from every other. Nearest-neighbour methods, kernel smoothers, clustering, and anomaly detection all rest on the premise that “near” is meaningfully different from “far”, a premise that evaporates.
conc <- function(d, n = 300, seed = 11) {
set.seed(seed + d)
X <- matrix(runif(n * d), n, d)
dd <- as.numeric(dist(X))
c(d = d, mean = mean(dd), sd = sd(dd),
relative_spread = sd(dd) / mean(dd),
contrast = (max(dd) - min(dd)) / min(dd))
}
tab <- do.call(rbind, lapply(c(2, 5, 10, 50, 100, 500, 1000), conc))
round(as.data.frame(tab), 4)#> d mean sd relative_spread contrast
#> 1 2 0.5339 0.2518 0.4716 397.9030
#> 2 5 0.8710 0.2482 0.2849 23.5588
#> 3 10 1.2580 0.2423 0.1926 5.9845
#> 4 50 2.8892 0.2427 0.0840 1.1748
#> 5 100 4.0926 0.2458 0.0601 0.6443
#> 6 500 9.1033 0.2430 0.0267 0.2340
#> 7 1000 12.9049 0.2424 0.0188 0.1588as.data.frame(tab) |>
ggplot(aes(d, relative_spread)) +
geom_line(linewidth = 1, colour = "firebrick") +
geom_point(size = 2.4) +
scale_x_log10() +
labs(title = "Distances concentrate as dimension grows",
subtitle = "Ratio of the standard deviation to the mean of all pairwise distances, uniform data on the unit cube",
x = "Dimension d (log scale)", y = "sd(distance) / mean(distance)") +
theme_dspa()
The effect is a surface over \((n,d)\), so it is worth rotating:
ns <- c(50, 100, 200, 400, 800)
ds <- c(2, 5, 10, 25, 50, 100, 250, 500)
Zc <- outer(ns, ds, Vectorize(function(nn, dd) {
set.seed(nn + dd)
X <- matrix(runif(nn * dd), nn, dd)
v <- as.numeric(dist(X)); sd(v) / mean(v)
}))
plot_ly(x = ds, y = ns, z = Zc, type = "surface", colorscale = "Viridis") |>
layout(title = "Distance concentration as a function of sample size and dimension",
scene = list(xaxis = list(title = "Dimension d", type = "log"),
yaxis = list(title = "Sample size n"),
zaxis = list(title = "sd / mean of distances")))Notice that the surface falls steeply along \(d\) and is nearly flat along \(n\). More data does not fix high dimension. Only reducing \(d\), or exploiting structure within it, does.
5.2 Intrinsic dimension
The saving grace is that real data rarely fill their ambient space. A \(28\times28\) image lives in \(\mathbb{R}^{784}\), but the set of images that look like handwritten digits is a far lower-dimensional subset, a manifold parameterized by stroke thickness, slant, curvature, and a handful of other factors.
The intrinsic dimension \(d_{\text{int}}\) is the dimension of that underlying set. A cheap estimator counts how the number of neighbours within radius \(r\) grows: on a \(d_{\text{int}}\)-dimensional manifold, \(N(r)\propto r^{d_{\text{int}}}\), so
\[\hat d_{\text{int}}=\frac{d\log N(r)}{d\log r}.\]
corr_dim <- function(X, q = c(0.02, 0.15)) {
dd <- as.numeric(dist(X))
rs <- quantile(dd, seq(q[1], q[2], length.out = 25))
Nr <- sapply(rs, \(r) mean(dd < r))
keep <- Nr > 0
unname(coef(lm(log(Nr[keep]) ~ log(rs[keep])))[2])
}
set.seed(19)
# A 2-D plane and a 1-D helix, both embedded in R^10
plane10 <- cbind(matrix(runif(500 * 2), 500, 2), matrix(0, 500, 8)) +
matrix(rnorm(500 * 10, sd = 0.01), 500, 10)
t_par <- runif(500, 0, 6 * pi)
helix10 <- cbind(cos(t_par), sin(t_par), t_par / 6, matrix(0, 500, 7)) +
matrix(rnorm(500 * 10, sd = 0.01), 500, 10)
cube10 <- matrix(runif(500 * 10), 500, 10)
c(plane_in_R10 = corr_dim(plane10),
helix_in_R10 = corr_dim(helix10),
uniform_cube_R10 = corr_dim(cube10))#> plane_in_R10 helix_in_R10 uniform_cube_R10
#> 1.97766 1.12203 6.91864All three datasets sit in \(\mathbb{R}^{10}\); their intrinsic dimensions are approximately 2, 1, and 10. Dimensionality reduction works exactly when \(d_{\text{int}}\ll d\), and its job is to find coordinates on that lower-dimensional set.
5.3 How far can you compress? Johnson–Lindenstrauss
There is a theorem that answers this, and its answer is surprising: the achievable target dimension does not depend on \(d\) at all.
Johnson–Lindenstrauss lemma (1984). For any \(0<\epsilon<1\) and any set of \(n\) points in \(\mathbb{R}^d\), if \[k\;\ge\;\frac{8\ln n}{\epsilon^2}\] then there exists a linear map \(f:\mathbb{R}^d\to\mathbb{R}^k\) such that for all pairs \(i,j\), \[(1-\epsilon)\|x_i-x_j\|^2\;\le\;\|f(x_i)-f(x_j)\|^2\;\le\;(1+\epsilon)\|x_i-x_j\|^2 .\]
Two consequences worth absorbing. The bound is logarithmic in \(n\) and independent of \(d\): a million points in a million dimensions can be compressed to a few thousand with 10% distortion. And the map is easy to construct, a random Gaussian matrix scaled by \(1/\sqrt k\) works with high probability. No data-dependent optimization required.
set.seed(23)
n_jl <- 500; d_jl <- 2000
Xjl <- matrix(rnorm(n_jl * d_jl), n_jl, d_jl)
D_orig <- as.numeric(dist(Xjl))
jl_test <- function(k) {
set.seed(100 + k)
R <- matrix(rnorm(d_jl * k), d_jl, k) / sqrt(k) # random projection
ratio <- (as.numeric(dist(Xjl %*% R)) / D_orig)^2
c(k = k, theoretical_eps = sqrt(8 * log(n_jl) / k),
observed_max_distortion = max(abs(ratio - 1)),
mean_ratio = mean(ratio))
}
do.call(rbind, lapply(c(50, 100, 250, 500, 1000), jl_test)) |> round(4)#> k theoretical_eps observed_max_distortion mean_ratio
#> [1,] 50 0.9972 1.2391 1.0067
#> [2,] 100 0.7051 0.7025 1.0046
#> [3,] 250 0.4459 0.4826 0.9991
#> [4,] 500 0.3153 0.2883 0.9991
#> [5,] 1000 0.2230 0.2117 1.0010set.seed(150)
Rk <- matrix(rnorm(d_jl * 150), d_jl, 150) / sqrt(150)
jl_df <- data.frame(original = D_orig, projected = as.numeric(dist(Xjl %*% Rk)))
ggplot(jl_df[sample(nrow(jl_df), 6000), ], aes(original, projected)) +
geom_point(alpha = 0.10, size = 0.6, colour = "steelblue") +
geom_abline(slope = 1, intercept = 0, colour = "firebrick", linewidth = 1) +
labs(title = "A random projection from 2,000 dimensions to 150",
subtitle = "No data was examined to build this map, yet distances are nearly preserved",
x = "Distance in R^2000", y = "Distance in R^150") +
theme_dspa()
The observed distortion is well inside the theoretical bound, because the lemma is a worst-case guarantee. Random projection is the cheapest reduction available, \(O(ndk)\) with no decomposition, and it is the right first move when \(d\) is enormous and you only need distances. What it does not give you is interpretable axes, which is why the rest of this chapter exists.
6 PART II: LINEAR METHODS
7 Principal component analysis
PCA finds an orthogonal rotation of the coordinate axes such that the first new axis carries the most variance, the second the most of what remains subject to being orthogonal to the first, and so on. Two apparently different objectives produce the same answer.
7.1 Derivation 1: maximize variance
Let \(X\) be \(n\times p\) with column-centred data, and let \(S=\frac{1}{n-1}X^\top X\) be the sample covariance. A linear combination \(\mathbf{z}=X\mathbf{a}\) has sample variance
\[\operatorname{Var}(\mathbf{z})=\mathbf{a}^\top S\,\mathbf{a}.\]
This is unbounded unless we constrain \(\mathbf{a}\), so require \(\mathbf{a}^\top\mathbf{a}=1\) and maximize the Lagrangian
\[\mathcal{L}(\mathbf{a},\lambda)=\mathbf{a}^\top S\mathbf{a}-\lambda\big(\mathbf{a}^\top\mathbf{a}-1\big).\]
Setting \(\nabla_{\mathbf{a}}\mathcal{L}=2S\mathbf{a}-2\lambda\mathbf{a}=0\) gives
\[\boxed{\;S\mathbf{a}=\lambda\mathbf{a}\;}\]
so the optimal direction is an eigenvector of \(S\), and the variance it achieves is \(\mathbf{a}^\top S\mathbf{a}=\lambda\mathbf{a}^\top\mathbf{a}=\lambda\), the corresponding eigenvalue. To maximize, take the eigenvector of the largest eigenvalue.
The \(k\)-th principal component is then \(\mathbf{z}_k=X\mathbf{a}_k\) where \(\mathbf{a}_k\) is the eigenvector of the \(k\)-th largest eigenvalue, subject to:
- \(\operatorname{Var}(\mathbf{z}_k)\) is maximal,
- \(\operatorname{Cov}(\mathbf{z}_k,\mathbf{z}_\ell)=0\) for all \(\ell<k\),
- \(\mathbf{a}_k^\top\mathbf{a}_k=1\).
Condition 2 is automatic: \(S\) is symmetric, so by the spectral theorem (Chapter 3, §3.6.4) its eigenvectors are orthogonal, and orthogonal loadings give uncorrelated scores.
Since \(\sum_k\lambda_k=\operatorname{tr}(S)=\sum_j\operatorname{Var}(X_j)\), the proportion of variance explained by the first \(k\) components is \(\sum_{j\le k}\lambda_j/\sum_j\lambda_j\).
7.2 Derivation 2: minimize reconstruction error
Ask instead for the \(k\)-dimensional subspace that best approximates the data. With orthonormal basis \(V_k\) (\(p\times k\)), the projection of row \(\mathbf{x}_i\) is \(V_kV_k^\top\mathbf{x}_i\), and we minimize
\[J(V_k)=\sum_{i=1}^{n}\big\|\mathbf{x}_i-V_kV_k^\top\mathbf{x}_i\big\|^2 =\|X-XV_kV_k^\top\|_F^2 .\]
Expanding, and using orthonormality,
\[J(V_k)=\operatorname{tr}(X^\top X)-\operatorname{tr}\!\big(V_k^\top X^\top XV_k\big) =\underbrace{(n-1)\operatorname{tr}(S)}_{\text{fixed}}-(n-1)\operatorname{tr}\!\big(V_k^\top SV_k\big).\]
Minimizing \(J\) is therefore maximizing \(\operatorname{tr}(V_k^\top SV_k)\) — the total variance retained. The two objectives are the same problem, and the solution is the top-\(k\) eigenvectors of \(S\).
This second view connects PCA directly to Eckart–Young–Mirsky (Chapter 3, §3.6.5): the optimal rank-\(k\) approximation of \(X\) in Frobenius norm is the truncated SVD, and its residual is
\[\|X-X_k\|_F^2=\sum_{j>k}\sigma_j^2=(n-1)\sum_{j>k}\lambda_j .\]
No other \(k\)-dimensional linear reduction does better. That optimality is PCA’s entire claim, and it is a theorem, not a heuristic.
7.3 PCA in practice
The PPMI subset records six clinical and imaging measures on 33 Parkinson’s disease participants.
library(rvest)
pd_sub <- read_html("https://wiki.socr.umich.edu/index.php/SMHS_PCA_ICA_FA") |>
html_nodes("table") |> _[[1]] |> html_table()
pd_sub <- pd_sub[, -1] # drop the patient ID column
dim(pd_sub); names(pd_sub)#> [1] 33 6#> [1] "Top_of_SN_Voxel_Intensity_Ratio" "Side_of_SN_Voxel_Intensity_Ratio"
#> [3] "Part_IA" "Part_IB"
#> [5] "Part_II" "Part_III"#> Top_of_SN_Voxel_Intensity_Ratio Side_of_SN_Voxel_Intensity_Ratio
#> Min. :1.06 Min. :0.931
#> 1st Qu.:1.33 1st Qu.:0.996
#> Median :1.48 Median :1.111
#> Mean :1.53 Mean :1.106
#> 3rd Qu.:1.75 3rd Qu.:1.198
#> Max. :2.15 Max. :1.381
#> Part_IA Part_IB Part_II Part_III
#> Min. :0.00 Min. : 0.00 Min. : 0.00 Min. : 0.0
#> 1st Qu.:0.00 1st Qu.: 2.00 1st Qu.: 0.00 1st Qu.: 2.0
#> Median :1.00 Median : 5.00 Median : 2.00 Median :12.0
#> Mean :1.24 Mean : 4.91 Mean : 4.09 Mean :13.4
#> 3rd Qu.:2.00 3rd Qu.: 7.00 3rd Qu.: 6.00 3rd Qu.:20.0
#> Max. :6.00 Max. :13.00 Max. :17.00 Max. :36.0Centering means setting each column’s mean to zero, not subtracting a single number from the whole matrix:
X <- as.matrix(pd_sub)
Xc <- sweep(X, 2, colMeans(X)) # equivalently scale(X, scale = FALSE)
round(colMeans(Xc), 12) # every column mean is now zero#> Top_of_SN_Voxel_Intensity_Ratio Side_of_SN_Voxel_Intensity_Ratio
#> 0 0
#> Part_IA Part_IB
#> 0 0
#> Part_II Part_III
#> 0 0#> [1] 131.5073 11.7834 6.0969 1.4244 0.0609 0.0080#> [,1] [,2] [,3]
#> [1,] -0.0075 -0.0182 0.0169
#> [2,] -0.0058 0.0006 0.0042
#> [3,] 0.0808 -0.0600 -0.0274
#> [4,] 0.2297 -0.2818 -0.9295
#> [5,] 0.2821 -0.8926 0.3445
#> [6,] 0.9279 0.3462 0.1279pca_cov <- prcomp(X, center = TRUE, scale. = FALSE)
# The loadings are the eigenvectors, up to an arbitrary sign per column
round(pca_cov$rotation[, 1:3], 4)#> PC1 PC2 PC3
#> Top_of_SN_Voxel_Intensity_Ratio 0.0075 0.0182 -0.0169
#> Side_of_SN_Voxel_Intensity_Ratio 0.0058 -0.0006 -0.0042
#> Part_IA -0.0808 0.0600 0.0274
#> Part_IB -0.2297 0.2818 0.9295
#> Part_II -0.2821 0.8926 -0.3445
#> Part_III -0.9279 -0.3462 -0.1279#> eigenvalues_match
#> TRUECommon misconception: “the loadings are the eigenvectors times \(-1\).” The sign of every principal component is arbitrary. If \(\mathbf{a}\) is a unit eigenvector then so is \(-\mathbf{a}\), with the same eigenvalue and the same subspace. Which one a routine returns depends on the LAPACK build, the algorithm (
prcompuses SVD,eigenuses symmetric QR), and even the platform. Never interpret the sign of a loading in isolation, only the relative signs within a component carry meaning, and if you need a convention, fix one yourself (for example, force the largest-magnitude loading in each column to be positive).
fix_signs <- function(rot) {
s <- sign(rot[cbind(apply(abs(rot), 2, which.max), seq_len(ncol(rot)))])
sweep(rot, 2, s, "*")
}
round(fix_signs(pca_cov$rotation)[, 1:3], 4)#> PC1 PC2 PC3
#> Top_of_SN_Voxel_Intensity_Ratio -0.0075 0.0182 -0.0169
#> Side_of_SN_Voxel_Intensity_Ratio -0.0058 -0.0006 -0.0042
#> Part_IA 0.0808 0.0600 0.0274
#> Part_IB 0.2297 0.2818 0.9295
#> Part_II 0.2821 0.8926 -0.3445
#> Part_III 0.9279 -0.3462 -0.12797.4 Scale dependence: the pitfall that matters most
\(S\) depends on the units of every column. A variable measured in cubic millimetres has a variance \(10^9\) times larger than the same quantity in cubic centimetres, and PCA will hand it the first component almost regardless of its scientific relevance.
data.frame(variable = colnames(X),
sd = round(apply(X, 2, sd), 4),
variance = round(apply(X, 2, var), 4),
share_of_total_variance = round(apply(X, 2, var) / sum(apply(X, 2, var)), 4))#> variable sd
#> Top_of_SN_Voxel_Intensity_Ratio Top_of_SN_Voxel_Intensity_Ratio 0.2676
#> Side_of_SN_Voxel_Intensity_Ratio Side_of_SN_Voxel_Intensity_Ratio 0.1256
#> Part_IA Part_IA 1.5213
#> Part_IB Part_IB 3.6260
#> Part_II Part_II 4.5371
#> Part_III Part_III 10.7120
#> variance share_of_total_variance
#> Top_of_SN_Voxel_Intensity_Ratio 0.0716 0.0005
#> Side_of_SN_Voxel_Intensity_Ratio 0.0158 0.0001
#> Part_IA 2.3144 0.0153
#> Part_IB 13.1477 0.0871
#> Part_II 20.5852 0.1364
#> Part_III 114.7462 0.7605pca_cor <- prcomp(X, center = TRUE, scale. = TRUE) # correlation-based
comparison <- data.frame(
PC = paste0("PC", 1:ncol(X)),
covariance_var_pct = round(100 * pca_cov$sdev^2 / sum(pca_cov$sdev^2), 2),
correlation_var_pct = round(100 * pca_cor$sdev^2 / sum(pca_cor$sdev^2), 2))
comparison#> PC covariance_var_pct correlation_var_pct
#> 1 PC1 87.16 53.27
#> 2 PC2 7.81 20.36
#> 3 PC3 4.04 9.50
#> 4 PC4 0.94 6.99
#> 5 PC5 0.04 5.39
#> 6 PC6 0.01 4.48data.frame(loading_on_PC1 = colnames(X),
covariance_PCA = round(fix_signs(pca_cov$rotation)[, 1], 4),
correlation_PCA = round(fix_signs(pca_cor$rotation)[, 1], 4))#> loading_on_PC1
#> Top_of_SN_Voxel_Intensity_Ratio Top_of_SN_Voxel_Intensity_Ratio
#> Side_of_SN_Voxel_Intensity_Ratio Side_of_SN_Voxel_Intensity_Ratio
#> Part_IA Part_IA
#> Part_IB Part_IB
#> Part_II Part_II
#> Part_III Part_III
#> covariance_PCA correlation_PCA
#> Top_of_SN_Voxel_Intensity_Ratio -0.0075 -0.2555
#> Side_of_SN_Voxel_Intensity_Ratio -0.0058 -0.3855
#> Part_IA 0.0808 0.3825
#> Part_IB 0.2297 0.4597
#> Part_II 0.2821 0.4251
#> Part_III 0.9279 0.4977The two analyses disagree substantially, and neither is “wrong”, they
answer different questions. Use the covariance matrix
when all variables share meaningful units and their relative variances
are informative (repeated measures of the same quantity, spectra, pixel
intensities). Use the correlation matrix, that is,
scale. = TRUE, when units are arbitrary or incomparable,
which covers most clinical, survey, and multi-instrument data. In
practice, standardize unless you have a specific reason not
to.
scree <- data.frame(
PC = seq_along(pca_cor$sdev),
variance = pca_cor$sdev^2,
proportion = pca_cor$sdev^2 / sum(pca_cor$sdev^2)) |>
mutate(cumulative = cumsum(proportion))
p1 <- ggplot(scree, aes(PC, variance)) +
geom_col(fill = "steelblue") +
geom_hline(yintercept = 1, linetype = "dashed", colour = "firebrick") +
scale_x_continuous(breaks = scree$PC) +
labs(title = "Scree plot (correlation PCA)",
subtitle = "Dashed line: eigenvalue 1, the Kaiser criterion",
x = "Component", y = "Eigenvalue") + theme_dspa(10)
p2 <- ggplot(scree, aes(PC, cumulative)) +
geom_line(linewidth = 0.9, colour = "firebrick") + geom_point(size = 2.2) +
geom_hline(yintercept = c(0.8, 0.9), linetype = "dotted", colour = "grey40") +
scale_x_continuous(breaks = scree$PC) +
scale_y_continuous(labels = scales::percent, limits = c(0, 1)) +
labs(title = "Cumulative variance explained", x = "Component", y = NULL) +
theme_dspa(10)
p1 | p2
7.5 Choosing the number of components
An “elbow” in the scree plot is a starting point, not a method, real scree plots are frequently smooth, and different readers see different elbows. Four defensible alternatives:
Proportion of variance. Retain enough components to reach a preset threshold (80% or 90%). Simple, arbitrary, and honest about being arbitrary.
Kaiser criterion. On the correlation matrix, retain components with \(\lambda_j>1\), those explaining more than a single standardized variable. Widely used and widely criticized for over-retaining.
Broken stick. Under a null of no structure, the expected share of the \(j\)-th ordered eigenvalue is \[b_j=\frac{1}{p}\sum_{m=j}^{p}\frac{1}{m}.\] Retain components whose observed share exceeds \(b_j\).
Parallel analysis (Horn, 1965). Generate many datasets of the same shape with independent columns, compute their eigenvalues, and retain components exceeding an upper percentile of the null distribution. This is the most defensible of the four because it is a calibrated comparison rather than a fixed cutoff.
Cross-validated reconstruction error. Hold out entries, reconstruct from \(k\) components, and choose the \(k\) minimizing held-out error. It answers the question you usually actually care about, predictive adequacy, and it does not assume anything about the eigenvalue distribution.
p_dim <- ncol(X)
obs_prop <- pca_cor$sdev^2 / sum(pca_cor$sdev^2)
broken_stick <- sapply(seq_len(p_dim), \(j) sum(1 / (j:p_dim)) / p_dim)
set.seed(31)
null_eigen <- replicate(1000, {
Xn <- apply(X, 2, \(col) sample(col)) # permute within columns: kills structure
prcomp(Xn, center = TRUE, scale. = TRUE)$sdev^2
})
pa_95 <- apply(null_eigen, 1, quantile, 0.95)
data.frame(
PC = 1:p_dim,
eigenvalue = round(pca_cor$sdev^2, 3),
proportion = round(obs_prop, 3),
broken_stick = round(broken_stick, 3),
keep_broken_stick = obs_prop > broken_stick,
parallel_95 = round(pa_95, 3),
keep_parallel = pca_cor$sdev^2 > pa_95,
keep_kaiser = pca_cor$sdev^2 > 1)#> PC eigenvalue proportion broken_stick keep_broken_stick parallel_95
#> 1 1 3.196 0.533 0.408 TRUE 1.904
#> 2 2 1.222 0.204 0.242 FALSE 1.486
#> 3 3 0.570 0.095 0.158 FALSE 1.196
#> 4 4 0.419 0.070 0.103 FALSE 0.996
#> 5 5 0.324 0.054 0.061 FALSE 0.817
#> 6 6 0.269 0.045 0.028 TRUE 0.634
#> keep_parallel keep_kaiser
#> 1 TRUE TRUE
#> 2 FALSE TRUE
#> 3 FALSE FALSE
#> 4 FALSE FALSE
#> 5 FALSE FALSE
#> 6 FALSE FALSEdata.frame(PC = 1:p_dim, Observed = pca_cor$sdev^2,
`Parallel analysis (95th pct)` = pa_95, check.names = FALSE) |>
pivot_longer(-PC, names_to = "series", values_to = "eigenvalue") |>
ggplot(aes(PC, eigenvalue, colour = series, linetype = series)) +
geom_line(linewidth = 1) + geom_point(size = 2.2) +
scale_x_continuous(breaks = 1:p_dim) +
scale_colour_manual(values = c("Observed" = "steelblue",
"Parallel analysis (95th pct)" = "firebrick")) +
labs(title = "Parallel analysis: keep components above the null band",
subtitle = "Null eigenvalues from 1,000 within-column permutations",
x = "Component", y = "Eigenvalue", colour = NULL, linetype = NULL) +
theme_dspa()
# Cross-validated reconstruction error: hold out a random 10% of ENTRIES,
# impute with the column mean, fit rank-k, and score only the held-out cells.
cv_rank <- function(X, k_max = ncol(X) - 1, folds = 10, seed = 41) {
set.seed(seed)
Xs <- scale(X)
n_cell <- length(Xs)
fold_id <- sample(rep(seq_len(folds), length.out = n_cell))
err <- matrix(NA_real_, folds, k_max)
for (f in seq_len(folds)) {
hold <- which(fold_id == f)
Xtr <- Xs; Xtr[hold] <- NA
cm <- colMeans(Xtr, na.rm = TRUE)
for (j in seq_len(ncol(Xtr))) Xtr[is.na(Xtr[, j]), j] <- cm[j]
sv <- svd(Xtr)
for (k in seq_len(k_max)) {
Xk <- sv$u[, 1:k, drop = FALSE] %*% diag(sv$d[1:k], k, k) %*%
t(sv$v[, 1:k, drop = FALSE])
err[f, k] <- mean((Xs[hold] - Xk[hold])^2)
}
}
data.frame(k = seq_len(k_max), cv_mse = colMeans(err), se = apply(err, 2, sd) / sqrt(folds))
}
cvres <- cv_rank(X)
cvres |> mutate(across(where(is.numeric), \(z) round(z, 4)))#> k cv_mse se
#> 1 1 0.6662 0.0834
#> 2 2 0.6661 0.0711
#> 3 3 0.8021 0.0972
#> 4 4 0.9055 0.0902
#> 5 5 0.9771 0.0819#> CV-optimal rank: 2ggplot(cvres, aes(k, cv_mse)) +
geom_ribbon(aes(ymin = cv_mse - se, ymax = cv_mse + se), fill = "grey85") +
geom_line(linewidth = 1, colour = "steelblue") +
geom_point(data = cvres[which.min(cvres$cv_mse), ], colour = "firebrick", size = 3.4) +
scale_x_continuous(breaks = cvres$k) +
labs(title = "Cross-validated reconstruction error selects the rank",
subtitle = "Held-out cells only; band is +/- 1 standard error",
x = "Number of components k", y = "Held-out MSE") +
theme_dspa()
7.6 Bootstrap uncertainty on variance explained
The proportion of variance captured is itself a statistic with sampling variability. Bootstrap it, taking care that the resampled analysis and the point estimate use the same PCA model.
set.seed(12)
B <- 1000
k_keep <- 3
boot_prop <- replicate(B, {
Xb <- X[sample(nrow(X), replace = TRUE), , drop = FALSE]
sdv <- prcomp(Xb, center = TRUE, scale. = TRUE)$sdev # SAME model as the estimate
sum(sdv[1:k_keep]^2) / sum(sdv^2)
})
point_est <- sum(pca_cor$sdev[1:k_keep]^2) / sum(pca_cor$sdev^2)
ci <- quantile(boot_prop, c(0.025, 0.975))
c(point_estimate = point_est, lower = ci[[1]], upper = ci[[2]])#> point_estimate lower upper
#> 0.831374 0.812461 0.898532ggplot(data.frame(p = boot_prop), aes(p)) +
geom_histogram(bins = 40, fill = "steelblue", colour = "white") +
geom_vline(xintercept = point_est, colour = "firebrick", linewidth = 1.1) +
geom_vline(xintercept = ci, colour = "grey30", linetype = "dashed") +
scale_x_continuous(labels = scales::percent) +
labs(title = sprintf("Bootstrap distribution: variance captured by the first %d components", k_keep),
subtitle = "Red: point estimate. Dashed: 95% percentile interval",
x = "Proportion of variance", y = "Count") +
theme_dspa()