Unsupervised Learning

Learn structure from unlabeled data. These methods are building blocks for EDA, pretraining, recommenders, and anomaly detection.

Table of contents

1. Taxonomy of unsupervised learning
2. Clustering (additional methods)
3. PCA (full derivation)
4. SVD and its applications
5. Linear Discriminant Analysis (LDA)
6. Independent Component Analysis (ICA)
7. t-SNE (visualization)
8. UMAP
9. Matrix factorization for recommenders
10. Anomaly detection

1. Taxonomy of unsupervised learning

2. Clustering (additional methods)

2.1 Hierarchical clustering

Builds a dendrogram (tree of merges/splits).

Agglomerative (bottom-up):

1. Start: each point is its own cluster.
2. At each step: merge the two closest clusters.
3. Stop when desired number of clusters reached (cut dendrogram).

Linkage criteria (how to measure distance between clusters):

Divisive (top-down): start with one cluster, recursively split. Less common.

Advantage: no need to specify $k$ in advance. Dendrogram shows all possible clusterings. Disadvantage: $O(n^2 \log n)$ or $O(n^3)$ depending on implementation. Not scalable.

from sklearn.cluster import AgglomerativeClustering
from scipy.cluster.hierarchy import dendrogram, linkage
import matplotlib.pyplot as plt

Z = linkage(X, method='ward')
plt.figure(figsize=(10,4))
dendrogram(Z, truncate_mode='lastp', p=30)
plt.show()

model = AgglomerativeClustering(n_clusters=4, linkage='ward')
labels = model.fit_predict(X)

2.2 DBSCAN (Density-Based Spatial Clustering)

Idea: clusters are dense regions separated by low-density regions.

Parameters: eps (neighborhood radius), min_samples (minimum neighbors to be a core point).

Definitions:

• Core point: $\geq$ min_samples points within radius eps.
• Border point: within eps of a core point but not core itself.
• Noise point: not within eps of any core point.

Algorithm:

1. For each unvisited point:

• If core point: start new cluster, expand to all density-reachable points.
• Otherwise: mark as noise (may later become border).

Advantages:

• Finds arbitrarily shaped clusters.
• Automatically identifies outliers (noise points).
• No need to specify $k$.

Disadvantages:

• Two parameters (eps, min_samples) to tune.
• Struggles when clusters have varying densities.
• Poor in high dimensions (all distances become similar).

from sklearn.cluster import DBSCAN
db = DBSCAN(eps=0.5, min_samples=5)
labels = db.fit_predict(X)
# -1 = noise/outlier
n_clusters = len(set(labels)) - (1 if -1 in labels else 0)

2.3 Cluster evaluation metrics (no ground truth)

Silhouette coefficient:

where $a(i)$ = mean distance to same-cluster, $b(i)$ = mean distance to nearest other cluster. Range: $[-1, 1]$. Higher = better. Average over all points.

Davies-Bouldin Index: lower = better. Ratio of within-cluster to between-cluster distances.

Calinski-Harabasz Index: ratio of between-cluster to within-cluster variance. Higher = better.

from sklearn.metrics import silhouette_score, davies_bouldin_score, calinski_harabasz_score

3. PCA (full derivation)

3.1 Goal

Find orthonormal directions $\mathbf{v}_1, \mathbf{v}_2, \ldots$ in $\mathbb{R}^d$ such that projected data $z_i = \mathbf{v}^\top\mathbf{x}_i$ has maximum variance.

Center data: $\mathbf{x}_i \leftarrow \mathbf{x}_i - \bar{\mathbf{x}}$.

3.2 First principal component

where $\mathbf{S} = \frac{1}{n-1}\mathbf{X}^\top\mathbf{X}$ is the sample covariance.

Lagrangian: $\mathcal{L} = \mathbf{v}^\top\mathbf{S}\mathbf{v} - \lambda(\mathbf{v}^\top\mathbf{v} - 1)$.

Stationarity: $\nabla_\mathbf{v}\mathcal{L} = 2\mathbf{S}\mathbf{v} - 2\lambda\mathbf{v} = \mathbf{0} \Rightarrow \mathbf{S}\mathbf{v} = \lambda\mathbf{v}$.

So $\mathbf{v}_1$ is an eigenvector of $\mathbf{S}$. The objective value is $\mathbf{v}^\top\mathbf{S}\mathbf{v} = \lambda$. Maximize → choose eigenvector for largest eigenvalue.

3.3 Subsequent components

$k$-th PC: eigenvector of $\mathbf{S}$ for $k$-th largest eigenvalue, subject to orthogonality with $\mathbf{v}_1,\ldots,\mathbf{v}_{k-1}$.

Spectral decomposition: $\mathbf{S} = \mathbf{V}\boldsymbol{\Lambda}\mathbf{V}^\top$, $\mathbf{V} = [\mathbf{v}_1,\ldots,\mathbf{v}_d]$, $\boldsymbol{\Lambda} = \text{diag}(\lambda_1 \geq \cdots \geq \lambda_d \geq 0)$.

3.4 Projection and reconstruction

Encode (reduce to $k$ dims): $\mathbf{z}_i = \mathbf{V}_k^\top\mathbf{x}_i \in \mathbb{R}^k$. Decode (reconstruct): $\hat{\mathbf{x}}_i = \mathbf{V}_k\mathbf{z}_i = \mathbf{V}_k\mathbf{V}_k^\top\mathbf{x}_i$. Reconstruction error: $\|\mathbf{x}_i - \hat{\mathbf{x}}_i\|^2 = \sum_{j>k}\lambda_j$. Minimized by choosing top $k$ eigenvectors (Eckart-Young for covariance matrix).

3.5 Explained variance ratio

Cumulative EVR: choose $k$ to retain 95% (or 99%) of variance.

3.6 PCA via SVD (computational approach)

Directly decompose data matrix (more numerically stable than computing covariance):

Then: eigenvectors of $\mathbf{S}$ = right singular vectors $\mathbf{V}$, eigenvalues of $\mathbf{S}$ = $\sigma_i^2$.

sklearn uses SVD by default (no need to form $\mathbf{X}^\top\mathbf{X}$).

3.7 When PCA helps (and doesn't)

Helps: correlated features, dimensionality reduction for downstream model, visualization, noise reduction.

Doesn't help: if relationship to target is nonlinear with respect to principal components; if features are already uncorrelated; if you need interpretable features (PCs are linear combinations of originals).

Supervised alternative: see Linear Discriminant Analysis (§5).

from sklearn.decomposition import PCA

pca = PCA(n_components=0.95)  # retain 95% variance
X_reduced = pca.fit_transform(X_train)
print(pca.n_components_)
print(pca.explained_variance_ratio_.cumsum())

# Visualize 2D
pca2 = PCA(n_components=2)
X_2d = pca2.fit_transform(X)

4. SVD and its applications

4.1 Full decomposition

$\mathbf{A} = \mathbf{U}\boldsymbol{\Sigma}\mathbf{V}^\top$ (see Note 01 for full details).

Truncated SVD (rank-$k$ approximation):

Minimizes $\|\mathbf{A}-\mathbf{B}\|_F$ over all rank-$k$ matrices $\mathbf{B}$ (Eckart-Young theorem).

4.2 Applications

5. Linear Discriminant Analysis (LDA)

5.1 Goal

Supervised dimensionality reduction: find projection that maximizes class separation.

Two scatter matrices:

• Within-class scatter: $\mathbf{S}_W = \sum_k \sum_{i:y_i=k}(\mathbf{x}_i-\boldsymbol\mu_k)(\mathbf{x}_i-\boldsymbol\mu_k)^\top$.
• Between-class scatter: $\mathbf{S}_B = \sum_k n_k(\boldsymbol\mu_k-\boldsymbol\mu)(\boldsymbol\mu_k-\boldsymbol\mu)^\top$.

5.2 Fisher criterion

Find projection $\mathbf{W}$ that maximizes:

Solution: generalized eigenvalue problem $\mathbf{S}_W^{-1}\mathbf{S}_B\mathbf{w} = \lambda\mathbf{w}$.

Projection onto top $\min(K-1, d)$ eigenvectors.

LDA as generative classifier: assumes each class is Gaussian with shared covariance $\boldsymbol\Sigma$. Under this model, Bayes-optimal boundary is linear.

QDA (Quadratic Discriminant Analysis): allows class-specific covariances → quadratic boundary.

6. Independent Component Analysis (ICA)

6.1 Problem (cocktail party)

Observe $\mathbf{x} = \mathbf{A}\mathbf{s}$, where $\mathbf{s}$ are independent non-Gaussian sources and $\mathbf{A}$ is unknown mixing matrix. Recover $\mathbf{s} = \mathbf{W}\mathbf{x}$.

Key difference from PCA: PCA finds uncorrelated components (second-order statistics); ICA finds statistically independent components (higher-order statistics).

Why non-Gaussian? Central limit theorem: sums of independent variables are more Gaussian. So a Gaussian component is maximally "mixed" — unmixing must find the most non-Gaussian projections.

6.2 FastICA

Maximize non-Gaussianity of $\mathbf{w}^\top\mathbf{x}$, measured by:

• Kurtosis: $\mathbb{E}[z^4] - 3(\mathbb{E}[z^2])^2$. Zero for Gaussian; non-zero for super/sub-Gaussian.
• Negentropy: $J(y) = H(y_\text{Gaussian}) - H(y)$. Always $\geq 0$, maximized for non-Gaussian.

Fixed-point iteration (FastICA) converges cubically (faster than gradient methods).

from sklearn.decomposition import FastICA
ica = FastICA(n_components=3, random_state=0)
S_estimated = ica.fit_transform(X)

Limitations: cannot determine order or scale of components; sources must be non-Gaussian (except at most one can be Gaussian).

7. t-SNE (visualization)

7.1 Purpose

t-distributed Stochastic Neighbor Embedding (van der Maaten & Hinton, 2008): map high-dimensional data to 2D/3D for visualization, preserving local structure (nearby points in high-D → nearby in low-D).

7.2 Algorithm

High-dimensional affinities: for each pair, define conditional probability:

$\sigma_i$ is set per-point based on perplexity (target number of effective neighbors).

Symmetrize: $p_{ij} = (p_{j|i}+p_{i|j})/(2n)$.

Low-dimensional affinities (heavy-tailed Student-t):

Why Student-t? Alleviates crowding problem: in high-D, moderate distances become equal; t-distribution has heavier tail to spread points in 2D.

Objective (minimize KL divergence):

Optimized via gradient descent on low-D positions $\mathbf{y}_i$.

7.3 Critical notes about t-SNE

• Not for general dimensionality reduction (distances not preserved globally; non-deterministic; not invertible).
• Perplexity (5–50): controls neighborhood size. Try multiple values.
• Cluster sizes and distances between clusters are not meaningful in t-SNE plots.
• Cannot be applied to new data (need to re-run on full dataset).
• Non-convex; results depend on random seed.

from sklearn.manifold import TSNE
tsne = TSNE(n_components=2, perplexity=30, learning_rate='auto', init='pca', random_state=42)
X_2d = tsne.fit_transform(X)

8. UMAP

8.1 Key differences from t-SNE

8.2 Brief theory

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

1. Construct a fuzzy topological representation of the high-dimensional data using k-nearest neighbors and Riemannian metric (local radius normalized).
2. Find a low-dimensional embedding that preserves this fuzzy structure, optimizing cross-entropy:

where $q_{ij} = (1+a\|\mathbf{y}_i-\mathbf{y}_j\|^{2b})^{-1}$ (learnable parameters $a,b$).

import umap
reducer = umap.UMAP(n_components=2, n_neighbors=15, min_dist=0.1, metric='euclidean', random_state=42)
X_2d = reducer.fit_transform(X)
# Transform new points:
X_new_2d = reducer.transform(X_new)

9. Matrix factorization for recommenders

9.1 Problem setup

Explicit feedback: user-item rating matrix $\mathbf{R} \in \mathbb{R}^{U\times I}$, mostly missing.

Goal: fill in missing entries to predict ratings/preferences.

9.2 Matrix factorization (SVD-based)

Approximate: $\mathbf{R} \approx \mathbf{P}\mathbf{Q}^\top$, where:

• $\mathbf{P} \in \mathbb{R}^{U\times k}$: user embeddings.
• $\mathbf{Q} \in \mathbb{R}^{I\times k}$: item embeddings.
• $k$ = latent factors (dimensionality).

Objective (observed entries only):

Alternating Least Squares (ALS): fix $\mathbf{Q}$, solve for each $\mathbf{p}_u$ (closed form); fix $\mathbf{P}$, solve for each $\mathbf{q}_i$. Parallelizable.

SGD: sample observed ratings, gradient update for $\mathbf{p}_u, \mathbf{q}_i$.

9.3 Non-negative Matrix Factorization (NMF)

$\mathbf{V} \approx \mathbf{W}\mathbf{H}$ with $\mathbf{W},\mathbf{H} \geq 0$ (element-wise).

Why non-negative? Parts-based representation: in images, each component corresponds to a part (not the signed PCA components). In NLP, topics are non-negative weighted combinations of words.

Multiplicative updates (Lee & Seung):

from sklearn.decomposition import NMF
model = NMF(n_components=10, init='nndsvda', max_iter=300, random_state=0)
W = model.fit_transform(X)  # user/document representation
H = model.components_       # topic/feature representation

9.4 Implicit feedback

Most real systems have implicit feedback (clicks, views, purchases) rather than explicit ratings.

Weighting: treat all interactions as positive; unobserved as potentially negative with lower confidence.

ALS for implicit (Hu et al., 2008): define confidence $c_{ui} = 1 + \alpha r_{ui}$, minimize:

where $p_{ui} = 1$ if interaction, 0 otherwise. implicit library implements this.

10. Anomaly detection

10.1 Problem types

• Novelty detection: training data is clean; find new data that differs from it.
• Outlier detection: training data may contain outliers; find them.

10.2 Isolation Forest

Idea: anomalies are rare and different → isolated quickly by random splits.

Build random trees that recursively split on random features at random thresholds. Anomaly score = average depth at which point is isolated (shorter path = more anomalous).

where $h(\mathbf{x})$ = path length, $c(n) = 2H(n-1) - 2(n-1)/n$ = expected path length for $n$ samples.

from sklearn.ensemble import IsolationForest
iso = IsolationForest(contamination=0.05, random_state=0)
labels = iso.fit_predict(X)  # -1 = outlier, 1 = inlier
scores = iso.score_samples(X)  # lower = more anomalous

10.3 Local Outlier Factor (LOF)

Compare local density of point to its neighbors. Anomalies are in regions of much lower density than neighbors.

where $\text{lrd}_k(i)$ = local reachability density.

LOF > 1 = outlier (lower density than neighbors).

10.4 One-class SVM

Fit a hypersphere (or hyperplane in kernel space) around the training data. Points outside are anomalies.

10.5 Statistical methods

Z-score: flag if $|x-\mu|/\sigma > 3$. Assumes Gaussian. IQR: flag if outside $[Q1-1.5\text{IQR}, Q3+1.5\text{IQR}]$. Non-parametric. Mahalanobis distance: multivariate generalization of Z-score; accounts for correlations.

*File: notes/09_unsupervised_learning.md — next: notes/10_practical_ml.md*