Principal component analysis and the SVD
15 July 2024, daily mean temperature (ERA5)
Twenty-two cities from Athens to Aleppo give twenty-two temperatures every day. Write down one number per day instead of twenty-two. How close can you get to all of them, and which number should it be?
By the end you can
Why fewer numbers at all: with many dimensions, data is sparse (most of the space is empty), distances between points become alike, and a model needs far more examples to cover the space. This is the curse of dimensionality.
Stack the days as rows and the cities as columns: a matrix of 5,844 rows and 22 columns. ERA5 is a reanalysis: a weather model kept close to observations, read at the grid cell nearest each city (Hersbach and colleagues, 2020). The data (CSV)
Beirut and Damascus are 85 km apart and both are warm in August. Gaziantep and Alexandria are not alike, but they too are warm in August. The 22 numbers of a day are far from independent.
Representing data with fewer numbers than it has columns, while keeping most of what is in it, is dimensionality reduction. To make the question precise we need two cities first.
Part 2 of 5
Which single line through the data keeps the most?
| x₁ | x₂ |
|---|---|
| 14.8 | 8.5 |
| 17.9 | 10.6 |
| 16.2 | 10.7 |
| 15.4 | 11.4 |
| ... | ... |
Each day \(i\) is a point \(x^{(i)} = (x_1^{(i)}, x_2^{(i)})\): Beirut's temperature and Damascus's, in °C. The first rows are the first four days of 2010.
The points lie along a diagonal (correlation 0.963): a warm day in Beirut is a warm day in Damascus. Damascus swings more: 8.3 °C standard deviation against 5.9 °C.
Principal component 1, \(v_1\), is the direction along which the days spread the most. Here it points 55° up from the Beirut axis: \(v_1 = \) (0.57, 0.82).
Principal component 2, \(v_2\), is the direction of most spread among those at a right angle to \(v_1\).
Principal components are orthogonal unit vectors, so together they are a new set of axes: a change of coordinates. Principal component analysis (PCA) finds them.
Mean, then centre the data
Variance, then normalise the data
All cities are in °C, so this deck centres and does not normalise: Gaziantep (9.4 °C) should count more than Alexandria (5.0 °C). Normalise when the columns have different units.
Keeping one number per day, \(z_1^{(i)}\), is a projection onto \(v_1\). The error of that projection is each day's distance in the direction of \(v_2\): the reconstruction error. How do we find \(v_1\) and \(v_2\)?
Pythagoras, for every day (the formula for the length along the line comes with the exercise in part 3). Summed over the days, the left side is fixed, 103.8 °C², so the line with the most variance along it is the line with the least squared distance to it. At \(v_1\) the error is 1.7 °C², a typical distance of 1.3 °C.
Part 3 of 5
How do we find that line, and how many numbers does a day need?
| \(n\) | number of examples: days (5,844), or digits (20,000) |
| \(d\) | numbers per example: cities (22), or pixels (256) |
| \(x^{(i)} \in \mathbb{R}^d\) | example \(i\), centred: a day minus the mean day \(\mu\) |
| \(X \in \mathbb{R}^{n \times d}\) | the centred data, one example per row |
| \(\Sigma = \tfrac{1}{n} X^\top X \in \mathbb{R}^{d \times d}\) | the covariance matrix |
| \(v_j \in \mathbb{R}^d,\ \lambda_j\) | principal directions (unit length, orthogonal) and their variances, \(\lambda_1 \ge \lambda_2 \ge \dots\) |
| \(z^{(i)} \in \mathbb{R}^k\) | the \(k\) numbers kept for example \(i\): \(z_j^{(i)} = v_j^\top x^{(i)}\) |
| \(X = U\,\Xi\,V^\top\) | the singular value decomposition, singular values \(\sigma_j\) on the diagonal of \(\Xi\) |
The F25 slides called the directions \(u\) and switched to \(v\) at the SVD; here they are \(v\) from the start, because they are the columns of \(V\), not of \(U\).
You know \(x^{(i)}\) (yellow) and the unit vector \(v\) (pink). Which point on the line is closest to \(x^{(i)}\), and how far from the origin is it? One dot product.
Drop the perpendicular (red dashes): the closest point on the line is its foot. What is lost is the red length.
\(z^{(i)}\) is one number, the signed distance along \(v\). The closest point is \(z^{(i)}\, v\) (green): \(z^{(i)}\) copies of the unit vector.
Over all days, the variance of the projections is the mean of \(\big(z^{(i)}\big)^2\) (centred data: the \(z^{(i)}\) average to zero). The next slide maximises it.
So the largest eigenvalue gives the first principal component. The next one, maximised among directions orthogonal to \(v_1\), is the eigenvector with the second largest eigenvalue, and so on.
Beirut and Damascus: \(\Sigma = \begin{bmatrix} 34.9 & 47.2 \\ 47.2 & 69.0 \end{bmatrix}\) °C², \(v_1 = \) (0.57, 0.82), \(\lambda_1 = \) 102.1, \(\lambda_2 = \) 1.7; the trace 103.8 is \(\lambda_1 + \lambda_2\), the total variance.
The circle of unit vectors becomes an ellipse whose axes are the eigenvectors and whose half-lengths are the eigenvalues: a symmetric matrix never shears, it only stretches along \(v_1\) and \(v_2\) (3Blue1Brown, Essence of linear algebra, chapter 14). The weather's \(\Sigma\) is so lopsided (98.4% along \(v_1\)) that it turns almost every \(u\) onto \(v_1\).
With \(k = d\) nothing is lost, since \(V V^\top = I\). For Beirut and Damascus, \(k = 1\) keeps 98.4% of the variance: one number per day, its position along \(v_1\). If the data were standardised as well (each feature divided by its spread \(s_c\), the centring slide), PCA runs on \(x_c / s_c\): the directions change, and the reconstruction scales back, \(\hat x = \mu + s \odot V_k z\).
eigh for a symmetric matrix, and sortimport numpy as np, pandas as pd
df = pd.read_csv('levant_era5_daily_mean_temp_2010_2025.csv')
X = df.drop(columns='date').to_numpy() # n = 5844 days by d = 22 cities, deg C
X_centered = X - X.mean(axis=0) # n by d
Sigma = X_centered.T @ X_centered / len(X) # d by d covariance
lam, V = np.linalg.eigh(Sigma) # eigh: symmetric, real, ascending order
lam, V = lam[::-1], V[:, ::-1] # largest first; column j of V is v_j
k = 3
Z = X_centered @ V[:, :k] # n by k: the k numbers per day
print(np.round(lam[:3] / lam.sum(), 3)) # shares 95.5%, 1.5%, 0.9%
print(Z.shape, np.round(Z[0], 1)) # (5844, 3) and the first day's three numbers
print(np.allclose(V.T @ V, np.eye(22))) # the directions are orthonormal
np.linalg.eig returns eigenvalues in no particular order and, through rounding, sometimes with tiny imaginary parts. eigh is for symmetric matrices: real results, sorted from smallest to largest, so reverse them. Without the \(1/n\) the eigenvectors are the same and every eigenvalue is \(n\) times too large.
Write the shape of each, for this cloud and for the weather (22 cities over 5,844 days), then turn the slide:
How many eigenvalues are large, and why? What happens to \(\lambda_3\), and to \(v_3\), as the cloud thickens?
The cloud. Each point is \(x^{(i)} \in \mathbb{R}^3\), so \(\Sigma\) is \(3 \times 3\) and there are three directions \(v_1, v_2, v_3 \in \mathbb{R}^3\). The points lie close to a plane: two eigenvalues are large, the third is small, and \(v_3\) is the normal to the plane. Thicken the cloud and \(\lambda_3\) grows; at some thickness there is no plane to find.
The total, \(\sum_j \lambda_j\), is the trace of \(\Sigma\): the sum of the 22 cities' variances, 1003 °C². The first eigenvalue is 957.8 of it, 95.5%. Rules of thumb for \(k\): the elbow of this plot, a share such as 95%, or, when the numbers feed a model, the \(k\) that predicts best on held-out data.
This day is off by (rms over the cities); over all days, .
One number per day keeps 95.5% of the variance and leaves 1.43 °C of error per city, against 6.75 °C for the mean day alone. Three numbers bring it under a degree (0.97 °C).
Yellow: \(z_1\) for every day of 2024, the 22 centred temperatures weighted by \(v_1\) and summed. Blue: its monthly means over all 16 years. Every weight is positive: \(z_1\) is how warm the whole region is, low in January (-42) and high in August (41). The weights follow each city's yearly swing (correlation 0.995): largest at Gaziantep (0.30), smallest at Alexandria (0.16).
\(U^\top U = I\) and \(V^\top V = I\): orthonormal columns
Columns of \(V\): the principal directions (\(v_1\) is Vt[0]). \(\sigma_j\): the singular values, largest first. The numbers kept: \(Z = XV = U\Xi\). Why this \(V\) is the eigenvalue problem's: next slide.
For the weather, \(\sigma_1 = \) 2366 and \(n = \) 5,844, so \(\lambda_1 = \sigma_1^2/n = \) 957.8 °C², the first bar of the scree plot. Computing \(X^\top X\) squares the condition number and loses digits of accuracy; the SVD works on \(X\) itself, which is why it is the standard route (Brunton and Kutz, 2022).
import numpy as np, pandas as pd
df = pd.read_csv('levant_era5_daily_mean_temp_2010_2025.csv')
Xc = df.drop(columns='date').to_numpy(); Xc = Xc - Xc.mean(axis=0) # 5844 by 22, centred
U, s, Vt = np.linalg.svd(Xc, full_matrices=False) # economy: U is 5844 by 22, s has 22 values, Vt is 22 by 22
for k in [1, 2, 3, 22]:
Xk = (U[:, :k] * s[:k]) @ Vt[:k] # k columns of U, k singular values, k rows of Vt
print(k, round(((Xc - Xk)**2).sum()), round((s[k:]**2).sum())) # the two numbers agree
Eckart and Young (1936): among all matrices of rank \(k\), \(X_k = U_k \Xi_k V_k^\top\) is the closest to \(X\) in squared error, and the error left is \(\|X - X_k\|_F^2 = \sum_{j > k} \sigma_j^2\). The PCA reconstruction with \(k\) components is exactly \(X_k\); dividing by \(n\) turns the dropped \(\sigma_j^2\) into the dropped eigenvalues of the scree plot. The three forms, in the shape of \(U\): full \(n \times n\), economy \(n \times d\) (full_matrices=False), truncated \(n \times k\).
import numpy as np, pandas as pd
df = pd.read_csv('levant_era5_daily_mean_temp_2010_2025.csv')
X = df.drop(columns='date').to_numpy() # 5844 days by 22 cities, deg C
Xc = X - X.mean(axis=0) # centre each city
U, s, Vt = np.linalg.svd(Xc, full_matrices=False) # rows of Vt are v_1, v_2, ...
lam = s**2 / len(X) # eigenvalues of the covariance
print(np.round(lam[:3] / lam.sum(), 3)) # 95.5%, 1.5%, 0.9%
for k in [0, 1, 2, 3]:
Xk = (U[:, :k] * s[:k]) @ Vt[:k] # rank-k reconstruction
print(k, round(np.sqrt(((Xc - Xk)**2).mean()), 2)) # rms error, deg C
v1 = Vt[0] * np.sign(Vt[0].sum()) # a direction has no sign: make v_1 positive
print(df.columns[1:][v1.argmax()], df.columns[1:][v1.argmin()]) # largest and smallest weight
Expected: shares 95.5%, 1.5%, 0.9%; errors 6.75 °C, 1.43 °C, 1.17 °C, 0.97 °C; then Gaziantep and Alexandria. A direction is defined only up to its sign, so libraries can return \(-v_j\); fix a convention before you read the weights.
A dataset has \(d = 6\) features. The eigenvalues of its covariance matrix are on the left.
One example, one row of \(X\), is:
Then the shapes (rows × columns). \(n\) and \(d\) are in the text; count, do not guess.
Over a dataset, the variance kept by \(v\) is the mean of \(z^2\) and the variance lost is the mean of (c): the two bars of the best-axis slide.
Part 4 of 5
What does PCA find in pictures of handwriting?
Introduction to Machine Learning
Here \(n = \) 20,000 digits and \(d = \) 256 pixels, values from 0 (black) to 1 (white). Each \(v_j\) has 256 entries, one per pixel, so it can be drawn as a picture: an eigen-digit. The F25 slide did this with faces (eigenfaces, Turk and Pentland, 1991); these are MNIST digits (LeCun and colleagues, 1998), shrunk to 16 by 16.
Pixels are all in the same unit, so the digits are centred and not normalised; 14 border pixels are black on every digit, have no spread, and could not be divided by it anyway. The first direction keeps 11.3% of the variance, far less than the weather's 95.5%: a digit has more independent parts than a day.
the digit, 256 numbers
Here each \(v_j\) is an eigen-digit, a picture of 256 numbers, and \(z_j = v_j^\top (x - \mu)\) is one number: how much of eigen-digit \(j\) this digit contains.
\(\hat x = \mu + \sum_{j=1}^k z_j v_j\): the mean digit plus \(k\) eigen-digits, each weighted by this digit's own score, the reconstruction of part 3 with \(d = 256\). These 20 digits were held out from the reader that judges the rebuilt digits; the eigen-digits themselves were found from all 20,000.
Two criteria for \(k\). Blue is the usual one, the share of the variance kept: 51 components for 90%, 125 for 99%. Yellow is the task: a reader is shown digits rebuilt from \(k\) components and we count how often it names them right. At \(k = 30\) only 81% of the variance is kept, yet 95% of the rebuilt digits are read correctly, against 96.6% for the originals.
How the numbers were made. Computed for this deck (scripts/verify_pca_svd.py in the course repository) on the MNIST digits of LeCun, Bottou, Bengio and Haffner (1998), shrunk to 16 by 16: PCA on 20,000 digits; a 3-nearest-neighbour reader fitted on 16,000 full digits; the other 4,000 rebuilt from \(k\) components and shown to it. Nothing here is quoted from a paper.
The variance thrown away is mostly the small wobbles of each person's pen, not the shape that makes a 4 a 4. Variance kept measures pixels, not meaning: for a task, choose \(k\) by how well the task goes.
Zeros gather at one end of \(z_1\) (their mean is +2.0) and ones at the other (−2.0); the other eight digits overlap in the middle, every mean within 0.4 of zero. Two numbers are enough to see some structure and too few to separate ten digits.
PCA keeps the directions of largest spread, which is not the same as the directions that separate classes. Methods that bend the map, such as t-SNE and UMAP, come with autoencoders.
PCA is one linear map from pictures to a few numbers and one linear map back. An autoencoder replaces both maps with networks: the autoencoder in the architectures deck is the place to see it next.
Part 5 of 5
What do the patterns mean, and what can they not do?
ML for Science
This pattern keeps .
The average year explains 87.9% of the variance, less than PCA's first number (95.5%), which also follows whole months that run warm or cool everywhere. On the anomalies (previous slide, second button) pattern 1 has one sign everywhere and keeps 71.9%: a warm or cold spell over the region, strongest inland (Zahle) and weakest at Alexandria. Pattern 2 is the same west against east seesaw (8.6%); pattern 3 is Cairo +0.40 against Gaziantep −0.42 (5.2%). 5 patterns keep 90% of the anomalies.
Athens's anomaly matches Damascus's two days later better (0.50) than on the same day (0.34): weather systems cross the region from west to east. In PCA's terms, \(z_2\) today predicts \(z_1\) two days later (correlation 0.37), although on the same day the two are uncorrelated by construction.
Each PCA term is a fixed map times a number that changes in time, \(z_j(t)\, v_j\). A travelling pattern needs several such terms, and PCA ranks them by variance, not by how they move. Dynamic mode decomposition looks for patterns that each move in time in one simple way.
| this deck | dynamic mode decomposition | for the weather |
|---|---|---|
| \(d\), numbers per example | \(n\), the state size | 22 cities |
| \(n\), examples | \(m\), snapshots in time | 5,844 days |
| \(v_j\), columns of \(V\) | \(\mathbf{u}_j\), columns of \(U\): spatial modes | the maps of the last slides |
| \(\Xi\), singular values \(\sigma_j\) | \(\Sigma\), the same \(\sigma_j\) | \(\sigma_1\) is the season |
| \(\sigma_j u_j\), the scores \(z_j\) | \(\sigma_j \mathbf{v}_j^\top\), rows of \(\Sigma V^\top\): temporal modes | \(z_1\) through the year |
There \(\Sigma\) holds the singular values; it is not a covariance. Climate scientists call the spatial modes empirical orthogonal functions, fluid dynamicists POD modes: both are this SVD.
np.linalg.svd(X, full_matrices=False) and plot \(\sigma_j^2 / \sum \sigma^2\). How many patterns keep 90%?Shuffle the 5,844 days and PCA returns the same patterns: it never used their order. Dynamic mode decomposition uses the order, by fitting the linear map that takes each day to the next, \(\mathbf{x}_{t+1} \approx A\,\mathbf{x}_t\), and starts from this SVD to do it.
That is the next deck, linear dynamical systems and DMD, on the wake behind a cylinder, where the SVD's spatial modes \(\mathbf{u}_j\) and temporal modes \(\sigma_j \mathbf{v}_j^\top\) appear in the notation of the last slide.