Principal Component Analysis
When data has many dimensions, the first question is: which directions actually carry the variation? Principal Component Analysis (PCA) answers it. The first principal component is the single direction along which the data spreads out the most — the long axis of the cloud.
We'll do it in 2D where you can see it. Generate a tilted, noisy cloud of points; PCA finds the arrow pointing along its length. The recipe: center the data, build its covariance matrix, and find that matrix's top eigenvector — here by a few rounds of power iteration, all by hand and all live in WebAssembly.
begin
using PlutoUI, WasmMakie
endcloud tilt =
scatter (off-axis noise) =
number of points =
seed =
let
# generate the cloud and accumulate its covariance in one pass
tilt = Float64(tilti) / 10.0
noise = Float64(noisei) / 10.0
n = npoints
s = (seed * 2654435761 + 12345) % 2147483647
if s == 0
s = 1
end
dx = Vector{Float64}(undef, n)
dy = Vector{Float64}(undef, n)
mx = 0.0
my = 0.0
for i in 1:n
s = (s * 16807) % 2147483647
u = Float64(s) / 2147483647.0
s = (s * 16807) % 2147483647
v = Float64(s) / 2147483647.0
x = 4.0 * (u - 0.5)
y = tilt * x + noise * (v - 0.5) * 2.0
dx[i] = x
dy[i] = y
mx = mx + x
my = my + y
end
nn = Float64(n)
mx = mx / nn
my = my / nn
cxx = 0.0
cyy = 0.0
cxy = 0.0
for i in 1:n
ax_ = dx[i] - mx
ay_ = dy[i] - my
cxx = cxx + ax_ * ax_
cyy = cyy + ay_ * ay_
cxy = cxy + ax_ * ay_
end
cxx = cxx / nn
cyy = cyy / nn
cxy = cxy / nn
# power iteration: repeatedly apply the covariance matrix to a vector and it lines
# up with the top eigenvector (the first principal component)
vx = 1.0
vy = 0.0
for it in 1:24
nx = cxx * vx + cxy * vy
ny = cxy * vx + cyy * vy
len = sqrt(nx * nx + ny * ny)
if len < 0.000001
len = 1.0
end
vx = nx / len
vy = ny / len
end
lam = vx * (cxx * vx + cxy * vy) + vy * (cxy * vx + cyy * vy)
half = 2.0 * sqrt(lam)
fig = Figure(size = (520, 520))
ax = Axis(fig[1, 1])
for i in 1:n
lines!(ax, [dx[i], dx[i]], [dy[i] - 0.06, dy[i] + 0.06]) # a data point
end
# the first principal component, drawn through the centroid
lines!(ax, [mx - half * vx, mx + half * vx], [my - half * vy, my + half * vy])
fig
endpca_stats = let
# compute the principal component in a SEPARATE bond-dependent cell so the markdown
# below interpolates it live (values inside a markdown cell's own `let` bake to the
# slider defaults).
tilt = Float64(tilti) / 10.0
noise = Float64(noisei) / 10.0
n = npoints
s = (seed * 2654435761 + 12345) % 2147483647
if s == 0
s = 1
end
mx = 0.0
my = 0.0
dx = Vector{Float64}(undef, n)
dy = Vector{Float64}(undef, n)
for i in 1:n
s = (s * 16807) % 2147483647
u = Float64(s) / 2147483647.0
s = (s * 16807) % 2147483647
v = Float64(s) / 2147483647.0
x = 4.0 * (u - 0.5)
y = tilt * x + noise * (v - 0.5) * 2.0
dx[i] = x
dy[i] = y
mx = mx + x
my = my + y
end
nn = Float64(n)
mx = mx / nn
my = my / nn
cxx = 0.0
cyy = 0.0
cxy = 0.0
for i in 1:n
cxx = cxx + (dx[i] - mx) * (dx[i] - mx)
cyy = cyy + (dy[i] - my) * (dy[i] - my)
cxy = cxy + (dx[i] - mx) * (dy[i] - my)
end
cxx = cxx / nn
cyy = cyy / nn
cxy = cxy / nn
vx = 1.0
vy = 0.0
for it in 1:24
nx = cxx * vx + cxy * vy
ny = cxy * vx + cyy * vy
len = sqrt(nx * nx + ny * ny)
if len < 0.000001
len = 1.0
end
vx = nx / len
vy = ny / len
end
lam1 = vx * (cxx * vx + cxy * vy) + vy * (cxy * vx + cyy * vy)
total = cxx + cyy
frac = 0.0
if total > 0.000001
frac = lam1 / total
end
(floor(vx * 100.0) / 100.0, floor(vy * 100.0) / 100.0, floor(frac * 1000.0) / 10.0)
end;First principal component direction: about (0.74, 0.66). It captures 97.5% of the cloud's total variance. Tighten the scatter and that number climbs toward 100% (the data becomes nearly one-dimensional); loosen it and the cloud rounds out so no single direction dominates.
Why this is everywhere
PCA is the workhorse of dimensionality reduction. Faces, gene expression, word counts — data living in thousands of dimensions often really varies along just a handful. Keep the top few principal components and you can compress, denoise, and visualize it, throwing away the directions that were mostly noise.
The same eigenvector idea — find the directions a matrix stretches the most — turns up again and again: in the SVD that compresses images (Module 1), in PageRank, in the normal modes of a vibrating structure. Variance points the way.
Appendix
The MIT lecture computes PCA with LinearAlgebra eigensolvers over Images/Plots. WebAssembly can't run that stack in the browser, so we build the 2x2 covariance by hand and extract its top eigenvector with a few rounds of power iteration, drawing with WasmMakie. For two dimensions this converges in a blink and matches the exact answer.