| 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)
August 26, 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.
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()
# --- Interactive equivalent ------------------------------------------------
plot_ly(x = boot_prop, type = "histogram", name = "Bootstrap") |>
layout(bargap = 0.1,
title = "Percent of variability captured by the first 3 PCs",
xaxis = list(title = "Proportion of variability"),
yaxis = list(title = "Frequency"))7.7 The biplot and the loadings triplot
A biplot overlays the scores (cases in PC space) with the loadings (original variables as arrows). Arrow direction shows which components a variable drives; arrow length shows how strongly.
scores <- as.data.frame(pca_cor$x)
loads <- as.data.frame(fix_signs(pca_cor$rotation))
loads$variable <- rownames(loads)
arrow_scale <- 0.9 * max(abs(scores[, 1:2])) / max(abs(as.matrix(loads[, 1:2])))
ggplot(scores, aes(PC1, PC2)) +
geom_point(alpha = 0.65, size = 2, colour = "grey30") +
geom_segment(data = loads,
aes(x = 0, y = 0, xend = PC1 * arrow_scale, yend = PC2 * arrow_scale),
arrow = arrow(length = unit(0.2, "cm")), colour = "firebrick",
linewidth = 0.7, inherit.aes = FALSE) +
geom_text(data = loads,
aes(x = PC1 * arrow_scale * 1.12, y = PC2 * arrow_scale * 1.12,
label = variable),
colour = "firebrick", size = 3, inherit.aes = FALSE) +
labs(title = "PCA biplot (correlation matrix)",
subtitle = sprintf("PC1 %.1f%% of variance, PC2 %.1f%%",
100 * obs_prop[1], 100 * obs_prop[2]),
x = sprintf("PC1 (%.1f%%)", 100 * obs_prop[1]),
y = sprintf("PC2 (%.1f%%)", 100 * obs_prop[2])) +
theme_dspa()
The three-component version is genuinely three-dimensional, and the loadings are directions in that space, so this one is interactive:
sc <- pca_cor$x
ld <- fix_signs(pca_cor$rotation)
scale_load <- 0.85 * max(abs(sc[, 1:3])) / max(abs(ld[, 1:3]))
p <- plot_ly() |>
add_trace(x = sc[, 1], y = sc[, 2], z = sc[, 3], type = "scatter3d",
mode = "markers", name = "Participants",
marker = list(size = 4, opacity = 0.75, color = sc[, 1],
colorscale = "Viridis", showscale = FALSE))
for (k in seq_len(ncol(ld))) {
p <- add_trace(p,
x = c(0, ld[k, 1]) * scale_load,
y = c(0, ld[k, 2]) * scale_load,
z = c(0, ld[k, 3]) * scale_load,
type = "scatter3d", mode = "lines",
name = rownames(ld)[k], line = list(width = 7))
}
p |> layout(title = "Scores on the first three PCs with variable loadings",
legend = list(orientation = "h"),
scene = list(xaxis = list(title = sprintf("PC1 (%.0f%%)", 100 * obs_prop[1])),
yaxis = list(title = sprintf("PC2 (%.0f%%)", 100 * obs_prop[2])),
zaxis = list(title = sprintf("PC3 (%.0f%%)", 100 * obs_prop[3]))))Each line radiates from the origin along one original variable’s loading vector. Variables pointing in similar directions are correlated; variables at right angles are approximately independent in this projection; a short line means that variable is poorly represented by the first three components.
7.8 Nonlinear PCA
Classical PCA assumes the informative structure is a
linear subspace and that the variables are
quantitative. Two generalizations relax this. Nonlinear
(categorical) PCA, optimal scaling, assigns numeric values to
categories so as to maximize the association among the quantified
variables, handling nominal and ordinal data directly; see Linting’s
treatment and the pcaMethods
package. Kernel PCA (§4.10) performs
linear PCA in an implicit high-dimensional feature space, which curves
the components in the original space.
8 The first principal component is not the regression line
This distinction is easy to state, easy to get wrong, and consequential.
Given bivariate data, four different lines can be fitted:
| Line | Minimizes | Slope |
|---|---|---|
| OLS \(y\) on \(x\) | \(\sum(y_i-\hat y_i)^2\) — vertical | \(s_{xy}/s_x^2\) |
| OLS \(x\) on \(y\) | \(\sum(x_i-\hat x_i)^2\) — horizontal | \(s_y^2/s_{xy}\) |
| PC1 | \(\sum \operatorname{dist}^2(p_i,\ell)\) — perpendicular | see below |
| Reduced major axis | Area of triangles | \(\operatorname{sign}(s_{xy})\,s_y/s_x\) |
The PC1 slope follows from the leading eigenvector of the \(2\times2\) covariance:
\[m_{\mathrm{PC1}}=\frac{s_y^2-s_x^2+\sqrt{\left(s_y^2-s_x^2\right)^2+4s_{xy}^2}}{2s_{xy}} .\]
All four coincide when \(|r|=1\). They diverge as the scatter grows.
Common misconception: “PC1 is the regression line.” OLS treats \(x\) as fixed and error-free and measures error only in \(y\); PC1 treats both coordinates symmetrically and measures error perpendicular to the line. The practical consequence is attenuation: as measurement noise in \(x\) increases, the OLS slope is biased toward zero, while PC1 (total least squares) is not. If both variables are measured with error, two instruments, two raters, two assays, the regression line is the wrong summary of the relationship, and which variable you happen to put on the \(y\) axis changes the answer.
ex <- data.frame(x = c(1, 3, 5, 6, 10, 16, 50),
y = c(4, 6, 5, 7, 10, 13, 12))
sx2 <- var(ex$x); sy2 <- var(ex$y); sxy <- cov(ex$x, ex$y)
xbar <- mean(ex$x); ybar <- mean(ex$y)
m_ols_yx <- sxy / sx2
m_ols_xy <- sy2 / sxy # slope in the y-on-x frame
m_pc1 <- (sy2 - sx2 + sqrt((sy2 - sx2)^2 + 4 * sxy^2)) / (2 * sxy)
m_rma <- sign(sxy) * sqrt(sy2 / sx2)
# Verify against prcomp: the PC1 direction is the leading loading vector
v1 <- prcomp(as.matrix(ex), center = TRUE, scale. = FALSE)$rotation[, 1]
c(closed_form = m_pc1, from_prcomp = unname(v1[2] / v1[1]))#> closed_form from_prcomp
#> 0.149925 0.149925lines_df <- data.frame(
method = factor(c("OLS y on x", "OLS x on y", "PC1 (total least squares)",
"Reduced major axis"),
levels = c("OLS y on x", "OLS x on y",
"PC1 (total least squares)", "Reduced major axis")),
slope = c(m_ols_yx, m_ols_xy, m_pc1, m_rma))
lines_df$intercept <- ybar - lines_df$slope * xbar
lines_df#> method slope intercept
#> 1 OLS y on x 0.146789 6.23460
#> 2 OLS x on y 0.292411 4.34152
#> 3 PC1 (total least squares) 0.149925 6.19384
#> 4 Reduced major axis 0.207178 5.44954ggplot(ex, aes(x, y)) +
geom_point(size = 3, colour = "grey25") +
geom_point(aes(x = xbar, y = ybar), colour = "gold", size = 5, shape = 21,
fill = "gold", stroke = 1.2) +
geom_abline(data = lines_df, aes(slope = slope, intercept = intercept,
colour = method), linewidth = 1) +
scale_colour_brewer(palette = "Set1") +
coord_fixed() +
labs(title = "Four lines through the same seven points",
subtitle = "All pass through the centroid (gold); only PC1 minimizes perpendicular distance",
x = "X", y = "Y", colour = NULL) +
theme_dspa()
# Draw the residual segments each method minimizes
seg <- function(slope, intercept, kind) {
if (kind == "vertical") {
data.frame(x = ex$x, y = ex$y, xend = ex$x, yend = intercept + slope * ex$x)
} else {
t0 <- (ex$x + slope * (ex$y - intercept)) / (1 + slope^2)
data.frame(x = ex$x, y = ex$y, xend = t0, yend = intercept + slope * t0)
}
}
s_ols <- seg(m_ols_yx, ybar - m_ols_yx * xbar, "vertical")
s_pc1 <- seg(m_pc1, ybar - m_pc1 * xbar, "perp")
p_a <- ggplot(ex, aes(x, y)) +
geom_abline(slope = m_ols_yx, intercept = ybar - m_ols_yx * xbar,
colour = "steelblue", linewidth = 1) +
geom_segment(data = s_ols, aes(x = x, y = y, xend = xend, yend = yend),
colour = "firebrick", linetype = "dashed") +
geom_point(size = 2.6) + coord_fixed() +
labs(title = "OLS minimizes VERTICAL residuals", x = "X", y = "Y") + theme_dspa(10)
p_b <- ggplot(ex, aes(x, y)) +
geom_abline(slope = m_pc1, intercept = ybar - m_pc1 * xbar,
colour = "darkgreen", linewidth = 1) +
geom_segment(data = s_pc1, aes(x = x, y = y, xend = xend, yend = yend),
colour = "firebrick", linetype = "dashed") +
geom_point(size = 2.6) + coord_fixed() +
labs(title = "PC1 minimizes PERPENDICULAR residuals", x = "X", y = "Y") + theme_dspa(10)
p_a | p_b
# --- Interactive equivalent ------------------------------------------------
plot_ly(ex) |>
add_markers(x = ~x, y = ~y, name = "Data") |>
add_lines(x = ~x, y = ~(ybar - m_ols_yx * xbar) + m_ols_yx * x,
name = "OLS y on x", line = list(width = 4)) |>
add_lines(x = ~x, y = ~(ybar - m_pc1 * xbar) + m_pc1 * x,
name = "PC1 (total least squares)", line = list(width = 4)) |>
add_markers(x = xbar, y = ybar, name = "Centroid",
marker = list(size = 20, color = "gold",
line = list(color = "black", width = 2))) |>
layout(title = "OLS versus the first principal component",
xaxis = list(title = "X", scaleanchor = "y"),
yaxis = list(title = "Y"), legend = list(orientation = "h"))The attenuation is measurable. Simulate a true relationship \(y=2x\) and add increasing measurement error to \(x\) only:
attenuate <- function(sd_x, reps = 400, n = 200, seed = 53) {
set.seed(seed)
out <- replicate(reps, {
x_true <- rnorm(n)
y <- 2 * x_true + rnorm(n, sd = 0.5)
x <- x_true + rnorm(n, sd = sd_x) # x is measured with error
sx2 <- var(x); sy2 <- var(y); sxy <- cov(x, y)
c(ols = sxy / sx2,
pc1 = (sy2 - sx2 + sqrt((sy2 - sx2)^2 + 4 * sxy^2)) / (2 * sxy))
})
data.frame(sd_x = sd_x, OLS = mean(out["ols", ]), PC1 = mean(out["pc1", ]))
}
att <- do.call(rbind, lapply(c(0, 0.1, 0.25, 0.5, 0.75, 1), attenuate))
att |> mutate(across(where(is.numeric), \(z) round(z, 3)))#> sd_x OLS PC1
#> 1 0.00 1.999 2.101
#> 2 0.10 1.980 2.099
#> 3 0.25 1.883 2.080
#> 4 0.50 1.600 2.008
#> 5 0.75 1.279 1.889
#> 6 1.00 0.998 1.725att |> pivot_longer(-sd_x, names_to = "method", values_to = "slope") |>
ggplot(aes(sd_x, slope, colour = method)) +
geom_hline(yintercept = 2, linetype = "dashed", colour = "grey40") +
geom_line(linewidth = 1) + geom_point(size = 2.4) +
scale_colour_manual(values = c(OLS = "steelblue", PC1 = "darkgreen")) +
labs(title = "Measurement error in x attenuates the OLS slope",
subtitle = "True slope is 2 (dashed). PC1 stays close; OLS drifts toward zero",
x = "Standard deviation of measurement error in x",
y = "Average estimated slope", colour = NULL) +
theme_dspa()
9 SVD and its exact relation to PCA
Chapter 3 (§3.6.5) established that every \(X\in\mathbb{R}^{n\times p}\) factors as \(X=UDV^\top\). Applied to a centred data matrix, this is PCA, no separate algorithm is needed.
9.1 The correspondences
Let \(X\) be \(n\times p\), column-centred, with thin SVD \(X=UDV^\top\) where \(U\) is \(n\times r\) with orthonormal columns (\(U^\top U=I_r\)), \(D=\operatorname{diag}(d_1\ge\cdots\ge d_r>0)\), and \(V\) is \(p\times r\) with \(V^\top V=I_r\), \(r=\operatorname{rank}(X)\).
A note on terminology. In the full SVD, \(U\) is \(n\times n\) and genuinely orthogonal; in the thin SVD used in practice, \(U\) has orthonormal columns but is rectangular, so \(U^\top U=I_r\) while \(UU^\top\ne I_n\). Only the square version deserves the name “orthogonal” (or “unitary” over \(\mathbb{C}\)). Note also that \(|\det|=1\) does not characterize orthogonality: a shear \(\left(\begin{smallmatrix}1&c\\0&1\end{smallmatrix}\right)\) has determinant 1 and distorts every angle. Orthogonality constrains angles; the determinant constrains volume.
Then:
\[S=\frac{X^\top X}{n-1}=\frac{VD^2V^\top}{n-1} \quad\Longrightarrow\quad \boxed{\;\lambda_j=\frac{d_j^2}{n-1}\;}\]
| PCA object | SVD expression |
|---|---|
| Principal directions (loadings) | columns of \(V\) |
| Principal component scores | columns of \(UD\) |
| Variance of PC \(j\) | \(\lambda_j=d_j^2/(n-1)\) |
| Standardized scores | columns of \(\sqrt{n-1}\,U\) |
| Loadings scaled to reproduce \(X\) | columns of \(VD/\sqrt{n-1}\) |
| Rank-\(k\) reconstruction | \(X_k=U_kD_kV_k^\top\) |
Note the \(\sqrt{n-1}\) in the fifth row. It is forced by requiring \(X=(\text{scores})(\text{loadings})^\top\): \[X=UDV^\top=\big(\sqrt{n-1}\,U\big)\left(\frac{VD}{\sqrt{n-1}}\right)^{\!\top}.\]
To work with the correlation matrix instead, standardize the columns first (centre and scale); everything else is unchanged.
Xs <- scale(X) # centre and scale -> correlation PCA
sv <- svd(Xs)
n_ <- nrow(Xs)
c(eigenvalues_agree = isTRUE(all.equal(sv$d^2 / (n_ - 1), pca_cor$sdev^2)),
scores_agree = isTRUE(all.equal(abs(sv$u %*% diag(sv$d)), abs(pca_cor$x),
check.attributes = FALSE)),
loadings_agree = isTRUE(all.equal(abs(sv$v), abs(pca_cor$rotation),
check.attributes = FALSE)))#> eigenvalues_agree scores_agree loadings_agree
#> TRUE TRUE TRUE#> [1] 3.1963 1.2219 0.5701 0.4193 0.3236 0.2688#> [,1] [,2] [,3]
#> [1,] -0.2555 0.7126 -0.3732
#> [2,] -0.3855 0.4721 0.3567
#> [3,] 0.3825 0.3729 0.7099
#> [4,] 0.4597 0.0980 -0.1117
#> [5,] 0.4251 0.3417 -0.4642
#> [6,] 0.4977 0.0626 0.0387Absolute values are compared because of the sign indeterminacy noted in §4.4.3.
9.2 Rank-\(k\) truncation error
Eckart–Young gives the reconstruction error exactly:
\[\|X-X_k\|_F^2=\sum_{j>k}d_j^2, \qquad \|X-X_k\|_2=d_{k+1}.\]
trunc_err <- function(k) {
Xk <- sv$u[, 1:k, drop = FALSE] %*% diag(sv$d[1:k], k, k) %*% t(sv$v[, 1:k, drop = FALSE])
c(k = k,
frobenius_observed = norm(Xs - Xk, "F"),
frobenius_theory = sqrt(sum(sv$d[-(1:k)]^2)),
spectral_observed = norm(Xs - Xk, "2"),
spectral_theory = sv$d[k + 1],
variance_retained = sum(sv$d[1:k]^2) / sum(sv$d^2))
}
do.call(rbind, lapply(1:5, trunc_err)) |> round(5)#> k frobenius_observed frobenius_theory spectral_observed spectral_theory
#> [1,] 1 9.47203 9.47203 6.25298 6.25298
#> [2,] 2 7.11474 7.11474 4.27122 4.27122
#> [3,] 3 5.69001 5.69001 3.66320 3.66320
#> [4,] 4 4.35399 4.35399 3.21804 3.21804
#> [5,] 5 2.93282 2.93282 2.93282 2.93282
#> variance_retained
#> [1,] 0.53271
#> [2,] 0.73636
#> [3,] 0.83137
#> [4,] 0.90126
#> [5,] 0.95520Observed and predicted agree to five decimals. This is not an approximation — it is the theorem.
9.3 Computing PCA when \(p\) is large
Three routes, with very different costs:
| Route | Cost | Use when |
|---|---|---|
| Eigen of \(S=X^\top X/(n-1)\) | \(O(np^2+p^3)\) | \(p\) small; but squares the condition number |
| Full SVD of \(X\) | \(O(np\min(n,p))\) | Default; prcomp() does this |
| Randomized SVD for top \(k\) | \(O(npk)\), or \(O(np\log k)\) with structured projections | \(k\ll\min(n,p)\) |
Forming \(X^\top X\) carries exactly
the numerical penalty from Chapter 3, §3.7.1: \(\kappa_2(X^\top X)=\kappa_2(X)^2\).
prcomp() uses the SVD of \(X\) for this reason and is preferred over
princomp(), which eigendecomposes the covariance matrix
(and divides by \(n\) rather than \(n-1\)).
Randomized SVD (Halko, Martinsson & Tropp, 2011) computes the top \(k\) components without touching the rest of the spectrum: sketch the column space with a random Gaussian matrix, orthonormalize, and take a small SVD in the sketched basis.
rsvd_simple <- function(A, k, oversample = 10, q = 2) {
n <- nrow(A); p <- ncol(A); l <- min(k + oversample, p)
Om <- matrix(rnorm(p * l), p, l)
Ymat <- A %*% Om
for (i in seq_len(q)) Ymat <- A %*% (t(A) %*% Ymat) # power iteration sharpens decay
Q <- qr.Q(qr(Ymat))
B <- t(Q) %*% A
s <- svd(B)
list(d = s$d[1:k], u = (Q %*% s$u)[, 1:k, drop = FALSE],
v = s$v[, 1:k, drop = FALSE])
}
set.seed(59)
n_big <- 2000; p_big <- 800; k_big <- 10
Utrue <- qr.Q(qr(matrix(rnorm(n_big * k_big), n_big)))
Vtrue <- qr.Q(qr(matrix(rnorm(p_big * k_big), p_big)))
Abig <- Utrue %*% diag(seq(50, 5, length.out = k_big)) %*% t(Vtrue) +
matrix(rnorm(n_big * p_big, sd = 0.05), n_big)
t_full <- system.time(sfull <- svd(Abig))[["elapsed"]]
t_rand <- system.time(srand <- rsvd_simple(Abig, k_big))[["elapsed"]]
c(full_svd_sec = round(t_full, 3),
randomized_sec = round(t_rand, 3),
speedup = round(t_full / t_rand, 1),
max_rel_error_top10 = max(abs(srand$d - sfull$d[1:k_big]) / sfull$d[1:k_big]))#> full_svd_sec randomized_sec speedup max_rel_error_top10
#> 2.0400000 0.1000000 20.4000000 0.0206828The top ten singular values are recovered to near machine precision
at a fraction of the cost. For gene-expression matrices, image
collections, or document–term matrices, anywhere \(\min(n,p)\) runs to thousands but you want
only a handful of components, this is the right tool. Production
implementations live in rsvd::rsvd() and
irlba::irlba().
10 Classical multidimensional scaling
PCA needs the data matrix. Sometimes all you have is a matrix of distances — from a survey of perceived similarity, a sequence-alignment score, or a proprietary metric. Classical MDS recovers coordinates from distances alone, and its relationship to PCA is an identity, not an analogy.
Given squared distances \(D^{(2)}=(d_{ij}^2)\), double-centre them:
\[B=-\tfrac12\,J\,D^{(2)}\,J,\qquad J=I-\tfrac1n\mathbf{1}\mathbf{1}^\top .\]
If the distances are Euclidean, then \(B=\tilde X\tilde X^\top\) where \(\tilde X\) is the centred configuration, the Gram matrix. Eigendecomposing \(B=Q\Lambda Q^\top\) and taking \(\tilde X=Q_k\Lambda_k^{1/2}\) recovers coordinates in \(k\) dimensions.
Classical MDS on Euclidean distances is exactly PCA. \(B=\tilde X\tilde X^\top\) and \(S\propto\tilde X^\top\tilde X\) share their nonzero eigenvalues, and the MDS coordinates equal the PCA scores up to sign. The difference is purely computational: PCA eigendecomposes a \(p\times p\) matrix, MDS an \(n\times n\) one. Choose whichever is smaller.
Dm <- dist(Xs) # Euclidean distances among the standardized cases
mds <- cmdscale(Dm, k = 3, eig = TRUE)
# Same eigenvalues (up to the n-1 scaling) and same coordinates up to sign
c(mds_eigen = round(mds$eig[1:3], 4),
pca_scaled = round(pca_cor$sdev[1:3]^2 * (n_ - 1), 4))#> mds_eigen1 mds_eigen2 mds_eigen3 pca_scaled1 pca_scaled2 pca_scaled3
#> 102.2807 39.0997 18.2433 102.2807 39.0997 18.2433cor_abs <- sapply(1:3, \(j) abs(cor(mds$points[, j], pca_cor$x[, j])))
c(coordinate_correlations = round(cor_abs, 8))#> coordinate_correlations1 coordinate_correlations2 coordinate_correlations3
#> 1 1 1# Double-centering from first principles
D2 <- as.matrix(Dm)^2
J <- diag(n_) - matrix(1 / n_, n_, n_)
B <- -0.5 * J %*% D2 %*% J
eB <- eigen(B, symmetric = TRUE)
c(manual_eigen = round(eB$values[1:3], 4),
cmdscale_eigen = round(mds$eig[1:3], 4),
gram_check = isTRUE(all.equal(B, tcrossprod(sweep(Xs, 2, colMeans(Xs))),
check.attributes = FALSE)))#> manual_eigen1 manual_eigen2 manual_eigen3 cmdscale_eigen1 cmdscale_eigen2
#> 102.2807 39.0997 18.2433 102.2807 39.0997
#> cmdscale_eigen3 gram_check
#> 18.2433 1.0000Non-metric MDS (MASS::isoMDS,
vegan::metaMDS) relaxes the requirement that distances be
Euclidean or even metric, preserving only their rank order by
minimizing a stress function. That is the first genuinely nonlinear
method in this chapter, and it foreshadows Part
III.
Cost: \(O(n^3)\) for the eigendecomposition and \(O(n^2)\) memory, so classical MDS becomes impractical beyond \(n\approx10^4\), the same wall hierarchical clustering hits (Chapter 3, §3.21).
11 Independent component analysis
11.1 The model
ICA assumes the observed variables are linear mixtures of statistically independent latent sources:
\[\mathbf{x}=A\mathbf{s},\qquad x_i=a_{i1}s_1+\cdots+a_{in}s_n,\]
where \(A\) is the unknown mixing matrix and the components of \(\mathbf{s}\) are mutually statistically independent. The task is to estimate an unmixing matrix \(W\approx A^{-1}\) so that \(\hat{\mathbf{s}}=W\mathbf{x}\) recovers the sources.
The canonical picture is the cocktail-party problem: \(n\) microphones record \(n\) speakers, each microphone capturing a different linear mixture, and ICA separates the voices.
Common misconception: “ICA does not assume independence.” Independence of the sources is ICA’s defining assumption, it is what the method is named for, and it is strictly stronger than the uncorrelatedness PCA imposes. Uncorrelatedness constrains second moments only; independence constrains all moments and the full joint density, \(p(s_1,\dots,s_n)=\prod_i p_i(s_i)\). What ICA drops relative to PCA is the assumption of Gaussianity, and that omission is essential rather than incidental, as the next paragraph shows.
11.2 Identifiability: what ICA can and cannot recover
Comon’s theorem (1994). The model \(\mathbf{x}=A\mathbf{s}\) with independent sources is identifiable up to permutation and scaling provided at most one source is Gaussian.
Three consequences you must internalize before interpreting any ICA output.
Scale (and sign) is unidentifiable. For any nonzero \(\alpha_i\), \(A\mathbf{s}=\big(A\Lambda^{-1}\big)\big(\Lambda\mathbf{s}\big)\) with \(\Lambda=\operatorname{diag}(\alpha_i)\). Multiplying a source by 2 and dividing its mixing column by 2 gives identical data. Implementations therefore fix the scale by convention (unit variance), and the sign is arbitrary.
Order is unidentifiable. For any permutation \(P\), \(A\mathbf{s}=(AP^{-1})(P\mathbf{s})\). There is no “first” independent component the way there is a first principal component, ICA components carry no natural ranking, because there is no variance to rank them by.
At most one Gaussian source. If two sources are Gaussian and independent, their joint density is spherically symmetric, so any orthogonal rotation of them is equally independent and equally Gaussian. The mixing matrix cannot be recovered even in principle.
set.seed(61)
n_ica <- 5000
# Two Gaussian sources: rotation is unidentifiable
Sg <- cbind(rnorm(n_ica), rnorm(n_ica))
theta <- pi / 6
Rot <- matrix(c(cos(theta), -sin(theta), sin(theta), cos(theta)), 2)
Sg_rot <- Sg %*% t(Rot)
c(both_independent = cor(Sg[, 1], Sg[, 2]),
rotation_also_independent = cor(Sg_rot[, 1], Sg_rot[, 2]),
same_joint_density = "yes -- spherical symmetry makes them indistinguishable")#> both_independent
#> "0.0083468889361991"
#> rotation_also_independent
#> "0.000750303399613017"
#> same_joint_density
#> "yes -- spherical symmetry makes them indistinguishable"11.3 Whitening and the non-Gaussianity objective
ICA proceeds in two stages.
Stage 1, whitening. Transform \(\mathbf{x}\) so its covariance is the identity: \(\tilde{\mathbf{x}}=ED^{-1/2}E^\top\mathbf{x}\) where \(\operatorname{Cov}(\mathbf{x})=EDE^\top\). This is PCA plus rescaling, and it reduces the problem from finding an arbitrary \(W\) to finding an orthogonal one, halving the free parameters from \(n^2\) to \(n(n-1)/2\).
Stage 2, maximize non-Gaussianity. By the central limit theorem, a sum of independent variables is more Gaussian than its summands. So a mixture \(\mathbf{w}^\top\tilde{\mathbf{x}}\) is most non-Gaussian precisely when it equals a single source. Maximizing non-Gaussianity therefore recovers the sources. Two standard measures:
\[\text{excess kurtosis:}\quad \operatorname{kurt}(y)=E[y^4]-3\big(E[y^2]\big)^2,\] \[\text{negentropy:}\quad J(y)=H(y_{\text{gauss}})-H(y)\;\ge\;0,\]
where \(H\) is differential entropy and \(y_{\text{gauss}}\) is Gaussian with the same variance. Negentropy is zero iff \(y\) is Gaussian, which makes it the theoretically ideal contrast, but it requires the unknown density, so FastICA uses the approximation
\[J(y)\;\propto\;\big(E[G(y)]-E[G(\nu)]\big)^2,\qquad \nu\sim N(0,1),\]
with \(G(u)=\frac1{a}\log\cosh(au)\)
(the logcosh contrast, robust) or \(G(u)=-e^{-u^2/2}\) (the exp
contrast, better for very heavy tails).
The FastICA fixed-point iteration for one component is
\[\mathbf{w}^{+}\leftarrow E\big[\tilde{\mathbf{x}}\,g(\mathbf{w}^\top\tilde{\mathbf{x}})\big]-E\big[g'(\mathbf{w}^\top\tilde{\mathbf{x}})\big]\mathbf{w}, \qquad \mathbf{w}\leftarrow \frac{\mathbf{w}^{+}}{\|\mathbf{w}^{+}\|},\]
with \(g=G'\). It is an approximate Newton step and converges cubically. Multiple components are extracted either in parallel (with symmetric orthogonalization each sweep) or by deflation (projecting out already-found directions).
Cost: \(O(np^2)\) for whitening, then \(O(\text{iterations}\times npk)\).
11.4 ICA in practice
set.seed(67)
n_s <- 5000
S_true <- cbind(runif(n_s), runif(n_s)) # uniform sources: NOT Gaussian
A_mix <- matrix(c(1, -1, 1, 3), 2, 2, byrow = TRUE)
Xmix <- S_true %*% A_mix
c(source_correlation = cor(S_true[, 1], S_true[, 2]),
mixed_correlation = cor(Xmix[, 1], Xmix[, 2]),
theoretical_mixed = -2 / sqrt(20))#> source_correlation mixed_correlation theoretical_mixed
#> -0.0033195 0.4421323 -0.4472136Independent uniform sources become strongly correlated after mixing, that correlation is the signature ICA undoes.
library(fastICA)
set.seed(71)
ica_fit <- fastICA(Xmix, n.comp = 2, alg.typ = "parallel", fun = "logcosh",
alpha = 1, row.norm = FALSE, maxit = 300, tol = 1e-6)
round(cor(ica_fit$S), 6) # recovered components: near-diagonal#> [,1] [,2]
#> [1,] 1 0
#> [2,] 0 1#> [,1] [,2]
#> [1,] 0.0024 1e+00
#> [2,] 1.0000 9e-04Read the second matrix carefully: each recovered component correlates near 1 with exactly one true source, but which one is arbitrary, that is the permutation ambiguity in action.
ica_df <- bind_rows(
data.frame(x = S_true[, 1], y = S_true[, 2], stage = "1. True sources (independent)"),
data.frame(x = Xmix[, 1], y = Xmix[, 2], stage = "2. Observed mixtures (correlated)"),
data.frame(x = ica_fit$X[, 1], y = ica_fit$X[, 2], stage = "3. Whitened"),
data.frame(x = ica_fit$S[, 1], y = ica_fit$S[, 2], stage = "4. ICA components (recovered)"))
ggplot(ica_df[sample(nrow(ica_df), 8000), ], aes(x, y)) +
geom_point(alpha = 0.10, size = 0.5, colour = "steelblue") +
facet_wrap(~ stage, scales = "free", nrow = 2) +
labs(title = "ICA recovers the parallelogram structure that mixing destroyed",
subtitle = "Whitening makes the cloud spherical; the final rotation aligns the axes with the sources",
x = NULL, y = NULL) +
theme_dspa(10)
# --- Interactive equivalent ------------------------------------------------
plot_ly() |>
add_markers(x = ica_fit$X[, 1], y = ica_fit$X[, 2], name = "Pre-processed data",
marker = list(color = "green", opacity = 0.3, symbol = 105)) |>
add_markers(x = ica_fit$S[, 1], y = ica_fit$S[, 2], name = "ICA components",
marker = list(color = "blue", opacity = 0.3, symbol = 5)) |>
layout(title = "Original (pre-processed) data and the corresponding ICA transform",
xaxis = list(title = "Component 1", scaleanchor = "y"),
yaxis = list(title = "Component 2"),
legend = list(orientation = "h"))Applied to real data:
#> Top_of_SN_Voxel_Intensity_Ratio
#> Top_of_SN_Voxel_Intensity_Ratio 1.000
#> Side_of_SN_Voxel_Intensity_Ratio 0.547
#> Part_IA -0.101
#> Part_IB -0.270
#> Part_II -0.044
#> Part_III -0.339
#> Side_of_SN_Voxel_Intensity_Ratio Part_IA
#> Top_of_SN_Voxel_Intensity_Ratio 0.547 -0.101
#> Side_of_SN_Voxel_Intensity_Ratio 1.000 -0.216
#> Part_IA -0.216 1.000
#> Part_IB -0.444 0.491
#> Part_II -0.377 0.504
#> Part_III -0.523 0.585
#> Part_IB Part_II Part_III
#> Top_of_SN_Voxel_Intensity_Ratio -0.270 -0.044 -0.339
#> Side_of_SN_Voxel_Intensity_Ratio -0.444 -0.377 -0.523
#> Part_IA 0.491 0.504 0.585
#> Part_IB 1.000 0.580 0.674
#> Part_II 0.580 1.000 0.639
#> Part_III 0.674 0.639 1.000set.seed(73)
ica_pd <- fastICA(as.matrix(pd_sub), n.comp = ncol(pd_sub), alg.typ = "parallel",
fun = "logcosh", alpha = 1, row.norm = FALSE,
maxit = 300, tol = 1e-6)
round(cor(ica_pd$S), 4)#> [,1] [,2] [,3] [,4] [,5] [,6]
#> [1,] 1 0 0 0 0 0
#> [2,] 0 1 0 0 0 0
#> [3,] 0 0 1 0 0 0
#> [4,] 0 0 0 1 0 0
#> [5,] 0 0 0 0 1 0
#> [6,] 0 0 0 0 0 1#> max_offdiagonal_correlation
#> 1.15574e-14The recovered components are essentially uncorrelated. Note that uncorrelatedness is necessary but not sufficient for independence, that is precisely why ICA optimizes a higher-order criterion rather than stopping at whitening.
12 Factor analysis
12.1 The model
Factor analysis is a generative latent-variable model, not a rotation. For a \(p\)-vector \(\mathbf{x}\) with mean \(\boldsymbol\mu\) and \(k<p\) latent factors:
\[\mathbf{x}-\boldsymbol\mu=L\mathbf{f}+\boldsymbol\varepsilon,\]
with \(L\) the \(p\times k\) loading matrix, \(\mathbf{f}\sim N(\mathbf{0},I_k)\) the common factors, and \(\boldsymbol\varepsilon\sim N(\mathbf{0},\Psi)\) with \(\Psi=\operatorname{diag}(\psi_1,\dots,\psi_p)\) the unique variances, independent of \(\mathbf{f}\). This implies a structured covariance:
\[\boxed{\;\Sigma=LL^\top+\Psi\;}\]
For each variable \(j\),
\[\underbrace{\sigma_{jj}}_{\text{total variance}}=\underbrace{\sum_{m=1}^{k}\ell_{jm}^2}_{\text{communality }h_j^2}+\underbrace{\psi_j}_{\text{uniqueness}}, \qquad h_j^2+\psi_j=1\ \text{ when standardized.}\]
The communality is the share of a variable’s variance explained by the common factors; the uniqueness is what is left, including measurement error.
Common misconception: “factor analysis is a generalization of PCA.” They are different objects answering different questions.
PCA Factor analysis Nature Deterministic rotation Generative probability model Error term None \(\boldsymbol\varepsilon\) with diagonal \(\Psi\) Direction Components are functions of the variables Variables are functions of the factors Fit test None available Likelihood-ratio \(\chi^2\) Solution Unique (up to sign) Unique only up to an orthogonal rotation of \(L\) Estimation Eigen/SVD, closed form Iterative maximum likelihood If anything the containment runs the other way: probabilistic PCA is the special case \(\Psi=\sigma^2I\), so PCA is a restricted factor model.
Rotational indeterminacy is intrinsic: for any orthogonal \(R\), \((LR)(LR)^\top=LRR^\top L^\top=LL^\top\), so \(L\) and \(LR\) fit identically. This is not a defect but a licence, since all rotations fit equally, choose the one that is most interpretable. Varimax maximizes the variance of squared loadings within columns, driving each toward a few large and many near-zero entries. Promax and oblimin allow correlated factors.
Heywood cases, estimated uniquenesses at or below
zero, signal an over-specified model, too few observations, or a
variable that is nearly a pure factor. factanal() warns
about them and they should never be ignored.
12.2 Choosing the number of factors
library(psych)
set.seed(79)
fa_par <- fa.parallel(pd_sub, fa = "fa", n.iter = 200, plot = FALSE,
show.legend = FALSE)#> Parallel analysis suggests that the number of factors = 1 and the number of components = NA#> suggested_factors
#> 1parallel_df <- data.frame(
factor_number = seq_along(fa_par$fa.values),
observed = fa_par$fa.values,
simulated = fa_par$fa.sim)
parallel_df |>
pivot_longer(-factor_number, names_to = "series", values_to = "eigenvalue") |>
ggplot(aes(factor_number, eigenvalue, colour = series)) +
geom_line(linewidth = 1) + geom_point(size = 2.4) +
geom_hline(yintercept = 0, linetype = "dashed", colour = "grey50") +
scale_x_continuous(breaks = parallel_df$factor_number) +
scale_colour_manual(values = c(observed = "steelblue", simulated = "firebrick")) +
labs(title = "Parallel analysis for factor retention",
subtitle = "Retain factors whose observed eigenvalue exceeds the simulated null",
x = "Factor", y = "Eigenvalue", colour = NULL) +
theme_dspa()
# --- Interactive equivalent ------------------------------------------------
plot_ly() |>
add_trace(y = fa_par$fa.values, type = "scatter", mode = "lines+markers",
name = "Observed eigenvalues",
marker = list(size = 12, symbol = 5)) |>
add_trace(y = fa_par$fa.sim, type = "scatter", mode = "lines+markers",
name = "Parallel analysis (simulated)",
marker = list(size = 12, symbol = 2)) |>
layout(title = "Scree plot with parallel analysis",
xaxis = list(title = "Factor"), yaxis = list(title = "Eigenvalue"),
legend = list(orientation = "h"))A caution on the likelihood-ratio test.
factanal()reports a \(\chi^2\) test of \(H_0\): “\(k\) factors are sufficient”. A large \(p\)-value means the data are consistent with \(k\) factors, it is not evidence that \(k\) is correct, and with small \(n\) the test has little power to detect misfit. Incrementing \(k\) until the test stops rejecting is model selection by repeated testing: the final \(p\)-value no longer has its nominal meaning and \(k\) is biased upward. Prefer parallel analysis, information criteria, or cross-validation, and use the LR test as a diagnostic rather than a selector.
#>
#> Call:
#> factanal(x = pd_sub, factors = 2, rotation = "varimax")
#>
#> Uniquenesses:
#> Top_of_SN_Voxel_Intensity_Ratio Side_of_SN_Voxel_Intensity_Ratio
#> 0.018 0.534
#> Part_IA Part_IB
#> 0.571 0.410
#> Part_II Part_III
#> 0.392 0.218
#>
#> Loadings:
#> Factor1 Factor2
#> Top_of_SN_Voxel_Intensity_Ratio 0.991
#> Side_of_SN_Voxel_Intensity_Ratio -0.417 0.540
#> Part_IA 0.650
#> Part_IB 0.726 -0.251
#> Part_II 0.779
#> Part_III 0.825 -0.318
#>
#> Factor1 Factor2
#> SS loadings 2.412 1.445
#> Proportion Var 0.402 0.241
#> Cumulative Var 0.402 0.643
#>
#> Test of the hypothesis that 2 factors are sufficient.
#> The chi square statistic is 1.35 on 4 degrees of freedom.
#> The p-value is 0.854# Communalities and uniquenesses must sum to 1 on standardized variables
data.frame(variable = rownames(fa_fit$loadings),
communality = round(1 - fa_fit$uniquenesses, 4),
uniqueness = round(fa_fit$uniquenesses, 4),
sum = round(1 - fa_fit$uniquenesses + fa_fit$uniquenesses, 4))#> variable communality
#> Top_of_SN_Voxel_Intensity_Ratio Top_of_SN_Voxel_Intensity_Ratio 0.9824
#> Side_of_SN_Voxel_Intensity_Ratio Side_of_SN_Voxel_Intensity_Ratio 0.4660
#> Part_IA Part_IA 0.4286
#> Part_IB Part_IB 0.5903
#> Part_II Part_II 0.6078
#> Part_III Part_III 0.7817
#> uniqueness sum
#> Top_of_SN_Voxel_Intensity_Ratio 0.0176 1
#> Side_of_SN_Voxel_Intensity_Ratio 0.5340 1
#> Part_IA 0.5714 1
#> Part_IB 0.4097 1
#> Part_II 0.3922 1
#> Part_III 0.2183 1#> any_heywood
#> FALSE# Rotations fit identically; they differ only in interpretability
fa_none <- factanal(pd_sub, factors = 2, rotation = "none")
fa_varimax <- factanal(pd_sub, factors = 2, rotation = "varimax")
fa_promax <- factanal(pd_sub, factors = 2, rotation = "promax")
c(loglik_none = fa_none$criteria[["objective"]],
loglik_varimax = fa_varimax$criteria[["objective"]],
loglik_promax = fa_promax$criteria[["objective"]],
identical_fit = isTRUE(all.equal(fa_none$criteria[["objective"]],
fa_varimax$criteria[["objective"]])))#> loglik_none loglik_varimax loglik_promax identical_fit
#> 0.0483586 0.0483586 0.0483586 1.0000000The objective is identical across rotations, confirming that rotation changes only the coordinate system, never the fit.
load_df <- as.data.frame(unclass(fa_fit$loadings))
load_df$variable <- rownames(load_df)
ggplot(load_df, aes(Factor1, Factor2)) +
geom_hline(yintercept = 0, colour = "grey70") +
geom_vline(xintercept = 0, colour = "grey70") +
geom_segment(aes(x = 0, y = 0, xend = Factor1, yend = Factor2),
arrow = arrow(length = unit(0.2, "cm")), colour = "grey45") +
geom_point(aes(colour = variable), size = 4) +
ggrepel::geom_text_repel(aes(label = variable), size = 3.2) +
scale_colour_manual(values = rainbow(nrow(load_df))) +
coord_fixed() +
labs(title = "Factor loadings after varimax rotation",
subtitle = "Variables far from the origin are well explained; direction shows which factor drives them",
x = "Factor 1", y = "Factor 2") +
theme_dspa() + theme(legend.position = "none")
# --- Interactive equivalent ------------------------------------------------
cols <- rainbow(nrow(load_df)) # note: rainbow() directly, NOT palette()
plot_ly(load_df, x = ~Factor1, y = ~Factor2, text = ~variable) |>
add_markers(marker = list(size = 18, color = cols)) |>
add_text(textfont = list(size = 14, color = cols), textposition = "top right") |>
layout(title = "Two-factor solution",
xaxis = list(title = "Factor 1", zeroline = TRUE),
yaxis = list(title = "Factor 2", zeroline = TRUE),
showlegend = FALSE)12.3 PCA, ICA, and FA side by side
| PCA | ICA | FA | |
|---|---|---|---|
| Assumes | Nothing distributional; informative structure is linear | Sources mutually independent; at most one Gaussian | \(\mathbf{x}=L\mathbf{f}+\boldsymbol\varepsilon\); \(\mathbf{f},\boldsymbol\varepsilon\) Gaussian; \(\Psi\) diagonal |
| Optimizes | Variance along orthogonal directions (2nd moments) | Non-Gaussianity / minimum mutual information (higher moments) | Likelihood of the structured covariance \(LL^\top+\Psi\) |
| Components are | Orthogonal, ranked by variance | Independent, not orthogonal, unranked | Latent causes, identified only up to rotation |
| Uniqueness | Unique up to sign | Up to permutation and scale | Up to any orthogonal rotation |
| Noise model | None | None (noise-free mixing) | Explicit, per-variable (\(\Psi\)) |
| Cost | \(O(np\min(n,p))\) | \(O(np^2 + Tnpk)\) | \(O(T p^3)\) |
| Use for | Compression, decorrelation, visualization | Source separation: EEG/MEG, fMRI, audio, artifact removal | Testing a theory of latent constructs |
Choose PCA to compress and decorrelate; ICA to separate physically distinct generating processes; FA to test a hypothesis about unobserved constructs.
13 PART III: NONLINEAR METHODS
14 When linear projection fails
Every method so far seeks a linear subspace. That works when the data lie near a flat sheet. When the manifold is curved, no linear projection can unroll it, and the failure is not subtle.
The canonical counterexample is the Swiss roll: a two-dimensional sheet rolled into a spiral in \(\mathbb{R}^3\). Its intrinsic dimension is 2, but the two dimensions are related to the ambient coordinates nonlinearly.
make_swiss_roll <- function(n = 1500, noise = 0.05, seed = 83) {
set.seed(seed)
t_par <- 1.5 * pi * (1 + 2 * runif(n)) # position ALONG the roll
h <- 21 * runif(n) # position ACROSS the roll
X <- cbind(t_par * cos(t_par), h, t_par * sin(t_par)) +
matrix(rnorm(3 * n, sd = noise), n, 3)
list(X = X, colour = t_par, height = h)
}
sr <- make_swiss_roll()
dim(sr$X)#> [1] 1500 3The roll must be rotated to be understood, so this figure is interactive:
plot_ly(x = sr$X[, 1], y = sr$X[, 2], z = sr$X[, 3],
type = "scatter3d", mode = "markers",
marker = list(size = 2.6, color = sr$colour, colorscale = "Viridis",
showscale = TRUE,
colorbar = list(title = "Position along roll"))) |>
layout(title = "Swiss roll: a 2-D sheet curved through 3-D space",
scene = list(xaxis = list(title = "X"), yaxis = list(title = "Y"),
zaxis = list(title = "Z"), aspectmode = "data"))Colour encodes position along the roll. A successful 2-D embedding should produce a smooth colour gradient with no mixing, points far apart along the sheet should stay far apart. PCA cannot do this:
pca_sr <- prcomp(sr$X, center = TRUE, scale. = FALSE)
sr_pca <- data.frame(PC1 = pca_sr$x[, 1], PC2 = pca_sr$x[, 2],
colour = sr$colour)
ggplot(sr_pca, aes(PC1, PC2, colour = colour)) +
geom_point(size = 1.1, alpha = 0.8) +
scale_colour_viridis_c(name = "Position\nalong roll") +
coord_fixed() +
labs(title = "PCA cannot unroll the Swiss roll",
subtitle = sprintf("PC1 and PC2 capture %.0f%% of variance, yet distant parts of the sheet overlap",
100 * sum(pca_sr$sdev[1:2]^2) / sum(pca_sr$sdev^2)),
x = "PC1", y = "PC2") +
theme_dspa()
The colours interleave: dark-purple and bright-yellow points, opposite ends of the sheet, land on top of each other. PCA has projected the roll onto its shadow, and a shadow cannot distinguish the near side of a spiral from the far side. The problem is not that PCA is badly tuned; it is that no linear map can do this.
The reason is that PCA preserves Euclidean distance, while the meaningful distance on a curved manifold is geodesic, measured along the surface. Two points separated by 1 unit through the air may be 30 units apart along the sheet.
# Euclidean vs. approximate geodesic distance for two points on opposite layers
i <- which.min(sr$colour); j <- which.max(sr$colour)
euclid_ij <- sqrt(sum((sr$X[i, ] - sr$X[j, ])^2))
# Geodesic approximated by shortest path in a k-NN graph
library(igraph)
knn_graph <- function(X, k = 10) {
Dm <- as.matrix(dist(X)); n <- nrow(Dm)
A <- matrix(0, n, n)
for (r in seq_len(n)) {
nb <- order(Dm[r, ])[2:(k + 1)]
A[r, nb] <- Dm[r, nb]; A[nb, r] <- Dm[nb, r]
}
graph_from_adjacency_matrix(A, mode = "undirected", weighted = TRUE)
}
g <- knn_graph(sr$X, k = 10)
geo_ij <- distances(g, v = i, to = j)[1, 1]
c(euclidean = round(euclid_ij, 2),
geodesic = round(geo_ij, 2),
ratio = round(geo_ij / euclid_ij, 2))#> euclidean geodesic ratio
#> 19.52 91.51 4.69Every nonlinear method below is, at bottom, a strategy for respecting geodesic rather than Euclidean structure.
15 Kernel PCA
Kernel PCA performs ordinary PCA in a high-dimensional feature space \(\mathcal{F}\) reached by a nonlinear map \(\phi:\mathbb{R}^p\to\mathcal{F}\) — without ever computing \(\phi\).
The kernel trick: PCA depends on the data only through inner products, so replace \(\langle\mathbf{x}_i,\mathbf{x}_j\rangle\) with a kernel \(K_{ij}=k(\mathbf{x}_i,\mathbf{x}_j)=\langle\phi(\mathbf{x}_i),\phi(\mathbf{x}_j)\rangle\). Common choices:
\[k_{\text{RBF}}(\mathbf{x},\mathbf{y})=\exp\!\big(-\sigma\|\mathbf{x}-\mathbf{y}\|^2\big), \qquad k_{\text{poly}}(\mathbf{x},\mathbf{y})=\big(\langle\mathbf{x},\mathbf{y}\rangle+c\big)^d .\]
Feature-space data must be centred, which cannot be done directly on \(\phi\). It can be done on \(K\):
\[\tilde K=K-\mathbf{1}_nK-K\mathbf{1}_n+\mathbf{1}_nK\mathbf{1}_n, \qquad \mathbf{1}_n=\tfrac1n\mathbf{1}\mathbf{1}^\top .\]
Eigendecomposing \(\tilde K=\alpha\Lambda\alpha^\top\), the projection of point \(i\) onto component \(m\) is \(\sqrt{\lambda_m}\,\alpha_{im}\).
Cost: \(O(n^2)\) memory for \(K\) and \(O(n^3)\) for the eigendecomposition, the same wall as classical MDS. Kernel PCA does not scale past \(n\approx10^4\) without approximation (Nyström, random Fourier features).
library(kernlab)
set.seed(89)
sub <- sample(nrow(sr$X), 800) # subsample: kPCA is O(n^3)
Xsr <- sr$X[sub, ]; col_sr <- sr$colour[sub]
kpca_fit <- kpca(Xsr, kernel = "rbfdot", kpar = list(sigma = 0.008), features = 3)
kp <- data.frame(KPC1 = rotated(kpca_fit)[, 1], KPC2 = rotated(kpca_fit)[, 2],
colour = col_sr)
ggplot(kp, aes(KPC1, KPC2, colour = colour)) +
geom_point(size = 1.2, alpha = 0.85) +
scale_colour_viridis_c(name = "Position\nalong roll") +
labs(title = "Kernel PCA with an RBF kernel",
subtitle = "Curved components separate the layers that linear PCA superimposed",
x = "Kernel PC1", y = "Kernel PC2") +
theme_dspa()
# --- Interactive equivalent (three kernel components) ----------------------
plot_ly(x = rotated(kpca_fit)[, 1], y = rotated(kpca_fit)[, 2],
z = rotated(kpca_fit)[, 3], type = "scatter3d", mode = "markers",
marker = list(size = 3, color = col_sr, colorscale = "Viridis")) |>
layout(title = "Kernel PCA embedding in 3-D",
scene = list(xaxis = list(title = "KPC1"), yaxis = list(title = "KPC2"),
zaxis = list(title = "KPC3")))Kernel PCA is the conceptual bridge: still an eigenproblem, still linear somewhere, but nonlinear in the original coordinates. The methods that follow abandon the eigenproblem entirely in favour of explicit optimization of a neighbourhood-preservation objective.
16 t-distributed stochastic neighbor embedding
t-SNE (van der Maaten & Hinton, 2008) embeds high-dimensional points into 2-D or 3-D by matching two probability distributions over pairs: one describing neighbourhoods in the original space, one in the embedding.
16.1 Step 1: similarities in the original space
For each point \(x_i\), define a conditional probability that \(x_i\) would pick \(x_j\) as its neighbour, under a Gaussian centred at \(x_i\):
\[p_{j\mid i}=\frac{\exp\!\big(-\|x_i-x_j\|^2/2\sigma_i^2\big)}{\sum_{k\ne i}\exp\!\big(-\|x_i-x_k\|^2/2\sigma_i^2\big)}, \qquad p_{i\mid i}=0 .\]
Symmetrize into a joint distribution over the \(N\) points:
\[p_{ij}=\frac{p_{j\mid i}+p_{i\mid j}}{2N}, \qquad \sum_{i\ne j}p_{ij}=1 .\]
Perplexity sets each bandwidth \(\sigma_i\). For a discrete distribution \(p\),
\[\mathrm{Perp}(p_i)=2^{H(p_i)},\qquad H(p_i)=-\sum_j p_{j\mid i}\log_2 p_{j\mid i},\]
and t-SNE binary-searches \(\sigma_i\) so that \(\mathrm{Perp}(p_i)\) equals a user-specified target, a smooth measure of the effective number of neighbours, typically 5–50. Because \(\sigma_i\) adapts per point, dense regions get narrow kernels and sparse regions wide ones.
16.2 Step 2: similarities in the embedding
In the low-dimensional map, use a Student-\(t\) with one degree of freedom — the Cauchy density \(f(z)\propto(1+z^2)^{-1}\), normalized over all pairs:
\[q_{ij}=\frac{\big(1+\|y_i-y_j\|^2\big)^{-1}}{\sum_{k\ne\ell}\big(1+\|y_k-y_\ell\|^2\big)^{-1}} .\]
Why a heavy tail? The volume of an \(N\)-ball of radius \(r\) is \(V_N(r)=\frac{\pi^{N/2}}{\Gamma(N/2+1)}r^N\). In high dimensions almost all of that volume sits near the surface: half of it lies outside radius \(2^{-1/N}r\approx r\). So a high-dimensional neighbourhood has far more “room” at moderate distance than a 2-D one does. Using a Gaussian in the map would force moderately-separated points implausibly close together, the crowding problem. The \(t\)-distribution’s heavier tail lets a given similarity be realized at a larger map distance, relieving the crush.
16.3 Step 3: minimize the divergence
Match \(Q\) to \(P\) by minimizing the Kullback–Leibler divergence
\[C=\mathrm{KL}(P\,\|\,Q)=\sum_{i\ne j}p_{ij}\log\frac{p_{ij}}{q_{ij}} .\]
The asymmetry matters. A large \(p_{ij}\) modelled by a small \(q_{ij}\) costs a lot; a small \(p_{ij}\) modelled by a large \(q_{ij}\) costs little. So t-SNE works hard to keep near points near and is comparatively indifferent to where far points go, which is exactly why between-cluster distances in a t-SNE plot are not interpretable.
The gradient with respect to the embedding coordinate \(y_i\) is
\[\boxed{\;\frac{\partial C}{\partial y_i}=4\sum_{j\ne i}\big(p_{ij}-q_{ij}\big)\big(1+\|y_i-y_j\|^2\big)^{-1}\big(y_i-y_j\big)\;}\]
Every term is a spring between \(y_i\) and \(y_j\): attractive when \(p_{ij}>q_{ij}\), repulsive otherwise, with stiffness falling off as the map distance grows. Note the argument is \(\|y_i-y_j\|\), the map distance. The input distances enter only through the fixed \(p_{ij}\).
Early exaggeration multiplies \(p_{ij}\) by a factor (typically 12) for the first ~250 iterations, forcing tight clusters to form before fine structure is resolved. Without it the optimization frequently stalls in poor local minima.
Cost. The exact gradient is \(O(N^2)\) per iteration. Barnes–Hut t-SNE
approximates the repulsive term with a quadtree, giving \(O(N\log N)\); FIt-SNE uses interpolation
and an FFT for \(O(N)\).
Rtsne() defaults to Barnes–Hut with \(\theta=0.5\).
16.4 Perplexity in practice
library(Rtsne)
set.seed(97)
sr_sub <- sample(nrow(sr$X), 1200)
Xt <- sr$X[sr_sub, ]; col_t <- sr$colour[sr_sub]
perps <- c(5, 30, 100)
tsne_runs <- lapply(perps, function(pp) {
set.seed(101)
r <- Rtsne(Xt, dims = 2, perplexity = pp, max_iter = 750, verbose = FALSE,
check_duplicates = FALSE)
data.frame(x = r$Y[, 1], y = r$Y[, 2], colour = col_t,
perplexity = paste("Perplexity =", pp))
})
bind_rows(tsne_runs) |>
ggplot(aes(x, y, colour = colour)) +
geom_point(size = 0.9, alpha = 0.85) +
scale_colour_viridis_c(name = "Position\nalong roll") +
facet_wrap(~ perplexity, nrow = 1, scales = "free") +
labs(title = "t-SNE unrolls the Swiss roll, and perplexity controls the scale",
subtitle = "Low perplexity fragments the sheet; high perplexity blurs local structure",
x = NULL, y = NULL) +
theme_dspa(10)
# --- Interactive equivalent (3-D embedding) --------------------------------
set.seed(101)
tsne3 <- Rtsne(Xt, dims = 3, perplexity = 30, max_iter = 750,
check_duplicates = FALSE)
plot_ly(x = tsne3$Y[, 1], y = tsne3$Y[, 2], z = tsne3$Y[, 3],
type = "scatter3d", mode = "markers",
marker = list(size = 3, color = col_t, colorscale = "Viridis")) |>
layout(title = "t-SNE 3-D embedding",
scene = list(xaxis = list(title = ""), yaxis = list(title = ""),
zaxis = list(title = "")))16.5 What a t-SNE plot does and does not tell you
Common misconception: “the clusters in my t-SNE plot are far apart, so they are very different.” Between-cluster distances in a t-SNE embedding are not meaningful. The KL objective penalizes misplacing near neighbours heavily and misplacing distant points barely at all, so the global layout is essentially arbitrary. Cluster sizes are also uninformative: the adaptive bandwidth \(\sigma_i\) expands sparse regions and compresses dense ones, so a visually large cluster may be a sparse one, not a numerous one.
Two further cautions. t-SNE will produce apparent clusters in data that has none, run it on uniform noise and see. And it is stochastic: different seeds give different maps of the same data, so set a seed and check stability across several.
set.seed(103)
noise <- matrix(rnorm(1000 * 30), 1000, 30) # pure noise, no structure at all
tsne_noise <- Rtsne(noise, dims = 2, perplexity = 30, max_iter = 750,
check_duplicates = FALSE)
ggplot(data.frame(x = tsne_noise$Y[, 1], y = tsne_noise$Y[, 2]), aes(x, y)) +
geom_point(size = 0.9, alpha = 0.5, colour = "steelblue") +
labs(title = "t-SNE applied to 30-dimensional Gaussian noise",
subtitle = "There is no structure in this data. The apparent clumping is an artifact of the algorithm",
x = NULL, y = NULL) +
theme_dspa()
t-SNE has no out-of-sample extension. The embedding
is a set of optimized coordinates, not a function, so a new observation
cannot be projected without re-running the whole optimization. This
rules t-SNE out as a preprocessing step inside a supervised pipeline,
you cannot fit on training data and transform a test set. UMAP, by
contrast, does supply a predict() method.
17 Uniform manifold approximation and projection
UMAP (McInnes, Healy & Melville, 2018) builds on Riemannian geometry and algebraic topology. Practically it resembles t-SNE, construct a weighted neighbourhood graph, then lay it out in low dimensions, but the theoretical framing and the objective differ, and it is substantially faster.
17.1 Simplicial complexes
UMAP represents the data’s topology with simplicial complexes. A \(k\)-simplex is the convex hull of \(k+1\) affinely independent points: a 0-simplex is a point, a 1-simplex an edge, a 2-simplex a triangle, a 3-simplex a tetrahedron.
plot_ly(type = "mesh3d",
x = c(0, 1/sqrt(3), -1/(2*sqrt(3)), -1/(2*sqrt(3))),
y = c(sqrt(2/3), 0, 0, 0),
z = c(0, 0, -1/2, 1/2),
i = c(0, 0, 0, 1), j = c(1, 2, 3, 2), k = c(2, 3, 1, 3),
facecolor = toRGB(viridisLite::viridis(4)),
opacity = 0.8, showscale = FALSE) |>
add_trace(x = c(-1/(2*sqrt(3)), -1/(2*sqrt(3))), y = c(0, 0), z = c(-1/2, 1/2),
type = "scatter3d", mode = "lines", showlegend = FALSE,
line = list(color = "rgb(20,20,20)", width = 14)) |>
add_trace(x = -1/(2*sqrt(3)), y = 0, z = 1/2, type = "scatter3d",
mode = "markers", showlegend = FALSE,
marker = list(size = 12, color = "blue")) |>
layout(title = "A 3-simplex (tetrahedron) with a highlighted 1-cell and 0-cell",
scene = list(xaxis = list(title = "X"), yaxis = list(title = "Y"),
zaxis = list(title = "Z")))Abstractly, a simplicial complex is a collection of sets \(X=\{X_i\}_{i\ge0}\) where each element of \(X^n\) is an \((n+1)\)-element set all of whose \(n\)-element subsets lie in \(X^{n-1}\), the faces of a simplex are also in the complex. A square in \(\mathbb{R}^2\) decomposes as
\[S_0=\{\{a\},\{b\},\{c\},\{d\}\},\quad S_1=\{\{a,b\},\{a,c\},\{a,d\},\{b,c\},\{c,d\}\},\quad S_2=\{\{a,b,c\},\{a,c,d\}\}.\]
Ordering matters for face maps: the \(i\)-th face map deletes the \(i\)-th vertex,
\[f_i:[x_0,\dots,x_i,\dots,x_n]\ \longrightarrow\ [x_0,\dots,x_{i-1},x_{i+1},\dots,x_n].\]
A fuzzy set replaces binary membership with \(\mu:A\to[0,1]\). Fuzzy simplicial sets require that a face’s membership strength be at least that of the simplex containing it, a coherence condition that makes the union operation below well defined.
The Vietoris–Rips construction builds a complex from a point cloud by connecting every pair within radius \(\epsilon\) and filling in every clique. Varying \(\epsilon\) generates a filtration, the object studied in persistent homology.
# --- Vietoris-Rips complex from a point cloud (requires reticulate + matplotlib)
# Run interactively; not evaluated during knitting so the chapter has no
# Python dependency.
library(reticulate)
py_run_string("
import numpy as np, matplotlib.pyplot as plt
def euclid(a, b): return np.linalg.norm(a - b)
def build_graph(X, eps=3.1, metric=euclid):
nodes = list(range(X.shape[0])); edges = []; weights = []
for i in range(X.shape[0]):
for j in range(i + 1, X.shape[0]):
d = metric(X[i], X[j])
if d <= eps:
edges.append({i, j}); weights.append([len(edges) - 1, d])
return nodes, edges, weights
def lower_nbrs(nodes, edges, v):
return {x for x in nodes if {x, v} in edges and v > x}
def rips(graph, k):
nodes, edges = graph[0:2]
VR = [{v} for v in nodes] + list(edges)
for i in range(k):
for simplex in [s for s in VR if len(s) == i + 2]:
for nb in set.intersection(*[lower_nbrs(nodes, edges, z) for z in simplex]):
VR.append(set.union(simplex, {nb}))
return VR
")
set.seed(107)
n_pt <- 30; theta <- seq(0, pi, length.out = n_pt); r <- 5
cloud <- cbind(r * cos(theta) + rnorm(n_pt, 0, 0.4),
r * sin(theta) + rnorm(n_pt, 0, 0.4))
py$X <- r_to_py(cloud)
py_run_string("VR = rips(build_graph(np.array(X), eps=3.1), 2); print(len(VR))")17.2 The UMAP objective
UMAP approximates the manifold locally. For each point \(x_i\), let \(\rho_i\) be the distance to its nearest neighbour and \(\sigma_i\) a local scale. The unsymmetrized edge weight is
\[v_{j\mid i}=\exp\!\left(-\frac{\max(0,\ r_{ij}-\rho_i)}{\sigma_i}\right),\]
with \(\sigma_i\) solved from the
constraint \(\sum_j v_{j\mid
i}=\log_2(2k)\) for \(k\) =
n_neighbors. Subtracting \(\rho_i\) enforces local
connectivity: every point is connected to its nearest neighbour
with weight 1, so no point is left isolated.
Symmetrize by fuzzy set union:
\[v_{ij}=v_{j\mid i}+v_{i\mid j}-v_{j\mid i}v_{i\mid j}, \qquad\text{in matrix form}\qquad V_{\text{sym}}=V+V^\top-V\circ V^\top,\]
where \(\circ\) is the Hadamard product.
In the embedding, weights come from a smooth family
\[w_{ij}=\frac{1}{1+a\,d_{ij}^{2b}},\]
with \(a,b\) fitted by nonlinear
least squares to approximate a piecewise function governed by
min_dist and spread; both typically lie in
\([0.5,5]\).
The objective is cross-entropy between the two fuzzy sets:
\[C_{\mathrm{UMAP}}=\sum_{i,j}\Bigg[\underbrace{v_{ij}\log\frac{v_{ij}}{w_{ij}}}_{\text{attractive}} +\underbrace{(1-v_{ij})\log\frac{1-v_{ij}}{1-w_{ij}}}_{\text{repulsive}}\Bigg].\]
The second term is the key structural difference from t-SNE. KL divergence contains only the first; cross-entropy adds an explicit repulsive term penalizing pairs that are close in the map but distant in the data. That is why UMAP tends to preserve more global structure.
The stochastic-gradient updates are
\[\frac{\partial C^{+}}{\partial y_i}=\frac{-2ab\,d_{ij}^{2(b-1)}}{1+a\,d_{ij}^{2b}}\,(y_i-y_j), \qquad \frac{\partial C^{-}}{\partial y_i}=\frac{2b}{\big(\epsilon+d_{ij}^2\big)\big(1+a\,d_{ij}^{2b}\big)}\,(y_i-y_j),\]
with \(\epsilon\approx10^{-3}\) preventing division by zero. Repulsion is applied by negative sampling, a handful of random non-neighbours per edge per epoch, rather than over all pairs, which is where most of UMAP’s speed advantage comes from.
Cost. Approximate \(k\)-NN via NN-descent is roughly \(O(n^{1.14})\) in practice; layout optimization is \(O(n\cdot\text{epochs}\cdot k)\). UMAP routinely handles \(10^6\) points where exact t-SNE would not finish.
library(umap)
run_umap2 <- function(X, n_neighbors, min_dist, seed = 109) {
cfg <- umap.defaults
cfg$n_neighbors <- n_neighbors
cfg$min_dist <- min_dist
cfg$n_components <- 2
cfg$random_state <- seed
umap(X, config = cfg)
}
settings <- list(c(5, 0.1), c(15, 0.1), c(50, 0.5))
umap_runs <- lapply(settings, function(s) {
u <- run_umap2(Xt, s[1], s[2])
data.frame(x = u$layout[, 1], y = u$layout[, 2], colour = col_t,
setting = sprintf("n_neighbors = %d, min_dist = %.1f", s[1], s[2]))
})
bind_rows(umap_runs) |>
ggplot(aes(x, y, colour = colour)) +
geom_point(size = 0.9, alpha = 0.85) +
scale_colour_viridis_c(name = "Position\nalong roll") +
facet_wrap(~ setting, nrow = 1, scales = "free") +
labs(title = "UMAP on the Swiss roll",
subtitle = "n_neighbors trades local against global structure; min_dist controls how tightly points pack",
x = NULL, y = NULL) +
theme_dspa(10)
# --- Interactive equivalent (3-D embedding) --------------------------------
cfg3 <- umap.defaults; cfg3$n_components <- 3; cfg3$n_neighbors <- 15
cfg3$random_state <- 109
u3 <- umap(Xt, config = cfg3)
plot_ly(x = u3$layout[, 1], y = u3$layout[, 2], z = u3$layout[, 3],
type = "scatter3d", mode = "markers",
marker = list(size = 3, color = col_t, colorscale = "Viridis")) |>
layout(title = "UMAP 3-D embedding",
scene = list(xaxis = list(title = ""), yaxis = list(title = ""),
zaxis = list(title = "")))17.3 UMAP hyperparameters
The umap() configuration exposes about twenty
parameters. Five matter most:
| Parameter | Effect |
|---|---|
n_neighbors |
Size of the local neighbourhood. Small → local detail, fragmented; large → global structure, blurred detail |
min_dist |
Minimum separation in the embedding. Small → tight clumps; large → evenly spread |
n_components |
Target dimension |
metric |
Distance in the original space: euclidean,
manhattan, cosine, pearson. Some
violate the triangle inequality, degrading the \(k\)-NN search |
random_state |
Seed — always set it |
Others: n_epochs (optimization iterations),
init ("spectral", the default, initializes
from graph-Laplacian eigenvectors, generally better than
"random"), set_op_mix_ratio (union
vs. intersection in the fuzzy symmetrization),
local_connectivity, bandwidth,
alpha and gamma (learning rate),
negative_sample_rate, a and b
(estimated from min_dist/spread when
NA), transform_state (seed for
predict()), knn_repeats, and
umap_learn_args for the Python backend.
Reproducibility versus replicability. Reproducibility is the strong condition: the same data, code, seeds, and software versions yield identical results. Replicability is weaker: an independent study of the same phenomenon, with its own data, reaches consistent conclusions. UMAP and t-SNE are stochastic, so reproducibility requires an explicit seed (
random_state), and replicability requires checking that conclusions survive re-running with different seeds and hyperparameters. An embedding that changes qualitatively under reseeding is not a finding.
set.seed(113)
seeds <- c(1, 2, 3)
stab <- lapply(seeds, \(s) run_umap2(Xt, 15, 0.1, seed = s)$layout)
# Compare embeddings by the correlation of their pairwise distance matrices
dcor <- function(a, b) cor(as.numeric(dist(a)), as.numeric(dist(b)))
data.frame(pair = c("seed 1 vs 2", "seed 1 vs 3", "seed 2 vs 3"),
distance_correlation = round(c(dcor(stab[[1]], stab[[2]]),
dcor(stab[[1]], stab[[3]]),
dcor(stab[[2]], stab[[3]])), 4))#> pair distance_correlation
#> 1 seed 1 vs 2 0.8944
#> 2 seed 1 vs 3 0.8862
#> 3 seed 2 vs 3 0.8952High correlation across seeds means the geometry is stable even though the raw coordinates differ, that is the check worth reporting.
17.4 What must not be read off a UMAP plot
- Hyperparameters change the picture. There is no universally optimal setting; scan a range and report stability.
- Cluster sizes are not meaningful. UMAP normalizes local density, so a large blob is not necessarily a numerous one.
- Between-cluster distances are not meaningful, less badly than t-SNE, thanks to the repulsive term and spectral initialization, but still not to be measured.
- Noise can cluster. Random variation is not preserved as random variation; spurious groupings appear.
- UMAP does support out-of-sample projection via
predict(), unlike t-SNE, a decisive practical advantage in supervised pipelines.
18 Measuring embedding quality
“Does it look good” is not a criterion. Two standard rank-based measures turn neighbourhood preservation into numbers.
Let \(r(i,j)\) be the rank of \(j\) among \(i\)’s neighbours in the original space, \(\hat r(i,j)\) the rank in the embedding, \(U_k(i)\) the points in \(i\)’s \(k\)-nearest set in the embedding but not the original, and \(V_k(i)\) the reverse.
\[T(k)=1-\frac{2}{Nk(2N-3k-1)}\sum_{i=1}^{N}\sum_{j\in U_k(i)}\big(r(i,j)-k\big) \qquad\textbf{(trustworthiness)}\]
\[C(k)=1-\frac{2}{Nk(2N-3k-1)}\sum_{i=1}^{N}\sum_{j\in V_k(i)}\big(\hat r(i,j)-k\big) \qquad\textbf{(continuity)}\]
Both lie in \([0,1]\) with 1 best. Trustworthiness penalizes points that appear close in the map but are not, false neighbours, the failure that misleads a reader. Continuity penalizes true neighbours torn apart. They trade off, so report both.
trust_cont <- function(X_hi, Y_lo, k = 12) {
n <- nrow(X_hi)
Rhi <- t(apply(as.matrix(dist(X_hi)), 1, \(r) rank(r, ties.method = "first")))
Rlo <- t(apply(as.matrix(dist(Y_lo)), 1, \(r) rank(r, ties.method = "first")))
norm_c <- 2 / (n * k * (2 * n - 3 * k - 1))
tr <- co <- 0
for (i in seq_len(n)) {
knn_hi <- which(Rhi[i, ] > 1 & Rhi[i, ] <= k + 1)
knn_lo <- which(Rlo[i, ] > 1 & Rlo[i, ] <= k + 1)
U <- setdiff(knn_lo, knn_hi); V <- setdiff(knn_hi, knn_lo)
if (length(U)) tr <- tr + sum(Rhi[i, U] - 1 - k)
if (length(V)) co <- co + sum(Rlo[i, V] - 1 - k)
}
c(trustworthiness = 1 - norm_c * tr, continuity = 1 - norm_c * co)
}
set.seed(127)
qsub <- sample(nrow(Xt), 500)
Xq <- Xt[qsub, ]
emb_pca <- prcomp(Xq)$x[, 1:2]
emb_tsne <- Rtsne(Xq, dims = 2, perplexity = 30, max_iter = 750,
check_duplicates = FALSE)$Y
emb_umap <- run_umap2(Xq, 15, 0.1)$layout
emb_mds <- cmdscale(dist(Xq), k = 2)
quality <- rbind(PCA = trust_cont(Xq, emb_pca),
`Classical MDS` = trust_cont(Xq, emb_mds),
`t-SNE` = trust_cont(Xq, emb_tsne),
UMAP = trust_cont(Xq, emb_umap))
round(quality, 4)#> trustworthiness continuity
#> PCA 0.9321 0.9810
#> Classical MDS 0.9321 0.9810
#> t-SNE 0.9978 0.9785
#> UMAP 0.9959 0.9888as.data.frame(quality) |>
tibble::rownames_to_column("method") |>
pivot_longer(-method, names_to = "metric", values_to = "value") |>
ggplot(aes(reorder(method, value), value, fill = metric)) +
geom_col(position = position_dodge(width = 0.8), width = 0.7) +
coord_flip(ylim = c(0.5, 1)) +
scale_fill_manual(values = c(trustworthiness = "steelblue",
continuity = "firebrick")) +
labs(title = "Neighbourhood preservation on the Swiss roll (k = 12)",
subtitle = "Trustworthiness penalizes false neighbours; continuity penalizes torn ones",
x = NULL, y = "Score (1 is best)", fill = NULL) +
theme_dspa()
Note that classical MDS scores the same as PCA, as §4.7 established, on Euclidean distances they are the same method. The nonlinear methods win on trustworthiness, which is the metric that matters when the plot will be used to claim two points are similar.
19 PART IV: APPLICATIONS
20 Handwritten digit recognition
The MNIST collection contains \(28\times28\) grayscale images of handwritten digits, each image a point in \(\mathbb{R}^{784}\). Its intrinsic dimension is far lower: the set of images that look like a handwritten “3” is parameterized by a handful of factors (slant, thickness, loop size, stroke curvature).
zip_path <- dspa_download(
"https://www.socr.umich.edu/people/dinov/2017/Spring/DSPA_HS650/data/DigitRecognizer_TrainingData.zip",
"DigitRecognizer_TrainingData.zip")
csv_path <- unzip(zip_path, exdir = dspa_cache_dir())[1]
mnist <- read.csv(csv_path)
dim(mnist)#> [1] 42000 785#>
#> 0 1 2 3 4 5 6 7 8 9
#> 4132 4684 4177 4351 4072 3795 4137 4401 4063 4188digit_labels <- mnist$label # 0-9
pixels <- as.matrix(mnist[, -1]) / 255 # rows = images, cols = pixels
# A ten-colour scale keyed to the DIGIT, not to a rank. Note rainbow() is used
# directly: palette() would set the global palette and return the OLD one.
digit_palette <- setNames(rainbow(10), as.character(0:9))
digit_cols <- digit_palette[as.character(digit_labels)]
c(images = nrow(pixels), pixels_per_image = ncol(pixels),
colours = length(digit_palette))#> images pixels_per_image colours
#> 42000 784 10Orientation, fixed once. MNIST stores each image
row-major in a 784-vector. Reshaping with matrix(v, 28, 28)
fills column-major, producing the transpose. One helper applies the
correction, and everything downstream uses it.
as_image <- function(v) matrix(v, nrow = 28, ncol = 28, byrow = TRUE)
image_df <- function(idx) {
m <- as_image(pixels[idx, ])
expand.grid(row = 1:28, col = 1:28) |>
mutate(intensity = as.vector(m),
label = paste0("digit ", digit_labels[idx], " (case ", idx, ")"))
}bind_rows(lapply(1:8, image_df)) |>
ggplot(aes(col, 29 - row, fill = intensity)) +
geom_raster() +
scale_fill_gradient(low = "white", high = "black", guide = "none") +
facet_wrap(~ label, nrow = 2) +
coord_fixed() +
labs(title = "The first eight training images", x = NULL, y = NULL) +
theme_void(base_size = 9) +
theme(strip.text = element_text(size = 8), plot.title = element_text(face = "bold"))
# --- Interactive equivalent (single image as a heatmap) --------------------
m1 <- as_image(pixels[1, ])
plot_ly(z = ~m1[28:1, ], type = "heatmap", showscale = FALSE,
colorscale = list(c(0, "white"), c(1, "black"))) |>
layout(title = sprintf("Digit %d", digit_labels[1]),
xaxis = list(title = "", scaleanchor = "y"),
yaxis = list(title = ""))20.1 PCA on MNIST
set.seed(131)
n_use <- 10000
use_idx <- sample(nrow(pixels), n_use)
Xm <- pixels[use_idx, ]
lab_m <- digit_labels[use_idx]
col_m <- digit_cols[use_idx]
# Drop always-blank border pixels (zero variance breaks scaling)
keep_px <- apply(Xm, 2, sd) > 0
c(pixels_kept = sum(keep_px), pixels_dropped = sum(!keep_px))#> pixels_kept pixels_dropped
#> 691 93t_pca <- system.time(pca_m <- prcomp(Xm[, keep_px], center = TRUE, scale. = FALSE))
cumvar <- cumsum(pca_m$sdev^2) / sum(pca_m$sdev^2)
c(seconds = round(t_pca[["elapsed"]], 2),
components_for_50pct = which(cumvar >= 0.50)[1],
components_for_90pct = which(cumvar >= 0.90)[1],
components_for_95pct = which(cumvar >= 0.95)[1],
ambient_dimension = sum(keep_px))#> seconds components_for_50pct components_for_90pct
#> 9.83 11.00 86.00
#> components_for_95pct ambient_dimension
#> 152.00 691.00Ninety percent of the variance in a 700-dimensional pixel space is captured by roughly 90 components, an eightfold compression before any nonlinear method is applied.
data.frame(k = 1:150, cumulative = cumvar[1:150]) |>
ggplot(aes(k, cumulative)) +
geom_line(linewidth = 1, colour = "steelblue") +
geom_hline(yintercept = c(0.5, 0.9, 0.95), linetype = "dashed",
colour = "grey45") +
scale_y_continuous(labels = scales::percent) +
labs(title = "Cumulative variance explained by MNIST principal components",
subtitle = "Dashed lines at 50%, 90%, and 95%",
x = "Number of components", y = NULL) +
theme_dspa()
The leading components are interpretable as images, eigendigits, the orthogonal basis of stroke patterns from which every digit is built:
eigen_img <- function(k) {
v <- numeric(ncol(pixels)); v[keep_px] <- pca_m$rotation[, k]
m <- as_image(v)
expand.grid(row = 1:28, col = 1:28) |>
mutate(loading = as.vector(m),
comp = sprintf("PC%d (%.1f%%)", k,
100 * pca_m$sdev[k]^2 / sum(pca_m$sdev^2)))
}
bind_rows(lapply(1:8, eigen_img)) |>
mutate(comp = factor(comp, levels = unique(comp))) |>
ggplot(aes(col, 29 - row, fill = loading)) +
geom_raster() +
scale_fill_gradient2(low = "firebrick", mid = "white", high = "steelblue",
midpoint = 0, guide = "none") +
facet_wrap(~ comp, nrow = 2) + coord_fixed() +
labs(title = "The first eight eigendigits",
subtitle = "Blue and red are opposite-signed loadings; each is a stroke pattern",
x = NULL, y = NULL) +
theme_void(base_size = 9) +
theme(strip.text = element_text(size = 8), plot.title = element_text(face = "bold"))
Rank-\(k\) reconstruction makes Eckart–Young visible on real data:
recon <- function(i, k) {
z <- pca_m$x[i, 1:k, drop = FALSE] %*% t(pca_m$rotation[, 1:k, drop = FALSE])
v <- numeric(ncol(pixels))
v[keep_px] <- as.vector(z) + pca_m$center
v
}
target <- 3
ranks <- c(2, 5, 10, 25, 50, 100)
bind_rows(lapply(ranks, function(k) {
m <- as_image(recon(target, k))
expand.grid(row = 1:28, col = 1:28) |>
mutate(intensity = as.vector(m), panel = sprintf("rank %d", k))
})) |>
bind_rows(expand.grid(row = 1:28, col = 1:28) |>
mutate(intensity = as.vector(as_image(Xm[target, ])),
panel = "original")) |>
mutate(panel = factor(panel, levels = c(sprintf("rank %d", ranks), "original"))) |>
ggplot(aes(col, 29 - row, fill = intensity)) +
geom_raster() +
scale_fill_gradient(low = "white", high = "black", guide = "none") +
facet_wrap(~ panel, nrow = 1) + coord_fixed() +
labs(title = sprintf("Rank-k reconstruction of a handwritten %d", lab_m[target]),
subtitle = "Each panel is the optimal rank-k approximation, by Eckart-Young",
x = NULL, y = NULL) +
theme_void(base_size = 9) +
theme(strip.text = element_text(size = 8), plot.title = element_text(face = "bold"))
Stacking those reconstructions as surfaces makes the convergence tangible, so this one is interactive:
offset <- 1.35
p <- plot_ly()
for (idx in seq_along(c(2, 10, 50))) {
k <- c(2, 10, 50)[idx]
p <- add_surface(p, z = as_image(recon(target, k)) + (idx - 1) * offset,
showscale = FALSE, opacity = 0.95, colorscale = "Greys",
name = paste("rank", k))
}
p <- add_surface(p, z = as_image(Xm[target, ]) + 3 * offset, showscale = FALSE,
opacity = 0.95, colorscale = "Greys", name = "original")
p |> layout(title = "Rank 2, 10, 50, and the original, stacked (rotate to compare)",
scene = list(xaxis = list(title = "column"),
yaxis = list(title = "row"),
zaxis = list(title = "intensity + offset")))pca_plot <- data.frame(PC1 = pca_m$x[, 1], PC2 = pca_m$x[, 2],
digit = factor(lab_m))
ggplot(pca_plot[1:3000, ], aes(PC1, PC2, colour = digit)) +
geom_point(size = 0.7, alpha = 0.6) +
scale_colour_manual(values = digit_palette) +
guides(colour = guide_legend(override.aes = list(size = 3), nrow = 1)) +
labs(title = "MNIST in the first two principal components",
subtitle = "A linear projection separates 0 and 1 but leaves 4, 7, and 9 entangled",
x = sprintf("PC1 (%.1f%%)", 100 * pca_m$sdev[1]^2 / sum(pca_m$sdev^2)),
y = sprintf("PC2 (%.1f%%)", 100 * pca_m$sdev[2]^2 / sum(pca_m$sdev^2))) +
theme_dspa()
20.2 t-SNE and UMAP on MNIST
Both are run on the leading 50 principal components rather than raw pixels. This is standard practice: it removes noise, drops the cost of the \(k\)-NN search from 700 dimensions to 50, and typically improves the embedding.
library(Rtsne)
Xpca50 <- pca_m$x[, 1:50]
set.seed(137)
t_tsne2 <- system.time(
tsne_2d <- Rtsne(Xpca50, dims = 2, perplexity = 30, max_iter = 750,
pca = FALSE, check_duplicates = FALSE, verbose = FALSE))
set.seed(137)
t_tsne3 <- system.time(
tsne_3d <- Rtsne(Xpca50, dims = 3, perplexity = 30, max_iter = 750,
pca = FALSE, check_duplicates = FALSE, verbose = FALSE))
c(tsne_2d_seconds = round(t_tsne2[["elapsed"]], 1),
tsne_3d_seconds = round(t_tsne3[["elapsed"]], 1))#> tsne_2d_seconds tsne_3d_seconds
#> 18.1 42.8tsne_df <- data.frame(x = tsne_2d$Y[, 1], y = tsne_2d$Y[, 2],
digit = factor(lab_m))
ggplot(tsne_df, aes(x, y, colour = digit)) +
geom_point(size = 0.45, alpha = 0.7) +
scale_colour_manual(values = digit_palette) +
guides(colour = guide_legend(override.aes = list(size = 3), nrow = 1)) +
labs(title = "t-SNE embedding of MNIST (784D to 2D, via 50 PCs)",
subtitle = "Ten well-separated clusters. Their relative positions and sizes are NOT interpretable",
x = NULL, y = NULL) +
theme_dspa()
show_n <- 3000
plot_ly(x = tsne_3d$Y[1:show_n, 1], y = tsne_3d$Y[1:show_n, 2],
z = tsne_3d$Y[1:show_n, 3],
type = "scatter3d", mode = "markers",
color = factor(lab_m[1:show_n]), colors = unname(digit_palette),
text = paste("digit", lab_m[1:show_n]),
marker = list(size = 2.2, opacity = 0.75)) |>
layout(title = "t-SNE 3-D embedding of MNIST",
scene = list(xaxis = list(title = ""), yaxis = list(title = ""),
zaxis = list(title = "")))library(umap)
cfg2 <- umap.defaults
cfg2$n_neighbors <- 15; cfg2$min_dist <- 0.1
cfg2$n_components <- 2; cfg2$random_state <- 139
cfg3 <- cfg2; cfg3$n_components <- 3
t_umap2 <- system.time(umap_2d <- umap(Xpca50, config = cfg2))
t_umap3 <- system.time(umap_3d <- umap(Xpca50, config = cfg3))
c(umap_2d_seconds = round(t_umap2[["elapsed"]], 1),
umap_3d_seconds = round(t_umap3[["elapsed"]], 1),
tsne_2d_seconds = round(t_tsne2[["elapsed"]], 1))#> umap_2d_seconds umap_3d_seconds tsne_2d_seconds
#> 46.3 62.7 18.1Note that the 2-D and 3-D embeddings come from separate
runs with matching n_components. A 2-D layout is a
different optimization from a 3-D one, it is not the first two
coordinates of the 3-D result.
umap_df <- data.frame(x = umap_2d$layout[, 1], y = umap_2d$layout[, 2],
digit = factor(lab_m))
ggplot(umap_df, aes(x, y, colour = digit)) +
geom_point(size = 0.45, alpha = 0.7) +
scale_colour_manual(values = digit_palette) +
guides(colour = guide_legend(override.aes = list(size = 3), nrow = 1)) +
labs(title = "UMAP embedding of MNIST (784D to 2D, via 50 PCs)",
subtitle = "Tighter clusters than t-SNE and more global structure: 4/7/9 sit adjacent, as they should",
x = NULL, y = NULL) +
theme_dspa()
# --- Interactive equivalent: labelled text markers -------------------------
sel <- 1:1000
plot_ly(x = umap_2d$layout[sel, 1], y = umap_2d$layout[sel, 2], mode = "text") |>
add_text(text = as.character(lab_m[sel]),
textfont = list(color = col_m[sel], size = 13)) |>
layout(title = "UMAP (784D to 2D) embedding",
xaxis = list(title = ""), yaxis = list(title = ""))plot_ly(x = umap_3d$layout[1:show_n, 1], y = umap_3d$layout[1:show_n, 2],
z = umap_3d$layout[1:show_n, 3],
type = "scatter3d", mode = "markers",
color = factor(lab_m[1:show_n]), colors = unname(digit_palette),
text = paste("digit", lab_m[1:show_n]),
marker = list(size = 2.2, opacity = 0.75)) |>
layout(title = "UMAP 3-D embedding of MNIST",
scene = list(xaxis = list(title = ""), yaxis = list(title = ""),
zaxis = list(title = "")))20.3 Out-of-sample projection
UMAP learns a mapping that can be applied to new data. t-SNE cannot, its output is a set of optimized coordinates, not a function.
set.seed(149)
new_idx <- setdiff(seq_len(nrow(pixels)), use_idx)[1:1500]
X_new <- pixels[new_idx, keep_px]
lab_new <- digit_labels[new_idx]
# Project new images into the SAME PC space, then through the fitted UMAP
X_new_pc <- scale(X_new, center = pca_m$center, scale = FALSE) %*%
pca_m$rotation[, 1:50]
proj_new <- predict(umap_2d, X_new_pc)
bind_rows(
data.frame(x = umap_2d$layout[, 1], y = umap_2d$layout[, 2],
digit = factor(lab_m), set = "Training (fitted)"),
data.frame(x = proj_new[, 1], y = proj_new[, 2],
digit = factor(lab_new), set = "New images (projected)")) |>
ggplot(aes(x, y, colour = digit)) +
geom_point(aes(alpha = set, size = set)) +
scale_colour_manual(values = digit_palette) +
scale_alpha_manual(values = c("Training (fitted)" = 0.16,
"New images (projected)" = 0.9)) +
scale_size_manual(values = c("Training (fitted)" = 0.4,
"New images (projected)" = 1.1)) +
guides(colour = guide_legend(override.aes = list(size = 3, alpha = 1), nrow = 1),
alpha = "none", size = "none") +
labs(title = "1,500 unseen images projected into the fitted UMAP embedding",
subtitle = "Faint points are the training embedding; bright points are new data landing in the right clusters",
x = NULL, y = NULL) +
theme_dspa()
# A 1-NN classifier in the 2-D embedding: how well does the projection preserve class?
library(class)
pred <- knn(train = umap_2d$layout, test = proj_new, cl = factor(lab_m), k = 5)
c(accuracy_5nn_in_2D_embedding = round(mean(pred == factor(lab_new)), 4),
chance = 0.1)#> accuracy_5nn_in_2D_embedding chance
#> 0.936 0.100A 5-nearest-neighbour classifier operating in two dimensions recovers the digit with high accuracy, the embedding has retained nearly all the class-relevant information from 784 dimensions. The SOCR TensorBoard activity offers an interactive version of this exploration on UK Biobank data.
21 Case study: Parkinson’s disease
21.1 Collecting the data
The SOCR PD dataset combines clinical, genetic, and neuroimaging measures.
pd_raw <- read_html("https://wiki.socr.umich.edu/index.php/SOCR_Data_PD_BiomedBigMetadata") |>
html_nodes("table") |> _[[1]] |> html_table()
dim(pd_raw); names(pd_raw)#> [1] 1128 33#> [1] "Cases" "L_caudate_ComputeArea"
#> [3] "L_caudate_Volume" "R_caudate_ComputeArea"
#> [5] "R_caudate_Volume" "L_putamen_ComputeArea"
#> [7] "L_putamen_Volume" "R_putamen_ComputeArea"
#> [9] "R_putamen_Volume" "L_hippocampus_ComputeArea"
#> [11] "L_hippocampus_Volume" "R_hippocampus_ComputeArea"
#> [13] "R_hippocampus_Volume" "cerebellum_ComputeArea"
#> [15] "cerebellum_Volume" "L_lingual_gyrus_ComputeArea"
#> [17] "L_lingual_gyrus_Volume" "R_lingual_gyrus_ComputeArea"
#> [19] "R_lingual_gyrus_Volume" "L_fusiform_gyrus_ComputeArea"
#> [21] "L_fusiform_gyrus_Volume" "R_fusiform_gyrus_ComputeArea"
#> [23] "R_fusiform_gyrus_Volume" "Sex"
#> [25] "Weight" "Age"
#> [27] "Dx" "chr12_rs34637584_GT"
#> [29] "chr17_rs11868035_GT" "UPDRS_part_I"
#> [31] "UPDRS_part_II" "UPDRS_part_III"
#> [33] "Time"#>
#> HC PD SWEDD
#> 400 400 32821.2 Preparing the data
Three decisions, made explicitly rather than by index arithmetic.
How to treat SWEDD. Scans without evidence of
dopaminergic deficit denotes participants carrying a clinical PD
diagnosis whose DaTscan is normal. Recoding them as healthy controls is
a clinical judgement, not data cleaning, and it changes
what any subsequent separation analysis means. The default here keeps
three classes; the SWEDD_CODING flag in the setup chunk
switches to the binary coding, and both are examined below.
id_cols <- intersect(c("Cases", "Case", "ID", "Time"), names(pd_raw))
label_col <- "Dx"
pd_labels <- if (SWEDD_CODING == "binary") {
factor(ifelse(pd_raw[[label_col]] == "PD", "PD", "Non-PD"),
levels = c("Non-PD", "PD"))
} else {
factor(pd_raw[[label_col]], levels = c("HC", "PD", "SWEDD"))
}
table(pd_labels)#> pd_labels
#> HC PD SWEDD
#> 400 400 328# Numeric feature block: drop identifiers and the label, keep everything numeric
pd_features <- pd_raw |>
dplyr::select(-any_of(c(id_cols, label_col))) |>
dplyr::select(where(is.numeric)) |>
dplyr::select(where(~ sd(.x, na.rm = TRUE) > 0)) # drop constant columns
pd_features <- pd_features[stats::complete.cases(pd_features), , drop = FALSE]
pd_labels <- pd_labels[stats::complete.cases(
pd_raw |> dplyr::select(-any_of(c(id_cols, label_col))) |> dplyr::select(where(is.numeric)) |>
dplyr::select(where(~ sd(.x, na.rm = TRUE) > 0)))]
c(cases = nrow(pd_features), features = ncol(pd_features))#> cases features
#> 1128 30The label is never an input. It is used only to colour the results. Feeding a diagnosis into an unsupervised embedding and then observing that the diagnosis separates is circular.
Duplicates. The table holds repeated measures, so
identical feature rows occur. Rtsne() rejects them by
default; we deduplicate once and report how many rows were removed.
#> duplicate_rows_removed
#> 846pd_X <- as.matrix(pd_features[!dup, , drop = FALSE])
pd_y <- pd_labels[!dup]
pd_cols <- setNames(c("HC" = "#3B7DD8", "PD" = "#D8433B", "SWEDD" = "#E0B33B",
"Non-PD" = "#3B7DD8")[levels(pd_y)], levels(pd_y))
c(final_cases = nrow(pd_X), final_features = ncol(pd_X))#> final_cases final_features
#> 282 3021.3 PCA
Variables span imaging volumes in the thousands and clinical scores in the tens, so the correlation matrix is the right choice.
pca_pd <- prcomp(pd_X, center = TRUE, scale. = TRUE)
pd_cumvar <- cumsum(pca_pd$sdev^2) / sum(pca_pd$sdev^2)
data.frame(PC = 1:10,
variance_pct = round(100 * pca_pd$sdev[1:10]^2 / sum(pca_pd$sdev^2), 2),
cumulative_pct = round(100 * pd_cumvar[1:10], 2))#> PC variance_pct cumulative_pct
#> 1 1 5.72 5.72
#> 2 2 5.47 11.19
#> 3 3 5.17 16.35
#> 4 4 4.77 21.13
#> 5 5 4.64 25.77
#> 6 6 4.31 30.08
#> 7 7 4.27 34.35
#> 8 8 4.08 38.43
#> 9 9 3.86 42.29
#> 10 10 3.79 46.08c(components_for_60pct = which(pd_cumvar >= 0.60)[1],
components_for_80pct = which(pd_cumvar >= 0.80)[1],
total_features = ncol(pd_X))#> components_for_60pct components_for_80pct total_features
#> 14 21 30set.seed(151)
pd_null <- replicate(300, {
Xn <- apply(pd_X, 2, sample)
prcomp(Xn, center = TRUE, scale. = TRUE)$sdev^2
})
pd_pa95 <- apply(pd_null, 1, quantile, 0.95)
n_keep_pa <- sum(pca_pd$sdev^2 > pd_pa95)
data.frame(PC = seq_along(pca_pd$sdev),
Observed = pca_pd$sdev^2,
`Parallel analysis (95th pct)` = pd_pa95, check.names = FALSE) |>
pivot_longer(-PC, names_to = "series", values_to = "eigenvalue") |>
ggplot(aes(PC, eigenvalue, colour = series)) +
geom_line(linewidth = 0.9) + geom_point(size = 1.5) +
scale_colour_manual(values = c("Observed" = "steelblue",
"Parallel analysis (95th pct)" = "firebrick")) +
labs(title = "PD data: scree plot with a calibrated null",
subtitle = sprintf("Parallel analysis retains %d components; note the absence of an elbow",
n_keep_pa),
x = "Component", y = "Eigenvalue", colour = NULL) +
theme_dspa()
There is no elbow. The eigenvalues decay gradually, each component explaining a comparable slice, the signature of genuinely high-dimensional data with no dominant low-rank structure. Parallel analysis still gives a defensible cutoff where visual inspection would not.
pd_scores <- data.frame(PC1 = pca_pd$x[, 1], PC2 = pca_pd$x[, 2],
PC3 = pca_pd$x[, 3], Dx = pd_y)
ggplot(pd_scores, aes(PC1, PC2, colour = Dx)) +
geom_point(size = 1.6, alpha = 0.6) +
stat_ellipse(level = 0.68, linewidth = 0.7) +
scale_colour_manual(values = pd_cols) +
labs(title = "PD cohorts in the first two principal components",
subtitle = sprintf("PC1 %.1f%%, PC2 %.1f%% of variance; ellipses are 68%% contours",
100 * pca_pd$sdev[1]^2 / sum(pca_pd$sdev^2),
100 * pca_pd$sdev[2]^2 / sum(pca_pd$sdev^2)),
x = "PC1", y = "PC2", colour = "Diagnosis") +
theme_dspa()
ld_pd <- fix_signs(pca_pd$rotation)
top_vars <- order(rowSums(ld_pd[, 1:3]^2), decreasing = TRUE)[1:8]
sc_pd <- pca_pd$x
scale_pd <- 0.8 * max(abs(sc_pd[, 1:3])) / max(abs(ld_pd[top_vars, 1:3]))
p <- plot_ly()
for (lev in levels(pd_y)) {
sel <- pd_y == lev
p <- add_trace(p, x = sc_pd[sel, 1], y = sc_pd[sel, 2], z = sc_pd[sel, 3],
type = "scatter3d", mode = "markers", name = lev,
marker = list(size = 3, opacity = 0.7,
color = unname(pd_cols[lev])))
}
for (k in top_vars) {
p <- add_trace(p, x = c(0, ld_pd[k, 1]) * scale_pd,
y = c(0, ld_pd[k, 2]) * scale_pd,
z = c(0, ld_pd[k, 3]) * scale_pd,
type = "scatter3d", mode = "lines",
name = rownames(ld_pd)[k], line = list(width = 5),
showlegend = TRUE)
}
p |> layout(title = "PD data on the first three PCs, with the eight strongest loadings",
legend = list(orientation = "h", font = list(size = 9)),
scene = list(xaxis = list(title = "PC1"), yaxis = list(title = "PC2"),
zaxis = list(title = "PC3")))The cohorts overlap heavily, separating only slightly along PC2. That is the honest finding: a linear projection of these features does not distinguish the clinical groups.
21.4 Factor analysis
set.seed(157)
pd_par <- psych::fa.parallel(pd_X, fa = "fa", n.iter = 100, plot = FALSE,
show.legend = FALSE)#> Parallel analysis suggests that the number of factors = 0 and the number of components = NA#> parallel_analysis_suggests
#> 0# Maximum factors estimable from p variables (df must remain non-negative)
p_pd <- ncol(pd_X)
k_max <- floor((2 * p_pd + 1 - sqrt(8 * p_pd + 1)) / 2)
c(variables = p_pd, max_identifiable_factors = k_max)#> variables max_identifiable_factors
#> 30 22k_fa <- pd_par$nfact
# fa_pd <- factanal(pd_X, factors = k_fa, rotation = "varimax")
if (k_fa < 1) {
message("Factor analysis is not suitable: Parallel analysis suggests 0 latent factors.")
message(sprintf("Bartlett's test p-value = %.3f; KMO overall MSA = %.2f",
psych::cortest.bartlett(cor(pd_X), n = nrow(pd_X))$p.value,
psych::KMO(pd_X)$MSA))
message(sprintf("These results are scientifically sound! \n Factor Analysis looks for unobserved latent constructs that explain why variables move together. \n Look at the 22 imaging (L_caudate_ComputeArea, ...) and the 8 non-imaging variables, including 3 deographics (Sex, Weight, Age), 2 genetic SNPs (chr12..., chr17...), and 3 clinical scales (UPDRS_part_I, II, III). Biologically, there is no shared underlying factor that drives a patient's neuroanatomical, sex, genetic mutations, and their weight simultaneously. Because these variables are fundamentally independent, the parallel analysis correctly identifies that no latent factors exist!"))
} else {
fa_pd <- factanal(pd_X, factors = min(k_fa, k_max), rotation = "varimax")
c(factors = k_fa,
chisq = round(fa_pd$STATISTIC, 2),
df = fa_pd$dof,
p_value = signif(fa_pd$PVAL, 4),
cumulative_variance = round(sum(colSums(unclass(fa_pd$loadings)^2)) / p_pd, 4))
head(sort(1 - fa_pd$uniquenesses, decreasing = TRUE), 8) |> round(3)
}Check what is the share of the total variance explained by the factors chosen by parallel analysis, and the likelihood-ratio test rejects. Rather than incrementing \(k\) until the test stops rejecting, model selection by repeated testing, which invalidates the reported \(p\)-value, we report the misfit as a finding, a low-dimensional common-factor structure does not describe these data.
Note that the factor analysis failed as it identified \(0\) factors, which is expected and apprpriate in this case. Demographics (Age, Sex), Genetics (SNPs), Clinical scores (UPDRS) and computer derived neuroimaging morphopmetry measures do not share a single common biological mechanism. Hence, a latent factor model cannot be mathematically justified for this mixed dataset.
21.5 ICA
Rpd <- cor(pd_X)
off <- Rpd[upper.tri(Rpd)]
c(median_abs_correlation = round(median(abs(off)), 4),
pct_above_0.3 = round(100 * mean(abs(off) > 0.3), 2),
max_abs = round(max(abs(off)), 4))#> median_abs_correlation pct_above_0.3 max_abs
#> 0.0432 0.0000 0.1764Most pairwise correlations are near zero. ICA seeks a rotation making already-uncorrelated variables independent; when the variables are close to uncorrelated to begin with and no strong non-Gaussian mixing structure is present, there is little for it to undo. This is a data property, not a software limitation, and it is worth reporting rather than forcing a fit.
pd_first5 <- as.data.frame(pd_X[, 1:5])
pd_first5$Dx <- pd_y
GGally::ggpairs(
pd_first5,
columns = 1:5,
aes(colour = Dx, alpha = 0.5),
title = "Parkinson's disease data: first five features",
upper = list(continuous = GGally::wrap("cor", size = 2.6)),
lower = list(continuous = GGally::wrap("points", size = 0.4))
) +
scale_colour_manual(values = pd_cols) +
scale_fill_manual(values = pd_cols) +
theme_dspa(8)
# --- Interactive equivalent (brushable) ------------------------------------
dims <- lapply(1:5, \(j) list(label = colnames(pd_X)[j], values = pd_X[, j]))
plot_ly() |>
add_trace(type = "splom", dimensions = dims, text = as.character(pd_y),
marker = list(color = as.integer(pd_y), size = 4, opacity = 0.6,
line = list(width = 0.5, color = "rgb(230,230,230)"))) |>
layout(title = "Parkinson's disease data pairs plot",
hovermode = "closest", dragmode = "select",
plot_bgcolor = "rgba(240,240,240,0.95)")21.6 t-SNE and UMAP
perp_pd <- min(30, floor((nrow(pd_X) - 1) / 3)) # perplexity must satisfy 3*perp < n
set.seed(163)
tsne_pd2 <- Rtsne(scale(pd_X), dims = 2, perplexity = perp_pd, max_iter = 1000,
check_duplicates = FALSE, verbose = FALSE)
set.seed(163)
tsne_pd3 <- Rtsne(scale(pd_X), dims = 3, perplexity = perp_pd, max_iter = 1000,
check_duplicates = FALSE, verbose = FALSE)
c(perplexity_used = perp_pd, n = nrow(pd_X))#> perplexity_used n
#> 30 282ggplot(data.frame(x = tsne_pd2$Y[, 1], y = tsne_pd2$Y[, 2], Dx = pd_y),
aes(x, y, colour = Dx)) +
geom_point(size = 1.6, alpha = 0.7) +
scale_colour_manual(values = pd_cols) +
labs(title = "t-SNE embedding of the PD data",
subtitle = "Colour is the diagnosis, which was NOT supplied to the algorithm",
x = NULL, y = NULL, colour = "Diagnosis") +
theme_dspa()
p <- plot_ly()
for (lev in levels(pd_y)) {
sel <- pd_y == lev
p <- add_trace(p, x = tsne_pd3$Y[sel, 1], y = tsne_pd3$Y[sel, 2],
z = tsne_pd3$Y[sel, 3], type = "scatter3d", mode = "markers",
name = lev, marker = list(size = 3.2, opacity = 0.8,
color = unname(pd_cols[lev])))
}
p |> layout(title = "PD t-SNE 3-D embedding",
scene = list(xaxis = list(title = ""), yaxis = list(title = ""),
zaxis = list(title = "")))cfg_pd <- umap.defaults
cfg_pd$n_neighbors <- min(15, nrow(pd_X) - 1)
cfg_pd$min_dist <- 0.1
cfg_pd$n_components <- 2
cfg_pd$random_state <- 167
umap_pd <- umap(scale(pd_X), config = cfg_pd)
ggplot(data.frame(x = umap_pd$layout[, 1], y = umap_pd$layout[, 2], Dx = pd_y),
aes(x, y, colour = Dx)) +
geom_point(size = 1.6, alpha = 0.7) +
scale_colour_manual(values = pd_cols) +
labs(title = "UMAP embedding of the PD data",
subtitle = "Same conclusion as t-SNE: the cohorts do not separate in two dimensions",
x = NULL, y = NULL, colour = "Diagnosis") +
theme_dspa()
# --- Interactive equivalent ------------------------------------------------
plot_ly(x = umap_pd$layout[, 1], y = umap_pd$layout[, 2], mode = "text") |>
add_text(text = as.character(pd_y),
textfont = list(color = unname(pd_cols[as.character(pd_y)]))) |>
layout(title = "UMAP PD embedding",
xaxis = list(title = ""), yaxis = list(title = ""))21.7 Reading the result honestly
# Quantify separation: how well does a 5-NN classifier do IN each embedding?
sep_score <- function(emb, y, k = 5) {
Dm <- as.matrix(dist(emb))
diag(Dm) <- Inf
pred <- apply(Dm, 1, \(r) {
nb <- order(r)[1:k]
names(which.max(table(y[nb])))
})
mean(pred == as.character(y))
}
baseline <- max(table(pd_y)) / length(pd_y)
data.frame(
representation = c("Majority-class baseline", "PCA (2-D)", "t-SNE (2-D)",
"UMAP (2-D)", "Full feature space"),
knn5_accuracy = round(c(baseline,
sep_score(pca_pd$x[, 1:2], pd_y),
sep_score(tsne_pd2$Y, pd_y),
sep_score(umap_pd$layout, pd_y),
sep_score(scale(pd_X), pd_y)), 4))#> representation knn5_accuracy
#> 1 Majority-class baseline 0.3546
#> 2 PCA (2-D) 0.3759
#> 3 t-SNE (2-D) 0.3936
#> 4 UMAP (2-D) 0.3369
#> 5 Full feature space 0.4326None of the embeddings improves on the majority-class baseline by much, and the full feature space is no better. The PD cohorts are not separable by these features in any low-dimensional representation, linear or nonlinear.
That is a substantive finding, not a failure of technique. It says the data are intrinsically high-dimensional, consistent with the flat scree plot, the factor-analysis misfit, and the near-zero correlation structure. Three independent diagnostics agreeing is worth more than any single embedding that happened to look tidy.
The methodological lesson generalizes. An embedding that fails to separate known classes is informative. An embedding that separates them is only informative if the labels were withheld from the algorithm. Supervised methods that use the labels explicitly, the subject of Chapter 5 onward, may still succeed where unsupervised projection does not. Further PD analyses appear on the SOCR publications site.
22 Computational complexity summary
\(n\) = cases, \(p\) = features, \(k\) = target dimension, \(T\) = iterations.
| Method | Time | Memory | Scales to | Notes |
|---|---|---|---|---|
| Random projection | \(O(npk)\) | \(O(pk)\) | \(10^7\)+ | No decomposition; distance-preserving by Johnson–Lindenstrauss |
| PCA via SVD of \(X\) | \(O(np\min(n,p))\) | \(O(np)\) | \(10^5\) | prcomp(); numerically preferred |
| PCA via eigen of \(X^\top X\) | \(O(np^2+p^3)\) | \(O(p^2)\) | \(10^5\) | princomp(); squares the condition number |
| Randomized SVD (top \(k\)) | \(O(npk)\) | \(O(nk)\) | \(10^7\) | rsvd, irlba; the right tool when \(k\ll\min(n,p)\) |
| Classical MDS | \(O(n^3)\) | \(\mathbf{O(n^2)}\) | \(10^4\) | Equivalent to PCA on Euclidean distances |
| Non-metric MDS | \(O(Tn^2)\) | \(O(n^2)\) | \(10^4\) | Preserves rank order only |
| ICA (FastICA) | \(O(np^2+Tnpk)\) | \(O(np)\) | \(10^5\) | Whitening dominates when \(p\) is large |
| Factor analysis (ML) | \(O(Tp^3)\) | \(O(p^2)\) | \(p\lesssim 10^3\) | Iterative; watch for Heywood cases |
| Kernel PCA | \(O(n^3)\) | \(\mathbf{O(n^2)}\) | \(10^4\) | Nyström or random features beyond that |
| t-SNE (exact) | \(O(Tn^2)\) | \(O(n^2)\) | \(5\times10^3\) | Every pair, every iteration |
| t-SNE (Barnes–Hut) | \(O(Tn\log n)\) | \(O(n)\) | \(10^5\) | Rtsne() default, \(\theta=0.5\) |
| t-SNE (FIt-SNE) | \(O(Tn)\) | \(O(n)\) | \(10^6\) | Interpolation + FFT |
| UMAP | \(\approx O(n^{1.14})\) + \(O(Tnk)\) | \(O(nk)\) | \(10^6\) | NN-descent for approximate \(k\)-NN; negative sampling |
| Trustworthiness / continuity | \(O(n^2\log n)\) | \(O(n^2)\) | \(5\times10^3\) | Full rank matrices in both spaces |
Three practical rules follow.
Reduce before you embed. Running t-SNE or UMAP on the leading 50 PCs rather than raw pixels cuts the \(k\)-NN cost by an order of magnitude and usually improves the embedding, because PCA strips high-frequency noise.
Memory, not time, is the usual wall. Kernel PCA,
classical MDS, and exact t-SNE all need an \(n\times n\) matrix. At \(n=10^5\) that is 80 GB in double precision,
the same constraint that stops hclust() (Chapter 3,
§3.21).
Match the algorithm to \(k\). If you want 10 components out of 10,000, computing all 10,000 and discarding 9,990 is waste. Randomized and truncated SVD exist for exactly this.
23 Common pitfalls
| # | Pitfall | Consequence | Fix |
|---|---|---|---|
| 1 | PCA on unscaled data with incompatible units | The largest-variance variable owns PC1 by construction | scale. = TRUE unless units are comparable |
| 2 | Interpreting the sign of a loading | Sign is arbitrary and platform-dependent | Only relative signs within a component mean anything |
| 3 | Subtracting a scalar and calling it centering | Column means are not zeroed | sweep(X, 2, colMeans(X)) or
scale(X, scale = FALSE) |
| 4 | Treating PC1 as the regression line | Attenuation when \(x\) has measurement error | PC1 is total least squares — perpendicular, not vertical |
| 5 | Choosing \(k\) by eyeballing an elbow | Real scree plots are often smooth | Parallel analysis, broken stick, or cross-validated error |
| 6 | Bootstrapping one model, plotting another’s estimate | The interval describes a different quantity | Same PCA specification throughout |
| 7 | palette(rainbow(k)) |
Returns the previous palette and mutates global state | rainbow(k) directly, or
scale_colour_manual() |
| 8 | Forgetting ICA’s independence assumption | Misreading what the method delivers | Independence is assumed; Gaussianity is what is dropped |
| 9 | Interpreting ICA component order or scale | Both are unidentifiable | Fix a convention yourself; never compare “IC1” across runs |
| 10 | Running ICA with two or more Gaussian sources | Mixing matrix is unrecoverable in principle | Comon’s theorem: at most one Gaussian |
| 11 | Calling FA “generalized PCA” | Different models, different questions | FA has an error term, a fit test, and a rotational indeterminacy |
| 12 | Incrementing factors until \(p>0.05\) | Invalidates the reported \(p\)-value; biases \(k\) upward | Parallel analysis or information criteria |
| 13 | Ignoring Heywood cases | Negative variance estimates signal misspecification | Reduce \(k\), or re-examine the variable |
| 14 | Measuring distances between t-SNE clusters | The KL objective barely constrains global layout | Report cluster membership, not geometry |
| 15 | Reading cluster size off t-SNE or UMAP | Density is normalized away | Count the points |
| 16 | Using one t-SNE seed | Stochastic; different seeds give different maps | Set a seed and verify stability across several |
| 17 | Believing apparent clusters in an embedding | t-SNE clusters pure noise | Compare against a permuted or noise null |
| 18 | Expecting t-SNE to project new data | It learns coordinates, not a function | UMAP’s predict(), or refit |
| 19 | Feeding the class label into an unsupervised embedding | Circular: separation is an artifact | Labels colour the result; they are never inputs |
| 20 | Embedding without deduplicating | Rtsne() errors on duplicate rows |
Deduplicate and report how many were dropped |
| 21 | Using princomp() by habit |
Eigen of covariance; divides by \(n\); less stable | prcomp() |
| 22 | Forming an \(n\times n\) kernel or distance matrix at scale | Memory blowup long before time becomes an issue | Nyström, landmarks, or approximate \(k\)-NN |
24 Practice problems
24.1 Problem 1: Prove the two PCA objectives coincide
Show that minimizing reconstruction error \(\|X-XV_kV_k^\top\|_F^2\) is equivalent to maximizing retained variance \(\operatorname{tr}(V_k^\top SV_k)\), and verify numerically.
Solution
Since \(V_k^\top V_k=I_k\), the projector \(P=V_kV_k^\top\) is symmetric idempotent (Chapter 3, §3.4.2). Then \[\|X-XP\|_F^2=\operatorname{tr}\big((X-XP)^\top(X-XP)\big) =\operatorname{tr}(X^\top X)-2\operatorname{tr}(PX^\top X)+\operatorname{tr}(PX^\top XP).\] Idempotence and cyclicity give \(\operatorname{tr}(PX^\top XP)=\operatorname{tr}(P^2X^\top X)=\operatorname{tr}(PX^\top X)\), so \[\|X-XP\|_F^2=\operatorname{tr}(X^\top X)-\operatorname{tr}\big(V_k^\top X^\top XV_k\big).\] The first term does not depend on \(V_k\). \(\blacksquare\)
set.seed(211)
Xp <- scale(matrix(rnorm(200 * 8), 200, 8) %*% matrix(rnorm(64), 8, 8))
pc <- prcomp(Xp, center = TRUE, scale. = FALSE)
Sp <- cov(Xp)
check <- function(k) {
Vk <- pc$rotation[, 1:k, drop = FALSE]
recon_err <- norm(Xp - Xp %*% Vk %*% t(Vk), "F")^2
variance_kept <- sum(diag(t(Vk) %*% Sp %*% Vk))
c(k = k, recon_error = recon_err,
total_minus_kept = sum(diag(crossprod(Xp))) - variance_kept * (nrow(Xp) - 1),
variance_kept = variance_kept,
sum_top_k_eigen = sum(pc$sdev[1:k]^2))
}
do.call(rbind, lapply(1:4, check)) |> round(6)#> k recon_error total_minus_kept variance_kept sum_top_k_eigen
#> [1,] 1 1059.640 1059.640 2.67518 2.67518
#> [2,] 2 694.608 694.608 4.50951 4.50951
#> [3,] 3 430.644 430.644 5.83596 5.83596
#> [4,] 4 238.836 238.836 6.79982 6.7998224.2 Problem 2: Random directions beat PC1 for reconstruction?
Compare the reconstruction error of the top-\(k\) PCA subspace against 1,000 random \(k\)-dimensional subspaces. How often does a random subspace win?
Solution
set.seed(223)
k <- 3
pca_err <- norm(Xp - Xp %*% pc$rotation[, 1:k] %*% t(pc$rotation[, 1:k]), "F")^2
rand_err <- replicate(1000, {
V <- qr.Q(qr(matrix(rnorm(ncol(Xp) * k), ncol(Xp), k)))
norm(Xp - Xp %*% V %*% t(V), "F")^2
})
c(pca_error = round(pca_err, 3),
best_random = round(min(rand_err), 3),
median_random = round(median(rand_err), 3),
times_random_beat_pca = sum(rand_err < pca_err))#> pca_error best_random median_random
#> 430.644 689.372 1004.043
#> times_random_beat_pca
#> 0.000ggplot(data.frame(e = rand_err), aes(e)) +
geom_histogram(bins = 40, fill = "grey80", colour = "white") +
geom_vline(xintercept = pca_err, colour = "firebrick", linewidth = 1.2) +
labs(title = "Reconstruction error: 1,000 random subspaces vs. PCA",
subtitle = "The red line is PCA. Eckart-Young guarantees nothing lies to its left",
x = "Squared Frobenius reconstruction error", y = "Count") +
theme_dspa()
24.3 Problem 3: Verify Johnson–Lindenstrauss numerically
For fixed \(\epsilon=0.2\), find empirically the smallest \(k\) at which a random projection keeps all pairwise distortions under \(\epsilon\), and compare with the theoretical bound \(8\ln n/\epsilon^2\).
Solution
set.seed(227)
n3 <- 300; d3 <- 1000; eps <- 0.2
X3 <- matrix(rnorm(n3 * d3), n3, d3)
D3 <- as.numeric(dist(X3))
max_distortion <- function(k, reps = 5) {
mean(replicate(reps, {
R <- matrix(rnorm(d3 * k), d3, k) / sqrt(k)
max(abs((as.numeric(dist(X3 %*% R)) / D3)^2 - 1))
}))
}
ks <- c(25, 50, 100, 200, 400, 800)
res3 <- data.frame(k = ks, mean_max_distortion = sapply(ks, max_distortion))
res3$within_eps <- res3$mean_max_distortion < eps
res3 |> mutate(mean_max_distortion = round(mean_max_distortion, 4))#> k mean_max_distortion within_eps
#> 1 25 1.7874 FALSE
#> 2 50 1.0203 FALSE
#> 3 100 0.7115 FALSE
#> 4 200 0.4666 FALSE
#> 5 400 0.3277 FALSE
#> 6 800 0.2361 FALSEc(theoretical_bound = ceiling(8 * log(n3) / eps^2),
empirical_sufficient_k = min(res3$k[res3$within_eps]))#> theoretical_bound empirical_sufficient_k
#> 1141 Inf24.4 Problem 4: When do OLS and PC1 agree?
Derive the condition under which the OLS slope equals the PC1 slope, and verify by simulation.
Solution
OLS gives \(m_{\text{OLS}}=s_{xy}/s_x^2\). Setting the PC1 expression equal: \[\frac{s_y^2-s_x^2+\sqrt{(s_y^2-s_x^2)^2+4s_{xy}^2}}{2s_{xy}}=\frac{s_{xy}}{s_x^2}.\] Cross-multiplying and simplifying yields \(s_x^2s_y^2=s_{xy}^2\), i.e. \[r^2=\frac{s_{xy}^2}{s_x^2s_y^2}=1 .\] The two lines coincide only when the correlation is \(\pm1\), perfect collinearity, where all four lines of §4.5 collapse.
slopes <- function(rho, n = 4000, seed = 229) {
set.seed(seed)
L <- t(chol(matrix(c(1, rho, rho, 1), 2)))
Z <- t(L %*% matrix(rnorm(2 * n), 2))
sx2 <- var(Z[, 1]); sy2 <- var(Z[, 2]); sxy <- cov(Z[, 1], Z[, 2])
c(rho = rho, ols = sxy / sx2,
pc1 = (sy2 - sx2 + sqrt((sy2 - sx2)^2 + 4 * sxy^2)) / (2 * sxy))
}
sl <- as.data.frame(do.call(rbind, lapply(c(0.2, 0.5, 0.8, 0.95, 0.99, 0.999), slopes)))
sl$gap <- abs(sl$pc1 - sl$ols)
round(sl, 4)#> rho ols pc1 gap
#> 1 0.200 0.2227 1.0346 0.8119
#> 2 0.500 0.5201 1.0242 0.5041
#> 3 0.800 0.8139 1.0152 0.2013
#> 4 0.950 0.9572 1.0075 0.0503
#> 5 0.990 0.9933 1.0033 0.0101
#> 6 0.999 1.0000 1.0010 0.0010ggplot(sl, aes(rho, gap)) +
geom_line(linewidth = 1, colour = "firebrick") + geom_point(size = 2.4) +
labs(title = "The OLS and PC1 slopes converge only as |r| approaches 1",
x = "Correlation", y = "|PC1 slope - OLS slope|") +
theme_dspa()
24.5 Problem 5: Two Gaussian sources defeat ICA
Show empirically that ICA recovers non-Gaussian sources but fails on two Gaussian ones.
Solution
recover <- function(S, seed = 233) {
set.seed(seed)
A <- matrix(c(2, 1, 1, 3), 2, 2)
Xm <- S %*% A
fit <- fastICA::fastICA(Xm, 2, alg.typ = "parallel", fun = "logcosh",
maxit = 400, tol = 1e-7)
M <- abs(cor(fit$S, S)) # 2x2 recovery matrix
# Best matching over the two possible permutations
max(mean(c(M[1, 1], M[2, 2])), mean(c(M[1, 2], M[2, 1])))
}
set.seed(239); n5 <- 5000
c(uniform_sources = round(recover(cbind(runif(n5), runif(n5))), 4),
laplace_sources = round(recover(cbind(rexp(n5) * sample(c(-1, 1), n5, TRUE),
rexp(n5) * sample(c(-1, 1), n5, TRUE))), 4),
one_gaussian = round(recover(cbind(rnorm(n5), runif(n5))), 4),
two_gaussian_FAILS = round(recover(cbind(rnorm(n5), rnorm(n5))), 4))#> uniform_sources laplace_sources one_gaussian two_gaussian_FAILS
#> 0.9999 0.9999 0.9986 0.902724.6 Problem 6: Perplexity, distortion, and neighbourhood scale
Embed a dataset with known cluster structure at several perplexities and measure trustworthiness at several \(k\). Explain the pattern.
Solution
set.seed(241)
n_per <- 150
clusters <- do.call(rbind, lapply(1:4, \(g)
matrix(rnorm(n_per * 10, mean = g * 4), n_per, 10)))
truth <- rep(1:4, each = n_per)
grid6 <- expand.grid(perp = c(5, 20, 50), k = c(5, 20, 50))
grid6$trust <- mapply(function(pp, kk) {
set.seed(251)
emb <- Rtsne(clusters, dims = 2, perplexity = pp, max_iter = 600,
check_duplicates = FALSE)$Y
trust_cont(clusters, emb, k = kk)[["trustworthiness"]]
}, grid6$perp, grid6$k)
grid6 |> mutate(trust = round(trust, 4))#> perp k trust
#> 1 5 5 0.9741
#> 2 20 5 0.9785
#> 3 50 5 0.9669
#> 4 5 20 0.9463
#> 5 20 20 0.9605
#> 6 50 20 0.9626
#> 7 5 50 0.9495
#> 8 20 50 0.9655
#> 9 50 50 0.9708ggplot(grid6, aes(factor(perp), trust, fill = factor(k))) +
geom_col(position = position_dodge(width = 0.8), width = 0.7) +
coord_cartesian(ylim = c(0.9, 1)) +
scale_fill_brewer(palette = "Set2", name = "Evaluation k") +
labs(title = "Trustworthiness depends on both perplexity and the scale you evaluate at",
x = "t-SNE perplexity", y = "Trustworthiness") +
theme_dspa()
24.7 Problem 7: Cross-validated rank selection on real data
Apply the entry-holdout cross-validation of §4.4.6 to the MNIST subset and compare with the 90%-variance rule.
Solution
set.seed(257)
Xm_small <- Xm[sample(nrow(Xm), 600), keep_px]
cv_mnist <- function(X, ks, folds = 5, seed = 263) {
set.seed(seed)
Xs <- scale(X, center = TRUE, scale = FALSE)
fold_id <- sample(rep(seq_len(folds), length.out = length(Xs)))
sapply(ks, function(k) {
mean(sapply(seq_len(folds), function(f) {
hold <- which(fold_id == f)
Xtr <- Xs; Xtr[hold] <- 0 # centred, so 0 is the mean
sv <- svd(Xtr, nu = k, nv = k)
Xk <- sv$u %*% diag(sv$d[1:k], k, k) %*% t(sv$v)
mean((Xs[hold] - Xk[hold])^2)
}))
})
}
ks7 <- c(2, 5, 10, 20, 40, 60, 90, 130, 200)
cvm <- data.frame(k = ks7, cv_mse = cv_mnist(Xm_small, ks7))
cvm |> mutate(cv_mse = round(cv_mse, 6))#> k cv_mse
#> 1 2 0.065125
#> 2 5 0.053792
#> 3 10 0.044290
#> 4 20 0.035354
#> 5 40 0.029918
#> 6 60 0.030580
#> 7 90 0.036826
#> 8 130 0.047525
#> 9 200 0.062661#> cv_optimal_k variance_rule_90pct
#> 40 86ggplot(cvm, aes(k, cv_mse)) +
geom_line(linewidth = 1, colour = "steelblue") + geom_point(size = 2.2) +
geom_point(data = cvm[which.min(cvm$cv_mse), ], colour = "firebrick", size = 3.6) +
scale_x_log10() +
labs(title = "Cross-validated rank selection on MNIST",
subtitle = "Held-out entries only; the minimum is the predictively optimal rank",
x = "Rank k (log scale)", y = "Held-out MSE") +
theme_dspa()
24.8 Problem 8: Build classical MDS from scratch
Implement double-centering, recover coordinates from a distance matrix alone, and confirm the result matches PCA up to sign and rotation.
Solution
my_cmdscale <- function(D, k = 2) {
D2 <- as.matrix(D)^2
n <- nrow(D2)
J <- diag(n) - matrix(1 / n, n, n)
B <- -0.5 * J %*% D2 %*% J
e <- eigen(B, symmetric = TRUE)
pos <- which(e$values > 1e-9)[1:k]
list(points = e$vectors[, pos, drop = FALSE] %*% diag(sqrt(e$values[pos]), k, k),
eig = e$values)
}
set.seed(269)
Xd <- scale(matrix(rnorm(120 * 6), 120, 6) %*% matrix(rnorm(36), 6, 6))
Dd <- dist(Xd)
mine <- my_cmdscale(Dd, k = 3)
base <- cmdscale(Dd, k = 3)
pc_d <- prcomp(Xd, center = TRUE, scale. = FALSE)$x[, 1:3]
c(matches_cmdscale = max(abs(abs(mine$points) - abs(base))),
matches_pca = max(abs(abs(mine$points) - abs(pc_d))),
distances_preserved = max(abs(as.numeric(dist(mine$points)) - as.numeric(Dd))))#> matches_cmdscale matches_pca distances_preserved
#> 1.60982e-14 8.88178e-15 3.22938e+0025 Checkpoint
- Two variables in your dataset are measured in cubic millimetres and
on a 0–4 clinical scale. You run
prcomp(X, center = TRUE). What have you actually computed, and what should you have done? - A colleague reports that reversing the sign of PC2 changed their biological interpretation. What do you tell them?
- Your t-SNE plot shows cluster A far from cluster B and adjacent to cluster C. What may you conclude about the relationships among A, B, and C?
- You have 4 microphones recording 4 speakers. Two speakers produce near-Gaussian output. What does ICA give you?
- A scree plot decays smoothly with no elbow across 30 components. Name two defensible ways to pick \(k\), and one thing the flat scree itself tells you.
- Why can UMAP project a new observation into an existing embedding while t-SNE cannot?
Answers
- You computed PCA on the covariance matrix, so the
variable with the larger raw variance, almost certainly the volume in
cubic millimetres — dominates PC1 for reasons of units, not biology. Use
scale. = TRUE(correlation PCA) whenever the variables have incomparable units. - The sign of every principal component is arbitrary: \(-\mathbf{a}\) is as valid an eigenvector as \(\mathbf{a}\), with the same eigenvalue and the same subspace, and which one is returned depends on the LAPACK build. Only the relative signs of loadings within a single component are interpretable. Any conclusion that flips when the sign flips is not a conclusion.
- Essentially nothing about the relationships. The KL objective heavily penalizes misplacing near neighbours and barely penalizes misplacing distant ones, so global layout is not constrained. Cluster sizes are also uninformative because the adaptive bandwidth normalizes density. You may conclude that the points within A are mutually similar, nothing more.
- Nothing usable for the two Gaussian speakers. Comon’s theorem requires at most one Gaussian source: two independent Gaussians have a spherically symmetric joint density, so any rotation of them is equally independent and equally Gaussian, and the mixing matrix cannot be recovered. The two non-Gaussian speakers may still be separated; the Gaussian pair will emerge as an arbitrary rotation of each other.
- Defensible choices: parallel analysis (retain components exceeding the 95th percentile of a permutation null) and cross-validated reconstruction error (hold out entries and minimize held-out MSE); broken-stick is a third. The flat scree itself is a finding: it says the data have no dominant low-rank structure and are intrinsically high-dimensional, so no two-dimensional picture will summarize them honestly.
- t-SNE optimizes a set of coordinates \(\{y_i\}\) directly, there is no function to
apply to a new point, so adding an observation requires re-running the
optimization and gives a different map. UMAP fits a parameterized
fuzzy-graph model whose embedding can be extended: a new point is placed
by finding its neighbours in the training graph and running a short
optimization against the frozen existing layout. That
is what
predict()does, and it is what makes UMAP usable inside a supervised pipeline where t-SNE is not.
26 Summary
Why reduce
- Distances concentrate in high dimensions: the ratio of their spread to their mean tends to zero, and more data does not help. Only reducing \(d\) does.
- Real data have low intrinsic dimension, that is what makes reduction possible rather than merely lossy.
- Johnson–Lindenstrauss bounds the achievable compression at \(k\ge 8\ln n/\epsilon^2\), logarithmic in \(n\) and independent of \(d\), achieved by a random matrix requiring no data.
Linear methods
- PCA maximizes variance and minimizes reconstruction error, the same problem, solved by the top eigenvectors of the covariance, with truncation error given exactly by Eckart–Young.
- PCA is scale dependent. Standardize unless the units are comparable.
- PC1 is not the regression line. OLS minimizes vertical residuals and attenuates under measurement error in \(x\); PC1 minimizes perpendicular ones and does not.
- SVD is PCA on centred data, with \(\lambda_j=d_j^2/(n-1)\). Randomized SVD gets the top \(k\) for \(O(npk)\).
- Classical MDS on Euclidean distances is PCA computed from the other margin.
- ICA assumes independence and drops Gaussianity; it is identifiable up to permutation and scale, and only if at most one source is Gaussian.
- FA is a generative model, not generalized PCA: it has an error term \(\Psi\), a testable fit, and a rotational indeterminacy that PCA lacks.
Nonlinear methods
- No linear projection can unroll a curved manifold, because the meaningful distance is geodesic, not Euclidean.
- Kernel PCA performs linear PCA in an implicit feature space, at \(O(n^3)\).
- t-SNE matches neighbourhood probability distributions by minimizing KL, with a heavy-tailed map kernel to relieve crowding. It preserves local structure, not global geometry, and has no out-of-sample extension.
- UMAP optimizes cross-entropy between fuzzy simplicial sets. The
explicit repulsive term and spectral initialization preserve more global
structure, it is far faster, and it supports
predict(). - Quantify embeddings with trustworthiness and continuity rather than judging them by eye.
Where these threads continue
| Thread | Continues in |
|---|---|
| Reduced features as classifier inputs | Supervised classification |
| Kernels: from kernel PCA to SVM | Black-box methods |
| Choosing \(k\) by cross-validation | Model assessment |
| Mixture models and cluster structure | Unsupervised clustering |
| PCA versus sparse and regularized selection | Feature selection |
| Autoencoders as learned nonlinear reduction | Deep learning |
| Gradient descent behind t-SNE and UMAP | Function optimization |
27 Chapter roadmap
- Chapter 1:
Foundations. R toolchain, reproducibility conventions,
dspa_read(), simulation. - Chapter 2: Data quality and exploratory visual analytics. Centrality, dispersion, density estimation, missingness.
- Chapter 3: Linear algebra, matrix computing, and regression. Rank, projection, SVD, conditioning, least squares.
- Supervised classification. kNN, naive Bayes, decision trees, logistic regression.
- Black-box methods. Neural networks, SVM and the kernel trick, random forests, gradient boosting.
- Model assessment, validation, improvement. Cross-validation, ROC/AUC, calibration, hyperparameter tuning.
- Unsupervised clustering. k-means, hierarchical, spectral, Gaussian mixture models.
- Variable importance and feature selection. Ridge, LASSO, elastic net, stability selection, FDR control.
- Function optimization. Gradient descent, Newton and quasi-Newton methods, EM, Bayesian optimization.
- Deep learning. Autoencoders, convolutional and recurrent architectures, representation learning.
28 Session information
#> R version 4.3.3 (2024-02-29 ucrt)
#> Platform: x86_64-w64-mingw32/x64 (64-bit)
#> Running under: Windows 11 x64 (build 26200)
#>
#> Matrix products: default
#>
#>
#> locale:
#> [1] LC_COLLATE=English_United States.utf8
#> [2] LC_CTYPE=English_United States.utf8
#> [3] LC_MONETARY=English_United States.utf8
#> [4] LC_NUMERIC=C
#> [5] LC_TIME=English_United States.utf8
#>
#> time zone: America/New_York
#> tzcode source: internal
#>
#> attached base packages:
#> [1] stats graphics grDevices utils datasets methods base
#>
#> other attached packages:
#> [1] class_7.3-22 umap_0.2.10.0 Rtsne_0.17 kernlab_0.9-32
#> [5] igraph_2.0.3 psych_2.4.6.26 fastICA_1.2-4 rvest_1.0.4
#> [9] MASS_7.3-60.0.1 plotly_4.12.0 patchwork_1.3.0 tidyr_1.3.1
#> [13] dplyr_1.1.4 ggplot2_4.0.1
#>
#> loaded via a namespace (and not attached):
#> [1] gtable_0.3.6 xfun_0.52 bslib_0.9.0 htmlwidgets_1.6.4
#> [5] websocket_1.4.1 processx_3.8.6 ggrepel_0.9.5 GGally_2.2.1
#> [9] lattice_0.22-6 vctrs_0.6.5 tools_4.3.3 crosstalk_1.2.1
#> [13] ps_1.9.0 generics_0.1.3 curl_6.2.0 parallel_4.3.3
#> [17] tibble_3.2.1 pkgconfig_2.0.3 Matrix_1.6-5 data.table_1.16.4
#> [21] RColorBrewer_1.1-3 S7_0.2.1 lifecycle_1.0.5 compiler_4.3.3
#> [25] farver_2.1.2 stringr_1.5.1 mnormt_2.1.1 chromote_0.4.0
#> [29] codetools_0.2-20 htmltools_0.5.8.1 sass_0.4.9 yaml_2.3.10
#> [33] lazyeval_0.2.2 later_1.4.1 pillar_1.10.1 jquerylib_0.1.4
#> [37] openssl_2.3.1 cachem_1.1.0 nlme_3.1-165 RSpectra_0.16-1
#> [41] ggstats_0.6.0 tidyselect_1.2.1 digest_0.6.37 stringi_1.8.4
#> [45] purrr_1.0.2 labeling_0.4.3 fastmap_1.2.0 grid_4.3.3
#> [49] cli_3.6.3 magrittr_2.0.3 withr_3.0.2 scales_1.4.0
#> [53] promises_1.3.2 rmarkdown_2.31 httr_1.4.7 otel_0.2.0
#> [57] reticulate_1.38.0 png_0.1-8 askpass_1.2.1 evaluate_1.0.3
#> [61] knitr_1.51 viridisLite_0.4.2 rlang_1.1.5 Rcpp_1.0.14
#> [65] glue_1.8.0 selectr_0.4-2 xml2_1.3.6 rstudioapi_0.18.0
#> [69] jsonlite_1.8.9 plyr_1.8.9 R6_2.6.1