Dimensionality Reduction

How many components should you keep? Explained variance ratios and scree plots help you decide the optimal number of dimensions to retain while preserving most of the information in your data.

Understanding Explained Variance

Each principal component captures a portion of the total variance (spread) in your data. The first component always captures the most (it's designed to!), the second captures the next most from what remains, and so on. The explained variance ratio tells you what fraction of total variance each component represents - like slices of a pie.

Why does this matter? If PC1 explains 90% of variance and PC2 explains 8%, keeping just these 2 components means you've retained 98% of your data's information! The remaining components (PC3, PC4, ...) might just be noise and can be safely discarded.

Computing Explained Variance

Fit PCA on all components first to see how variance is distributed:

import numpy as np
import matplotlib.pyplot as plt
from sklearn.decomposition import PCA
from sklearn.preprocessing import StandardScaler
from sklearn.datasets import load_iris

# Load and scale data
iris = load_iris()
scaler = StandardScaler()
X_scaled = scaler.fit_transform(iris.data)

# Fit PCA on ALL components
pca = PCA()  # No n_components = keep all
pca.fit(X_scaled)

# Get explained variance
variance_ratio = pca.explained_variance_ratio_
cumulative = np.cumsum(variance_ratio)

print("Individual variance:", variance_ratio.round(3))
print("Cumulative variance:", cumulative.round(3))

Creating a Scree Plot

A scree plot visualizes the explained variance for each component. Look for an "elbow" where adding more components yields diminishing returns:

# Create a scree plot
fig, ax1 = plt.subplots(figsize=(10, 6))

# Bar chart for individual variance
components = range(1, len(variance_ratio) + 1)
ax1.bar(components, variance_ratio, alpha=0.7, color='steelblue', label='Individual')
ax1.set_xlabel('Principal Component')
ax1.set_ylabel('Explained Variance Ratio', color='steelblue')
ax1.tick_params(axis='y', labelcolor='steelblue')

# Line for cumulative variance
ax2 = ax1.twinx()
ax2.plot(components, cumulative, 'ro-', label='Cumulative')
ax2.axhline(y=0.95, color='g', linestyle='--', label='95% threshold')
ax2.set_ylabel('Cumulative Explained Variance', color='red')
ax2.tick_params(axis='y', labelcolor='red')

plt.title('Scree Plot with Cumulative Variance')
fig.legend(loc='center right', bbox_to_anchor=(0.85, 0.5))
plt.tight_layout()
plt.show()

Automatic Component Selection

Instead of manually choosing, you can specify a variance threshold:

# Method 1: Specify variance to retain
pca_95 = PCA(n_components=0.95)  # Keep 95% of variance
X_reduced = pca_95.fit_transform(X_scaled)
print(f"Components for 95%: {pca_95.n_components_}")

# Method 2: Specify exact number
pca_2 = PCA(n_components=2)
X_2d = pca_2.fit_transform(X_scaled)
print(f"Variance with 2 components: {sum(pca_2.explained_variance_ratio_):.2%}")

Variance Thresholds Table

Common rules of thumb for choosing the number of components:

Reconstructing Data from Components

You can transform data back to the original space. The difference between original and reconstructed data shows what information was lost:

# Reduce to 2 components
pca = PCA(n_components=2)
X_reduced = pca.fit_transform(X_scaled)

# Reconstruct back to original dimensions
X_reconstructed = pca.inverse_transform(X_reduced)

# Calculate reconstruction error
mse = np.mean((X_scaled - X_reconstructed) ** 2)
print(f"Reconstruction MSE: {mse:.4f}")

# Compare original vs reconstructed for first sample
print("Original:", X_scaled[0].round(2))
print("Reconstructed:", X_reconstructed[0].round(2))

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

iris = load_iris()
X_scaled = StandardScaler().fit_transform(iris.data)

pca = PCA()  # Keep all components
pca.fit(X_scaled)

print("Variance ratio:", pca.explained_variance_ratio_.round(4))

import numpy as np
from sklearn.decomposition import PCA
from sklearn.preprocessing import StandardScaler
from sklearn.datasets import load_iris

iris = load_iris()
X_scaled = StandardScaler().fit_transform(iris.data)

pca = PCA()
pca.fit(X_scaled)

cumulative = np.cumsum(pca.explained_variance_ratio_)
print("Cumulative variance:", cumulative.round(4))

import numpy as np
import matplotlib.pyplot as plt
from sklearn.decomposition import PCA
from sklearn.preprocessing import StandardScaler
from sklearn.datasets import load_iris

iris = load_iris()
X_scaled = StandardScaler().fit_transform(iris.data)
pca = PCA().fit(X_scaled)

variance = pca.explained_variance_ratio_
cumulative = np.cumsum(variance)
components = range(1, len(variance) + 1)

fig, ax1 = plt.subplots(figsize=(8, 5))
ax1.bar(components, variance, alpha=0.7, color='blue')
ax1.set_xlabel('Component')
ax1.set_ylabel('Individual Variance')

ax2 = ax1.twinx()
ax2.plot(components, cumulative, 'r-o')
ax2.set_ylabel('Cumulative Variance')

plt.title('Scree Plot')
plt.show()

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

iris = load_iris()
X_scaled = StandardScaler().fit_transform(iris.data)

pca = PCA(n_components=0.90)
pca.fit(X_scaled)

print(f"Components for 90%: {pca.n_components_}")
print(f"Actual variance: {sum(pca.explained_variance_ratio_):.2%}")

t-SNE (t-distributed Stochastic Neighbor Embedding) excels at revealing clusters and patterns in high-dimensional data by creating stunning 2D or 3D visualizations that preserve local structure.

What is t-SNE?

Technical explanation: Unlike PCA which preserves global variance (overall spread of data), t-SNE focuses on preserving local neighborhoods. Points that are close together in high-dimensional space will remain close in the low-dimensional projection. This makes t-SNE excellent for discovering clusters that might not be visible with PCA.

Important caveat: Because t-SNE distorts global structure, the distances between clusters are NOT meaningful. Two clusters might appear far apart in a t-SNE plot but could actually be close in the original space. t-SNE is for visualization and cluster discovery ONLY, not for measuring distances!

Applying t-SNE in Python

Scikit-learn provides an easy-to-use t-SNE implementation:

from sklearn.manifold import TSNE
from sklearn.preprocessing import StandardScaler
from sklearn.datasets import load_iris
import matplotlib.pyplot as plt

# Load and scale data
iris = load_iris()
X = iris.data
y = iris.target

scaler = StandardScaler()
X_scaled = scaler.fit_transform(X)

# Apply t-SNE
tsne = TSNE(n_components=2, random_state=42, perplexity=30)
X_tsne = tsne.fit_transform(X_scaled)

print(f"Original shape: {X.shape}")  # (150, 4)
print(f"t-SNE shape: {X_tsne.shape}")  # (150, 2)

# Visualize t-SNE results
plt.figure(figsize=(10, 8))
scatter = plt.scatter(X_tsne[:, 0], X_tsne[:, 1], c=y, cmap='viridis', alpha=0.7, s=50)
plt.xlabel('t-SNE Dimension 1')
plt.ylabel('t-SNE Dimension 2')
plt.title('t-SNE Visualization of Iris Dataset')
plt.colorbar(scatter, label='Species')
plt.show()

The Perplexity Parameter

Perplexity is the most important hyperparameter. It controls how many neighbors each point considers:

# Compare different perplexity values
fig, axes = plt.subplots(1, 3, figsize=(15, 5))
perplexities = [5, 30, 50]

for ax, perp in zip(axes, perplexities):
    tsne = TSNE(n_components=2, perplexity=perp, random_state=42)
    X_tsne = tsne.fit_transform(X_scaled)
    ax.scatter(X_tsne[:, 0], X_tsne[:, 1], c=y, cmap='viridis', alpha=0.7)
    ax.set_title(f'Perplexity = {perp}')
    ax.set_xlabel('Dim 1')
    ax.set_ylabel('Dim 2')

plt.tight_layout()
plt.show()

t-SNE Best Practices

Speed Optimization with PCA

For datasets with many features, first reduce with PCA before applying t-SNE:

from sklearn.decomposition import PCA
from sklearn.manifold import TSNE
from sklearn.datasets import fetch_openml

# Load MNIST (784 features)
mnist = fetch_openml('mnist_784', version=1, as_frame=False)
X, y = mnist.data[:5000], mnist.target[:5000].astype(int)

# Step 1: Reduce to 50 dimensions with PCA
pca = PCA(n_components=50)
X_pca = pca.fit_transform(X)

# Step 2: Apply t-SNE on reduced data (much faster!)
tsne = TSNE(n_components=2, random_state=42)
X_tsne = tsne.fit_transform(X_pca)

print(f"784 dims -> 50 dims -> 2 dims")

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

iris = load_iris()
X_scaled = StandardScaler().fit_transform(iris.data)

tsne = TSNE(n_components=2, random_state=42)
X_tsne = tsne.fit_transform(X_scaled)

print(f"Shape: {X_tsne.shape}")

import matplotlib.pyplot as plt
from sklearn.manifold import TSNE
from sklearn.preprocessing import StandardScaler
from sklearn.datasets import load_iris

iris = load_iris()
X_scaled = StandardScaler().fit_transform(iris.data)
tsne = TSNE(n_components=2, random_state=42)
X_tsne = tsne.fit_transform(X_scaled)

plt.figure(figsize=(8, 6))
plt.scatter(X_tsne[:, 0], X_tsne[:, 1], c=iris.target, cmap='viridis')
plt.xlabel('t-SNE 1')
plt.ylabel('t-SNE 2')
plt.colorbar(label='Species')
plt.title('t-SNE of Iris')
plt.show()

import matplotlib.pyplot as plt
from sklearn.manifold import TSNE
from sklearn.preprocessing import StandardScaler
from sklearn.datasets import load_iris

iris = load_iris()
X_scaled = StandardScaler().fit_transform(iris.data)

fig, axes = plt.subplots(1, 3, figsize=(15, 5))
for ax, perp in zip(axes, [5, 30, 50]):
    tsne = TSNE(n_components=2, perplexity=perp, random_state=42)
    X_tsne = tsne.fit_transform(X_scaled)
    ax.scatter(X_tsne[:, 0], X_tsne[:, 1], c=iris.target, cmap='viridis')
    ax.set_title(f'Perplexity = {perp}')
plt.tight_layout()
plt.show()

from sklearn.decomposition import PCA
from sklearn.manifold import TSNE
from sklearn.datasets import load_digits

digits = load_digits()
X = digits.data

# Step 1: PCA to 30 dimensions
pca = PCA(n_components=30)
X_pca = pca.fit_transform(X)

# Step 2: t-SNE on reduced data
tsne = TSNE(n_components=2, random_state=42)
X_tsne = tsne.fit_transform(X_pca)

print(f"Original: {X.shape} -> PCA: {X_pca.shape} -> t-SNE: {X_tsne.shape}")

UMAP (Uniform Manifold Approximation and Projection) is a modern technique that preserves both local and global structure, runs faster than t-SNE, and works well for both visualization and as a preprocessing step for machine learning.

Why UMAP? (The Modern Champion)

UMAP (Uniform Manifold Approximation and Projection) was developed in 2018 and has quickly become the preferred choice for many data scientists and machine learning practitioners. Why? Because it combines the best aspects of PCA (speed, usability as features for ML) with the visualization quality of t-SNE.

The UMAP advantage: Unlike t-SNE which only preserves local structure, UMAP preserves both local AND global structure. This means cluster separations and relative positions are more meaningful. Plus, it scales to millions of data points efficiently!

Installing UMAP

UMAP is not included in scikit-learn. Install it separately:

# Install UMAP (run once)
# pip install umap-learn

# Note: the package is 'umap-learn' but you import 'umap'
import umap

print("UMAP installed successfully!")

Applying UMAP in Python

import umap
from sklearn.preprocessing import StandardScaler
from sklearn.datasets import load_iris
import matplotlib.pyplot as plt

# Load and scale data
iris = load_iris()
X = iris.data
y = iris.target

scaler = StandardScaler()
X_scaled = scaler.fit_transform(X)

# Apply UMAP
reducer = umap.UMAP(n_components=2, random_state=42)
X_umap = reducer.fit_transform(X_scaled)

print(f"Original shape: {X.shape}")  # (150, 4)
print(f"UMAP shape: {X_umap.shape}")  # (150, 2)

# Visualize UMAP results
plt.figure(figsize=(10, 8))
scatter = plt.scatter(X_umap[:, 0], X_umap[:, 1], c=y, cmap='viridis', alpha=0.7, s=50)
plt.xlabel('UMAP Dimension 1')
plt.ylabel('UMAP Dimension 2')
plt.title('UMAP Visualization of Iris Dataset')
plt.colorbar(scatter, label='Species')
plt.show()

Key UMAP Parameters

UMAP has two main hyperparameters that control the embedding:

# Experiment with parameters
fig, axes = plt.subplots(2, 2, figsize=(12, 12))

params = [
    {'n_neighbors': 5, 'min_dist': 0.1},
    {'n_neighbors': 50, 'min_dist': 0.1},
    {'n_neighbors': 15, 'min_dist': 0.01},
    {'n_neighbors': 15, 'min_dist': 0.5},
]

for ax, p in zip(axes.flat, params):
    reducer = umap.UMAP(**p, random_state=42)
    X_umap = reducer.fit_transform(X_scaled)
    ax.scatter(X_umap[:, 0], X_umap[:, 1], c=y, cmap='viridis', alpha=0.7)
    ax.set_title(f"n_neighbors={p['n_neighbors']}, min_dist={p['min_dist']}")

plt.tight_layout()
plt.show()

UMAP vs t-SNE Comparison

Using UMAP for ML Preprocessing

Unlike t-SNE, UMAP can be used as a preprocessing step for machine learning:

import umap
from sklearn.model_selection import train_test_split
from sklearn.ensemble import RandomForestClassifier
from sklearn.datasets import load_digits

# Load data
digits = load_digits()
X, y = digits.data, digits.target

# Split data
X_train, X_test, y_train, y_test = train_test_split(X, y, test_size=0.2, random_state=42)

# Fit UMAP on training data
reducer = umap.UMAP(n_components=10, random_state=42)
X_train_umap = reducer.fit_transform(X_train)

# Transform test data using fitted UMAP
X_test_umap = reducer.transform(X_test)

# Train classifier on UMAP features
clf = RandomForestClassifier(random_state=42)
clf.fit(X_train_umap, y_train)
print(f"Accuracy: {clf.score(X_test_umap, y_test):.2%}")

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

iris = load_iris()
X_scaled = StandardScaler().fit_transform(iris.data)

reducer = umap.UMAP(n_components=2, random_state=42)
X_umap = reducer.fit_transform(X_scaled)

print(f"Shape: {X_umap.shape}")

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

iris = load_iris()
X_scaled = StandardScaler().fit_transform(iris.data)

reducer = umap.UMAP(n_neighbors=30, min_dist=0.05, random_state=42)
X_umap = reducer.fit_transform(X_scaled)

print(f"Shape: {X_umap.shape}")

import umap
from sklearn.model_selection import train_test_split
from sklearn.preprocessing import StandardScaler
from sklearn.datasets import load_iris

iris = load_iris()
X_train, X_test, y_train, y_test = train_test_split(
    iris.data, iris.target, test_size=0.2, random_state=42
)

scaler = StandardScaler()
X_train_scaled = scaler.fit_transform(X_train)
X_test_scaled = scaler.transform(X_test)

reducer = umap.UMAP(n_components=2, random_state=42)
X_train_umap = reducer.fit_transform(X_train_scaled)
X_test_umap = reducer.transform(X_test_scaled)

print(f"Train: {X_train_umap.shape}, Test: {X_test_umap.shape}")

import umap
import matplotlib.pyplot as plt
from sklearn.preprocessing import StandardScaler
from sklearn.datasets import load_iris

iris = load_iris()
X_scaled = StandardScaler().fit_transform(iris.data)

reducer = umap.UMAP(n_components=2, random_state=42)
X_umap = reducer.fit_transform(X_scaled)

plt.figure(figsize=(8, 6))
plt.scatter(X_umap[:, 0], X_umap[:, 1], c=iris.target, cmap='viridis')
plt.xlabel('UMAP 1')
plt.ylabel('UMAP 2')
plt.colorbar(label='Species')
plt.title('UMAP of Iris Dataset')
plt.show()

Confused about when to use PCA, t-SNE, or UMAP? Follow this visual guide to choose the right dimensionality reduction technique for your specific situation!