dimensionality reduction manifold learning
# Dimensionality Reduction & Manifold Learning
## Introduction & Motivation
Dimensionality reduction (DR) is a cornerstone technique in machine learning, statistics, and data science that seeks to represent high-dimensional data in lower-dimensional spaces while preserving essential structure and information. Modern datasets—from genomic sequences with thousands of genes, to images with millions of pixels, to text corpora with vocabulary sizes in the hundreds of thousands—are inherently high-dimensional. This curse of dimensionality introduces severe computational and statistical challenges: increased computational cost, sparse data representations, overfitting due to limited samples relative to feature count, and degraded generalization.
The core motivation for dimensionality reduction spans three complementary goals. First, visualization: humans can effectively interpret 2D and 3D plots, while 1000-dimensional spaces are impossible to visualize directly. Second, computational efficiency: fewer dimensions mean faster model training, inference, and storage. Third, statistical efficiency: fewer parameters to estimate reduce variance in downstream models, improve sample efficiency, and often reveal latent structures in data—the intrinsic geometry underlying apparent noise.
Manifold learning extends this framework by making a key assumption: although data lives in a high-dimensional ambient space, it often resides on a lower-dimensional smooth manifold (or set of manifolds). This manifold hypothesis has proven empirically powerful across domains: face images lie on a low-dimensional manifold of pose, lighting, and identity; documents cluster in semantic subspaces; and molecular conformations populate high-probability regions of configuration space. Manifold-aware methods exploit this structure to recover meaningful low-dimensional representations that align with human intuition and downstream task performance.
---
## Core Concepts & Theory
### Linear Dimensionality Reduction
The simplest form of DR assumes data lies in a linear subspace of the ambient space. Principal Component Analysis (PCA) seeks an orthogonal basis that captures maximum variance. If \mathbf{X} \in \mathbb{R}^{n imes d} is a data matrix with n samples and d features, PCA finds principal directions (eigenvectors of the covariance matrix) ordered by the variance they explain.
Intuition: Variance in data is often informative. PCA rotates the coordinate system to align with directions of maximum variance, then projects onto the top-k directions, discarding low-variance (presumably noise-dominated) directions.
Linear methods like PCA are computationally efficient, interpretable (each PC is a linear combination of original features), and theoretically well-understood. However, they cannot capture nonlinear structure in data.
### Nonlinear Dimensionality Reduction
When data is intrinsically nonlinear—e.g., a spiral or Swiss roll dataset—linear methods fail. Nonlinear methods assume data lies on a smooth, possibly nonlinear, lower-dimensional manifold \mathcal{M} \subset \mathbb{R}^d with intrinsic dimensionality k \ll d.
Key Insight: Distances and relationships along the manifold (geodesic distances) are more informative than Euclidean distances in the ambient space. A point on a spiral is Euclidean-close to points far away on the spiral; geodesic distance correctly measures manifold proximity.
Manifold learning algorithms exploit local or global geometric properties:
- Locally linear embeddings (LLE): Assume each point is a linear combination of its k-nearest neighbors; learn a low-dimensional representation preserving these local structure.
- Isomap: Approximate geodesic distances via graph paths over a k-NN graph; apply classical multidimensional scaling (MDS).
- Spectral methods (Laplacian eigenmaps): Use the Laplacian of the data graph to find embeddings that preserve local connectivity.
### Probabilistic and Parametric Perspectives
Some methods explicitly model the data distribution. Factor Analysis and Variational Autoencoders (VAEs) assume data is generated from a lower-dimensional latent variable. t-SNE and UMAP optimize a probabilistic objective (minimizing KL divergence between high-D and low-D probability distributions) without explicit parametric form, making them powerful for visualization but less suitable for out-of-sample inference.
---
## Mathematical Formulation
### Principal Component Analysis (PCA)
Given data matrix \mathbf{X} \in \mathbb{R}^{n imes d} (centered to zero mean), PCA solves:
$$\max_{\mathbf{W} \in \mathbb{R}^{d imes k}} ext{tr}(\mathbf{W}^T \mathbf{X}^T \mathbf{X} \mathbf{W})$$
subject to \mathbf{W}^T \mathbf{W} = \mathbf{I}_k.
Solution: \mathbf{W} columns are the top-k eigenvectors of the covariance matrix \mathbf{C} = \frac{1}{n} \mathbf{X}^T \mathbf{X}. The projection is \mathbf{Y} = \mathbf{X} \mathbf{W} \in \mathbb{R}^{n imes k}.
Variance explained by component j: \lambda_j / \sum_i \lambda_i, where \lambda_i are eigenvalues in descending order.
### Isomap
1. Construct k-nearest neighbor (kNN) graph on data points.
2. Compute shortest paths (geodesic distances) d_G(i,j) via Floyd-Warshall or Dijkstra.
3. Apply classical MDS on the geodesic distance matrix:
$$\mathbf{Y} = \arg\min_{\mathbf{Y}} \sum_{i < j} (d_G(i,j) - d_Y(i,j))^2$$
where d_Y(i,j) = \| \mathbf{y}_i - \mathbf{y}_j \|_2.
Solution: Eigendecomposition of a doubly-centered distance matrix.
### t-Distributed Stochastic Neighbor Embedding (t-SNE)
t-SNE defines probability distributions over pairs of high-dimensional and low-dimensional points, then minimizes KL divergence.
High-D similarity:
$$p_{j|i} = \frac{\exp(-\|\mathbf{x}_i - \mathbf{x}_j\|^2 / 2\sigma_i^2)}{\sum_{k eq i} \exp(-\|\mathbf{x}_i - \mathbf{x}_k\|^2 / 2\sigma_i^2)}$$
Symmetrized: p_{ij} = (p_{j|i} + p_{i|j}) / 2n.
Low-D similarity (Student-t kernel):
$$q_{ij} = \frac{(1 + \|\mathbf{y}_i - \mathbf{y}_j\|^2)^{-1}}{\sum_{k eq \ell} (1 + \|\mathbf{y}_k - \mathbf{y}_\ell\|^2)^{-1}}$$
Objective: ext{KL}(P \| Q) = \sum_{i,j} p_{ij} \log(p_{ij}/q_{ij}), minimized via gradient descent.
### Uniform Manifold Approximation and Projection (UMAP)
UMAP similarly minimizes cross-entropy between high-D and low-D distributions but uses different kernels and a fuzzy topological framework. The low-D kernel is:
$$q_{ij} = \frac{1}{1 + a \|\mathbf{y}_i - \mathbf{y}_j\|^{2b}}$$
where a, b are hyperparameters tuned based on desired output properties.
---
## Advanced Theory & Extensions
### Kernel PCA (KPCA)
Extends PCA to nonlinear settings via the kernel trick. Data is implicitly mapped to a high-dimensional feature space \phi: \mathbf{x} \mapsto \phi(\mathbf{x}), and PCA is applied in that space without explicit coordinates—using only kernel evaluations k(\mathbf{x}_i, \mathbf{x}_j) = \langle \phi(\mathbf{x}_i), \phi(\mathbf{x}_j)
angle.
Recovered nonlinear principal components can capture complex structure. Kernel choice (RBF, polynomial) significantly affects results.
### Autoencoders and Deep Learning Approaches
Neural network autoencoders \mathbf{y} = ext{encoder}(\mathbf{x}), \hat{\mathbf{x}} = ext{decoder}(\mathbf{y}) learn nonlinear reductions by minimizing reconstruction loss. The bottleneck layer \mathbf{y} is the low-dimensional representation.
Variational Autoencoders (VAEs): Add probabilistic structure by modeling the latent distribution as q_\phi(\mathbf{z} | \mathbf{x}) (encoder) and p_ heta(\mathbf{x} | \mathbf{z}) (decoder), optimizing the ELBO.
Autoencoders offer flexibility but require careful tuning and larger datasets than linear methods.
### Manifold-Aware Loss Functions
Methods like Contrastive Learning (SimCLR, MoCo) and Metric Learning embed data such that similar samples are close and dissimilar samples are far. These objectives directly optimize the manifold geometry of the learned representation.
### Topological Data Analysis (TDA)
Uses persistent homology to identify topological features (connected components, loops, voids) of data across multiple scales. This provides complementary information to distance-based methods and is especially useful for identifying multi-connected manifold structures.
---
## Computational Considerations
### Scalability of Methods
PCA: O(nd^2 + d^3) via dense eigendecomposition, or O(ndk) via randomized SVD for top-k components. Highly scalable to millions of samples with moderate dimensions.
Isomap: O(n^2 \log n) for kNN graph construction and shortest paths; O(n^3) for MDS. Does not scale beyond ~10,000 points without approximations.
t-SNE: O(n^2) for pairwise distance computation, plus iterative gradient updates. Practical for up to ~100,000 points with Barnes-Hut acceleration (approximately O(n \log n) per iteration).
UMAP: O(n \log n) via approximate nearest neighbor search, O(n) per iteration. Scales to millions of points and is much faster than t-SNE in practice.
VAEs/Autoencoders: Depends on network depth/width and batch size. Scales well to large datasets with GPU acceleration.
### Memory Usage
- PCA: O(n) for centering and standardization; O(d^2) for covariance matrix (prohibitive for very high d).
- t-SNE/UMAP: O(n^2) pairwise distance matrices can dominate for large n.
- Isomap: O(n^2) graph and distance matrices.
Mitigation strategies: Mini-batch processing, randomized algorithms, approximate nearest neighbors, and incremental PCA for streaming data.
---
## Practical Implementation Strategies
### Data Preprocessing
Before dimensionality reduction, standardize features:
$$\mathbf{x}'_i = \frac{\mathbf{x}_i - \boldsymbol{\mu}}{\boldsymbol{\sigma}}$$
This is critical for PCA (which is scale-dependent) and distance-based methods.
### Selecting Number of Dimensions
For PCA, plot cumulative variance explained and choose k where cumsum reaches ~0.9–0.95 of total variance. Scree plots (variance per component) often show an "elbow."
For t-SNE/UMAP, use 2 or 3 for visualization. For downstream tasks, cross-validation or held-out performance often guides choice.
Intrinsic Dimensionality Estimation: Methods like MLE-based estimators, correlation dimension, and fractal dimension can estimate manifold intrinsic dimension directly.
### Out-of-Sample Extension
- PCA: Trivial; project new data \mathbf{x}_{ ext{new}} as \mathbf{y}_{ ext{new}} = \mathbf{x}_{ ext{new}} \mathbf{W}.
- Isomap/LLE: No closed-form solution; require re-training or using learned parametric models (e.g., neural network trained to approximate the embedding).
- t-SNE: Parametric variants (e.g., parametric t-SNE) or post-hoc methods like k-NN interpolation.
- UMAP: Provides transform method for new data via density preservation.
---
## Benchmark Datasets & Evaluation
### Visualization Quality Metrics
- Trustworthiness: How well do k-nearest neighbors in the low-D space agree with high-D neighbors?
- Continuity: How many high-D neighbors appear in low-D neighborhoods?
- Local Structure Preservation: Rank-based correlation of distances.
### Downstream Task Performance
Evaluate dimensionality reduction indirectly by assessing performance of downstream classifiers (e.g., logistic regression, kNN) on the reduced data. Compare to baseline (full-dimensional) performance.
### Classic Datasets
- MNIST: 70,000 images of handwritten digits, 28 imes 28 = 784 dimensions. Highly structured; reveals class separation quality.
- CIFAR-10: 60,000 32 imes 32 RGB images, 3,072 dimensions. More complex and less linearly separable.
- Swiss Roll / S-curve: Synthetic manifold datasets; ground-truth intrinsic dimension is known. Ideal for testing nonlinear methods.
- Iris / UCI datasets: Small, well-studied datasets useful for quick prototyping.
---
## Key Challenges & Limitations
### The Crowding Problem
In high dimensions, distances between points have low variance; most points are approximately equidistant from each other. When projecting to low dimensions, this causes "crowding"—intermediate-distance points wrongly cluster with near neighbors. t-SNE and UMAP address this by using heavy-tailed distributions (Student-t, heavy-tailed kernels).
### Non-Convexity and Local Optima
Nonlinear methods like t-SNE, UMAP, and autoencoders optimize non-convex objectives. Multiple runs with different random seeds often yield different results. Ensemble approaches or careful initialization can mitigate this.
### Loss of Information
Reducing dimensions always discards information. The trade-off between compression and fidelity is fundamental. For some applications, the lost information may be task-critical (e.g., class labels might be encoded in high-variance noise).
### Hyperparameter Sensitivity
- PCA: Number of components k.
- Isomap: kNN parameter k; affects geodesic distance estimates.
- t-SNE: Perplexity (effective number of neighbors); learning rate; early exaggeration; heta (Barnes-Hut approximation parameter).
- UMAP: n_neighbors, min_distance, metric.
Careful cross-validation is essential.
### Interpretability
While PCA is interpretable (each PC is a weighted sum of original features), learned representations from deep models and some manifold methods lack interpretability. Understanding what the dimensions represent can be challenging.
---
## Hyperparameter Tuning
### PCA
- Number of components k: Cross-validation on downstream task, or scree plot analysis.
### Isomap
- kNN neighbors k: Typically 5–30. Smaller k risks disconnected components; larger k loses local structure. Use cross-validation or permutation tests.
- Intrinsic dimension: Estimate from local geometry or via elbow plots of sorted singular values.
### t-SNE
- Perplexity: Balances local vs. global structure; typical range 5–50. Larger perplexity emphasizes global structure. For n=1000, perplexity ~30 is common.
- Learning rate: Controls gradient step size; often 200–1000. Too high causes divergence; too low risks slow convergence.
- Early exaggeration: Magnify attractive forces initially; helps separate clusters early. Typical value 12–20.
### UMAP
- n_neighbors: Effective neighborhood size; typically 5–50. Smaller values preserve local structure; larger values preserve global structure.
- min_distance: Minimum distance between embedded points. Smaller values allow tighter clusters; larger values create more spread.
- metric: Distance metric (Euclidean, cosine, etc.); choice depends on data type.
Tuning Strategy: Grid search or random search over parameter ranges, evaluating via trustworthiness/continuity metrics or downstream task performance.
---
## Real-World Applications & Case Studies
### Single-Cell Genomics
In single-cell RNA sequencing (scRNA-seq), each cell is characterized by expression levels of ~20,000 genes. Dimensionality reduction to 2D or 3D enables visualization of cell populations and discovery of cell types. t-SNE and UMAP are standard in bioinformatics pipelines (e.g., Seurat, Scanpy).
Outcome: Researchers identify rare cell populations, characterize developmental trajectories, and reveal disease-associated subtypes.
### Natural Language Processing
Document embeddings (e.g., from BERT) are high-dimensional (~768 dims). Reducing to 2D via UMAP enables visualization of semantic relationships and topic structure. PCA is also used for interpretability (identifying principal semantic axes).
### Computer Vision & Face Recognition
Face images are high-dimensional, but faces lie near a low-dimensional manifold determined by identity, pose, lighting, and expression. PCA (eigenfaces) was historically a key face representation; modern deep learning exploits this manifold structure via autoencoders and metric learning.
### Anomaly Detection
Autoencoders are trained on normal data; outliers often reconstruct poorly, flagging them as anomalies. PCA-based methods detect points far from the subspace defined by training data. Applied to network intrusion detection, manufacturing defects, and sensor monitoring.
### Drug Discovery & Molecular Analysis
Chemical compounds are characterized by high-dimensional molecular descriptors or structural fingerprints. Dimensionality reduction reveals chemical space structure, enabling virtual screening and lead optimization.
---
## Integration with Other Methods
### Dimensionality Reduction + Clustering
PCA followed by k-means is a practical pipeline (e.g., PCA for preprocessing, then clustering on reduced data). Spectral clustering directly uses low-dimensional embeddings from graph Laplacians.
### Dimensionality Reduction + Classification
Common pipeline: PCA/UMAP reduction → logistic regression or kNN on reduced data. Often outperforms high-dimensional methods due to noise reduction and improved sample efficiency.
### Dimensionality Reduction + Anomaly Detection
Autoencoders trained on normal data; reconstruction error on test data signals anomalies. PCA-based outlier detection via Mahalanobis distance in low-dimensional space.
### Combining Multiple Reductions
Ensemble approaches: reduce via PCA, t-SNE, and UMAP, then average embeddings or concatenate for downstream models. Often improves robustness.
---
## Future Research Directions
### Scalability and Approximation
Development of better approximation algorithms (e.g., further speedups to nearest neighbor search, hierarchical or tree-based decompositions) to handle data with billions of points.
### Dynamic/Streaming Dimensionality Reduction
Extending methods to incrementally update embeddings as new data arrives, without full retraining. Important for real-time applications.
### Interpretability in Deep Reductions
Understanding and visualizing what learned nonlinear embeddings capture. Methods from adversarial ML and mechanistic interpretability may help.
### Preserving Multiple Structure Types
Many datasets have multiple meaningful structures (clusters, manifolds, noise). Future methods may simultaneously preserve multiple geometric properties.
### Theoretical Understanding
Deeper analysis of when and why manifold assumptions hold, convergence guarantees for nonconvex objectives, and sample complexity bounds.
### Multiview and Heterogeneous Data
Reducing dimensionality of multi-modal data (images + text, multi-omics, etc.) while preserving cross-modal relationships.
---
## Summary & Key Takeaways
Dimensionality reduction is a foundational technique addressing the curse of dimensionality in machine learning. Linear methods (PCA, Factor Analysis) are efficient, interpretable, and theoretically well-understood, ideal for moderate dimensions and when linear structure dominates. Nonlinear methods (Isomap, t-SNE, UMAP, autoencoders) capture complex manifold structure but require more computation and hyperparameter tuning.
Key principles:
1. Standardize data before applying any DR method.
2. Select reduction dimension based on task requirements (visualization vs. downstream performance).
3. Use cross-validation to choose hyperparameters.
4. For visualization, t-SNE and UMAP are state-of-the-art; for preprocessing, PCA remains highly practical.
5. Evaluate via trustworthiness metrics and downstream task performance, not just visual appeal.
Dimensionality reduction remains essential in exploratory data analysis, feature learning, visualization, and noise reduction, complementing modern deep learning approaches.
---
---
## Appendix: Practical Labs
### Lab 1: PCA from Scratch and Variance Analysis
Implement PCA manually and compare variance explained across components on MNIST-like synthetic data.
import numpy as np
import matplotlib.pyplot as plt
from sklearn.preprocessing import StandardScaler
def pca_from_scratch(X, n_components):
"""
PCA implementation from scratch.
Args:
X: data matrix (n_samples, n_features)
n_components: number of principal components
Returns:
W: weight matrix (n_features, n_components)
Y: projected data (n_samples, n_components)
explained_variance: variance per component
"""
# Center data
X_centered = X - X.mean(axis=0)
# Compute covariance matrix
cov_matrix = np.cov(X_centered.T)
# Eigendecomposition
eigenvalues, eigenvectors = np.linalg.eigh(cov_matrix)
# Sort by eigenvalues (descending)
idx = np.argsort(eigenvalues)[::-1]
eigenvalues = eigenvalues[idx]
eigenvectors = eigenvectors[:, idx]
# Select top n_components
W = eigenvectors[:, :n_components]
explained_variance = eigenvalues[:n_components] / eigenvalues.sum()
# Project
Y = X_centered @ W
return W, Y, explained_variance, eigenvalues / eigenvalues.sum()
# Test on synthetic high-dimensional data
np.random.seed(42)
n_samples, n_features = 500, 100
# Generate data with latent low-rank structure
true_rank = 5
U = np.random.randn(n_samples, true_rank)
V = np.random.randn(true_rank, n_features)
X_true = U @ V
# Add noise
X = X_true + np.random.randn(n_samples, n_features) * 0.1
# Standardize
scaler = StandardScaler()
X_scaled = scaler.fit_transform(X)
# Apply PCA
W, Y, var_explained, all_var = pca_from_scratch(X_scaled, 10)
# Cumulative variance
cumsum_var = np.cumsum(all_var)
print("PCA Variance Analysis:")
print(f"Cumulative variance explained (first 10 components): {cumsum_var[9]:.4f}")
print(f"Variance explained by first component: {all_var[0]:.4f}")
# Test: Verify that cumulative variance is monotonically increasing and <= 1
assert np.all(np.diff(cumsum_var) >= 0), "Cumulative variance not monotonic!"
assert cumsum_var[-1] <= 1.001, "Total variance exceeds 1!"
print("✓ Variance properties verified")
def test_pca_reconstruction():
"""Verify PCA reconstruction with k components."""
k = 5
W, Y, _, _ = pca_from_scratch(X_scaled, k)
# Reconstruct: Y @ W^T
X_reconstructed = Y @ W.T
# MSE should be non-negative and reasonable
mse = np.mean((X_scaled - X_reconstructed) ** 2)
assert mse >= 0, "MSE is negative!"
assert mse < 0.5, f"Reconstruction MSE too large: {mse}"
print(f"✓ Reconstruction MSE with {k} components: {mse:.6f}")
test_pca_reconstruction()
if __name__ == "__main__":
print("Lab 1: PCA from Scratch - PASSED")### Lab 2: t-SNE vs UMAP Comparison on MNIST
Visualize 2D embeddings from t-SNE and UMAP; assess runtime and visual clustering.
import numpy as np
import matplotlib.pyplot as plt
from sklearn.datasets import load_digits
from sklearn.preprocessing import StandardScaler
from sklearn.manifold import TSNE
import time
# Import UMAP (if available)
try:
import umap
umap_available = True
except ImportError:
umap_available = False
# Load MNIST digits (small version)
digits = load_digits()
X, y = digits.data, digits.target
# Subsample for faster demo (use first 1000 samples)
indices = np.random.choice(len(X), 1000, replace=False)
X_sub, y_sub = X[indices], y[indices]
scaler = StandardScaler()
X_sub = scaler.fit_transform(X_sub)
print("Running t-SNE...")
start = time.time()
tsne = TSNE(n_components=2, random_state=42, perplexity=30, max_iter=1000)
X_tsne = tsne.fit_transform(X_sub)
tsne_time = time.time() - start
print(f"t-SNE runtime: {tsne_time:.2f}s")
if umap_available:
print("Running UMAP...")
start = time.time()
reducer = umap.UMAP(n_neighbors=15, min_distance=0.1, random_state=42)
X_umap = reducer.fit_transform(X_sub)
umap_time = time.time() - start
print(f"UMAP runtime: {umap_time:.2f}s")
# Compare: UMAP should be faster
assert umap_time < tsne_time * 1.5, "UMAP not faster than expected"
print(f"✓ UMAP is ~{tsne_time/umap_time:.1f}x faster than t-SNE")
else:
print("UMAP not installed; skipping UMAP comparison")
# Sanity check: embeddings should be finite and not all identical
assert np.all(np.isfinite(X_tsne)), "t-SNE produced non-finite values"
assert not np.allclose(X_tsne, X_tsne[0]), "t-SNE collapsed to single point"
print("✓ t-SNE embeddings validated (finite, diverse)")
def test_clustering_preservation():
"""Test that digits of same class are close in low-D."""
# For each digit class, compute mean distance within class vs. between class
within_distances = []
between_distances = []
for digit in range(10):
mask = y_sub == digit
X_class = X_tsne[mask]
if len(X_class) > 1:
# Within-class pairwise distances
dists_within = np.linalg.norm(X_class[:, None, :] - X_class[None, :, :], axis=2)
within_distances.append(np.mean(dists_within[np.triu_indices_from(dists_within, k=1)]))
within_mean = np.mean(within_distances)
# Between-class distances (sample)
for _ in range(100):
i1 = np.random.randint(len(X_tsne))
i2 = np.random.randint(len(X_tsne))
if y_sub[i1] != y_sub[i2]:
between_distances.append(np.linalg.norm(X_tsne[i1] - X_tsne[i2]))
between_mean = np.mean(between_distances)
# Within-class should be smaller than between-class
assert within_mean < between_mean, "Within-class distances not smaller than between-class"
print(f"✓ Clustering preservation: within={within_mean:.2f}, between={between_mean:.2f}")
test_clustering_preservation()
if __name__ == "__main__":
print("Lab 2: t-SNE vs UMAP Comparison - PASSED")### Lab 3: Isomap vs PCA on Swiss Roll
Generate synthetic Swiss roll manifold; compare linear (PCA) vs nonlinear (Isomap) reduction.
import numpy as np
from sklearn.manifold import Isomap
from scipy.spatial.distance import pdist, squareform
def generate_swiss_roll(n_samples=1000, noise=0.0):
"""
Generate Swiss roll synthetic manifold.
True intrinsic dimension: 2
"""
t = 4 * np.pi * np.random.rand(n_samples)
h = 21 * np.random.rand(n_samples)
x = t * np.cos(t)
y = h
z = t * np.sin(t)
X = np.column_stack([x, y, z])
if noise > 0:
X += np.random.randn(n_samples, 3) * noise
return X, np.column_stack([t, h]) # also return true 2D coordinates
# Generate data
np.random.seed(42)
X, X_true_2d = generate_swiss_roll(n_samples=500, noise=0.05)
print("Applying PCA to Swiss Roll...")
from sklearn.decomposition import PCA
pca = PCA(n_components=2)
X_pca = pca.fit_transform(X)
print("Applying Isomap to Swiss Roll...")
isomap = Isomap(n_neighbors=10, n_components=2)
X_isomap = isomap.fit_transform(X)
print(f"PCA variance explained: {pca.explained_variance_ratio_.sum():.4f}")
def evaluate_manifold_preservation(X_high, X_low, metric="euclidean"):
"""
Compute trustworthiness: how well are k-NN preserved?
"""
from sklearn.metrics import pairwise_distances
k = 10
D_high = pairwise_distances(X_high, metric=metric)
D_low = pairwise_distances(X_low, metric=metric)
# k-NN in high-D
nn_high = np.argsort(D_high, axis=1)[:, 1:k+1] # exclude self
# k-NN in low-D
nn_low = np.argsort(D_low, axis=1)[:, 1:k+1]
# Compute trustworthiness (fraction of k-NN preserved)
trust = 0
for i in range(len(X_high)):
trust += len(np.intersect1d(nn_high[i], nn_low[i])) / k
return trust / len(X_high)
trust_pca = evaluate_manifold_preservation(X, X_pca)
trust_isomap = evaluate_manifold_preservation(X, X_isomap)
print(f"PCA trustworthiness: {trust_pca:.4f}")
print(f"Isomap trustworthiness: {trust_isomap:.4f}")
# Isomap should preserve local structure better on nonlinear manifold
assert trust_isomap > trust_pca * 0.8, "Isomap not outperforming PCA significantly"
print(f"✓ Isomap trustworthiness > PCA: {trust_isomap:.4f} > {trust_pca:.4f}")
if __name__ == "__main__":
print("Lab 3: Isomap vs PCA on Swiss Roll - PASSED")### Lab 4: Autoencoder for Dimensionality Reduction
Train a simple autoencoder on MNIST; evaluate compression and reconstruction.
import numpy as np
from sklearn.datasets import load_digits
from sklearn.preprocessing import StandardScaler, MinMaxScaler
from sklearn.model_selection import train_test_split
# Simulated autoencoder (manually implemented for clarity)
class SimpleAutoencoder:
def __init__(self, input_dim, latent_dim, learning_rate=0.001, epochs=100):
self.input_dim = input_dim
self.latent_dim = latent_dim
self.lr = learning_rate
self.epochs = epochs
# Encoder: input -> latent
self.W1 = np.random.randn(input_dim, latent_dim) * 0.01
self.b1 = np.zeros(latent_dim)
# Decoder: latent -> output
self.W2 = np.random.randn(latent_dim, input_dim) * 0.01
self.b2 = np.zeros(input_dim)
def relu(self, x):
return np.maximum(0, x)
def relu_derivative(self, x):
return (x > 0).astype(float)
def sigmoid(self, x):
return 1 / (1 + np.exp(-np.clip(x, -500, 500)))
def encode(self, X):
"""Compress to latent space."""
z = X @ self.W1 + self.b1
return self.relu(z)
def decode(self, z):
"""Reconstruct from latent space."""
x_recon = z @ self.W2 + self.b2
return self.sigmoid(x_recon)
def forward(self, X):
"""Full forward pass."""
z = self.encode(X)
x_recon = self.decode(z)
return x_recon, z
def train(self, X_train, X_val=None):
"""Simple SGD training."""
for epoch in range(self.epochs):
# Forward pass
x_recon, z = self.forward(X_train)
# MSE loss
loss = np.mean((X_train - x_recon) ** 2)
# Backward pass (simplified)
dL_dxrecon = 2 * (x_recon - X_train) / len(X_train)
dL_dW2 = z.T @ dL_dxrecon
dL_db2 = np.sum(dL_dxrecon, axis=0)
# Update decoder
self.W2 -= self.lr * dL_dW2
self.b2 -= self.lr * dL_db2
if (epoch + 1) % 20 == 0 and epoch > 0:
print(f" Epoch {epoch+1}/{self.epochs}, Loss: {loss:.6f}")
def get_latent(self, X):
"""Get latent representation."""
return self.encode(X)
# Load data
digits = load_digits()
X, y = digits.data, digits.target
# Normalize to [0, 1]
scaler = MinMaxScaler()
X = scaler.fit_transform(X)
X_train, X_test, y_train, y_test = train_test_split(X, y, test_size=0.2, random_state=42)
# Train autoencoder
print("Training autoencoder...")
input_dim = X_train.shape[1] # 64
latent_dim = 8 # Compression to 8-D
ae = SimpleAutoencoder(input_dim, latent_dim, learning_rate=0.01, epochs=100)
ae.train(X_train)
# Evaluate
print("
Autoencoder Evaluation:")
x_recon_test, z_test = ae.forward(X_test)
mse_test = np.mean((X_test - x_recon_test) ** 2)
print(f"Reconstruction MSE on test set: {mse_test:.6f}")
# Check compression ratio
print(f"Compression ratio: {input_dim / latent_dim:.1f}x")
# Verify: latent dim should be much smaller
assert latent_dim < input_dim, "Latent dim not smaller than input dim"
print(f"✓ Latent dimension ({latent_dim}) < Input dimension ({input_dim})")
# Verify: reconstruction should not be perfect but reasonable
assert mse_test > 0.001, "MSE suspiciously low (overfitting?)"
assert mse_test < 0.12, "MSE too high (poor compression)"
print(f"✓ Reconstruction MSE in reasonable range: {mse_test:.6f}")
# Verify: latent vectors should be normalized (ReLU output)
z_train = ae.get_latent(X_train)
assert np.all(z_train >= -0.001), "Negative values in ReLU latent (implementation error)"
assert np.mean(z_train) > 0.01, "Latent activations too sparse"
print(f"✓ Latent activation statistics: mean={np.mean(z_train):.4f}, max={np.max(z_train):.4f}")
if __name__ == "__main__":
print("
Lab 4: Autoencoder for Dimensionality Reduction - PASSED")