39 - Dimensionality reduction: PCA from scratch and t-SNE

A lamp projects onto the wall the shadow of a three-dimensional object. If the lamp is pointed from the wrong side, the shadow of a bottle is a circle and you cannot tell what it is; if you point it from the right side, the shadow keeps almost all of the shape. Dimensionality reduction looks for the right side from which to “project” data with many features onto a few dimensions, losing as little as possible. It is used to compress, to speed up other models, to remove redundancy and, with 2 or 3 dimensions, to look at the data.

Two tools, two philosophies. PCA is linear: it rotates the axes to put the most variance on the first ones. t-SNE is nonlinear and aims at a single goal, visualization: preserving neighbourhoods. We implement both by hand and also see the interpretation errors they induce.

The ingredients, one by one

• Variance: how much the values of a feature spread around their mean. A feature with zero variance is a constant and carries no information.
• Covariance: how much two features vary together. Positive if when one grows the other grows too, negative if one grows and the other falls, zero if they are uncorrelated. The covariance matrix C collects all the pairwise covariances; on the diagonal are the variances.
• Centring the data: subtracting from each column its mean, so that the cloud of points is centred on the origin.
• Eigenvector and eigenvalue: a matrix transforms vectors by rotating and stretching them. An eigenvector is a special direction that the matrix does not rotate, only lengthens or shortens; the eigenvalue is the stretching factor. For the covariance matrix, the eigenvectors are the natural axes of the data cloud and the eigenvalues say how “long” the cloud is along each of them.

The four items chain into a single procedure: you centre the data, compute the covariance matrix, look for its eigenvectors, and the eigenvalues tell you how much variance each axis carries. If you understand this sequence, PCA is just its application; the rest of the lesson checks it with numbers.

The drawing: a cloud elongated along the diagonal. The first principal component (PC1, red) is the direction along which the cloud is most extended; the second (PC2, green) is perpendicular. Reducing to one dimension means projecting each point onto PC1: only the dispersion along PC2 is lost, and it is small.

PCA: finding the direction of maximum variance

Given X with n rows and d columns, already centred (the mean of each column subtracted). We look for a direction, that is, a unit vector w (‖w‖ = 1), along which the projected data have maximum variance. The projection of a point is x·w and the variance of the projections is

Var(Xw) = (1/(n−1)) · (Xw)ᵀ(Xw) = wᵀ C w          with C = XᵀX / (n−1)   (covariance matrix)

The problem is: maximize wᵀCw subject to wᵀw = 1. With a Lagrange multiplier λ:

L(w, λ) = wᵀCw − λ (wᵀw − 1)
∂L/∂w = 2Cw − 2λw = 0     →     C w = λ w

The stationary directions are therefore the eigenvectors of C. And the variance along each is wᵀCw = wᵀ(λw) = λ: the eigenvalue is the variance. The first principal component is the eigenvector with the largest eigenvalue; the second is the next one (orthogonal to the first, because C is symmetric), and so on. The share of variance explained by component i is λ_i / Σλ.

To reduce to k dimensions: Z = X W_k, where W_k has the first k eigenvectors as columns. To reconstruct: X̂ = Z W_kᵀ (+ the mean). It can be proved that this choice also minimizes the squared reconstruction error: the two points of view (maximum variance, minimum error) coincide. Here we check it numerically instead of proving it.

Implementation by hand on handwritten digits

The Digits dataset contains 1797 8×8 images of digits, that is, 64 features per image.

import numpy as np
from sklearn.datasets import load_digits, load_wine
from sklearn.decomposition import PCA

X, y = load_digits(return_X_y=True)
print("shape:", X.shape)

mu = X.mean(axis=0)
Xc = X - mu                                         # 1) centre
C = Xc.T @ Xc / (len(X) - 1)                        # 2) covariance matrix (64 x 64)
valori, vettori = np.linalg.eigh(C)                 # 3) eigenvalues/eigenvectors (C symmetric)
ordine = np.argsort(valori)[::-1]                   #    largest first
valori, vettori = valori[ordine], vettori[:, ordine]
quota = valori / valori.sum()

sk = PCA().fit(X)
print("variance shares equal to sklearn:", np.allclose(quota, sk.explained_variance_ratio_))
print("maximum difference between components (in absolute value):", np.abs(np.abs(vettori[:, :5].T) - np.abs(sk.components_[:5])).max())
print("first 5 variance shares:", np.round(quota[:5], 3))

shape: (1797, 64)
variance shares equal to sklearn: True
maximum difference between components (in absolute value): 3.469446951953614e-15
first 5 variance shares: [0.149 0.136 0.118 0.084 0.058]

In the listing the data are centred (Xc), the 64×64 covariance matrix C is built with the formula XcᵀXc/(n−1), and the eigenvalue problem is solved with np.linalg.eigh, the variant for symmetric matrices (faster and more stable than eig, and it guarantees real eigenvalues and orthogonal eigenvectors). eigh returns the eigenvalues in increasing order, so they are reordered with argsort(...)[::-1]. The variance share of component i is valori[i] / valori.sum(). The three prints compare with sklearn.decomposition.PCA: the shares coincide and the components differ by about 10⁻¹⁵ in absolute value.

The comparison uses the absolute value because an eigenvector is defined up to its sign: w and −w are the same direction, and different libraries may return one or the other. It is not an error, but it can surprise you if you compare the numbers by eye.

How many components to keep

cum = np.cumsum(quota)
for soglia in (0.80, 0.90, 0.95, 0.99):
    print("components to explain %2d%% of the variance: %d out of 64" % (soglia * 100, np.searchsorted(cum, soglia) + 1))

for k in (2, 5, 10, 20, 40):
    W = vettori[:, :k]
    errore = ((Xc - Xc @ W @ W.T) ** 2).sum(axis=1).mean()
    teorico = (valori[k:].sum()) * (len(X) - 1) / len(X)
    print("k=%2d  mean reconstruction error %8.2f   discarded variance %8.2f" % (k, errore, teorico))

components to explain 80% of the variance: 13 out of 64
components to explain 90% of the variance: 21 out of 64
components to explain 95% of the variance: 29 out of 64
components to explain 99% of the variance: 41 out of 64
k= 2  mean reconstruction error   858.94   discarded variance   858.94
k= 5  mean reconstruction error   546.72   discarded variance   546.72
k=10  mean reconstruction error   314.51   discarded variance   314.51
k=20  mean reconstruction error   126.99   discarded variance   126.99
k=40  mean reconstruction error    14.17   discarded variance    14.17

With np.searchsorted(cum, soglia) you find the first position where the cumulative variance exceeds the threshold (hence the +1 to count the components). Then for each k, W is built with the first k eigenvectors (a 64×k matrix), the data are projected (Xc @ W: n×k), reconstructed (@ W.T: n×64 again) and the mean squared error per sample is measured; it is compared with the sum of the discarded eigenvalues, multiplied by (n−1)/n to align the normalization.

The reconstruction error coincides with the variance of the discarded components: it is the numerical proof that maximum variance and minimum error are the same thing. With 29 components out of 64, 95% of the variance is kept: the digits have a lot of redundancy between adjacent pixels.

S = [ 5  2 ]
    [ 2  2 ]

PCA as a step in a pipeline

Reducing the dimensions does not always make a classifier worse, and it can speed it up. Let us measure it with cross-validation, with PCA inside the pipeline (so that it is redone on each fold, without leakage):

from sklearn.pipeline import make_pipeline
from sklearn.preprocessing import StandardScaler
from sklearn.linear_model import LogisticRegression
from sklearn.model_selection import cross_val_score

for k in (2, 5, 10, 20, 40, 64):
    pipe = make_pipeline(StandardScaler(), PCA(n_components=k), LogisticRegression(max_iter=2000))
    print("components %2d  CV accuracy %.3f" % (k, cross_val_score(pipe, X, y, cv=5).mean()))

components  2  CV accuracy 0.535
components  5  CV accuracy 0.771
components 10  CV accuracy 0.840
components 20  CV accuracy 0.899
components 40  CV accuracy 0.914
components 64  CV accuracy 0.920

The Pipeline chains standardization, PCA with n_components=k and logistic regression (a linear classifier, described in logistic regression); cross_val_score with cv=5 trains and evaluates it five times on different splits, and .mean() reports the average. Every time the scaler and the PCA see only the training part of that fold.

With 2 components the accuracy is 0.535, with 10 it rises to 0.840, with 20 to 0.899 and with 40 to 0.914, against 0.920 with all 64: the last 24 components add half a point, the rest was largely redundancy. The number of components is a hyperparameter, to be chosen in validation.

The scale trap

PCA looks for variance, and variance depends on the unit of measure. A column in thousands of units has a gigantic variance and takes the first component, whatever it says. On the Wine dataset:

Xw, yw = load_wine(return_X_y=True)
print("without scaling, share of the first 3 components:", np.round(PCA().fit(Xw).explained_variance_ratio_[:3], 4))
print("with standardization:                            ", np.round(PCA().fit(StandardScaler().fit_transform(Xw)).explained_variance_ratio_[:3], 4))

without scaling, share of the first 3 components: [9.981e-01 1.700e-03 1.000e-04]
with standardization:                             [0.362  0.1921 0.1112]

The listing compares the explained-variance shares of the first three components with and without standardization. The difference is the sign that, without scaling, PCA is ordering the features by unit of measure, not by information content.

Without scaling the first component “explains” 99.8% of the variance, but it is just the column proline (values in the thousands) disguised as a component. After standardization the variance is spread among the components (36%, 19%, 11%). The rule: standardize when the features have different units; do not do it only when the units are homogeneous and the scale has a meaning (for example the pixels of an image).

t-SNE: preserving neighbourhoods to visualize

PCA is linear: if the data lie on a curved surface, the linear projection squashes it. t-SNE (t-distributed Stochastic Neighbor Embedding) builds a 2D map in which points that are close in the original space stay close. The idea has three steps.

1. Neighbourhoods in the original space. For each point i, the probability that j is its neighbour is a Gaussian centred on i:

p(j|i) = exp(−‖x_i − x_j‖² / 2σ_i²) / Σ_{k≠i} exp(−‖x_i − x_k‖² / 2σ_i²)
p_ij = ( p(j|i) + p(i|j) ) / 2n

The width σ_i is chosen for each point, so that the point has an “effective” number of neighbours equal to the perplexity: Perp = 2^H(P_i), with H the Shannon entropy of the distribution p(·|i). In dense zones σ_i is small, in sparse ones it is large. It is found with a binary search.

2. Neighbourhoods in the map. For the points y_i of the map a Student t distribution with one degree of freedom is used, with heavy tails:

q_ij = (1 + ‖y_i − y_j‖²)⁻¹ / Σ_{k≠l} (1 + ‖y_k − y_l‖²)⁻¹

The heavy tails solve the “crowding problem”: in 2D there is no room to put all the moderate neighbours of a point close to it; with the t, moderately distant points can sit very far apart in the map without paying much.

3. Optimization. y is moved to minimize the Kullback-Leibler divergence between P and Q, KL(P‖Q) = Σ p_ij log(p_ij/q_ij), with gradient descent. The gradient has a compact form:

∂KL/∂y_i = 4 Σ_j (p_ij − q_ij) (y_i − y_j) (1 + ‖y_i − y_j‖²)⁻¹

Each pair acts like a spring: it attracts the points if p_ij > q_ij (they are close in the data but far in the map), and repels them in the opposite case.

Implementation by hand (500 digits)

The listing has two functions: affinita builds the matrix P with the binary search of σ_i for each point, and tsne optimizes the 2D positions with the gradient written above. At the end it compares our map with the scikit-learn one and with plain two-component PCA, using the trustworthiness (how much the neighbours in the map are also neighbours in the original data).

from sklearn.manifold import TSNE, trustworthiness

def affinita(D, perplexity):
    n = len(D); P = np.zeros((n, n)); target = np.log(perplexity)
    for i in range(n):
        d = np.delete(D[i], i); d = d - d.min()                 # shifting does not change the distribution
        lo, hi, beta = 0.0, np.inf, 1.0 / d.mean()              # beta = 1 / (2 sigma^2)
        for _ in range(60):                                     # binary search on beta
            pr = np.exp(-beta * d); s = pr.sum()
            H = np.log(s) + beta * (d * pr).sum() / s           # entropy (in nats)
            if H > target: lo, beta = beta, (beta * 2 if np.isinf(hi) else (beta + hi) / 2)
            else:          hi, beta = beta, (beta + lo) / 2
        P[i, np.arange(n) != i] = np.exp(-beta * d) / np.exp(-beta * d).sum()
    return np.maximum((P + P.T) / (2 * n), 1e-12)

def tsne(X, perplexity=30, n_iter=500, lr=200.0, seed=0):
    n = len(X)
    D = ((X[:, None, :] - X[None, :, :]) ** 2).sum(axis=2)
    P = affinita(D, perplexity)
    rng = np.random.default_rng(seed)
    Y = rng.normal(0, 1e-4, (n, 2)); V = np.zeros_like(Y)
    for it in range(n_iter):
        esagerazione = 12.0 if it < 100 else 1.0                # first steps: P enlarged so that clusters form
        num = 1.0 / (1.0 + ((Y[:, None, :] - Y[None, :, :]) ** 2).sum(axis=2))
        np.fill_diagonal(num, 0.0)
        Q = np.maximum(num / num.sum(), 1e-12)
        PQ = (esagerazione * P - Q) * num
        grad = 4 * (PQ.sum(axis=1)[:, None] * Y - PQ @ Y)
        V = (0.5 if it < 250 else 0.8) * V - lr * grad          # descent with momentum
        Y = Y + V; Y = Y - Y.mean(axis=0)
    return Y, (P * np.log(P / Q)).sum()

idx = np.random.default_rng(0).choice(len(X), 500, replace=False)
Xs, ys = X[idx], y[idx]

Y_mia, kl = tsne(Xs)
Y_sk = TSNE(perplexity=30, init="pca", random_state=0).fit(Xs)
Y_pca = PCA(n_components=2).fit_transform(Xs)

print("final KL by hand: %.3f    sklearn: %.3f" % (kl, Y_sk.kl_divergence_))
for nome, Y in (("PCA 2D", Y_pca), ("t-SNE by hand", Y_mia), ("t-SNE sklearn", Y_sk.embedding_)):
    print("%-14s trustworthiness (10 neighbours): %.3f" % (nome, trustworthiness(Xs, Y, n_neighbors=10)))

final KL by hand: 0.550    sklearn: 0.519
PCA 2D         trustworthiness (10 neighbours): 0.830
t-SNE by hand  trustworthiness (10 neighbours): 0.981
t-SNE sklearn  trustworthiness (10 neighbours): 0.983

In the listing affinita builds the matrix P: for each point i it runs the binary search on beta (= 1/2σ²) until the entropy of the neighbour distribution equals the logarithm of the perplexity; then it symmetrizes and normalizes. tsne initializes the 2D points with tiny noise, and for 500 iterations computes Q with the t distribution, the gradient in the compact form, and updates the positions with gradient descent with momentum (the momentum accumulates a fraction of the previous step to cross flat zones).

In the first 100 iterations P is multiplied by 12 (early exaggeration): enlarging the attractions helps the clusters form before they are refined. The final comparison uses 500 randomly drawn digits.

The trustworthiness measures how much the neighbours in the map are also neighbours in the original space (1 = perfect). The linear projection of 2-dimensional PCA obtains 0.830, the by-hand t-SNE 0.981 and the scikit-learn one 0.983: preserving neighbourhoods is exactly the purpose of t-SNE. Our implementation comes close to the scikit-learn one (final KL 0.550 against 0.519), which initializes the points with PCA (init="pca"), an option that the scikit-learn guide recommends to better preserve the global structure. For the details of this experiment: 500 digits instead of 1797, to keep the by-hand computation times acceptable.

How NOT to read a t-SNE

The map is beautiful and deceptive. The most common reading errors:

• Distances between clusters mean nothing: t-SNE preserves local neighbourhoods, not global distances. Two clusters far apart in the map are not necessarily “more different” than two close ones.
• Cluster sizes mean nothing: the density normalization (σ_i) expands dense clusters and compresses sparse ones.
• The shape depends on the hyperparameters: the same t-SNE with a different perplexity draws maps of different appearance (clusters that split or merge; look at them by plotting the graphs with coloured labels).

The three errors have the same origin: the t-SNE map preserves only local neighbourhoods, and everything else (distances, sizes, shapes) is a side effect of the optimization and of the hyperparameters. A t-SNE figure is read for who is with whom, never for how far away it is. Of the effect of the hyperparameters we measure what can be measured in a table:

for perp in (2, 5, 30, 100):
    t = TSNE(perplexity=perp, init="pca", random_state=0).fit(Xs)
    print("perplexity %3d  KL %.3f   trustworthiness %.3f" % (perp, t.kl_divergence_, trustworthiness(Xs, t.embedding_, n_neighbors=10)))

perplexity   2  KL 0.583   trustworthiness 0.970
perplexity   5  KL 0.597   trustworthiness 0.980
perplexity  30  KL 0.519   trustworthiness 0.983
perplexity 100  KL 0.405   trustworthiness 0.984

On the plane of local neighbourhoods the trustworthiness stays high for all the perplexities (from 0.970 to 0.984): the “local” quality is similar, but the global arrangement of the islands can change, and the table does not show that. Beware of the KL value: when the perplexity changes, the distribution P itself changes, so the divergences are not comparable across perplexities. A lower KL with perplexity 100 does not mean that the map is “better”. Another rule: the perplexity must be smaller than the number of points, and typical values are between 5 and 50. The scikit-learn guide defines the perplexity as 2 raised to the Shannon entropy of the conditional distribution, that is, the effective number of neighbours considered.

The listing repeats scikit-learn's t-SNE with four perplexities (on the same 500 samples) and prints KL and trustworthiness: it is the table just commented. Three limits of use, not of reading, remain:

• No transform function: t-SNE directly optimizes the positions of the given points; there is no way to project a new point. PCA has a matrix W and projects any point.
• Do not use it as a step before a clustering or a classifier: it is meant for looking, not as feature extraction.
• It varies with the seed: repeat with several seeds before trusting a structure.

These three limits, together with the previous ones, explain why t-SNE is an exploration tool and not a pipeline step: the scikit-learn guide says it in different terms, recalling that its main purpose is the visualization of high-dimensional data, that the algorithm is stochastic and that the global structure is not explicitly preserved.

PCA versus t-SNE: which to choose

• PCA: fast, deterministic (up to sign), linear, invertible (you can reconstruct), projects new data, the components are interpretable through the weights (loadings). It does not capture curved structures.
• t-SNE: captures nonlinear structures, excellent for exploring; slow (O(n²) in the exact version), stochastic, no projection of new points, global distances not interpretable.
• Common practice: PCA to 30-50 components to remove noise and speed up, then t-SNE on those to visualize.

The choice comes down to a question: do you need a transformation that you also apply to new data and that you can invert, or a map to look at? In the first case PCA, in the second t-SNE. If you want both, PCA reduces the noise first and t-SNE draws afterwards, but the t-SNE map does not go into production.

Real-world traps

• PCA without scaling with different units (see above).
• PCA fitted on the whole dataset before the train/test split: leakage. Inside the pipeline.
• Mistaking “variance” for “information useful for classification”: the direction of maximum variance can be irrelevant to the class you want to predict. PCA does not look at the labels.
• Interpreting the axes of t-SNE: the axes have no meaning.
• Sign of the components: two runs can give components with opposite sign. If you save a model, save the matrix W too.

The five traps are of two kinds. Those about PCA (scale, leakage, variance versus usefulness, sign) falsify the result without signalling it; the one about the axes of t-SNE is a reading error. For all of them the rule of the lesson holds: measure before interpreting, and compare the pipeline with the one without reduction.

Three ways to compute PCA: eigenvalues, SVD and scikit-learn

So far PCA has been presented as the diagonalization of the covariance matrix. In practice the covariance is almost never computed explicitly: the singular value decomposition (SVD) of the centred data matrix X_c = U Σ Vᵀ is used. The link is direct: the covariance is X_cᵀ X_c / (n − 1) = V Σ² Vᵀ / (n − 1), so the columns of V (the “right singular vectors”) are the eigenvectors of the covariance, and the eigenvalues are σ_i² / (n − 1). The SVD is preferred because it is numerically more stable (it does not square the condition numbers) and does not require building the d × d matrix, which with thousands of features would be enormous. Let us verify that the three paths give the same thing on the digits.

import warnings; warnings.filterwarnings("ignore")
from sklearn.model_selection import train_test_split, cross_val_score
from sklearn.neighbors import KNeighborsClassifier
from sklearn.pipeline import make_pipeline
from sklearn.linear_model import LogisticRegression

X, y = load_digits(return_X_y=True)
Xc = X - X.mean(axis=0)                                                  # centred data

cov = Xc.T @ Xc / (len(X) - 1)
eig_values, eig_vectors = np.linalg.eigh(cov)                            # path 1: eigenvalues of the covariance
order = np.argsort(eig_values)[::-1]
eig_values, eig_vectors = eig_values[order], eig_vectors[:, order]

U, S, Vt = np.linalg.svd(Xc, full_matrices=False)                        # path 2: SVD of the data
svd_values = S ** 2 / (len(X) - 1)

pca = PCA().fit(X)                                                       # path 3: scikit-learn

print("maximum difference between the first 20 variances, eigenvalues vs SVD: %.1e" % np.abs(eig_values[:20] - svd_values[:20]).max())
print("maximum difference, scikit-learn vs SVD: %.1e" % np.abs(pca.explained_variance_[:20] - svd_values[:20]).max())
sign = np.sign((eig_vectors[:, :10] * Vt[:10].T).sum(axis=0))             # eigenvectors are defined up to sign
print("maximum difference between the first 10 vectors (after aligning the sign): %.1e" % np.abs(eig_vectors[:, :10] * sign - Vt[:10].T).max())
print("rank of the centred matrix:", np.linalg.matrix_rank(Xc), "  zero-variance pixels:", int((X.std(axis=0) == 0).sum()))

maximum difference between the first 20 variances, eigenvalues vs SVD: 4.5e-13
maximum difference, scikit-learn vs SVD: 2.6e-13
maximum difference between the first 10 vectors (after aligning the sign): 1.2e-14
rank of the centred matrix: 61   zero-variance pixels: 3

The listing computes the same decomposition in three ways: the eigenvalues of the explicit covariance (eigh, which exploits the symmetry of the matrix), the SVD of the centred data, and scikit-learn's PCA class (which internally uses the SVD). The first 20 variances coincide up to 4.5·10⁻¹³ between eigenvalues and SVD and up to 2.6·10⁻¹³ between scikit-learn and SVD: floating-point rounding differences, not differences of substance. The eigenvectors coincide up to 1.2·10⁻¹⁴ after aligning the sign: each eigenvector v is as valid as −v, and different libraries may return it with opposite signs, so a test that compares the components without aligning the sign is a fragile test.

The last line is a useful reminder: the centred matrix has rank 61 and not 64, because 3 border pixels are always zero across the whole dataset (zero variance) and the last three components carry no information. If you did not know this fact, the whitening below would surprise you.

Removing noise with PCA and the whitening trap

PCA has a second use, besides dimension reduction for models: reconstruction. If you project an image onto the first k components and then go back with inverse_transform, you obtain the best rank-k approximation (in the least-squares sense). The noise, which is spread uniformly over all the directions, is largely discarded along with the low-variance components, while the signal, concentrated in the first components, remains. Let us measure it: we add Gaussian noise with deviation 4 (on pixels that range from 0 to 16) to the test images and reconstruct with k components computed on the training set alone.

X_train, X_test, y_train, y_test = train_test_split(X, y, test_size=0.3, random_state=0, stratify=y)
rng3 = np.random.default_rng(0)
noisy = X_test + rng3.normal(0, 4, X_test.shape)
rmse = lambda a, b: np.sqrt(((a - b) ** 2).mean())

print("error (RMSE) with respect to the clean images; the noise alone is worth %.2f" % (rmse(noisy, X_test)))
print("k    RMSE after reconstruction   RMSE if the image were already clean")
for k in (2, 5, 10, 20, 40, 64):
    pc = PCA(k).fit(X_train)
    denoised = pc.inverse_transform(pc.transform(noisy))
    clean_rec = pc.inverse_transform(pc.transform(X_test))
    print("%-4d  %-27.3f  %.3f" % (k, rmse(denoised, X_test), rmse(clean_rec, X_test)))

error (RMSE) with respect to the clean images; the noise alone is worth 4.00
k    RMSE after reconstruction   RMSE if the image were already clean
2     3.743                        3.676
5     3.160                        2.950
10    2.750                        2.247
20    2.668                        1.448
40    3.196                        0.485
64    3.996                        0.000

The noise alone is worth a root mean squared error of 4.00 (it is the deviation we chose). Then for each k it fits PCA on the training set, projects the noisy images, reconstructs them and measures the error with respect to the clean images. The error falls from 3.743 with k = 2 to a minimum of 2.668 with k = 20, then rises: 3.196 with k = 40 and 3.996 with k = 64, where the reconstruction is the identity and the noise comes back whole.

This is the bias-variance trade-off in the guise of compression: with a small k you throw away the signal too (the right-hand column, the error on already clean images, starts at 3.676 and falls to zero with 64 components), with a large k you keep the noise too. The optimum, k = 20, cuts the error from 4.00 to 2.67 (about a third less) using only 20 numbers per image. Note that the optimal k for noise removal is not the one that “keeps 95% of the variance”: it depends on the amount of noise, and it is another parameter to be chosen with validation data.

Whitening: useful and dangerous

Whitening divides each component by its standard deviation, so that all the components have variance 1 and are uncorrelated: the covariance of the transformed data becomes the identity matrix. Some models (those based on distance, or with isotropic regularization) benefit from it because no direction dominates. But there is a downside: dividing by a small deviation amplifies precisely the directions with the least variance, which are those dominated by noise, and if a deviation is practically zero the division explodes.

pw = PCA(20, whiten=True).fit(X_train)
cov_w = np.cov(pw.transform(X_train).T)
print("covariance after whitening: diagonal between %.3f and %.3f, largest off-diagonal element %.1e" % (cov_w.diagonal().min(), cov_w.diagonal().max(), np.abs(cov_w - np.diag(cov_w.diagonal())).max()))

print("\ntest accuracy, with and without whitening:")
print("k    KNN without   KNN with   logistic regr. without   logistic regr. with")
for k in (10, 20, 40, 64):
    row = []
    for model in (KNeighborsClassifier(5), LogisticRegression(max_iter=10000)):
        for whiten in (False, True):
            row.append(make_pipeline(PCA(k, whiten=whiten), model).fit(X_train, y_train).score(X_test, y_test))
    print("%-4d  %-12.3f  %-9.3f  %-24.3f  %.3f" % (k, row[0], row[1], row[2], row[3]))

covariance after whitening: diagonal between 1.000 and 1.000, largest off-diagonal element 1.6e-15

test accuracy, with and without whitening:
k    KNN without   KNN with   logistic regr. without   logistic regr. with
10    0.969         0.956      0.941                     0.939
20    0.974         0.972      0.948                     0.948
40    0.980         0.957      0.961                     0.965
64    0.980         0.237      0.967                     0.961

The first line confirms the definition: after whitening 20 components, the covariance of the transformed data has a diagonal of exactly 1.000 and the largest off-diagonal element is 1.6·10⁻¹⁵ (zero up to rounding). Then it compares four configurations with 10, 20, 40 and 64 components. With few components whitening changes almost nothing: with 20 components KNN goes from 0.974 to 0.972 and logistic regression stays at 0.948. With 40, KNN gets worse (0.980 → 0.957) while logistic regression improves slightly (0.961 → 0.965): small effects of different sign, which on a test of 540 images cannot be told apart from chance.

With 64 components the disaster happens: KNN with whitening collapses to 0.237, while without whitening it stays at 0.980. The cause is the announced one: the last components have zero or almost zero variance (remember the 3 always-zero pixels) and dividing them by their very small deviation turns them into enormous numbers, which dominate the Euclidean distance. No error, no warning: the code runs and returns a terrible model. The practical rule is that whitening should be used only together with a choice of k that excludes the negligible components, or with a regularization term in the division; and it should always be compared with the same pipeline without whitening.

How many components really: explained variance versus accuracy

The number of components to keep is usually chosen with an explained-variance threshold. Let us compare that threshold with a classifier's accuracy, to see whether the two measures say the same thing.

print("k    explained variance   CV accuracy (logistic regression, 5 folds)")
for k in (2, 5, 10, 20, 30, 64):
    scores = cross_val_score(make_pipeline(PCA(k), LogisticRegression(max_iter=10000)), X, y, cv=5)
    print("%-4d  %-18.3f  %.3f +- %.3f" % (k, pca.explained_variance_ratio_[:k].sum(), scores.mean(), scores.std()))

k    explained variance   CV accuracy (logistic regression, 5 folds)
2     0.285               0.574 +- 0.015
5     0.545               0.823 +- 0.040
10    0.738               0.888 +- 0.033
20    0.894               0.895 +- 0.029
30    0.959               0.910 +- 0.030
64    1.000               0.913 +- 0.024

For each k the listing puts PCA inside a pipeline with logistic regression and measures the accuracy with 5-fold cross-validation, which redoes the PCA on each training fold (so there is no information leak from the validation fold). It compares the cumulative explained variance (how much geometric information you keep) with the accuracy (how much the classifier really needs). With k = 2, 28.5% of the variance is kept and the accuracy is 0.574; with k = 5, 54.5% and 0.823; with k = 10, 73.8% and 0.888; with k = 20, 89.4% and 0.895; with k = 30, 95.9% and 0.910; with all 64, 100% and 0.913.

The accuracy flattens out much earlier than the variance: between 10 and 64 components it gains 2.5 points, while the standard deviation across folds is 0.024–0.033, so the differences between k = 20, 30 and 64 are of the same order as the measurement noise. The consequence is an operating rule: instead of fixing “95% of the variance” in advance, choose the smallest k whose accuracy is not distinguishable from the maximum within the variability across folds. Here it would be a value between 20 and 30, which halves or more the size of the input. Keep in mind that this result concerns the 8×8 digits, a very easy problem: on different data the curve will have another shape.

Three exercises to try

Exercise 1. Repeat the noise-removal experiment with noise of deviation 2 and of deviation 8: does the optimal k change? Why do you expect it to shift in that direction? Exercise 2. Whiten only the components with explained variance above a threshold and leave the others unchanged: do you recover the advantages of whitening without the collapse of KNN? Exercise 3. Choose k as the smallest value for which the CV accuracy stays within one standard deviation of the maximum, and compare it with the 95%-explained-variance rule: do they coincide?

To practice on the anthropometric measurements for helmet shells there is a Python exercise skeleton with three parts to complete: eigenvalues and eigenvectors of a 2x2 matrix in closed form, explained variance and reconstruction error, and the t-SNE conditional probabilities that do not change when all the distances are tripled. As long as the points to complete are empty, it prints which exercises are missing.

What you can do now

You can derive PCA as an eigenvalue problem and verify that maximum variance and minimum reconstruction error coincide. You can choose how many components to keep by looking at accuracy and not just variance, use PCA to remove noise and recognize the collapse caused by whitening. You can write t-SNE by hand and read it without attributing meaning to the distances and sizes of the groups. In the next lesson, linear regression, we start again from a straight line.

Sources

Previous lesson in the course: Clustering: K-Means and DBSCAN without labels. Other lessons in the series that link to this one: Linear Algebra 2: Eigenvalues, SVD and PCA, KNN: classifying by looking at neighbours (and the curse of dimensionality).