24 - Linear Algebra 2: Eigenvalues, SVD and PCA

Machine Learning Engineer course · Module 02 · Mathematics and statistics · Lesson 9

Goal. Find the special directions of a matrix (eigenvectors), understand that every matrix is "a rotation, a stretching along axes, another rotation" (the singular value decomposition, SVD), use it to compress data and images with the best low-rank approximation, and arrive at PCA (principal component analysis), deriving it from scratch and comparing it with scikit-learn.

Myth to dispel. "PCA chooses the most important variables." No: it does not select anything; it builds new variables, linear combinations of all the original ones, ordered by how much variance they carry. The difference changes how the result is interpreted.

Introduction

In the previous lesson we saw that a matrix is a transformation of space. Now we ask: are there directions in which the transformation behaves in a simple way? And how do we describe any rectangular matrix (for example an n × d dataset) so as to find its "main structure"? The answers are eigenvalues and the SVD, and their most famous application to data is PCA.

The path defines eigenvalues and eigenvectors, moves to symmetric matrices and the spectral theorem, builds the SVD and the low-rank approximation, and arrives at PCA by deriving it from scratch. It closes by saying what PCA can and cannot do. The prerequisites are the lesson on vectors and matrices (rank, conditioning, normal equations) and the covariance; dimensionality reduction applied to models is taken up again in PCA from scratch and t-SNE.

Eigenvalues and eigenvectors

The analogy: the rubber sheet

Imagine stretching a rubber sheet. Most points move in a direction different from the one in which they were with respect to the centre (they rotate). But there are special directions along which points merely move away from or towards the centre, without rotating: they stay on the same line. Those directions are the eigenvectors; the stretching factor along each is the corresponding eigenvalue.

Definition. v ≠ 0 is an eigenvector of A with eigenvalue λ if

A v = λ v.

The effect of A on v is only to multiply it by the number λ (stretching if |λ| > 1, compression if |λ| < 1, flipping if λ < 0).

How they are computed

A v = λ v is equivalent to (A − λI) v = 0, where I is the identity. This homogeneous system has a solution v ≠ 0 only if the matrix A − λI is singular (previous lesson: otherwise the only solution would be v = 0), that is if

det(A − λI) = 0 (the characteristic polynomial).

For A = [[a, b], [c, d]]: det(A − λI) = (a − λ)(d − λ) − bc = λ² − (a + d)λ + (ad − bc) = 0, hence

λ² − tr(A)·λ + det(A) = 0

with tr(A) = a + d the trace (sum of the elements on the diagonal). From the quadratic formula two facts can be read: the sum of the eigenvalues is the trace and the product is the determinant.

Example. A = [[2, 1], [1, 2]]: tr = 4, det = 3, hence λ² − 4λ + 3 = 0, λ = 3 or λ = 1. For λ = 3: (A − 3I)v = [[−1, 1], [1, −1]]v = 0 gives v ∝ (1, 1). For λ = 1: [[1, 1], [1, 1]]v = 0 gives v ∝ (1, −1).

import numpy as np

A = np.array([[2.0, 1.0],
              [1.0, 2.0]])
autovalori, autovettori = np.linalg.eig(A)                # eigenvectors in COLUMNS
print("eigenvalues:", autovalori)
print("eigenvectors (columns):\n", np.round(autovettori, 4))

for i in range(2):
    v = autovettori[:, i]
    print(f"A v{i+1} = {np.round(A @ v, 4)} | λ{i+1} v{i+1} = {np.round(autovalori[i] * v, 4)}")

print("sum of the eigenvalues =", autovalori.sum(), "= trace", np.trace(A))
print("product of the eigenvalues =", round(float(autovalori.prod()), 6), "= det", round(float(np.linalg.det(A)), 6))

Output:

eigenvalues: [3. 1.]
eigenvectors (columns):
 [[ 0.7071 -0.7071]
 [ 0.7071  0.7071]]
A v1 = [2.1213 2.1213] | λ1 v1 = [2.1213 2.1213]
A v2 = [-0.7071  0.7071] | λ2 v2 = [-0.7071  0.7071]
sum of the eigenvalues = 4.0 = trace 4.0
product of the eigenvalues = 3.0 = det 3.0

What the code does. np.linalg.eig returns two arrays: the eigenvalues and a matrix whose columns are the eigenvectors (normalised to length 1). The loop verifies the definition Av = λv for each pair: the two vectors coincide. The last lines confirm the sum (trace) and the product (determinant). Careful: eig may return the signs of the eigenvectors one way or the other (if v is an eigenvector so is −v).

What they are for

Powers of a matrix. Applying A many times: A·A·…·A v. If you write v = c₁v₁ + c₂v₂ (combination of eigenvectors) then Aᵏ v = c₁λ₁ᵏ v₁ + c₂λ₂ᵏ v₂: each component grows (or decays) like the power of its eigenvalue, and after many applications the eigenvector with the eigenvalue of largest modulus dominates. It is the power method.

Markov chains. A system that moves between states with fixed probabilities (the weather: sun or rain, with the probabilities of passing from one day to the next). The "steady-state" distribution is the eigenvector with eigenvalue 1 of the transition matrix.

def metodo_delle_potenze(M, iterazioni=50, seme=0):
    v = np.random.default_rng(seme).normal(size=M.shape[0])
    for _ in range(iterazioni):
        v = M @ v
        v /= np.linalg.norm(v)                 # renormalise to prevent it from exploding
    return v, float(v @ M @ v)                 # eigenvector and estimate of the eigenvalue (Rayleigh quotient)

v, lam = metodo_delle_potenze(A)
print("power method:", np.round(v, 4), "| estimated λ =", round(lam, 4), "(expected 3, direction (1,1)/√2 ≈ 0.7071)")

# weather: P[i, j] = probability of going from state j to state i (columns sum to 1). Invented numbers.
P = np.array([[0.8, 0.4],           # sun tomorrow | sun today, rain today
              [0.2, 0.6]])          # rain tomorrow
stato = np.array([0.0, 1.0])        # today it rains for certain
for giorno in range(1, 6):
    stato = P @ stato
    print(f"day {giorno}: P(sun) = {stato[0]:.4f}, P(rain) = {stato[1]:.4f}")

val, vec = np.linalg.eig(P)
regime = vec[:, np.argmax(val.real)].real
print("steady-state distribution (eigenvector of λ = 1, normalised to sum 1):", np.round(regime / regime.sum(), 4))

Output:

power method: [-0.7071 -0.7071] | estimated λ = 3.0 (expected 3, direction (1,1)/√2 ≈ 0.7071)
day 1: P(sun) = 0.4000, P(rain) = 0.6000
day 2: P(sun) = 0.5600, P(rain) = 0.4400
day 3: P(sun) = 0.6240, P(rain) = 0.3760
day 4: P(sun) = 0.6496, P(rain) = 0.3504
day 5: P(sun) = 0.6598, P(rain) = 0.3402
steady-state distribution (eigenvector of λ = 1, normalised to sum 1): [0.6667 0.3333]

What the code does. metodo_delle_potenze starts from a random vector and repeatedly applies M, renormalising each time (otherwise the numbers would grow like 3ᵏ); it converges to the dominant eigenvector (1, 1)/√2. The ratio vᵀMv (Rayleigh quotient) estimates the eigenvalue, close to 3. The second part is a Markov chain with invented probabilities: starting from "rain", the distribution approaches in a few days a stable value, which coincides with the eigenvector of λ = 1 of the transition matrix (2/3 sun, 1/3 rain). This same logic, on an enormous matrix of links between web pages, is at the basis of the PageRank algorithm (Brin and Page, 1998).

Symmetric matrices and the spectral theorem

A matrix is symmetric if A = Aᵀ (equal to its transpose). They are everywhere in data: the covariance matrix, XᵀX, distance matrices. For them a very strong result holds, the spectral theorem:

1. The eigenvalues are real numbers.
2. There exists a basis of eigenvectors pairwise perpendicular.
3. Hence A = Q Λ Qᵀ, with Q an orthogonal matrix (the columns are the eigenvectors of length 1, QᵀQ = I) and Λ diagonal with the eigenvalues.

Why the eigenvectors are perpendicular (for two different eigenvalues). Let Av₁ = λ₁v₁, Av₂ = λ₂v₂ with λ₁ ≠ λ₂, and A symmetric. We compute (Av₁)·v₂ in two ways: = λ₁ (v₁·v₂) and, by symmetry, (Av₁)·v₂ = v₁·(Aᵀv₂) = v₁·(Av₂) = λ₂ (v₁·v₂). Hence (λ₁ − λ₂)(v₁·v₂) = 0 and, since λ₁ ≠ λ₂, v₁·v₂ = 0. ∎

The geometric meaning: a symmetric matrix stretches space along perpendicular axes (the axes of an ellipse, as in the figure), each with its factor λ.

rng = np.random.default_rng(1)
M = rng.normal(size=(4, 4))
S = M + M.T                                   # any symmetric matrix
lam, Q = np.linalg.eigh(S)                    # eigh: specialised for symmetric matrices
print("eigenvalues (real, sorted):", np.round(lam, 4))
print("Q^T Q = I ?                  ", np.allclose(Q.T @ Q, np.eye(4)))
print("S = Q Λ Q^T ?                ", np.allclose(Q @ np.diag(lam) @ Q.T, S))
print("trace = sum of eigenvalues   :", np.isclose(np.trace(S), lam.sum()))

Output:

eigenvalues (real, sorted): [-2.1903  0.1975  1.4328  3.3985]
Q^T Q = I ?                   True
S = Q Λ Q^T ?                 True
trace = sum of eigenvalues   : True

What the code does. It builds a random symmetric matrix (M + Mᵀ always is) and computes its eigenvalues with eigh, the version for symmetric matrices (faster, more accurate, returns the eigenvalues already sorted in ascending order and orthonormal eigenvectors). It verifies the three statements of the theorem: QᵀQ = I, S = QΛQᵀ and the property of the trace.

Why it matters. The covariance matrix of the data is symmetric and positive semidefinite (all the eigenvalues are ≥ 0, because uᵀ C u is a variance, which cannot be negative). Its eigenvectors are the principal directions of the data, and the eigenvalues the variances along them: this is PCA, shortly.

The singular value decomposition (SVD)

Eigenvalues exist for square matrices (and are not always real). For any m × n matrix (a dataset!) there is a decomposition that always exists:

A = U Σ Vᵀ

where U (m × m) and V (n × n) are orthogonal matrices (rotations, possibly with a flip) and Σ is "diagonal" (m × n) with the singular values σ₁ ≥ σ₂ ≥ … ≥ 0 on the diagonal. In geometric terms: every linear transformation is the composition of a rotation (Vᵀ), a stretching along the axes with factors σᵢ (Σ) and another rotation (U). The SVD is the circle that becomes an ellipse, in general.

Where it comes from (construction from the eigenvalues)

The matrix AᵀA is symmetric and positive semidefinite (for every x, xᵀAᵀAx = ‖Ax‖² ≥ 0). By the spectral theorem it has orthonormal eigenvectorsv₁, …, vₙwith eigenvaluesλᵢ ≥ 0. We defineσᵢ = √λᵢ`. Then:

• ‖A vᵢ‖² = vᵢᵀAᵀA vᵢ = λᵢ = σᵢ²: the vector vᵢ is stretched by σᵢ.
• The vectors A vᵢ are perpendicular to each other: (A vᵢ)·(A vⱼ) = vᵢᵀAᵀA vⱼ = λⱼ (vᵢ·vⱼ) = 0 for i ≠ j.
• So, normalising, uᵢ = A vᵢ / σᵢ (for σᵢ > 0) are orthonormal vectors, and A vᵢ = σᵢ uᵢ: this is exactly A V = U Σ, that is A = U Σ Vᵀ. ∎

The eigenvalues of AᵀA are therefore the squares of the singular values: σᵢ² = λᵢ(AᵀA). The three items say in order what A does to each vᵢ (it stretches it by σᵢ), that the images remain perpendicular to each other and that, once normalised, they are the columns of U: the construction requires nothing beyond the spectral theorem, and it is the one that the listing below retraces in numbers.

rng = np.random.default_rng(2)
A = rng.normal(size=(6, 4))                                  # any rectangular matrix

U, s, Vt = np.linalg.svd(A, full_matrices=False)             # "economy" form
print("shapes: U", U.shape, "| s", s.shape, "| Vt", Vt.shape)
print("singular values:", np.round(s, 4))
print("A = U Σ V^T ?", np.allclose(U @ np.diag(s) @ Vt, A))

# construction from the eigenvalues of A^T A
lam, V = np.linalg.eigh(A.T @ A)
ordine = np.argsort(lam)[::-1]                               # descending order
lam, V = lam[ordine], V[:, ordine]
print("sqrt(eigenvalues of A^T A):", np.round(np.sqrt(np.clip(lam, 0, None)), 4), "= singular values")

U_costr = (A @ V) / np.sqrt(lam)                             # u_i = A v_i / σ_i
print("orthonormal u_i:", np.allclose(U_costr.T @ U_costr, np.eye(4)))
print("rank =", int((s > 1e-10).sum()), "| κ(A) = σ_max/σ_min =", round(float(s[0] / s[-1]), 4), "| np.linalg.cond:", round(float(np.linalg.cond(A)), 4))

Output:

shapes: U (6, 4) | s (4,) | Vt (4, 4)
singular values: [3.0945 2.1676 1.4749 0.365 ]
A = U Σ V^T ? True
sqrt(eigenvalues of A^T A): [3.0945 2.1676 1.4749 0.365 ] = singular values
orthonormal u_i: True
rank = 4 | κ(A) = σ_max/σ_min = 8.4786 | np.linalg.cond: 8.4786

What the code does. np.linalg.svd(..., full_matrices=False) returns U, the singular values s (in a vector, already sorted in descending order) and Vt (Vᵀ). The reconstruction A = UΣVᵀ is verified. Then it retraces the construction of the proof: eigenvalues of AᵀA, square root (obtaining the singular values), and uᵢ = Avᵢ/σᵢ, which turn out orthonormal. The last lines link the SVD to the concepts of the previous lesson: the rank is the number of non-zero singular values, and the condition number is σ_max/σ_min (the definition we had anticipated).

Low-rank approximation

Writing the row-by-column product, A = Σᵢ σᵢ uᵢ vᵢᵀ: a sum of rank-1 matrices, sorted by importance (σ₁ ≥ σ₂ ≥ …). Keeping only the first k terms one obtains

Aₖ = Σᵢ₌₁ᵏ σᵢ uᵢ vᵢᵀ.

The Eckart-Young theorem (1936) says that Aₖ is the best rank-k approximation of A: among all matrices of rank at most k, it is the one with the smallest error (in the Frobenius norm, the root of the sum of the squares of all the elements). The error is

‖A − Aₖ‖_F = √(σ²ₖ₊₁ + σ²ₖ₊₂ + …).

(We state it without proving it; see Strang. We will verify it numerically.) The space occupied falls: instead of m·n numbers one stores k·(m + n + 1).

Compression of a synthetic image

We build a 64 × 64 grayscale image with simple structures: a gradient, a circle and some stripes, plus a little noise.

m = n = 64
y, x = np.mgrid[0:m, 0:n]
immagine = 0.5 * (x / n) + 0.3 * (np.hypot(y - 32, x - 32) < 14) + 0.2 * (np.sin(x / 3.0) > 0.5)
immagine = immagine + np.random.default_rng(3).normal(0, 0.02, immagine.shape)

U, s, Vt = np.linalg.svd(immagine, full_matrices=False)
print("first 8 singular values:", np.round(s[:8], 3))
print("cumulative energy of the first k (sum of σ²):")
energia = np.cumsum(s ** 2) / (s ** 2).sum()
for k in (1, 2, 3, 5, 10, 20):
    Ak = (U[:, :k] * s[:k]) @ Vt[:k]
    errore = np.linalg.norm(immagine - Ak)
    teoria = np.sqrt((s[k:] ** 2).sum())
    memoria = k * (m + n + 1) / (m * n)
    print(f"k = {k:>2}: energy {energia[k-1]:.4f} | error {errore:.4f} = theory {teoria:.4f} | memory {memoria:6.1%} of the original")

Output:

first 8 singular values: [26.343  3.873  1.667  0.941  0.721  0.558  0.465  0.412]
cumulative energy of the first k (sum of σ²):
k =  1: energy 0.9703 | error 4.6069 = theory 4.6069 | memory   3.1% of the original
k =  2: energy 0.9913 | error 2.4950 = theory 2.4950 | memory   6.3% of the original
k =  3: energy 0.9952 | error 1.8563 = theory 1.8563 | memory   9.4% of the original
k =  5: energy 0.9971 | error 1.4281 = theory 1.4281 | memory  15.7% of the original
k = 10: energy 0.9984 | error 1.0649 = theory 1.0649 | memory  31.5% of the original
k = 20: energy 0.9992 | error 0.7609 = theory 0.7609 | memory  63.0% of the original

What the code does. np.mgrid generates the coordinates of each pixel; the image is the sum of a horizontal gradient, a disc and stripes, with small noise. The SVD gives its singular values. The loop reconstructs Aₖ with the first k terms ((U[:, :k] * s[:k]) @ Vt[:k] is Σ σᵢ uᵢ vᵢᵀ written with a matrix product), measures the error and compares it with the formula √(Σ_{i>k} σᵢ²) of the theorem: they coincide. The cumulative energy Σᵢ≤ₖ σᵢ²/Σσ² says which fraction of the "substance" of the image is captured. With very few terms (k = 5) the error is already small and the memory falls to 15.7% of the original (k·(m + n + 1)/(m·n)) and the energy captured is 99.7%: the image is "almost low rank" because it has simple structures. A real image rich in detail decays more slowly.

Where it is used. Compression, noise reduction (small σ are discarded, often noise), recommendation systems (the users × products matrix approximated with a few latent factors), latent semantic analysis of texts.

PCA

The analogy: the right shadow

You have a three-dimensional object and project it onto a wall with a light: the shadow is a 2D version. Depending on the orientation of the wall, the shadow keeps more or less information: a pencil seen tip-on is a dot, seen from the side it is a long segment. PCA looks for the orientation of the wall that maximises the dispersion of the shadow, that is keeps as much "shape" as possible. The best 1D wall is the direction of maximum variance (the first principal component), the second axis is the direction of maximum variance perpendicular to the first, and so on.

The derivation

Given n examples with d features, let X be the centred n × d matrix (from each column its mean has been subtracted). The covariance matrix is C = XᵀX/(n − 1) (d × d, symmetric, positive semidefinite).

The variance of the data projected on the direction u (a unit vector, ‖u‖ = 1) is the variance of the numbers Xu: (1/(n−1))·‖Xu‖² = uᵀ C u.

We want to maximise uᵀCu with the constraint uᵀu = 1. At the maximum the gradient of the function must be parallel to the gradient of the constraint (otherwise we could move along the constraint and increase the function: this is the method of Lagrange multipliers). The gradient of uᵀCu is 2Cu (C symmetric), that of uᵀu is 2u. Hence

C u = λ u:

u must be an eigenvector of the covariance matrix, and the variance obtained uᵀCu = λ uᵀu = λ is its eigenvalue. To maximise it one takes the eigenvector with the largest eigenvalue. The second component is the eigenvector of the second eigenvalue, perpendicular to the first by the spectral theorem; and so on. The principal components are therefore the eigenvectors of C, and the explained variances are its eigenvalues.

The link with the SVD. If X = UΣVᵀ, then XᵀX = VΣ²Vᵀ: the eigenvectors of C are the columns of V and the eigenvalues are σᵢ²/(n−1). Therefore PCA is computed in a numerically more stable way with the SVD of the centred data, without forming XᵀX (which would double the conditioning: previous lesson). It is what scikit-learn does.

rng = np.random.default_rng(12)
n = 500
dati = rng.normal(size=(n, 2)) @ np.array([[2.2, 1.4], [0.0, 0.6]])         # cloud elongated along the diagonal
Xc = dati - dati.mean(axis=0)                                               # centring

# (a) PCA with the eigenvalues of the covariance
C = Xc.T @ Xc / (n - 1)
lam, V = np.linalg.eigh(C)
ordine = np.argsort(lam)[::-1]
lam, V = lam[ordine], V[:, ordine]

# (b) PCA with the SVD of the centred data
U, s, Vt = np.linalg.svd(Xc, full_matrices=False)
var_svd = s ** 2 / (n - 1)

# (c) scikit-learn
from sklearn.decomposition import PCA
pca = PCA(n_components=2).fit(dati)

print("variances (eigenvalues of C)   :", np.round(lam, 4))
print("variances (σ²/(n-1) from SVD)  :", np.round(var_svd, 4))
print("scikit-learn variances         :", np.round(pca.explained_variance_, 4))
print("explained variance ratio       :", np.round(pca.explained_variance_ratio_, 4))
print("first component (eigh)   :", np.round(V[:, 0], 4))
print("first component (SVD)    :", np.round(Vt[0], 4))
print("first component (sklearn):", np.round(pca.components_[0], 4), "(the sign may be opposite)")

punteggi = Xc @ V                                                           # data in the new coordinates
print("variance of the scores per component:", np.round(punteggi.var(axis=0, ddof=1), 4))
print("covariance between the components (≈ 0):", round(float(np.cov(punteggi.T)[0, 1]), 10))

Output:

variances (eigenvalues of C)   : [5.8321 0.2833]
variances (σ²/(n-1) from SVD)  : [5.8321 0.2833]
scikit-learn variances         : [5.8321 0.2833]
explained variance ratio       : [0.9537 0.0463]
first component (eigh)   : [-0.8276 -0.5613]
first component (SVD)    : [-0.8276 -0.5613]
first component (sklearn): [0.8276 0.5613] (the sign may be opposite)
variance of the scores per component: [5.8321 0.2833]
covariance between the components (≈ 0): 0.0

What the code does. It generates a cloud of 500 points elongated along the diagonal and centres it (an indispensable step: PCA works on the dispersion around the mean). It computes the components in three ways: (a) eigenvalues and eigenvectors of the covariance with eigh (then sorted from the largest), (b) SVD of the centred data (the variances are σ²/(n−1)), (c) sklearn.decomposition.PCA. The three variances coincide and the components coincide up to the sign (an eigenvector and its opposite are both valid). Finally it projects the data onto the components (Xc @ V): the variances of the scores are precisely the eigenvalues, and the covariance between the components is practically zero: the new coordinates are uncorrelated.

A real example: the iris flowers

The flower dataset (150 examples, 4 measurements) is included in scikit-learn. Let us see how many dimensions are really needed, with and without standardisation.

from sklearn.datasets import load_iris
from sklearn.preprocessing import StandardScaler

Xi, yi = load_iris(return_X_y=True)

for etichetta, dati_i in (("unstandardised", Xi), ("standardised", StandardScaler().fit_transform(Xi))):
    p = PCA().fit(dati_i)
    print(f"{etichetta:>21}: explained variance per component = {np.round(p.explained_variance_ratio_, 4)} | first 2: {p.explained_variance_ratio_[:2].sum():.4f}")

p2 = PCA(n_components=2).fit(StandardScaler().fit_transform(Xi))
print("components (rows) on the 4 measurements, standardised:\n", np.round(p2.components_, 3))

# reconstruction from the first 2 components: how much is lost?
Z = StandardScaler().fit_transform(Xi)
rico = p2.inverse_transform(p2.transform(Z))
print("mean squared reconstruction error with 2 components:", round(float(((Z - rico) ** 2).mean()), 4))

Output:

       unstandardised: explained variance per component = [0.9246 0.0531 0.0171 0.0052] | first 2: 0.9777
         standardised: explained variance per component = [0.7296 0.2285 0.0367 0.0052] | first 2: 0.9581
components (rows) on the 4 measurements, standardised:
 [[ 0.521 -0.269  0.58   0.565]
 [ 0.377  0.923  0.024  0.067]]
mean squared reconstruction error with 2 components: 0.0419

What the code does. For the four measurements of the flowers it computes the share of variance explained by each component, both on the raw data and on the standardised data (mean 0, standard deviation 1: lesson 5). The two versions give different results: PCA depends on scale, because the variance of a variable measured in large units weighs more. Then it shows the weights (the "recipes") of the first two components on the four measurements and measures the reconstruction error by going back with inverse_transform: with two components (instead of four) little information is lost.

How to choose how many components to keep

• Cumulative variance: keep the components needed to reach a threshold (for example 95%), and declare the threshold.
• Scree plot: the variance explained per component; one looks for the "elbow" where it drops abruptly.
• Downstream purpose: if PCA precedes a model, choose k by the cross-validation of the final performance, inside the Pipeline (lesson 8).

The three criteria go from the most descriptive to the most operational. The cumulative threshold and the elbow look only at the data and are arbitrary, but communicable; cross-validation looks at the goal and is usually the one that counts, because the share of variance does not say whether the discarded information was needed for the prediction.

What PCA can and cannot do

• It is linear: it captures structures that are described with straight directions. On data arranged along a spiral or a circle it does not find the structure.
• It depends on scale: standardise (unless the units have a common meaning).
• It maximises variance, not "relevance": the direction that separates the classes may have little variance. For classification there are supervised methods.
• The components are combinations of all the variables: they are harder to interpret than the original variables.
• It is trained on the training set only: if you compute PCA on the whole dataset before splitting, you have the data leakage of lesson 8. Put PCA in a Pipeline.
• The signs are arbitrary: do not compare the signs of components computed at different times.

The limits come down to two questions to ask before using PCA: do the data really have a linear structure, and is variance what matters downstream? If the answer to either is no (spiral, classes separated by directions of little variance), PCA compresses well and informs badly, and a non-linear or supervised method is preferable.

Typical mistakes (symptom, cause, remedy)

• Not centring the data. The first component points towards the mean instead of along the dispersion. Remedy: Xc = X - X.mean(axis=0) (scikit-learn does it by itself).
• Not standardising variables on different scales. Remedy: StandardScaler.
• eig instead of eigh for symmetric matrices: it does not exploit the symmetry, does not sort the eigenvalues and, if the matrix is symmetric only up to rounding, may return complex ones. Remedy: eigh.
• Slightly negative eigenvalues (−1e-17) in a covariance. They are rounding errors: np.clip(lam, 0, None).
• PCA on the whole dataset before the split. Remedy: Pipeline.
• Interpreting the components as "the important variables". Remedy: look at the weights (components_), remembering that they mix everything.
• Comparing the signs of eigenvectors. Remedy: normalise the convention (for example largest component positive).
• Using the full SVD when the economy one is enough on large matrices: waste of memory. Remedy: full_matrices=False or truncated algorithms.

The eight mistakes fall into three groups. Data preparation: centring, standardising, doing PCA after the split. Choice of tool: eigh instead of eig, economy SVD. Reading of the results: negative eigenvalues from rounding, signs, components mistaken for variables. The first group changes the numbers, the other two change the conclusions.

How to work in practice

Start from the plot of the data and of the explained variance. Standardise, centre and only then apply PCA, inside a pipeline, computing it with the SVD without forming XᵀX. Always verify the reconstruction A ≈ UΣVᵀ and the orthogonality of the components, declare the variance threshold and the number of components chosen, and for symmetric matrices use eigh.

The vocabulary. An eigenvector is a direction that the matrix does not rotate and the eigenvalue the stretching factor along it; the trace is the sum of the diagonal elements, equal to the sum of the eigenvalues. The spectral theorem writes a symmetric matrix as QΛQᵀ; the singular values are the square roots of the eigenvalues of AᵀA and the SVD is A = UΣVᵀ, that is rotation, stretching, rotation. A principal component is an eigenvector of the covariance matrix, the explained variance the share of the total variance it carries, the Lagrange multiplier the technique for maximising with a constraint.

Mini-project

1. Apply PCA to 3 synthetic clouds: (a) elongated along the diagonal, (b) spherical, (c) two distant clusters. For each print the explained variances and comment on what it means for the effective dimension.
2. With scikit-learn's digits dataset (load_digits, 8×8 images, 64 features), find the number of components needed to explain 90% of the variance, and compare the accuracy of a logistic regression with 64 features and with those components, using a Pipeline(StandardScaler, PCA, LogisticRegression) and cross-validation.
3. Compress a grayscale image of your choice (or the synthetic one of the lesson, enlarged) with k = 1, 5, 20 and report error and memory.

The first point trains you to read the explained variance as a measure of the effective dimension: you expect a single dominant component in case (a), all equal in (b) and one dominant component that separates the groups in (c). The second is an honest comparison, because PCA is inside the cross-validation; the third takes up the compression of the lesson with an image of your own.

Challenge

• Easy. Compute by hand the eigenvalues and eigenvectors of [[4, 1], [2, 3]] and verify them with np.linalg.eig.
• Medium. Prove that the eigenvalues of Aᵏ are λᵏ and that a matrix is invertible if and only if no eigenvalue is zero (hint: det A = Π λᵢ).
• Hard. Prove that the first principal component minimises the sum of the squared distances of the points from the line through the mean: "maximising the projected variance" is equivalent to "minimising the reconstruction error" (hint: Pythagoras, for each point ‖x‖² = (u·x)² + ‖x − (u·x)u‖²).
• Expert. Implement PCA with the power method and deflation (remove the first component and repeat) and compare it with eigh.
• Enterprise. Take a dataset from your organisation with at least 15 numeric features. Apply PCA inside a pipeline, estimate how many components are needed for 90% of the variance, and verify whether a model trained on the components has performance comparable to that on the original variables.

The first two challenges are exercises on eigenvalues that can be checked with eig; the Hard one is the most instructive, because it proves that maximising the projected variance and minimising the reconstruction error are the same problem. The Expert one makes you write by hand what PCA hides.

What you can do now

You can compute eigenvalues and eigenvectors of a 2×2 matrix with the characteristic polynomial, prove that those of a symmetric matrix are perpendicular, build the SVD from the eigenvalues of AᵀA and use the Eckart-Young theorem to compress an image while measuring the error. You can derive PCA with Lagrange multipliers, compute it with the SVD of the centred data and compare it with scikit-learn, and you know when not to trust its results. You now have all the bricks of the mathematics module: the next step is to use them in a model, starting from The ML project process.

The same series also includes: 25 - Matrix Calculus and Automatic Differentiation: From the Jacobian to an Autograd Written From Scratch (the next lesson), 26 - Constrained Optimisation: Lagrange Multipliers, Duality and the KKT Conditions, 29 - Sampling and Monte Carlo Methods: Estimating What You Cannot Compute, 30 - Markov Chains: The Property That Forgets the Past, and the Distribution That Stops Changing (which refer back to this lesson).