Unsupervised Learning: PCA, Manifolds & Clustering

1. Motivation & Roadmap

In scientific research, we often collect vast amounts of data before we have clear, well-defined labels. Unsupervised learning provides the tools to explore this data, uncover its intrinsic structure, and generate new hypotheses. It primarily addresses two core questions:

1. Compression & Denoising: "What are the dominant factors of variation in my data?" This is the domain of dimensionality reduction.
2. Segmentation & Discovery: "Which of my samples behave similarly?" This is the domain of clustering.

1.1 The Scientific Discovery Workflow

Unsupervised learning is a key part of the broader scientific discovery process. This entire workflow can be visualized as a cycle of hypothesis and discovery.

Our focus is on the "Unsupervised Modeling" block, which itself is a detailed pipeline. Let's explore this with a concrete example.

1.2 A Practical Roadmap: From Raw Data to Insight

Imagine we have collected a dataset of 200 Raman spectra, each with 1024 data points (features). Our goal is to see if there are distinct material phases present in the samples. The path from this raw data (a 200x1024 matrix) to scientific insight follows a structured roadmap.

Step 1: Preprocessing - The Crucial First Step

Raw data is rarely ready for modeling. Preprocessing ensures that the patterns found by the model are scientifically meaningful, not artifacts of data collection or scaling. A good checklist includes:

• Handling Missing Values: Did an instrument fail for some data points? Strategies include removing the sample/feature, or imputing the missing value (e.g., with the mean or median of the feature).
• Outlier Detection: Are there obvious instrumental glitches or corrupted measurements? Robust methods like Robust PCA (discussed later) can help, or they can be removed manually.
• Feature Scaling: This is non-negotiable for many algorithms. Features with large variances can dominate distance calculations and PCA.
  • Standardization (StandardScaler): Rescales features to have a mean of 0 and a standard deviation of 1. This is the default choice for PCA and distance-based clustering.
  • Normalization (MinMaxScaler): Rescales features to a fixed range, typically [0, 1]. Useful when the absolute bounds are important.
• Transformations: If a feature has a highly skewed distribution (e.g., concentrations), a log transform can make its distribution more symmetric and easier to model.

Step 2: Dimensionality Reduction

Working with 1024 features directly is computationally expensive and prone to the "curse of dimensionality." We use techniques like PCA or UMAP to project the data into a low-dimensional space (e.g., 2 or 3 dimensions) that captures the most important variance. This allows for visualization and often improves clustering performance.

Step 3: Clustering

On these low-dimensional embeddings, we apply a clustering algorithm (like K-Means, HDBSCAN, or GMM) to find natural groups in our 200 spectra. The choice of algorithm depends on the expected geometry of the data (e.g., spherical vs. arbitrary shapes).

Step 4: Validation

How good are the clusters? We use internal metrics like the Silhouette score to quantify how dense and well-separated our clusters are. We can also assess their stability by seeing if the same clusters emerge when we re-run the analysis on bootstrapped samples of our data.

Step 5: Interpretation

This is where data science becomes science. If we found 3 stable clusters, we must go back to the original, high-dimensional data. We can plot the average spectrum for each of the 3 groups. Do they correspond to known material phases? Does one cluster show a new, unexpected peak? This is the step that generates new hypotheses to be tested in the lab.

2. Principal Component Analysis (PCA)

Principal Component Analysis (PCA) is a cornerstone of unsupervised learning. It performs a linear transformation to convert a set of observations of possibly correlated variables into a set of values of linearly uncorrelated variables called principal components. This new coordinate system is chosen such that the first principal component has the largest possible variance, and each succeeding component in turn has the highest variance possible under the constraint that it is orthogonal to the preceding components.

2.1 The Mathematics of PCA: Two Perspectives

Understanding PCA can be approached from two equivalent mathematical viewpoints: eigenvalue decomposition and singular value decomposition (SVD).

Perspective 1: Eigenvalue Decomposition of the Covariance Matrix

This is the traditional view. For a centered data matrix \(X \in \mathbb{R}^{n \times d}\) (where \(n\) is samples, \(d\) is features), PCA seeks to find the directions of maximum variance. This is framed as an eigenvalue problem on the covariance matrix \(C = \frac{1}{n-1}X^T X\):

\[C w_i = \lambda_i w_i\]

The eigenvectors \(w_i\) (the principal components) form the axes of the new feature space, and the corresponding eigenvalues \(\lambda_i\) represent the variance captured by each component. The eigenvectors are ordered by their eigenvalues, from largest to smallest.

Perspective 2: Singular Value Decomposition (SVD)

A more numerically stable and general approach is to use SVD on the data matrix \(X\). SVD decomposes \(X\) into three matrices:

\[X = U \Sigma V^T\]

• \(U \in \mathbb{R}^{n \times n}\): The left singular vectors (scores of samples in the new basis).
• \(\Sigma \in \mathbb{R}^{n \times d}\): A diagonal matrix of singular values \(\sigma_i\). The variance explained by each component is \(\lambda_i = \frac{\sigma_i^2}{n-1}\).
• \(V \in \mathbb{R}^{d \times d}\): The right singular vectors, which are the principal components (equivalent to the eigenvectors \(w_i\) from the first perspective).

This view is often preferred in practice due to its numerical stability. The computational cost of SVD-based PCA is typically \(\mathcal{O}(nd^2)\) or \(\mathcal{O}(n^2d)\), depending on which dimension is larger.

2.2 The Toolkit for PCA Interpretation

Performing PCA is easy, but interpreting the results to extract scientific insight is the real skill. Several tools are essential for this.

Scree Plot and Cumulative Variance

How many principal components should we keep? The scree plot helps us decide. It plots the variance (eigenvalue) explained by each principal component, in descending order. We look for an "elbow" in the plot, a point where the explained variance begins to level off. A cumulative variance plot is often overlaid to see the total percentage of variance captured by the first \(k\) components.

Loading Plot

The loadings are the coefficients of the linear combination of the original variables from which the principal components are constructed. A loading plot (often a bar chart) shows the contribution of each original feature to a given principal component. This is the key to understanding the physical meaning of the PCs.

Biplot

A biplot is a powerful visualization that overlays the scores (the projected data points) and the loadings (the original feature directions) in the same plot, typically for PC1 vs. PC2. It allows us to see the relationships between samples and how they relate to the original variables simultaneously.

# Python pseudocode for PCA interpretation
import numpy as np
import matplotlib.pyplot as plt
from sklearn.decomposition import PCA
from sklearn.preprocessing import StandardScaler

# Assume X is our n x d data matrix
X_scaled = StandardScaler().fit_transform(X)
pca = PCA(n_components=10)
X_pca = pca.fit_transform(X_scaled)

# 1. Scree Plot
plt.figure()
plt.bar(range(10), pca.explained_variance_ratio_)
plt.plot(range(10), np.cumsum(pca.explained_variance_ratio_), 'r-o')
plt.xlabel("Principal Component")
plt.ylabel("Explained Variance Ratio")
plt.title("Scree Plot")
plt.show()

# 2. Loading Plot for PC1
loadings_pc1 = pca.components_[0, :]
plt.figure()
plt.bar(range(len(loadings_pc1)), loadings_pc1)
plt.xlabel("Original Feature Index")
plt.ylabel("Loading Value")
plt.title("Loadings for PC1")
plt.show()

# 3. Biplot (conceptual)
# Requires more complex code to scale arrows, but the concept is:
plt.figure()
plt.scatter(X_pca[:, 0], X_pca[:, 1]) # Plot scores
for i, feature in enumerate(feature_names):
    plt.arrow(0, 0, pca.components_[0, i]*5, pca.components_[1, i]*5, color='r', head_width=0.1)
    plt.text(pca.components_[0, i]*5.5, pca.components_[1, i]*5.5, feature, color='r')
plt.xlabel("PC1")
plt.ylabel("PC2")
plt.title("Biplot")
plt.grid()
plt.show()

2.3 PCA Whitening

An optional but useful transformation in PCA is whitening. After projecting the data onto the principal components, this step rescales each component to have unit variance. The resulting data matrix has a covariance matrix that is the identity matrix, meaning the features are uncorrelated and all have the same variance. This can be a useful preprocessing step for other algorithms and is common in signal processing for decorrelating signals.

3. Advanced PCA: Kernel & Robust PCA

While powerful, standard PCA has two key limitations: it's linear and it's sensitive to outliers. Advanced PCA variants address these issues.

3.1 Kernel PCA

When the variance in data follows a non-linear structure (e.g., curves, spirals), standard PCA fails to capture it effectively. Kernel PCA overcomes this using the "kernel trick." Instead of explicitly mapping the data with a non-linear function, it replaces the dot products in the standard PCA algorithm with a non-linear kernel function (e.g., RBF, polynomial). This implicitly maps the data into a higher-dimensional feature space where, hopefully, the structure becomes linearly separable.

Choosing Kernel PCA Parameters

The performance of Kernel PCA is critically dependent on the choice of kernel and its parameters.

The gamma (\(\gamma\)) parameter for the RBF kernel is crucial. It controls the width of the kernel. If \(\gamma\) is too large, the model can overfit, treating each point as its own cluster. If it's too small, the model underfits and fails to capture the non-linear structure. A common heuristic for a starting value is 1 / (n_features * X.var()).

3.2 Robust PCA

Because it tries to maximize variance, standard PCA is extremely sensitive to large outliers. A single corrupted measurement can completely skew the direction of the principal components. Robust PCA addresses this by decomposing the data matrix \(X\) into two separate matrices: a low-rank matrix \(L\), representing the underlying structure, and a sparse matrix \(S\), representing the outliers or sparse noise.

\[ X = L + S \]

This decomposition is achieved by solving the optimization problem: \(\min_{L,S} \|L\|_* + \lambda\|S\|_1\), where \(\|L\|_*\) is the nuclear norm (sum of singular values) and \(\|S\|_1\) is the L1-norm (sum of absolute values of elements). This is like separating the background (low-rank) from the foreground (sparse).

Practical Applications and Parameters

Robust PCA is highly effective for tasks like removing background in video surveillance or removing instrumental glitches, like cosmic ray hits, from a set of spectra. The regularization parameter \(\lambda\) balances the trade-off between the low-rank and sparse components. A good rule of thumb for a starting value is \(\lambda = 1 / \sqrt{\max(n, d)}\), where \(n\) and \(d\) are the dimensions of the data matrix.

# Pseudocode for Robust PCA using the r_pca library
from r_pca import R_pca
import numpy as np

# X is our n x d data matrix (e.g., video frames or spectra)
# Set the lambda parameter
lambda_param = 1 / np.sqrt(max(X.shape))
rpca = R_pca(X, lmbda=lambda_param)
L, S = rpca.fit(max_iter=10000, iter_print=100)

# L is the clean background / underlying structure
# S is the sparse outliers / moving objects

4. Non-linear Manifold Learning

While Kernel PCA can find non-linear structures, a broader class of algorithms known as manifold learning is specifically designed for this purpose, primarily for visualization. The core idea is the "manifold hypothesis": that high-dimensional data often lies on or near a much lower-dimensional, non-linear manifold embedded within the high-dimensional space. These algorithms try to "unroll" this manifold to see its true geometry.

4.1 A Tour of Manifold Algorithms

Different algorithms make different assumptions about the manifold's structure.

• Isomap (Isometric Mapping): A global approach that approximates the geodesic distance (the shortest path along the manifold) between points and tries to preserve these distances in the low-dimensional embedding. It's good for simple, well-behaved global structures.
• LLE (Locally Linear Embedding): A local approach that assumes the manifold is locally linear. Each point is reconstructed as a linear combination of its neighbors, and LLE tries to preserve these reconstruction weights in the low-dimensional space.
• t-SNE (t-Distributed Stochastic Neighbor Embedding): Focuses intensely on preserving local neighborhood structure. It converts high-dimensional Euclidean distances into conditional probabilities representing similarities. It excels at visualizing well-separated clusters but often shatters the global structure of the manifold.
• UMAP (Uniform Manifold Approximation and Projection): A newer algorithm based on topological data analysis. It builds a fuzzy topological representation of the data and then optimizes a low-dimensional embedding to be as structurally similar as possible. It is often much faster than t-SNE and provides a better balance between preserving local detail and global structure.
• Others: Methods like Diffusion Maps (based on random walks on the data) and TriMAP also offer different trade-offs and are useful in specific contexts.

4.2 The Impact of Hyperparameters

The output of these algorithms is highly sensitive to their hyperparameters. Understanding them is key to producing meaningful visualizations.

• Perplexity (t-SNE): Loosely related to the number of nearest neighbors considered for each point. Typical values are between 5 and 50. Low perplexity focuses on very local structure, while high perplexity can reveal larger-scale patterns.
• n_neighbors (UMAP): Similar to perplexity, this determines the size of the local neighborhood used to learn the manifold structure. A small value will focus on very local structure, potentially missing the bigger picture, while a large value will push UMAP to focus more on the global structure.
• min_dist (UMAP): This parameter controls how tightly packed points are allowed to be in the low-dimensional embedding. Low values are good for clustering, creating tight clumps, while high values will result in a more uniform, spread-out embedding.

4.3 Evaluating Embedding Quality

Since these are visualization tools, visual inspection is the primary evaluation method. However, quantitative metrics can help assess how well the embedding preserves the original data's structure.

• Trustworthiness: Measures to what extent the local structure is preserved. It asks: "Of the points that are neighbors in the low-D embedding, how many were also neighbors in the high-D space?" A score close to 1 is good.
• Continuity: The converse of trustworthiness. It asks: "Of the points that were neighbors in the high-D space, how many are still neighbors in the low-D embedding?"

Because these algorithms often have a random component, it's good practice to assess their stability by running them multiple times with different random seeds and evaluating the average quality and variance of the results.

5. Clustering Algorithms

Clustering algorithms seek to partition data into groups, or "clusters," where samples within a group are more similar to each other than to those in other groups. The choice of algorithm is critical as it imposes a certain structure on the data. No single algorithm is best for all datasets.

5.1 A Map of Clustering Algorithms

We can broadly categorize clustering methods based on their underlying principles:

• Prototype-based (Centroid-based): These algorithms define clusters by a central prototype.
  • Examples: k-Means, k-Medoids.
• Density-based: These algorithms define clusters as dense regions of points separated by low-density regions. They are excellent for finding non-spherical clusters and handling noise.
  • Examples: DBSCAN, OPTICS, HDBSCAN.
• Hierarchical: These algorithms build a tree-like hierarchy of clusters, which can be either agglomerative (bottom-up) or divisive (top-down).
  • Example: Agglomerative Clustering.
• Graph-based: These algorithms model the data as a graph, where nodes are samples and edges represent similarity, and then partition the graph.
  • Example: Spectral Clustering.

5.2 Prototype-Based Clustering

k-Means

The goal of k-Means is to partition \(n\) observations into \(k\) clusters in which each observation belongs to the cluster with the nearest mean (cluster centroid). It assumes clusters are convex and spherical, which is a major limitation in scientific data.

5.3 Hierarchical Clustering

Agglomerative Clustering

This method starts with each point as its own cluster and progressively merges the closest pair of clusters until only one cluster (or a specified number) remains. The result is a dendrogram, a tree diagram that visualizes the merge hierarchy.

Interpreting a Dendrogram: The y-axis of the dendrogram represents the distance or dissimilarity between clusters. To obtain a specific number of clusters, you "cut" the dendrogram horizontally. The number of vertical lines your horizontal cut intersects is the number of clusters you get.

5.4 Density-Based Clustering

These methods are often a superior choice for complex scientific data as they make fewer assumptions about cluster shape.

DBSCAN (Density-Based Spatial Clustering of Applications with Noise)

DBSCAN groups together points that are closely packed, marking as outliers points that lie alone in low-density regions. It requires two key parameters:

• eps (\(\epsilon\)): The maximum distance between two samples for one to be considered as in the neighborhood of the other.
• min_samples: The number of samples in a neighborhood for a point to be considered as a core point.

Parameter Selection: Choosing eps is critical. A common heuristic is to use a k-distance plot. For a chosen min_samples (e.g., 5), calculate the distance to the 5th nearest neighbor for every point. Sort these distances and plot them. The "elbow" of this curve is a good estimate for an optimal eps value.

HDBSCAN and OPTICS

HDBSCAN improves upon DBSCAN by converting it into a hierarchical clustering algorithm, effectively allowing for variable-density clusters and eliminating the need to choose eps. OPTICS is another variant that can also find clusters in data of varying density. For most applications, HDBSCAN is a powerful and easy-to-use starting point.

5.5 Graph-Based Clustering

Spectral Clustering

This technique uses the connectivity information of the data to perform dimensionality reduction before clustering in a lower-dimensional space. It constructs a similarity graph and uses the eigenvectors of the graph Laplacian to create embeddings. K-Means is then applied to these embeddings. It is excellent for non-convex clusters like the "two moons" example but can be computationally expensive for large datasets.

# Pseudocode for key clustering algorithms
from sklearn.cluster import KMeans, DBSCAN, AgglomerativeClustering, SpectralClustering

# k-Means
kmeans = KMeans(n_clusters=3, random_state=0)
labels_kmeans = kmeans.fit_predict(X)

# Agglomerative Clustering
agg = AgglomerativeClustering(n_clusters=3, linkage='ward')
labels_agg = agg.fit_predict(X)

# DBSCAN
# Find eps from k-distance plot first
dbscan = DBSCAN(eps=0.5, min_samples=5)
labels_dbscan = dbscan.fit_predict(X)

# Spectral Clustering
spectral = SpectralClustering(n_clusters=2, affinity='nearest_neighbors', random_state=0)
labels_spectral = spectral.fit_predict(X)

6. Probabilistic Clustering: Gaussian Mixture Models

A Gaussian Mixture Model (GMM) is a powerful probabilistic algorithm that assumes the data is generated from a mixture of a finite number of Gaussian distributions (i.e., "elliptical clouds") with unknown parameters. It's a generalization of k-Means that offers more flexibility.

6.1 Core Concepts

• Soft Clustering: Unlike k-Means which assigns each point to a single cluster ("hard assignment"), a GMM provides a probability of each point belonging to each of the Gaussian components ("soft assignment"). This is useful for understanding uncertainty in the data.
• Expectation-Maximization (EM): GMMs are fit using the EM algorithm, which iteratively refines the parameters (mean, covariance, weight) of each Gaussian component to best fit the data.
• Flexible Cluster Shapes: By controlling the covariance matrix of each Gaussian, GMMs can model clusters that are not just spherical but also elliptical and oriented in different directions.

6.2 Model Selection: How many clusters?

The optimal number of Gaussian components (clusters) is often chosen by fitting GMMs with a range of components and selecting the one that minimizes an information criterion. These criteria balance model fit with model complexity.

• Akaike Information Criterion (AIC): A good general-purpose metric.
• Bayesian Information Criterion (BIC): Tends to prefer simpler models (fewer components) than AIC, as it penalizes complexity more heavily.

7. Cluster Validation & Stability

In unsupervised learning, there is no "ground truth" label to measure against. Therefore, we rely on internal validation metrics to assess the quality of the clustering results. A good clustering result is one where the clusters are dense (points are close to each other) and well-separated (points are far from other clusters).

7.1 Key Validation Metrics

No single metric is perfect; it's best to use a combination to get a complete picture.

7.2 The Importance of Stability Analysis

Perhaps more important than any single metric is the stability of the clustering solution. If the clusters are real, they should be found consistently even if the data is slightly perturbed.

Bootstrapping for Stability:

1. Create multiple "bootstrap" samples of your data by sampling with replacement.
2. Run your clustering algorithm on each bootstrap sample.
3. Compare the clustering results between the original data and each bootstrap sample using a metric like the Adjusted Rand Index (ARI). ARI measures the similarity between two clusterings, where 1 is a perfect match.

A solution that is consistently found across many bootstrap samples (e.g., average ARI > 0.85) is highly likely to be a meaningful and robust feature of the data.

8. The Full Unsupervised Workflow: A Decision Guide

The combination of these techniques forms a powerful, standard workflow for scientific data exploration. It's not a rigid pipeline but a series of decisions guided by your data and your scientific goals.

The Iterative Cycle of Discovery

This workflow is rarely linear. The final step, "Interpretation," is the most critical. If the discovered clusters don't make scientific sense or are not stable, you must loop back. This might involve:

• Revisiting preprocessing (e.g., trying a different scaling method).
• Choosing a different dimensionality reduction technique.
• Trying a more appropriate clustering algorithm for the observed data geometry.

This iterative process of modeling, validating, and interpreting is the core engine of data-driven scientific discovery.

9. Lab: Automatic Grouping of CV Curves

Goal: We have a large, unlabeled dataset of Cyclic Voltammograms. Can we automatically group them by their underlying electrochemical mechanism without any prior knowledge?

10. Practical Tips & Pitfalls

• Standardize Features: Always standardize features before PCA or distance-based clustering. It's a non-negotiable step.
• k-Means Assumptions: K-Means assumes clusters are convex and isotropic (spherical). If your data violates this (e.g., elongated or non-linear shapes), the results will be poor. Use Silhouette plots and visual inspection to verify.
• DBSCAN Parameter Tuning: The `eps` parameter in DBSCAN is sensitive. A good heuristic is to plot the k-distance graph (distance to the k-th nearest neighbor, where k = `min_samples`) and look for the "elbow" to find a suitable `eps` value.
• Interpreting PCA Loadings: The real scientific insight from PCA comes from plotting the loadings for each principal component. This tells you which original features (e.g., potential ranges) contribute most to the variance and helps you give a physical meaning to the axes.

11. Key Takeaways

• Unsupervised learning is a powerful tool for exploratory data analysis, hypothesis generation, and discovering hidden structure when labels are unavailable.
• PCA is an essential technique for compressing high-dimensional scientific data, making complex patterns easier to visualize and analyze.
• Clustering algorithms automatically find natural groupings in data, which can correspond to different physical mechanisms, sample types, or experimental outcomes.
• The combination of PCA for dimensionality reduction followed by clustering is a robust and highly effective workflow for modern scientific data analysis.