Dimensionality Reduction & Manifold Learning

In modern machine learning, datasets frequently span hundreds to millions of features—from genomic sequencing matrices (D > 20,000) to high-resolution imagery and LLM latent embeddings (D > 4,096). As the feature dimensionality D expands, Euclidean geometry exhibits pathological behaviors termed the Curse of Dimensionality (Bellman, 1961):

Machine learning practitioners must carefully distinguish between Feature Selection and Feature Reduction (Extraction):

Principal Component Analysis (PCA) is an unsupervised linear transformation technique that identifies orthogonal axes—termed Principal Components—along which the variance of the data is maximized. It constructs a new coordinate system where feature covariance is zero (complete decorrelation).

Step-by-Step Algorithmic Derivation

1. Step 1: Standardization & Mean-Centering: Because PCA is variance-driven, variables with large numerical magnitudes dominate principal axes. Data must be zero-centered (and standardized to unit variance):
   Formulation: Formulation: Xstd = (X - μ / σ) quad such that quad E[X] = 0, quad Var(x_j) = 1
2. Step 2: Covariance Matrix Computation: Compute the sample covariance matrix Σ in RD × D, capturing all pairwise linear relationships:
   Formulation: Formulation: Σ = (1 / N - 1) Xstd^T Xstd
3. Step 3: Spectral Eigenvalue Decomposition: Compute the eigenvalues λ_i and corresponding eigenvectors v_i satisfying the characteristic equation:
   Formulation: Formulation: Σ v_i = λ_i v_i quad ⇔ quad (Σ - λ_i I) v_i = 0
   Every eigenvector v_i in R^D represents an orthogonal principal direction; its eigenvalue λ_i equals the exact variance of the data projected onto that axis.
4. Step 4: Sorting & Projection Matrix Assembly: Sort the eigenvalues in descending order (λ_1 ≥ λ_2 ≥ dots ≥ λ_D ≥ 0). Select the top d eigenvectors to form the transformation matrix W_d = [v_1, v_2, dots, v_d] in RD × d.
5. Step 5: Latent Subspace Projection: Project the original N × D data matrix onto the d-dimensional subspace:
   Formulation: Formulation: Z = Xstd W_d in RN × d

The total variance of the standardized dataset equals the trace of the covariance matrix, which is the sum of all eigenvalues: Total Var = Sum i=1^D λ_i = D.

Heuristics for Retaining Components

• Variance Threshold Rule: Select the minimum d such that CEV_d ≥ 0.90 or 0.95 (retaining 90--95 % of total system information).
• Scree Plot Elbow Criterion (Cattell): Plot eigenvalues λ_i against component index i. Locate the point of inflection ("elbow") where the curve levels off into trivial random noise.
• Kaiser-Guttman Criterion: For standardized variables (Var=1), retain only components with λ_i > 1.0 (components that explain more variance than a single original variable).

While PCA is unsupervised (ignoring class labels and seeking maximum overall variance), Linear Discriminant Analysis (LDA) is a supervised dimensionality reduction technique. LDA maximizes class separability by identifying linear combinations of features that maximize between-class variance while simultaneously minimizing within-class variance.

Fisher's Criterion & Generalized Eigensystem

Fisher's criterion maximizes the Rayleigh quotient:

Taking the gradient with respect to w and setting it to zero yields the Generalized Eigenvalue Problem:

Linear methods (PCA, LDA) project data onto flat hyperplanes. When datasets inhabit curved, non-linear manifolds (e.g., Swiss rolls or high-dimensional single-cell RNA clusters), linear projections cause distant folds to overlap catastrophically. Non-linear manifold learning preserves local geometric neighborhoods.

1. t-SNE (t-Distributed Stochastic Neighbor Embedding)

Developed by van der Maaten and Hinton (2008), t-SNE converts Euclidean distances between data points into conditional probabilities that represent affinities:

• High-Dimensional Gaussian Affinity:
  Formulation:

  pji = ( exp(-x_i - x_j^2 / 2σ_i^2) / Sum k ≠ i exp(-x_i - x_k^2 / 2σ_i^2)), quad pij = (pji + pij / 2N)
  The bandwidth σ_i is determined dynamically via binary search to match a user-defined Perplexity (effective neighborhood size, typically 30--50).
• Low-Dimensional Student-t Distribution (The Crowding Problem Solution): In low dimensions (d=2), points are mapped to y_i, y_j. To prevent moderate-distance points from crowding together into an indistinguishable clump, t-SNE uses a heavy-tailed Cauchy (Student-t with 1 DOF) distribution:
  Formulation:

  qij = ((1 + y_i - y_j^2)-1 / Sum k Sum l ≠ k (1 + y_k - y_l^2)-1)
• KL Divergence Optimization: Points y_i are positioned by minimizing the Kullback-Leibler divergence between high-dimensional distribution P and low-dimensional distribution Q via gradient descent:
  Formulation:

  Lt-SNE = DKL(P parallel Q) = Sum _i Sum _j pij log ( (pij / qij) right)

2. UMAP (Uniform Manifold Approximation & Projection)

Introduced by McInnes et al. (2018), UMAP is grounded in Riemannian geometry and algebraic topology:

The following self-contained script synthesizes a high-dimensional multiclass dataset (N=1,200, D=30) and benchmarks PCA variance decomposition, supervised LDA class compression, and t-SNE non-linear 2D projection.

# =====================================================================
# DevAscent Academy: Dimensionality Reduction & Manifold Learning Pipeline
# Requires: numpy, pandas, scikit-learn
# =====================================================================

import numpy as np
import pandas as pd
from sklearn.datasets import make_classification
from sklearn.preprocessing import StandardScaler
from sklearn.decomposition import PCA
from sklearn.discriminant_analysis import LinearDiscriminantAnalysis as LDA
from sklearn.manifold import TSNE
from sklearn.linear_model import LogisticRegression
from sklearn.model_selection import train_test_split
from sklearn.metrics import accuracy_score

def execute_reduction_benchmark():
    # -----------------------------------------------------------------
    # 1. Synthesize High-Dimensional Multiclass Data (N=1200, D=30, C=3)
    # -----------------------------------------------------------------
    X_raw, y = make_classification(
        n_samples=1200,
        n_features=30,
        n_informative=12,
        n_redundant=10,
        n_repeated=0,
        n_classes=3,
        n_clusters_per_class=2,
        random_state=42
    )

    X_train, X_test, y_train, y_test = train_test_split(
        X_raw, y, test_size=0.25, random_state=42, stratify=y
    )

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

    print(f"[Dataset] Features: {X_train.shape[1]} | Samples: {X_train.shape[0]} train, {X_test.shape[0]} test")

    # Baseline Classifier on Raw 30 Dimensions
    clf_baseline = LogisticRegression(max_iter=1000, random_state=42)
    clf_baseline.fit(X_train_scaled, y_train)
    acc_baseline = accuracy_score(y_test, clf_baseline.predict(X_test_scaled))
    print(f"[Baseline] 30-D Raw Feature Space Accuracy: {acc_baseline * 100:.2f}%")

    # -----------------------------------------------------------------
    # 2. Principal Component Analysis (Spectral Variance Decomposition)
    # -----------------------------------------------------------------
    print("\n--- [Phase 1] Principal Component Analysis (PCA) ---")
    pca_full = PCA().fit(X_train_scaled)
    evr = pca_full.explained_variance_ratio_
    cev = np.cumsum(evr)

    # Determine dimensions for 85% and 95% variance
    d_85 = int(np.argmax(cev >= 0.85) + 1)
    d_95 = int(np.argmax(cev >= 0.95) + 1)

    print(f"Top 5 Individual Component Variances: {[round(float(v), 4) for v in evr[:5]]}")
    print(f"Components required for 85% Variance: {d_85} of 30")
    print(f"Components required for 95% Variance: {d_95} of 30")

    # Fit Optimal PCA (Retaining 85% variance)
    pca_optimal = PCA(n_components=d_85, random_state=42)
    X_train_pca = pca_optimal.fit_transform(X_train_scaled)
    X_test_pca = pca_optimal.transform(X_test_scaled)

    clf_pca = LogisticRegression(max_iter=1000, random_state=42)
    clf_pca.fit(X_train_pca, y_train)
    acc_pca = accuracy_score(y_test, clf_pca.predict(X_test_pca))
    print(f"PCA ({d_85}-D Subspace) Downstream Accuracy: {acc_pca * 100:.2f}%")

    # -----------------------------------------------------------------
    # 3. Linear Discriminant Analysis (Supervised Separation)
    # -----------------------------------------------------------------
    print("\n--- [Phase 2] Linear Discriminant Analysis (LDA) ---")
    # For C=3 classes, max components = C - 1 = 2
    lda = LDA(n_components=2)
    X_train_lda = lda.fit_transform(X_train_scaled, y_train)
    X_test_lda = lda.transform(X_test_scaled)

    clf_lda = LogisticRegression(max_iter=1000, random_state=42)
    clf_lda.fit(X_train_lda, y_train)
    acc_lda = accuracy_score(y_test, clf_lda.predict(X_test_lda))
    print(f"LDA (2-D Projection) Downstream Accuracy: {acc_lda * 100:.2f}%")
    print(f"LDA Explained Variance Ratio (C-1=2): {[round(float(v), 4) for v in lda.explained_variance_ratio_]}")

    # -----------------------------------------------------------------
    # 4. t-SNE Non-Linear Manifold Embedding (Visualization Subspace)
    # -----------------------------------------------------------------
    print("\n--- [Phase 3] t-SNE Manifold Visualization Embedding ---")
    tsne = TSNE(n_components=2, perplexity=35, random_state=42, max_iter=1000)
    X_tsne_2d = tsne.fit_transform(X_train_scaled[:400])
    print(f"Successfully projected 400 high-D points to 2D manifold embedding.")
    print(f"t-SNE Divergence Loss (Kullback-Leibler): {round(float(tsne.kl_divergence_), 4)}")

    # -----------------------------------------------------------------
    # 5. Reduction Summary Matrix
    # -----------------------------------------------------------------
    summary_df = pd.DataFrame([
        {'Method': 'Raw Uncompressed', 'Dimensions': 30, 'Test_Accuracy': f"{acc_baseline*100:.2f}%", 'Type': 'Baseline'},
        {'Method': f'PCA (85% Var)', 'Dimensions': d_85, 'Test_Accuracy': f"{acc_pca*100:.2f}%", 'Type': 'Unsupervised Linear'},
        {'Method': 'LDA (C - 1)', 'Dimensions': 2, 'Test_Accuracy': f"{acc_lda*100:.2f}%", 'Type': 'Supervised Linear'},
        {'Method': 't-SNE (2D)', 'Dimensions': 2, 'Test_Accuracy': 'Visual Only', 'Type': 'Non-linear Manifold'}
    ])
    print("\n" + summary_df.to_string(index=False))

if __name__ == '__main__':
    execute_reduction_benchmark()