Clustering and PCA

A developer-support team is preparing semantic search over billing, login, and deploy-policy messages. The embedding model has produced vectors, but nobody has reviewed their neighborhoods. When someone searches for “I was charged twice,” will the nearest message be another billing case or an unrelated deploy note?

Without reviewed query-document pairs, we can't yet score retrieval relevance. We can start by inspecting the vectors: k-means groups nearby points, while principal component analysis (PCA) gives them fewer coordinates for inspection or compression. Neither method reads the messages. This chapter builds their core update rules from scratch, connects them to coarse vector search indexing, and tests whether changing units, similarity metrics, or dimensions breaks a useful neighbor relationship.

Start with six messages, not six labels

We'll use six invented two-number embeddings, chosen to make the arithmetic readable. Real embeddings often contain hundreds or thousands of coordinates; here we can see every coordinate and check every distance. The examples use NumPy and assume you're comfortable with arrays and coordinate-wise averages.

Before looking for groups, identify the pairs you'd expect from reading the messages. We can then check whether distance agrees with that judgment.

Names in the first column serve our audit rather than feeding into the algorithm. In practice, the clustering logic receives only six coordinate pairs without ever seeing billing, login, or deploy.

Looking at the billing pair, the points differ by 0.2 horizontally and 0.8 vertically. Straight-line distance between them is $\sqrt{0.2^2+0.8^2}=\sqrt{0.68}\approx0.82$. For any two points $a=(a_x,a_y)$ and $b=(b_x,b_y)$, that Euclidean distance is:

$$d(a,b)=\sqrt{(a_x-b_x)^2 + (a_y-b_y)^2}$$

Evaluating this formula measures separation in the representation we currently have. Run it before clustering. A sensible local neighborhood gives the later groups something concrete to explain.

Our code computes all 36 pairwise distances. Broadcasting creates a (6, 6, 2) array of coordinate differences; axis=2 reduces each coordinate pair to one distance. Setting the diagonal to infinity prevents each message from selecting itself. This all-pairs approach is for a small audit, not a large search index.

Each message's closest neighbor belongs to the same readable theme. That's useful evidence about this representation, not a promise about future messages. Keep that distinction in view as we move from pairs to groups.

Partition unlabeled vectors with Lloyd's algorithm

Search doesn't require clusters. Here, clustering helps us select groups of messages to review together. On these six unlabeled vectors, k-means receives the coordinates and a requested number of groups, $k$. A centroid is a cluster's center, calculated by averaging its points.

The algorithm minimizes the within-cluster sum of squares (WCSS), commonly referred to as inertia:^{[1]Reference 1The Elements of Statistical Learning.https://hastie.su.domains/ElemStatLearn/[2]Reference 2Pattern Recognition and Machine Learning.https://www.microsoft.com/en-us/research/publication/pattern-recognition-machine-learning/}

$$J = \sum_{k=1}^{K} \sum_{i \in S_k} \left\lVert x_i - \mu_k \right\rVert_2^2$$

Here, $x_i$ is message vector $i$, $S_k$ is the set of points assigned to cluster $k$, and $\mu_k$ is the centroid of cluster $k$. Lloyd's algorithm optimizes this objective through coordinate descent, alternating between two steps until convergence:^{[3]Reference 3Clusteringhttps://scikit-learn.org/stable/modules/clustering.html}

1. Assignment step: Hold the centroids $\mu_k$ fixed. Assign each vector $x_i$ to its nearest centroid:

$$z_i = \arg\min_{k \in \{1,\dots,K\}} \left\lVert x_i - \mu_k \right\rVert_2^2$$

This minimizes $J$ with respect to the assignments $z$ while keeping centroids stationary.

1. Update step: Hold the assignments $S_k$ fixed. Recompute each centroid $\mu_k$ to minimize $J$:

$$\frac{\partial J}{\partial \mu_k} = -2 \sum_{i \in S_k} (x_i - \mu_k) = 0 \implies \mu_k = \frac{1}{|S_k|} \sum_{i \in S_k} x_i$$

Taking the arithmetic mean isn't an arbitrary choice. The derivative shows that the sample mean is the unique point that minimizes the sum of squared Euclidean distances to all members of the set.

Because both steps strictly decrease or preserve $J$, the sequence of objective values is monotonically non-increasing. Since there are only finitely many ways to partition $n$ points into $k$ clusters, Lloyd's algorithm is guaranteed to terminate. It settles on a local minimum, not necessarily the global optimum.

Start with $k=3$ and pick one seed from each visible area:

• $c_0=(1.0,1.0)$ from billing_A
• $c_1=(4.0,4.2)$ from login_A
• $c_2=(8.0,1.0)$ from deploy_A

Before any update, the initial inertia evaluated at these seeds is:

$$J_{\text{seed}} = (0^2 + 0^2) + (0.2^2 + 0.8^2) + (0^2 + 0^2) + (1.0^2 + (-0.4)^2) + (0^2 + 0^2) + (1.0^2 + 0.8^2) = 0.68 + 1.16 + 1.64 = 3.48$$

After assignment, each pair stays with its local seed. Averaging the coordinates in each pair produces:

Evaluating the objective around these updated means cuts inertia in half:

$$J_{\text{updated}} = 2 \cdot (0.1^2 + 0.4^2) + 2 \cdot (0.5^2 + 0.2^2) + 2 \cdot (0.5^2 + 0.4^2) = 0.34 + 0.58 + 0.82 = 1.74$$

Why can't you write cluster 0 = billing into a long-lived dashboard?

Answer

Cluster integers have no stable semantic meaning. Another initialization or retraining run can number the same groups differently. Persist a run ID, inspect current members, and attach a reviewed human label separately.

Build the k-means loop and prevent initialization traps

Production libraries provide careful seeding and convergence controls. Our scratch implementation accepts a matrix of vectors and initial centers. Iteration stops when cluster assignments repeat, which guarantees that the calculated means remain unchanged too. Halting on stable assignments avoids relying on arbitrary coordinate tolerances that might stop prematurely.

The initialization trap: why k-means++ matters

Uniform random initialization frequently picks multiple seeds that sit within the same dense cluster. When that happens, two seeds split a single coherent topic while another distant region is left without any seed at all.

Arthur and Vassilvitskii addressed this vulnerability with k-means++ initialization.^{[3]Reference 3Clusteringhttps://scikit-learn.org/stable/modules/clustering.html} Rather than sampling centroids uniformly, k-means++ spreads initial seeds across the feature space:

1. Select the first centroid $\mu_1$ uniformly at random from the dataset $X$.
2. For each point $x \in X$, compute $D(x)$, the shortest Euclidean distance from $x$ to any already chosen centroid:

$$D(x) = \min_{j \in \{1,\dots,m\}} \|x - \mu_j\|_2$$

1. Sample the next centroid from $X$ with probability proportional to the squared distance:

$$P(x) = \frac{D(x)^2}{\sum_{x' \in X} D(x')^2}$$

1. Repeat steps 2 and 3 until all $K$ centroids have been selected.

Points close to existing centroids have small $D(x)^2$, making duplicate seeding in the same neighborhood improbable. Points far away have large $D(x)^2$, drawing the next centroid toward unexplored clusters. This simple stochastic heuristic guarantees an expected bound of $O(\log K)$ times the optimal WCSS.

We can see this probability distribution in action on our developer messages. Suppose billing_A was chosen as the first seed:

The neighboring billing message has less than a 0.5% chance of being selected as the next seed. The distant deploy messages receive over 72% of the cumulative selection probability.

The loop also exposes a failure that a library call can hide. If a center gets no assigned vectors, mean is undefined. The guard raises instead of writing nan into the next iterate:

Don't average an empty slice. A production implementation needs an explicit repair policy, such as reseeding the empty cluster at the point farthest from any current center.

One embedding dimension ranges from 0 to 1, while another ranges from 0 to 10,000. What happens if Euclidean k-means uses them without scaling?

Answer

The large-range dimension can dominate distance and manufacture clusters around its units. Standardize or otherwise justify feature scales before interpreting the assignments.

Evaluate cluster geometry with elbow curves and silhouette scores

The toy map made $k=3$ feel natural because it contained three distinct pairs. A production corpus won't hand you its correct cluster count. Treat $k$ as an architectural decision to evaluate, not a fixed constant.

A cluster audit combines two quantitative diagnostic tools with human validation:

1. The elbow method: Plot inertia $J(k)$ as a function of $k$. Since adding more centroids grants extra degrees of freedom, the best achievable inertia can't increase as $k$ grows, reaching zero when $k=n$. Look for the inflection point (the "elbow") where the rate of improvement drops:

$$\Delta J(k) = J(k) - J(k+1)$$

1. Silhouette analysis: For each vector $x_i$, compute its intra-cluster cohesion $a(i)$ (the average distance to all other points in its assigned cluster):

$$a(i) = \frac{1}{|S_{z_i}| - 1} \sum_{j \in S_{z_i}, j \neq i} \left\lVert x_i - x_j \right\rVert_2$$

Next, find its nearest neighbor cluster and compute separation $b(i)$ (the smallest mean distance to any other cluster):

$$b(i) = \min_{k \neq z_i} \frac{1}{|S_k|} \sum_{j \in S_k} \left\lVert x_i - x_j \right\rVert_2$$

The silhouette score for point $i$ is:^{[4]Reference 4silhouette_scorehttps://scikit-learn.org/stable/modules/generated/sklearn.metrics.silhouette_score.html}

$$s(i) = \frac{b(i) - a(i)}{\max(a(i), b(i))}$$

Scores near $+1$ indicate that a point is well inside its cluster and far from neighbors. Scores near $0$ mean the point lies on a cluster boundary. Negative scores mean the point is closer to a neighboring cluster on average than to its assigned center.

Watch how adding a fourth cluster continues lowering inertia while degrading silhouette:

Moving from $k=2$ to $k=3$ creates a steep drop in inertia (from 20.06 down to 1.74) while driving silhouette up to 0.76. Moving to $k=4$ reduces inertia only marginally (from 1.74 down to 0.92) while splitting a natural pair, which pulls silhouette down to 0.51.

Equal scores can still tell different stories

An objective can be stable while the semantic interpretation remains ambiguous. Place four messages at the vertices of a square. Two orthogonal two-cluster partitions yield the exact same inertia: split left from right, or split bottom from top.

Both groupings optimize the objective equally well. If horizontal separation represented user tier and vertical separation represented issue severity, geometry alone wouldn't tell you which grouping your support team needs. Equal objective values don't make the product stories interchangeable.

Compress representations with principal component analysis

K-means assigns vectors to discrete categories. PCA addresses a complementary question: can we describe the cloud with fewer continuous coordinates while preserving as much variance as possible?

PCA has two equivalent mathematical foundations:^{[2]Reference 2Pattern Recognition and Machine Learning.https://www.microsoft.com/en-us/research/publication/pattern-recognition-machine-learning/}

1. Maximum variance formulation: Find an orthonormal unit vector $w$ ($\|w\|_2=1$) such that projecting the centered data onto $w$ maximizes the empirical variance of the projections:

$$\max_{\|w\|=1} \frac{1}{n-1} \sum_{i=1}^n (x_i^\top w)^2 = \max_{\|w\|=1} w^\top C w$$

where $C = \frac{1}{n-1} X_c^\top X_c$ is the sample covariance matrix.

1. Minimum reconstruction error formulation: Find an orthonormal basis $w$ that minimizes the sum of squared reconstruction errors when projecting points back to the original space:

$$\min_{\|w\|=1} \sum_{i=1}^n \left\lVert x_i - (x_i^\top w) w \right\rVert_2^2$$

These two views are identical because of the Pythagorean theorem for orthogonal projections. For any point $x_i$ and unit direction $w$:

$$\|x_i\|_2^2 = \|(x_i^\top w) w\|_2^2 + \|x_i - (x_i^\top w) w\|_2^2$$

Summing over all points in the centered dataset yields:

$$\text{Total Variance} = \text{Projected Variance} + \text{Reconstruction Error}$$

Because total variance is fixed for any given dataset, maximizing projected variance automatically minimizes reconstruction error.

Covariance eigendecomposition versus singular value decomposition

First center the matrix by subtracting the column-wise mean vector $\mu$:

$$X_c = X - \mathbf{1} \mu^\top$$

From $X_c$, two computational routes lead to the same principal components:^{[5]Reference 5sklearn.decomposition.PCAhttps://scikit-learn.org/stable/modules/generated/sklearn.decomposition.PCA.html}

• Route A: Eigendecomposition of the covariance matrix. Form $C = \frac{1}{n-1} X_c^\top X_c$, then solve the eigenvalue equation:

$$C v_j = \lambda_j v_j$$

The eigenvectors $v_j$ form the principal directions, and the eigenvalues $\lambda_j$ represent the variance captured along each axis.

• Route B: Singular Value Decomposition (SVD) directly on $X_c$. Factor the centered data matrix into:

$$X_c = U \Sigma V^\top$$

Plugging this factorization into the covariance definition reveals the exact equivalence:

$$X_c^\top X_c = (V \Sigma^\top U^\top)(U \Sigma V^\top) = V \Sigma^2 V^\top$$

Because $U^\top U = I$, the right singular vectors $V$ are precisely the eigenvectors of $X_c^\top X_c$, and the singular values $\sigma_j$ relate directly to the covariance eigenvalues:

$$\lambda_j = \frac{\sigma_j^2}{n-1}$$

Production libraries (such as scikit-learn and PyTorch) compute PCA via SVD rather than explicit covariance eigendecomposition for two numerical reasons:

1. Condition number squaring: Forming $X_c^\top X_c$ squares the condition number of the data matrix: $\kappa(X_c^\top X_c) = \kappa(X_c)^2$. If a feature has small singular values, squaring them can push values below machine precision, causing catastrophic numerical cancellation.
2. Memory and computation at scale: If you have 100 documents represented by 1536-dimensional embeddings ($n \ll d$), forming $X_c^\top X_c$ requires allocating and operating on a $1536 \times 1536$ covariance matrix. Truncated SVD computes the top $k$ principal components in $O(n \cdot d \cdot k)$ without ever materializing the $d \times d$ covariance matrix.

The explained variance ratio measures the proportion of total linear spread retained by each component:

$$\text{Ratio}_j = \frac{\lambda_j}{\sum_{k=1}^d \lambda_k} = \frac{\sigma_j^2}{\sum_{k=1}^d \sigma_k^2}$$

For our six developer messages, PC1 runs almost horizontally, accounting for 85.2% of the total geometric variance:

PCA produces a single coordinate, pc1, as a linear combination of original dimensions. Sign orientation is arbitrary: flipping the sign of an eigenvector reverses coordinates without affecting pairwise distances or reconstruction error.

⚠️ Common mistake: Cluster a 2D PCA plot and assume you've recovered the neighborhoods in the full embedding. PCA before clustering can help remove noise or reduce computation, but it changes the space being grouped. Compare it with clustering the original vectors and evaluate the result against the task. A clean plot alone doesn't justify the transform.

Measure what compression discards

Connecting explained variance to reconstruction error demonstrates the trade-off. Project the data onto one component, reconstruct it in the original coordinate system, and measure the squared residual:

$$\hat{X} = (X_c v_1) v_1^\top + \mathbf{1} \mu^\top$$

The boundary is concrete: explained variance measures geometry, not semantic importance. A subtle feature separating critical edge cases might carry low total variance, yet determine whether a query succeeds or fails.

Two principal components explain 95% of total variance. Does that prove they preserve retrieval quality?

Answer

No. PCA preserves high-variance directions, not task relevance. A low-variance direction may still carry the signal needed for a query, so measure downstream neighbors or retrieval metrics after projection.

Connect clustering to vector search: inverted file indexing

Clustering isn't just an exploratory diagnostic. In production vector databases (such as Faiss, Milvus, Qdrant, and pgvector), k-means forms the foundation of the Inverted File Index (IVF) for fast approximate nearest neighbor (ANN) search.^{[6]Reference 6Billion-scale similarity search with GPUs.https://arxiv.org/abs/1702.08734[7]Reference 7Faiss indexeshttps://github.com/facebookresearch/faiss/wiki/Faiss-indexes}

When serving millions of $d$-dimensional embeddings, exhaustive brute-force search (a flat index) compares the query against every stored vector, requiring $O(N \cdot d)$ operations per query. Under strict latency budgets (often under 10 milliseconds), brute-force scanning fails at scale.

IVF turns this linear search into a sublinear candidate lookup using k-means:

1. Training phase: Run k-means over the document embeddings to learn $K$ coarse centroids (often $K = 4096$ or $K \approx \sqrt{N}$), creating a Voronoi partition of the vector space.
2. Index building: Assign each document vector to its closest centroid. Instead of storing a flat array, the index maintains an inverted posting list for each centroid containing its assigned document IDs.
3. Query phase:
   • Compute distances from the incoming query vector $q$ to the $K$ coarse centroids.
   • Select the $nprobe$ closest centroids (typically $nprobe \ll K$, such as $nprobe = 16$).
   • Scan only the document vectors in the inverted posting lists belonging to those $nprobe$ cells.

This prunes $(1 - \frac{nprobe}{K})$ of the search space. If $K=4096$ and $nprobe=16$, the engine skips 99.6% of the database during retrieval.

We can simulate an IVF coarse index over our developer messages:

The $nprobe$ parameter governs the trade-off between search latency and recall:

• Low $nprobe$ ($nprobe=1$): Maximal query throughput. If a document lies close to a Voronoi boundary, the query might land in the adjacent cell, causing a boundary recall miss.
• High $nprobe$ ($nprobe \to K$): High recall approaching exact brute-force search, but at the cost of evaluating more candidates.

Failure test 1: let units manufacture the clusters

An embedding often sits beside tabular metadata such as message length, incident count, or timestamps. Mixing unlike measurements into one matrix causes k-means or PCA to track the largest raw numerical units rather than semantic meaning.

Here an invented topic coordinate separates billing and deploy, while message length varies across both. In raw units, a short billing note is closer to a short deploy note than to a long billing note. We'll run the same two-center loop twice, always seeding from bill_100 and dep_110. Only feature scaling changes.

The raw 800-character gap dwarfs the topic difference of 6: the two short messages group together, regardless of topic. With zero mean and unit variance, the loop recovers the intended topic split.

Don't standardize every embedding blindly. If a model's vectors come with a specified normalization contract, follow it. When fitting a scaler or PCA, learn its parameters on the development corpus and reuse that transform for held-out queries. As with supervised validation, fitting scalers on the evaluation set causes data leakage.

Failure test 2: let the similarity rule change neighbors

Text retrieval often compares vector directions with cosine similarity, while ordinary k-means minimizes Euclidean distance. Dot product rewards magnitude as well as angular alignment.

OpenAI's embedding models return vectors normalized to length 1.^{[8]Reference 8Vector embeddingshttps://developers.openai.com/api/docs/guides/embeddings} Under that contract, cosine similarity equals the dot product, and squared Euclidean distance is an inverted linear transform of cosine similarity:

$$\|u - v\|_2^2 = \|u\|_2^2 + \|v\|_2^2 - 2 (u \cdot v) = 2 - 2 \cos(u, v)$$

Because the relation is strictly monotonic, cosine ranking and Euclidean nearest-neighbor ranking produce identical orderings for unit vectors. Mix in unnormalized metadata or switch to an unnormalized model, and that shortcut fails.

The equivalence applies to nearest neighbors, not standard k-means centroids. Standard k-means computes the coordinate mean of assigned points. The mean of several unit vectors is generally not unit length, causing centroids to drift off the hypersphere:

If you cluster the same vectors you'll cosine-search later, align the algorithm with the serving contract:^{[9]Reference 9Faiss FAQ: k-means for the cosine similarityhttps://github.com/facebookresearch/faiss/wiki/FAQ#k-means-for-the-cosine-similarity}

• Run spherical k-means: Assign points by cosine similarity (maximum dot product on unit vectors), then project centroids back to the unit sphere ($\mu_k \leftarrow \mu_k / \|\mu_k\|_2$) after each update step.
• Or treat running Euclidean k-means on normalized vectors and normalizing centroids only at inference time as an approximation. Validate held-out recall against a spherical baseline before deploying it.

K-means also favors roughly spherical clusters of comparable size. Elongated manifolds or nested topics score poorly under inertia even when they represent meaningful semantic groupings.

Failure test 3: centering changes cosine geometry

Cosine similarity measures the angle between two vectors from the origin:

$$\text{cosine\_similarity}(u, v) = \frac{u \cdot v}{\|u\| \|v\|}$$

PCA introduces an origin trap. It centers the dataset by subtracting the column-wise mean from each vector, shifting the origin to the center of the data cloud.^{[5]Reference 5sklearn.decomposition.PCAhttps://scikit-learn.org/stable/modules/generated/sklearn.decomposition.PCA.html} If we then perform cosine retrieval on the centered vectors, we're measuring angles from an entirely different point.

Shifting the origin alters angles between vectors, so original cosine distances and neighbor rankings aren't preserved. Two vectors pointing in nearly identical directions from the original origin can point in opposite directions after centering.

Fit the mean on corpus documents, then apply it to the query:

After centering, doc2 sits more than 90 degrees from the query. The transformation altered the geometric question the angle answers. If you reduce dimensions for cosine retrieval, treat the compressed space as a new representation and validate retrieval metrics directly.

Failure test 4: PCA can hide the dimension your query needs

A retrieval system doesn't need a globally optimal low-rank approximation. It needs the correct document for a specific query. Suppose horizontal variation across billing and deploy policies dominates corpus variance, while the distinction between password-reset and general billing pages lies in a quieter vertical direction.

Project documents to one principal component, then execute nearest-document lookup for a login question:

After projection, billing_help, login_2fa, and runbook_note collapse to the exact same one-dimensional coordinate. PCA eliminated the quieter vertical dimension that separated the correct login document from tied alternatives. An argmin call now returns whichever candidate appears first in storage.

Reconstruction error and query relevance answer different questions. Measure retrieval metrics (like Recall@K and NDCG) on labeled query sets before deploying compressed embeddings.

Turn every plot into a representation audit

An unsupervised plot earns its keep when it generates testable hypotheses. Before claiming that an embedding space is organized, compile an audit record:

Use reviewed examples as an audit set. If you repeatedly tune your pipeline after inspecting them, treat them as development data, not an untouched evaluation benchmark.

On this teaching set, purity is 1.0 because the coordinates were designed to be clean. A production audit must evaluate boundary edge cases, multilingual queries, and ambiguous intents.

Practice: try to break the representation

Use the clean six-message run as a baseline, then perturb one decision at a time:

1. Empty center repair. Modify the initial seeds in kmeans-from-scratch.py so one cluster empties. Implement a reseed policy that restarts the orphaned centroid at the point having the highest current squared distance from any center.
2. K-means++ restart test. Implement the full k-means++ selection loop. Compare the variance of final inertia across 50 random restarts between uniform random seeding and k-means++ seeding.
3. Competing units. Add a metadata column (like character count) with values in $[100, 1000]$. Compare cluster assignments before and after z-score standardization.
4. IVF recall boundary. In ivf-coarse-quantization-search.py, move the query point to $[3.0, 3.0]$ (equidistant between billing and login). Measure whether $nprobe=1$ misses the true nearest neighbor and check if $nprobe=2$ recovers it.
5. Lost query signal. Retain both principal components in compression-can-change-retrieval.py. Confirm that the login query recovers login_2fa, and explain why high explained variance didn't prevent the rank-1 tie.

What to look for

• An empty mean(axis=0) produces an invalid center. Detect members.shape[0] == 0 and trigger an explicit reseed rather than propagating NaN.
• K-means++ yields consistently lower inertia and faster convergence across restarts by avoiding bunched seeds.
• When unscaled metadata dominates distance, clustering reflects units rather than semantics. Verify feature normalization before interpreting clusters.
• On Voronoi boundaries, setting $nprobe=1$ causes recall drops. Production vector indexes calibrate $nprobe$ against target Recall@K.
• PCA preserves global variance, not task labels. Evaluating downstream retrieval metrics on held-out pairs is mandatory before deploying compressed representations.