Dimensionality Reduction

Dimensionality Reduction

Three methods, three different objectives. PCA preserves global variance and is reversible. t-SNE preserves local neighbourhoods and distorts everything else. UMAP attempts both and fails differently. Knowing which distortion each one introduces is the whole content of this page.¹

1. Why the statistics push you to reduce

Beyond the four goals above there is a structural reason: in high dimensions, geometry stops behaving the way two- and three-dimensional intuition expects. Volume runs to the surface, the inscribed ball all but vanishes, random directions become almost orthogonal, and — the one that actually breaks things — pairwise distances concentrate, so "nearest neighbour" gradually stops meaning anything.

That last effect is why t-SNE and UMAP both recommend running PCA first when the input has hundreds or thousands of columns. It is not only speed: it is the quality of the neighbourhoods those methods are built on.

2. PCA: two derivations, one answer

Centre the data so \(\mu = 0\), and look for a unit direction \(\mathbf{w}\) that summarizes it well. There are two ways to say "well", and the elegance of PCA is that they coincide.

Try to find the best angle before pressing snap to PC1. The two bars always sum to the same number, at every angle: that is the identity above, applied to the whole cloud.

The second component is the eigenvector of the second largest eigenvalue, necessarily orthogonal to the first because \(\Sigma\) is symmetric. Stacking the first \(k\) into \(W\) gives the projection \(Z = XW\) and the reconstruction \(\hat X = ZW^\top\).

Choosing \(k\)

"""PCA on Iris: what each component costs, and what it is worth.

Four standardized measurements, four components. `explained` is the share of
total variance each one carries, `recon_mse` is what is lost by keeping only
the first k and projecting back, `trustworthiness` is how many of each point's
twelve nearest neighbours in four dimensions are still its neighbours in k, and
`knn_acc` is a 5-NN classifier trained on the k components.

The last column is the one to look at twice. Two components carry 95.8% of the
variance and reconstruct almost perfectly — and they classify *worse* than
three. PCA is unsupervised: it maximizes variance, and nothing guarantees that
the directions with the most variance are the directions that separate the
classes.

Printed as a markdown table, in identifiers only, so one artifact serves both
the English and the Portuguese page.
"""

import numpy as np
from sklearn.datasets import load_iris
from sklearn.decomposition import PCA
from sklearn.manifold import trustworthiness
from sklearn.model_selection import StratifiedKFold, cross_val_score
from sklearn.neighbors import KNeighborsClassifier
from sklearn.pipeline import make_pipeline
from sklearn.preprocessing import StandardScaler

X, y = load_iris(return_X_y=True)
X_std = StandardScaler().fit_transform(X)        # PCA maximizes variance, and variance has units
cv = StratifiedKFold(5, shuffle=True, random_state=0)

ratios = PCA().fit(X_std).explained_variance_ratio_

print("| `k` | `explained` | `cumulative` | `recon_mse` | `trustworthiness` | `knn_acc` |")
print("|---:|---:|---:|---:|---:|---:|")
for k in (1, 2, 3, 4):
    pca = PCA(k).fit(X_std)
    scores = pca.transform(X_std)
    recon = np.mean((X_std - pca.inverse_transform(scores)) ** 2)
    # the accuracy is cross-validated with PCA *inside* the pipeline, or it would leak
    acc = cross_val_score(
        make_pipeline(StandardScaler(), PCA(k), KNeighborsClassifier(5)), X, y, cv=cv
    ).mean()
    print(f"| {k} | {ratios[k - 1]:.4f} | {ratios[:k].sum():.4f} | {recon:.4f} "
          f"| {trustworthiness(X_std, scores, n_neighbors=12):.3f} | **{acc:.3f}** |")

Iris, standardized, four components. The first two carry 95.8% of the variance and reconstruct with an error of 0.042 against 0.270 for one. By the usual criteria — "keep 95% of the variance", "look for the elbow in the scree plot" — the answer is \(k = 2\) and you stop.

Now read the last column. \(k = 2\) classifies at 0.900 and \(k = 3\) at 0.960. The third component carries 3.7% of the variance and it is worth six points of accuracy.

from sklearn.decomposition import PCA
from sklearn.pipeline import make_pipeline
from sklearn.preprocessing import StandardScaler

pipe = make_pipeline(StandardScaler(), PCA(n_components=0.95), MLPClassifier())
# n_components as a float = "keep enough components for this share of the variance"
# and, inside the pipeline, it is re-fitted on every training fold

Doing it by hand, once

The eigendecomposition is worth writing out once, so that PCA(n_components=2) stops being a black box. Both scripts produce the same projection, up to the arbitrary sign of each component.

import matplotlib.pyplot as plt
import numpy as np
import pandas as pd
from io import StringIO
from sklearn.datasets import load_iris
from sklearn.preprocessing import StandardScaler

# Loading Iris dataset
iris = load_iris()

# Transform in dataframe
df = pd.DataFrame(
    data=iris.data,
    columns=['sepal_l', 'sepal_w', 'petal_l', 'petal_w']
)
df['class'] = iris.target_names[iris.target]

X = df.iloc[:,0:4].values
y = df.iloc[:,4].values

# Standardizing
X_std = StandardScaler().fit_transform(X)

# Covariance
cov_mat = np.cov(X_std.T)

# Calculate autovalues and autovectors
eig_vals, eig_vecs = np.linalg.eig(cov_mat)

print('Eigenvectors \n%s' %eig_vecs)
print('\nEigenvalues \n%s' %eig_vals)

# Make a list of (eigenvalue, eigenvector) tuples
eig_pairs = [(np.abs(eig_vals[i]), eig_vecs[:,i]) for i in range(len(eig_vals))]

# Sort the (eigenvalue, eigenvector) tuples from high to low
eig_pairs.sort(key=lambda x: x[0], reverse=True)

# Visually confirm that the list is correctly sorted by decreasing eigenvalues
print('Eigenvalues in descending order:')
for i in eig_pairs: print(i[0])

# Sum the cummulative of each eigen value
tot = sum(eig_vals)
var_exp = [(i / tot)*100 for i in sorted(eig_vals, reverse=True)]
cum_var_exp = np.cumsum(var_exp)

n_eigen = [1, 2, 3, 4]

# Plot the cumulative for each eign value
plt.figure(figsize=(6, 4))
plt.bar(n_eigen, var_exp, alpha=0.5, align='center',
    label='individual explained variance')
plt.step(n_eigen, cum_var_exp, where='mid',
    label='cumulative explained variance')
plt.ylabel('Explained variance ratio')
plt.xlabel('Principal components')
plt.legend(loc='best')
plt.tight_layout()

# Take the only the two firsts eigen values
matrix_w = np.hstack((eig_pairs[0][1].reshape(4,1),
                      eig_pairs[1][1].reshape(4,1)))

print('*' * 10)
print('Reduced to 2-D')
print('Matrix W:\n', matrix_w)

# Calculate the new Y for all samples
Y = X_std.dot(matrix_w)

# Plot the data for the 2 firsts principal components
plt.figure(figsize=(6, 4))
for lab, col in zip(('setosa', 'versicolor', 'virginica'), ('blue', 'red', 'green')):
    plt.scatter(Y[y==lab, 0],
                Y[y==lab, 1],
                label=lab,
                c=col)
plt.xlabel('Principal Component 1')
plt.ylabel('Principal Component 2')
plt.legend(loc='lower center')
plt.tight_layout()

# Para imprimir na página HTML
buffer = StringIO()
plt.savefig(buffer, format="svg", transparent=True)
print(buffer.getvalue())
plt.close()

import matplotlib.pyplot as plt
import pandas as pd
from io import StringIO
from sklearn.datasets import load_iris
from sklearn.decomposition import PCA as pca
from sklearn.preprocessing import StandardScaler

# Loading Iris dataset
iris = load_iris()

# Transform in dataframe
df = pd.DataFrame(
    data=iris.data,
    columns=['sepal_l', 'sepal_w', 'petal_l', 'petal_w']
)
df['class'] = iris.target_names[iris.target]

X = df.iloc[:,0:4].values
y = df.iloc[:,4].values

# Standardizing
X_std = StandardScaler().fit_transform(X)

sklearn_pca = pca(n_components=2)
Y = sklearn_pca.fit_transform(X_std)

# Plot the data for the 2 firsts principal components
plt.figure(figsize=(6, 4))
for lab, col in zip(('setosa', 'versicolor', 'virginica'), ('blue', 'red', 'green')):
    plt.scatter(Y[y==lab, 0],
                Y[y==lab, 1],
                label=lab,
                c=col)
plt.xlabel('Principal Component 1')
plt.ylabel('Principal Component 2')
plt.legend(loc='lower center')
plt.tight_layout()

# Para imprimir na página HTML
buffer = StringIO()
plt.savefig(buffer, format="svg", transparent=True)
print(buffer.getvalue())
plt.close()

3. Where PCA fails

PCA finds a linear subspace. When the data lies on a curved manifold — a spiral, a swiss roll, the set of all images of a rotating object — no plane summarizes it, and the first components will happily describe a direction that cuts straight through the curve.

The number under each projection counts false neighbours: of the ten nearest points in the map, how many are more than a tenth of the curve away along it. PCA scores 72% — three out of four of each point's apparent neighbours are somewhere else entirely. The starred pair is the extreme case: four turns apart on the spring, 0.007 apart after projection.

Two things are worth noticing before moving on.

First, this failure is invisible to trustworthiness, the metric section 7 recommends. It scores PCA at 0.956 here, because it compares the map against distances in the ambient 3-D space — and in 3-D, points on adjacent coils genuinely are close. What the projection destroys is the intrinsic structure, the distance along the curve, and you only see that by measuring against it. Which metric you pick decides what counts as a failure.

Second, look at what t-SNE did. It gets the neighbourhoods almost perfectly — 1% false — and it tears the curve into arcs and scatters them. That is not a defect to fix; it is the trade. The neighbourhoods are right and the arrangement of the pieces means nothing, which is exactly why the seven forbidden statements exist.

That is the opening for the non-linear methods, and they change the question. PCA asks "which directions concentrate variance?" — a global question with a linear, reversible answer. t-SNE and UMAP ask "who is near whom?" — a local question, answered by deliberately sacrificing the global geometry.

4. t-SNE: matching neighbourhood distributions²

In one sentence: turn distances into neighbourhood probabilities in both spaces, and move the points of the 2-D map until the two distributions agree.

Step 1 — neighbourhoods in the original space. For each point \(i\), a gaussian kernel over its neighbours:

Each point gets its own \(\sigma_i\), found by binary search until the perplexity \(2^{H(P_i)}\) hits the value you asked for. Perplexity ≈ the effective number of neighbours. Dense regions get a small \(\sigma\), sparse regions a large one — an automatic adaptation of scale. The result is symmetrized, \(p_{ij} = (p_{j|i} + p_{i|j})/2n\), so that every point contributes to the loss, including isolated ones.

Step 2 — neighbourhoods in the map, with a heavy tail.

A Student-\(t\) with one degree of freedom (a Cauchy), and this is the essential difference from the original SNE. Why the heavy tail — the crowding problem: the volume of a ball grows as \(r^d\), so a point in high dimensions has room for many neighbours at moderate distance, and in 2-D that room does not exist. With a gaussian on both sides, every moderate neighbour would be forced toward the centre and the map would collapse. The heavy tail lets moderately distant pairs sit very far apart in the map at low cost, which frees space to separate the groups.

The side effect is the thing you must remember: distances between groups in the picture stop being interpretable.

Step 3 — minimize the divergence.

The KL is asymmetric, and the asymmetry is the whole behaviour of the method. High \(p\) with low \(q\) — real neighbours drawn apart — is punished heavily. Low \(p\) with high \(q\) — distant points drawn together — is nearly free. Hence the rule:

What is together in a t-SNE map may not be together in reality.

Watch it happen. The panel runs the real algorithm — the perplexity search, the early exaggeration, the gradient descent — a few iterations per frame.

The dataset to spend the most time on is the last one. It is twelve columns of uniform noise, with nothing to find, and t-SNE will hand you tidy separated groups anyway.

5. UMAP: a fuzzy graph, drawn

UMAP reaches a similar result by a different route. Instead of distributions over all pairs, it builds an explicit \(k\)-neighbour graph with fuzzy weights and then draws it.

Step 1 — a local graph with guaranteed connectivity.

\(\rho_i\) is the distance to the nearest neighbour and \(\sigma_i\) is calibrated so that \(\sum_j w_{i \to j} = \log_2 k\). Subtracting \(\rho_i\) guarantees every point has at least one edge of weight 1, so nothing is left isolated. The two directed views are combined by a fuzzy union — the analogue of t-SNE's symmetrization.

Step 2 — draw the graph. The loss is a binary cross-entropy rather than a KL:

The second term is what t-SNE does not have: an explicit penalty for placing far things near each other. That is the origin of UMAP's slightly better global preservation. Optimization uses SGD with negative sampling — one edge is drawn and its endpoints attract, a few random points are drawn as negatives and repel — which avoids summing over all pairs and is why UMAP is faster.

# not a course dependency: pip install umap-learn
import umap
Z = umap.UMAP(n_neighbors=15, min_dist=0.1, random_state=0).fit_transform(X_std)

n_neighbors trades local detail against global structure; min_dist controls how tightly points may pack in the map, and it is purely aesthetic — it changes the picture without changing what was learned.

The panel below builds the graph and then draws it, with the same four datasets. Turn on graph to see the thing being optimized: the layout is that graph, untangled.

Run the coiled curve in both panels and the difference in the losses shows up directly — UMAP's second term, the explicit penalty for placing far things near each other, keeps the curve connected noticeably longer.

6. Choosing, and what you may say afterwards

flowchart TD
    A["what is the reduction for?"] --> B{"will a model consume<br/>the output?"}
    B -->|yes| P["<b>PCA</b><br/><small>honest fit/transform, reversible</small>"]
    B -->|"no — it is a picture"| C{"how many rows?"}
    C -->|"up to ~10k"| T["<b>t-SNE</b><br/><small>best local fidelity</small>"]
    C -->|"more, or new rows arrive"| U["<b>UMAP</b><br/><small>faster, approximate transform</small>"]
    T --> V["and always: PCA to ~50 first"]
    U --> V

    classDef ok fill:#e6f4ea,stroke:#3fb950,color:#14532d
    classDef q  fill:#eef2f7,stroke:#8b949e,color:#1f2937
    class P,T,U,V ok
    class A,B,C q

The last box is not optional. With hundreds or thousands of input columns, run PCA down to ~50 before t-SNE or UMAP — not only for speed, but because the neighbourhoods those methods consume are computed with the very distances that A3 above is about.

7. How to judge a reduction

A picture is not a result. If a reduction is going to support a claim, it needs a number. The three that are used:

• Reconstruction error — only meaningful for reversible methods (PCA, an autoencoder). It measures what was lost.
• Trustworthiness — of each point's \(k\) nearest neighbours in the map, how many were genuinely among its nearest neighbours in the original space. Applies to any method, including t-SNE and UMAP. It is in sklearn.manifold.trustworthiness.
• Downstream performance — the accuracy of a kNN trained in the reduced space against the original. The most honest test when the reduction is preprocessing.

"""What t-SNE is good at, and the thing it will do to you if you let it.

Two measurements, each answering one half of the question.

First: does a two-dimensional map keep the neighbourhoods of the original
space? `trustworthiness` counts, for every point, how many of its twelve
nearest neighbours in the map were genuinely among its nearest neighbours in
64 dimensions. This is what t-SNE optimizes and PCA does not, and the gap is
the reason t-SNE exists.

Second: what does t-SNE do to data with no structure at all? The input is
uniform noise in twenty dimensions — no clusters, by construction. A 5-means
silhouette near zero is the correct answer, and the original space gives it.
The t-SNE map does not: the optimization has to put the points somewhere, and
where it puts them looks like groups.

Printed as a markdown table, in identifiers only, so one artifact serves both
the English and the Portuguese page.
"""

import numpy as np
from sklearn.cluster import KMeans
from sklearn.datasets import load_digits
from sklearn.decomposition import PCA
from sklearn.manifold import TSNE, trustworthiness
from sklearn.metrics import silhouette_score
from sklearn.preprocessing import StandardScaler

NEIGHBOURS, PERPLEXITY = 12, 30

digits = StandardScaler().fit_transform(load_digits(return_X_y=True)[0])
embed = lambda X: TSNE(n_components=2, init="pca", perplexity=PERPLEXITY,
                       random_state=0).fit_transform(X)

print(f"| `digits_1797x64 -> 2d` | `trustworthiness` |")
print("|---|---:|")
for name, Z in (("PCA", PCA(2, random_state=0).fit_transform(digits)),
                ("t-SNE", embed(digits))):
    print(f"| `{name}` | **{trustworthiness(digits, Z, n_neighbors=NEIGHBOURS):.3f}** |")

print()
noise = np.random.default_rng(0).random((600, 20))       # no structure whatsoever
print("| `uniform_noise_600x20` | `silhouette_of_5_means` |")
print("|---|---:|")
for name, Z in (("original_20d", noise), ("t-SNE_map_2d", embed(noise))):
    labels = KMeans(5, n_init=10, random_state=0).fit_predict(Z)
    print(f"| `{name}` | **{silhouette_score(Z, labels):.3f}** |")