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 those update rules from scratch. It then uses the same groups as coarse cells in a vector index, and it tests whether a useful neighbor survives a change of units, similarity, or kept dimensions.

View larger

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 NumPy examples need only arrays and coordinate-wise averages to get started.

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:

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 supports a small audit; large search indexes need a scalable neighbor-search method.

import numpy as np

names = np.array(["billing_A", "billing_B", "login_A", "login_B", "deploy_A", "deploy_B"])
vectors = np.array([
    [1.0, 1.0],
    [1.2, 1.8],
    [4.0, 4.2],
    [5.0, 3.8],
    [8.0, 1.0],
    [9.0, 1.8],
])

distances = np.linalg.norm(vectors[:, None, :] - vectors[None, :, :], axis=2)
np.fill_diagonal(distances, np.inf)

for row, name in enumerate(names):
    neighbor = distances[row].argmin()
    print(f"{name:10} -> {names[neighbor]:10} distance={distances[row, neighbor]:.2f}")

billing_A  -> billing_B  distance=0.82
billing_B  -> billing_A  distance=0.82
login_A    -> login_B    distance=1.08
login_B    -> login_A    distance=1.08
deploy_A   -> deploy_B   distance=1.28
deploy_B   -> deploy_A   distance=1.28

Each message's closest neighbor belongs to the same readable theme. This sample shows the representation grouping these messages consistently. Check more messages before relying on that pattern for future data. If the pairs already agree with the text, what does a cluster add that the neighbor list does not?

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 objective is the within-cluster sum of squares (WCSS), commonly called inertia:^{[1]Reference 1Pattern Recognition and Machine Learning.https://www.microsoft.com/en-us/research/publication/pattern-recognition-machine-learning/[2]Reference 2KMeanshttps://scikit-learn.org/stable/modules/generated/sklearn.cluster.KMeans.html}

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:

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$:

For a nonempty cluster, the sample mean uniquely minimizes the sum of squared Euclidean distances to its members. An empty cluster has no defined mean; our implementation handles that case explicitly.

Each step decreases or preserves $J$. With exact updates, consistent tie handling, and nonempty clusters, Lloyd's algorithm reaches a fixed assignment in finitely many steps.^{[4]Reference 4k-means++: The Advantages of Careful Seedinghttps://theory.stanford.edu/~sergei/papers/kMeansPP-soda.pdf} Changing tie decisions or repeatedly reseeding centers needs separate care; a non-increasing objective alone doesn't rule out a cycle. A converged assignment needn't achieve the global minimum, and a library's iteration limit or tolerance can stop a run before that fixed point.

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 moving the centers, let's measure how well these three seeds fit the six points. Their initial inertia is:

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:

import numpy as np

vectors = np.array([
    [1.0, 1.0], [1.2, 1.8],
    [4.0, 4.2], [5.0, 3.8],
    [8.0, 1.0], [9.0, 1.8],
])
centroids = vectors[[0, 2, 4]].copy()

squared_distance = ((vectors[:, None, :] - centroids[None, :, :]) ** 2).sum(axis=2)
cluster = squared_distance.argmin(axis=1)
updated = np.vstack([vectors[cluster == index].mean(axis=0) for index in range(3)])
inertia = ((vectors - updated[cluster]) ** 2).sum()

print("assignment:", cluster.tolist())
print("updated centroids:", np.round(updated, 2).tolist())
print(f"inertia after update: {inertia:.2f}")

assignment: [0, 0, 1, 1, 2, 2]
updated centroids: [[1.1, 1.4], [4.5, 4.0], [8.5, 1.4]]
inertia after update: 1.74

View larger

Build the k-means loop and prevent initialization traps

Production libraries provide seeding and convergence controls. Our scratch implementation accepts a matrix of vectors and initial centers. It stops when consecutive assignments repeat, so their means remain unchanged too. NumPy's argmin consistently resolves a distance tie to the first center. A 20-step cap limits this teaching run; a production implementation needs an explicit convergence policy.

Check the feature axis before calculating distances: vectors must have shape (n, d) and centers (k, d). A mistaken center matrix of shape (k, 1) can silently broadcast its one coordinate across all d features instead of raising an error.^{[5]Reference 5Broadcasting.https://numpy.org/doc/stable/user/basics.broadcasting.html} Reject that mismatch and non-finite coordinates before assignment. Only finite distances from correctly shaped inputs should determine cluster membership.

import numpy as np

names = np.array(["billing_A", "billing_B", "login_A", "login_B", "deploy_A", "deploy_B"])
vectors = np.array([
    [1.0, 1.0], [1.2, 1.8],
    [4.0, 4.2], [5.0, 3.8],
    [8.0, 1.0], [9.0, 1.8],
])

def kmeans(x: np.ndarray, initial: np.ndarray, max_steps: int = 20):
    x = np.asarray(x, dtype=float)
    centers = np.asarray(initial, dtype=float).copy()
    if x.ndim != 2 or 0 in x.shape:
        raise ValueError("vectors must be a nonempty (n, d) matrix")
    if centers.ndim != 2 or len(centers) == 0 or centers.shape[1] != x.shape[1]:
        raise ValueError("centers must be a nonempty (k, d) matrix with matching features")
    if not np.isfinite(x).all() or not np.isfinite(centers).all():
        raise ValueError("vectors and centers must contain finite coordinates")
    previous_labels = None
    for step in range(1, max_steps + 1):
        squared = ((x[:, None, :] - centers[None, :, :]) ** 2).sum(axis=2)
        labels = squared.argmin(axis=1)
        next_centers = []
        for index in range(len(centers)):
            members = x[labels == index]
            if members.shape[0] == 0:
                raise ValueError(f"cluster {index} is empty")
            next_centers.append(members.mean(axis=0))
        next_centers = np.vstack(next_centers)
        if np.array_equal(labels, previous_labels):
            inertia = ((x - next_centers[labels]) ** 2).sum()
            return labels, next_centers, inertia, step
        centers = next_centers
        previous_labels = labels.copy()
    raise RuntimeError("k-means did not converge")

labels, centers, inertia, steps = kmeans(vectors, vectors[[0, 2, 4]])
for index in range(3):
    members = names[labels == index].tolist()
    print(f"cluster {index}: {members} center={centers[index].round(2).tolist()}")
print(f"converged in {steps} steps with inertia={inertia:.2f}")

cluster 0: ['billing_A', 'billing_B'] center=[1.1, 1.4]
cluster 1: ['login_A', 'login_B'] center=[4.5, 4.0]
cluster 2: ['deploy_A', 'deploy_B'] center=[8.5, 1.4]
converged in 2 steps with inertia=1.74

The initialization trap: why k-means++ matters

Uniform random initialization can pick several seeds in one dense region and leave another region without a seed. Later updates may recover, or they may settle on a poor partition. Topic names still come from our message audit, not from the seeding rule.

Arthur and Vassilvitskii introduced k-means++ initialization in 2007.^{[4]Reference 4k-means++: The Advantages of Careful Seedinghttps://theory.stanford.edu/~sergei/papers/kMeansPP-soda.pdf} After the first uniform draw, the original algorithm uses squared-distance sampling:

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:

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

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

Farther points receive more probability. Already selected coordinates have probability zero, but another point in the same neighborhood can still be selected. If every remaining distance is zero, the formula's denominator vanishes: there are fewer distinct locations than requested centers.

For this original sampling rule, the initial objective satisfies $\mathbb{E}[J_{\mathrm{seed}}]\leq8(\ln K+2)J^*$, where $J^*$ is the globally optimal WCSS. Lloyd updates can only reduce it.^{[4]Reference 4k-means++: The Advantages of Careful Seedinghttps://theory.stanford.edu/~sergei/papers/kMeansPP-soda.pdf} The bound averages over random seeds. Individual runs can start above the optimum or converge at different speeds.

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

import numpy as np

vectors = np.array([
    [1.0, 1.0], [1.2, 1.8],
    [4.0, 4.2], [5.0, 3.8],
    [8.0, 1.0], [9.0, 1.8],
])
names = ["billing_A", "billing_B", "login_A", "login_B", "deploy_A", "deploy_B"]

# First seed picked at index 0 (billing_A)
seed_0 = vectors[0]
d2 = ((vectors - seed_0) ** 2).sum(axis=1)
prob = d2 / d2.sum()

for name, dist2, p in zip(names, d2, prob):
    print(f"{name:10} D(x)^2={dist2:5.2f}  P(selection)={p:5.3f}")

billing_A  D(x)^2= 0.00  P(selection)=0.000
billing_B  D(x)^2= 0.68  P(selection)=0.004
login_A    D(x)^2=19.24  P(selection)=0.122
login_B    D(x)^2=23.84  P(selection)=0.151
deploy_A   D(x)^2=49.00  P(selection)=0.311
deploy_B   D(x)^2=64.64  P(selection)=0.411

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.

Library detail, checked September 21, 2026: scikit-learn 1.9.1 uses greedy k-means++: it tries several candidates at a seeding step and keeps the best. Its n_init="auto" performs one run for init="k-means++". Set an explicit restart count when you want several runs. The theorem above describes the original sampling distribution; this greedy variant has different sampling behavior.^{[2]Reference 2KMeanshttps://scikit-learn.org/stable/modules/generated/sklearn.cluster.KMeans.html}

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:

import numpy as np

vectors = np.array([
    [1.0, 1.0], [1.2, 1.8],
    [4.0, 4.2], [5.0, 3.8],
    [8.0, 1.0], [9.0, 1.8],
])
centers = np.array([[1.0, 1.0], [1.2, 1.8], [100.0, 100.0]])
squared = ((vectors[:, None, :] - centers[None, :, :]) ** 2).sum(axis=2)
labels = squared.argmin(axis=1)
for index in range(3):
    print(f"cluster {index}: {(labels == index).sum()} points")

cluster 0: 1 points
cluster 1: 5 points
cluster 2: 0 points

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.

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 against $k$. Globally optimal inertia can't increase as you add centers and reaches zero by $k=n$. Actual Lloyd runs can violate that trend because they find different local solutions. Compare restarts before interpreting a curve. An "elbow" is a visible slowdown in improvement, not necessarily a mathematical inflection or a uniquely correct topic count:

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):

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

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

Scores near $+1$ mean low within-cluster distances relative to separation. A score near $0$ means $a(i)$ and $b(i)$ are similar; it needn't locate a centroid's Voronoi boundary. A negative score means smaller average distance to another cluster's members, not necessarily to its center.

For a singleton cluster, $a(i)$'s displayed denominator is zero; the convention used below and in scikit-learn assigns silhouette 0. We also return 0 when both distances vanish. Mean silhouette requires 2 through $n-1$ nonempty clusters.^{[6]Reference 6silhouette_scorehttps://scikit-learn.org/stable/modules/generated/sklearn.metrics.silhouette_score.html}

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

from itertools import combinations

import numpy as np

vectors = np.array([
    [1.0, 1.0], [1.2, 1.8],
    [4.0, 4.2], [5.0, 3.8],
    [8.0, 1.0], [9.0, 1.8],
])

def kmeans(x: np.ndarray, initial: np.ndarray, max_steps: int = 20):
    centers = initial.astype(float).copy()
    previous_labels = None
    for _ in range(max_steps):
        squared = ((x[:, None, :] - centers[None, :, :]) ** 2).sum(axis=2)
        labels = squared.argmin(axis=1)
        next_centers = []
        for index in range(len(centers)):
            members = x[labels == index]
            if members.shape[0] == 0:
                return None
            next_centers.append(members.mean(axis=0))
        next_centers = np.vstack(next_centers)
        if np.array_equal(labels, previous_labels):
            inertia = ((x - next_centers[labels]) ** 2).sum()
            return labels, inertia
        centers = next_centers
        previous_labels = labels.copy()
    return None

def mean_silhouette(x: np.ndarray, labels: np.ndarray) -> float:
    if not 2 <= len(np.unique(labels)) <= len(x) - 1:
        raise ValueError("silhouette needs 2 through n-1 nonempty clusters")
    scores = []
    for i, point in enumerate(x):
        same = np.where(labels == labels[i])[0]
        others = same[same != i]
        if others.size == 0:
            scores.append(0.0)
            continue
        cohesion = np.linalg.norm(x[others] - point, axis=1).mean()
        separation = min(
            np.linalg.norm(x[labels == cluster] - point, axis=1).mean()
            for cluster in np.unique(labels) if cluster != labels[i]
        )
        scale = max(cohesion, separation)
        scores.append((separation - cohesion) / scale if scale > 0 else 0.0)
    return float(np.mean(scores))

for k in (2, 3, 4):
    best = None
    for idx in combinations(range(len(vectors)), k):
        result = kmeans(vectors, vectors[list(idx)])
        if result is None:
            continue
        labels, inertia = result
        if best is None or inertia < best[0]:
            best = (inertia, labels)
    if best is None:
        raise RuntimeError(f"no converged nonempty partition for k={k}")
    inertia, labels = best
    print(f"k={k}: inertia={inertia:.2f}, silhouette={mean_silhouette(vectors, labels):.2f}")

k=2: inertia=20.06, silhouette=0.56
k=3: inertia=1.74, silhouette=0.76
k=4: inertia=0.92, silhouette=0.51

Moving from $k=2$ to $k=3$ cuts inertia from 20.06 to 1.74 and raises silhouette to 0.76. With $k=4$, inertia falls to 0.92, but splitting a reviewed pair lowers silhouette to 0.51. The code searches the tested combinations of data-point seeds and returns the best converged result; a global optimum would require a separate guarantee.

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.

import numpy as np

square = np.array([[0.0, 0.0], [0.0, 2.0], [2.0, 0.0], [2.0, 2.0]])
seeds = {
    "left/right": np.array([[0.0, 1.0], [2.0, 1.0]]),
    "bottom/top": np.array([[1.0, 0.0], [1.0, 2.0]]),
}

for name, seed in seeds.items():
    squared = ((square[:, None, :] - seed[None, :, :]) ** 2).sum(axis=2)
    labels = squared.argmin(axis=1)
    centers = np.vstack([square[labels == i].mean(axis=0) for i in range(2)])
    inertia = ((square - centers[labels]) ** 2).sum()
    groups = [np.where(labels == i)[0].tolist() for i in range(2)]
    print(f"{name:10} groups={groups} inertia={inertia:.1f}")

left/right groups=[[0, 1], [2, 3]] inertia=4.0
bottom/top groups=[[0, 2], [1, 3]] inertia=4.0

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 groups. PCA asks a complementary question: can we describe the cloud with fewer continuous coordinates while preserving as much variance as possible?

In the next two formulations, $x_i$ means a centered vector: subtract the same sample mean from every row of $X$ to form $X_c$. Assume $n>1$. PCA has two equivalent mathematical foundations:^{[1]Reference 1Pattern Recognition and Machine Learning.https://www.microsoft.com/en-us/research/publication/pattern-recognition-machine-learning/}

1. Maximum variance formulation: Find a unit vector $w$ ($\|w\|_2=1$) whose projections have maximum sample variance:

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

1. Minimum reconstruction error formulation: Find a unit direction $w$ that minimizes the squared residual sum in the centered space:

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

Summing over all points in the centered dataset yields:

Divide all three terms by $n-1$ to express the identity in variance units. The left side is fixed, so maximizing projected variance minimizes the residual sum. For $m$ components, choose $m$ mutually orthogonal principal directions; reconstruction adds the fitted mean back.

Covariance eigendecomposition versus singular value decomposition

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

From $X_c$, two computational routes lead to the same principal components:^{[7]Reference 7sklearn.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:

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$. Use the thin factorization, with $r=\min(n,d)$, $U\in\mathbb{R}^{n\times r}$, diagonal $\Sigma\in\mathbb{R}^{r\times r}$, and $V\in\mathbb{R}^{d\times r}$:

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

Because $U^\top U = I_r$, columns of $V$ are covariance eigenvectors. If $r<d$, additional orthogonal directions have zero covariance eigenvalue. Singular values and covariance eigenvalues satisfy:

These routes are equivalent in exact arithmetic, but their numerical and storage costs differ:

1. Small directions can disappear. On the nonzero singular spectrum, forming the Gram matrix squares the ratio $\sigma_{\max}/\sigma_{\min}$. A small squared singular value can be lost through rounding relative to larger entries. For full column rank, this is $\kappa_2(X_c^\top X_c)=\kappa_2(X_c)^2$. If $n\ll d$, the covariance matrix is rank deficient and its ordinary condition number is infinite; don't apply that finite-ratio formula without the qualification.
2. Shape matters. With 100 documents and 1536 features, a covariance approach materializes a $1536\times1536$ matrix. A direct SVD avoids that matrix. Randomized methods can approximate a small set of leading directions with matrix products; their cost and accuracy depend on oversampling and power iterations, so $O(ndm)$ isn't a universal exact-SVD cost.^{[8]Reference 8randomized_svdhttps://scikit-learn.org/stable/modules/generated/sklearn.utils.extmath.randomized_svd.html}

Solver choice, checked September 21, 2026: scikit-learn 1.9.1 PCA supports both SVD solvers and covariance_eigh. Its automatic policy selects covariance eigendecomposition for a sufficiently tall, narrow matrix. That can be efficient when $n\gg d$ and $d$ is small; it has a different precision trade-off from direct SVD.^{[7]Reference 7sklearn.decomposition.PCAhttps://scikit-learn.org/stable/modules/generated/sklearn.decomposition.PCA.html} Choose a solver for your data shape and accuracy needs rather than assuming every PCA implementation takes the same route.

import numpy as np

vectors = np.array([
    [1.0, 1.0], [1.2, 1.8],
    [4.0, 4.2], [5.0, 3.8],
    [8.0, 1.0], [9.0, 1.8],
])
n = len(vectors)
centered = vectors - vectors.mean(axis=0)

# Route A: SVD on centered data
_, s, vt = np.linalg.svd(centered, full_matrices=False)
eigenvalues_svd = (s ** 2) / (n - 1)

# Route B: Eigendecomposition of sample covariance matrix
cov = (centered.T @ centered) / (n - 1)
eigenvalues_cov = np.sort(np.linalg.eigvalsh(cov))[::-1]

print("SVD eigenvalues:       ", np.round(eigenvalues_svd, 4).tolist())
print("Covariance eigenvalues:", np.round(eigenvalues_cov, 4).tolist())
print("Eigenvalues match:     ", np.allclose(eigenvalues_svd, eigenvalues_cov))

SVD eigenvalues:        [11.1825, 1.9442]
Covariance eigenvalues: [11.1825, 1.9442]
Eigenvalues match:      True

The explained variance ratio measures the proportion of total linear spread retained by each component, provided total variance is nonzero. Include zero eigenvalues when extending a thin SVD to $d$ directions:

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

import numpy as np

names = np.array(["billing_A", "billing_B", "login_A", "login_B", "deploy_A", "deploy_B"])
vectors = np.array([
    [1.0, 1.0], [1.2, 1.8],
    [4.0, 4.2], [5.0, 3.8],
    [8.0, 1.0], [9.0, 1.8],
])

centered = vectors - vectors.mean(axis=0)
_, singular_values, vt = np.linalg.svd(centered, full_matrices=False)
pc1 = vt[0]
if pc1[0] < 0:
    pc1 = -pc1
projection = centered @ pc1
ratio = singular_values**2 / (singular_values**2).sum()

print("mean:", np.round(vectors.mean(axis=0), 2).tolist())
print("pc1:", np.round(pc1, 3).tolist())
print("explained variance:", np.round(ratio, 3).tolist())
for name, value in zip(names, projection):
    print(f"{name:10} pc1={value:5.2f}")

mean: [4.7, 2.27]
pc1: [1.0, -0.016]
explained variance: [0.852, 0.148]
billing_A  pc1=-3.68
billing_B  pc1=-3.49
login_A    pc1=-0.73
login_B    pc1= 0.27
deploy_A   pc1= 3.32
deploy_B   pc1= 4.31

Keeping PC1 produces one coordinate as a linear combination of centered dimensions. A direction's sign is arbitrary: reversing it changes the sign of its scores while preserving distances and reconstruction. Equal eigenvalues also allow rotations within the tied eigenspace; cutting through such a tie doesn't specify a unique retained subspace.

View larger

⚠️ 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:

import numpy as np

vectors = np.array([
    [1.0, 1.0], [1.2, 1.8],
    [4.0, 4.2], [5.0, 3.8],
    [8.0, 1.0], [9.0, 1.8],
])
mean = vectors.mean(axis=0)
centered = vectors - mean
_, singular_values, vt = np.linalg.svd(centered, full_matrices=False)

for components in (1, 2):
    basis = vt[:components]
    reconstructed = (centered @ basis.T) @ basis + mean
    error = ((vectors - reconstructed) ** 2).sum()
    kept = (singular_values[:components] ** 2).sum() / (singular_values ** 2).sum()
    print(f"components={components}: variance_kept={kept:.3f}, squared_error={error:.2f}")

components=1: variance_kept=0.852, squared_error=9.72
components=2: variance_kept=1.000, squared_error=0.00

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.

Connect clustering to vector search: inverted file indexing

Explained variance still doesn't tell a query which documents to touch. An inverted file (IVF) index uses k-means for a coarser approximate nearest-neighbor (ANN) search: which few cells might hold the nearest document? Faiss provides IVF indexes. pgvector provides IVFFlat.^{[9]Reference 9Faiss indexeshttps://github.com/facebookresearch/faiss/wiki/Faiss-indexes[10]Reference 10pgvectorhttps://github.com/pgvector/pgvector} Vector databases choose among different index structures. Qdrant's current dense-vector index uses HNSW, a graph method.^{[11]Reference 11Indexing: Vector Indexhttps://qdrant.tech/documentation/manage-data/indexing/#vector-index}

Exhaustive search compares a query with all $N$ stored $d$-dimensional vectors, using $O(Nd)$ arithmetic. Whether that meets a latency budget depends on corpus size, hardware, batching, filters, and memory access. Keep a flat-search baseline; an approximate index isn't automatically faster for a small corpus.

A simple Euclidean IVFFlat index uses k-means to reduce the candidate scan:

1. Training phase: Learn $K$ coarse centroids from a representative corpus sample. Choose $K$ by measuring routing cost, candidate scans, and recall; a square-root-of-$N$ starting rule is a heuristic, not an optimum.^{[9]Reference 9Faiss indexeshttps://github.com/facebookresearch/faiss/wiki/Faiss-indexes}
2. Index building: Assign each document vector to its closest centroid. Each centroid's inverted list stores its assigned vectors and document IDs. IVFFlat retains full vectors; IVF variants can instead store compressed codes.
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.

With equal-size lists, scanning $nprobe$ of $K$ lists visits $N(nprobe/K)$ documents. Including the flat centroid search gives estimated work $O(Kd+N(nprobe/K)d)$. If $K=4096$ and $nprobe=16$, the candidate scan skips about 99.6% of documents under that balance assumption. Unequal list lengths change the fraction, and routing still costs work.^{[9]Reference 9Faiss indexeshttps://github.com/facebookresearch/faiss/wiki/Faiss-indexes} For fixed $K$ and $nprobe$, the candidate scan still grows linearly with $N$.

View larger

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

import numpy as np

vectors = np.array([
    [1.0, 1.0], [1.2, 1.8],  # billing
    [4.0, 4.2], [5.0, 3.8],  # login
    [8.0, 1.0], [9.0, 1.8],  # deploy
])
names = ["billing_A", "billing_B", "login_A", "login_B", "deploy_A", "deploy_B"]

# Centroids learned during index training (K=3)
centroids = np.array([
    [1.1, 1.4],  # C0: billing
    [4.5, 4.0],  # C1: login
    [8.5, 1.4],  # C2: deploy
])

# Build inverted lists
assignments = ((vectors[:, None, :] - centroids[None, :, :]) ** 2).sum(axis=2).argmin(axis=1)
inverted_lists = {cluster_id: [] for cluster_id in range(3)}
for doc_idx, cluster_id in enumerate(assignments):
    inverted_lists[cluster_id].append(doc_idx)

# Incoming user query: "Cannot authenticate with my 2FA token"
query = np.array([4.2, 4.0])

# Step 1: Coarse search over 3 centroids
centroid_distances = np.linalg.norm(centroids - query, axis=1)
nprobe = 1
probed_clusters = centroid_distances.argsort()[:nprobe]

# Step 2: Fine search only within probed inverted lists
candidate_indices = [idx for c in probed_clusters for idx in inverted_lists[c]]
candidate_distances = np.linalg.norm(vectors[candidate_indices] - query, axis=1)
best_match = candidate_indices[candidate_distances.argmin()]

print("probed clusters:", probed_clusters.tolist())
print("scanned docs:   ", [names[i] for i in candidate_indices])
print("nearest match:  ", names[best_match], f"(distance={candidate_distances.min():.2f})")
print(f"pruning ratio:   {len(names) - len(candidate_indices)} of {len(names)} vectors skipped")

probed clusters: [1]
scanned docs:    ['login_A', 'login_B']
nearest match:   login_A (distance=0.28)
pruning ratio:   4 of 6 vectors skipped

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

• Low $nprobe$: Fewer lists to scan. The nearest document may belong to an unprobed cell, especially near a boundary; throughput still depends on list sizes and execution costs.
• High $nprobe$: More candidates and usually better ANN recall. Scanning all lists with exact distances and consistent tie handling makes IVFFlat exhaustive. Compressed IVF codes can still change rankings even when all lists are scanned; candidate recovery and distance precision are separate issues.

Make the nearest centroid give the wrong answer

The login query above happens to succeed with one probe. Move the query to $[2.8,2.71]$, just on the C1 side of the C0/C1 boundary. Before running the code, compare the two questions: "Which center is closest?" and "Which document is closest?"

import numpy as np

names = np.array(["billing_A", "billing_B", "login_A", "login_B", "deploy_A", "deploy_B"])
vectors = np.array([
    [1.0, 1.0], [1.2, 1.8], [4.0, 4.2],
    [5.0, 3.8], [8.0, 1.0], [9.0, 1.8],
])
centroids = np.array([[1.1, 1.4], [4.5, 4.0], [8.5, 1.4]])
assignment = ((vectors[:, None, :] - centroids[None, :, :]) ** 2).sum(axis=2).argmin(axis=1)
query = np.array([2.8, 2.71])
coarse = np.linalg.norm(centroids - query, axis=1)
distance = np.linalg.norm(vectors - query, axis=1)
exact = int(distance.argmin())
print(f"C0 distance={coarse[0]:.3f}; C1 distance={coarse[1]:.3f}")
print(f"flat: {names[exact]} at {distance[exact]:.3f}")

for nprobe in (1, 2, 3):
    cells = coarse.argsort()[:nprobe]
    candidates = np.flatnonzero(np.isin(assignment, cells))
    fine = np.linalg.norm(vectors[candidates] - query, axis=1)
    winner = int(candidates[fine.argmin()])
    comparisons = len(centroids) + len(candidates)
    print(f"nprobe={nprobe}: cells={cells.tolist()} winner={names[winner]} "
          f"ANN recall@1={int(winner == exact)} distances={comparisons}")

C0 distance=2.146; C1 distance=2.134
flat: billing_B at 1.841
nprobe=1: cells=[1] winner=login_A ANN recall@1=0 distances=5
nprobe=2: cells=[1, 0] winner=billing_B ANN recall@1=1 distances=7
nprobe=3: cells=[1, 0, 2] winner=billing_B ANN recall@1=1 distances=9

The comparison counts describe one ANN search, excluding the six-distance flat baseline we ran for verification. ANN recall@1 asks whether the approximate winner matches the exact geometric winner. For top $r$, compare overlap with the exact top $r$, resolving ties consistently. This differs from semantic recall against human relevance labels. Even 100% ANN recall can reproduce poor embedding neighborhoods perfectly.

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.

import numpy as np

# Column 0 carries topic; column 1 is message length in characters.
names = ["bill_100", "bill_500", "bill_900", "dep_110", "dep_510", "dep_910"]
features = np.array([
    [-3.2, 100.0], [-3.0, 500.0], [-2.8, 900.0],
    [ 2.8, 110.0], [ 3.0, 510.0], [ 3.2, 910.0],
])

def cluster(x: np.ndarray) -> list[int]:
    centers = x[[0, 3]].copy()  # Same source rows in both runs.
    previous = None
    for _ in range(20):
        squared = ((x[:, None, :] - centers[None, :, :]) ** 2).sum(axis=2)
        labels = squared.argmin(axis=1)
        if np.array_equal(labels, previous):
            return labels.tolist()
        groups = [x[labels == i] for i in range(2)]
        if any(len(group) == 0 for group in groups):
            raise ValueError("empty cluster")
        centers = np.vstack([group.mean(axis=0) for group in groups])
        previous = labels.copy()
    raise RuntimeError("k-means did not converge")

print("raw dist bill_100 to dep_110:", round(float(np.linalg.norm(features[0] - features[3])), 1))
print("raw dist bill_100 to bill_900:", round(float(np.linalg.norm(features[0] - features[2])), 1))

raw_labels = cluster(features)
print("raw groups, seeds bill_100 and dep_110:")
for index in range(2):
    members = [names[i] for i, label in enumerate(raw_labels) if label == index]
    print(f"  cluster {index}: {members}")

scaled = (features - features.mean(axis=0)) / features.std(axis=0)
scaled_labels = cluster(scaled)
print("scaled groups, same seed rows:")
for index in range(2):
    members = [names[i] for i, label in enumerate(scaled_labels) if label == index]
    print(f"  cluster {index}: {members}")

raw dist bill_100 to dep_110: 11.7
raw dist bill_100 to bill_900: 800.0
raw groups, seeds bill_100 and dep_110:
  cluster 0: ['bill_100', 'dep_110']
  cluster 1: ['bill_500', 'bill_900', 'dep_510', 'dep_910']
scaled groups, same seed rows:
  cluster 0: ['bill_100', 'bill_500', 'bill_900']
  cluster 1: ['dep_110', 'dep_510', 'dep_910']

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. Scaling gives coordinates a chosen relative weight; it doesn't reveal semantic importance. If a model specifies normalization, follow that contract. Fit the scaler or PCA on permitted corpus/development data and reuse the transform for held-out queries. Corpus documents available at serving time can legitimately define an unsupervised index or transform. Fitting on held-out query data violates an evaluation that promises unseen queries; an explicitly transductive protocol makes a different promise. As with supervised validation, specify the boundary rather than declaring every use of unlabeled evaluation data equivalent.

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 documents that API embedding outputs are L2-normalized by default, including after shortening with dimensions (checked September 21, 2026).^{[12]Reference 12Embeddings Frequently Asked Questionshttps://help.openai.com/en/articles/6824809-embeddings-frequently-asked-questions} For unit vectors, cosine similarity equals dot product, and squared Euclidean distance satisfies:

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.

import numpy as np

query = np.array([1.0, 0.0])
names = ["duplicate_charge_policy", "very_long_deploy_page"]
documents = np.array([
    [1.0, 0.1],
    [8.0, 4.0],
])

dot = documents @ query
cosine = dot / (np.linalg.norm(documents, axis=1) * np.linalg.norm(query))

print("dot-product winner:", names[int(dot.argmax())])
print("cosine winner:     ", names[int(cosine.argmax())])
print("cosine scores:     ", np.round(cosine, 3).tolist())

dot-product winner: very_long_deploy_page
cosine winner:      duplicate_charge_policy
cosine scores:      [0.995, 0.894]

That ranking equivalence requires the compared vectors to have the same normalization contract. Standard k-means centers are coordinate means, generally not unit vectors:

import numpy as np

a = np.array([1.0, 0.0])
b = np.array([0.0, 1.0])
mean = (a + b) / 2
print("mean:", np.round(mean, 3).tolist())
print(f"mean L2 norm: {np.linalg.norm(mean):.3f}")

mean: [0.5, 0.5]
mean L2 norm: 0.707

For a unit query $x$, ordinary Euclidean assignment minimizes $1+\|\mu_k\|^2-2x^\top\mu_k$. Different center norms change the decision. Normalizing only during serving can therefore route documents differently from the assignment used to build lists. Two constructed groups below demonstrate the difference; they aren't presented as a converged clustering run.

import numpy as np

# C0 is the mean of two unit vectors; C1 contains one unit vector.
c0_members = np.array([[0.1, np.sqrt(0.99)], [0.1, -np.sqrt(0.99)]])
centers = np.vstack([c0_members.mean(axis=0), [0.8, 0.6]])
query = np.array([1.0, 0.0])
squared = ((centers - query) ** 2).sum(axis=1)
cosine = (centers @ query) / np.linalg.norm(centers, axis=1)
print("squared Euclidean distances:", np.round(squared, 2).tolist())
print("cosine similarities:", np.round(cosine, 2).tolist())
print("Euclidean routes to C" + str(int(squared.argmin())))
print("cosine routes to C" + str(int(cosine.argmax())))

squared Euclidean distances: [0.81, 0.4]
cosine similarities: [1.0, 0.8]
Euclidean routes to C1
cosine routes to C0

For a cosine clustering objective, use spherical k-means: normalize nonzero inputs, assign by largest dot product to unit centers, and normalize each cluster's vector sum after updates.^{[13]Reference 13Faiss FAQ: k-means for the cosine similarityhttps://github.com/facebookresearch/faiss/wiki/FAQ#k-means-for-the-cosine-similarity} Empty clusters and a zero sum, such as two opposite unit vectors, require an explicit fallback. If you choose Euclidean k-means instead, keep build and query routing consistent and measure the index's held-out recall.

Euclidean k-means favors compact groups under squared distance, but it doesn't impose a statistical equal-size or spherical distribution assumption. Its Voronoi assignment regions are convex. It can split an elongated semantic group or dedicate a center to an outlier; inspect whether the partition suits the task.

Failure test 3: centering changes cosine geometry

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

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.^{[7]Reference 7sklearn.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 can alter angles and rankings. By contrast, subtracting the same mean preserves Euclidean pairwise differences: $(u-\mu)-(v-\mu)=u-v$. A complete orthonormal change of basis also preserves those Euclidean distances; truncation can discard them.

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

import numpy as np

query = np.array([1.0, 2.0])
doc1 = np.array([1.1, 2.1])  # Close to query direction
doc2 = np.array([2.0, 1.0])  # Further away in direction

def cos_sim(u, v):
    denominator = np.linalg.norm(u) * np.linalg.norm(v)
    if denominator == 0:
        raise ValueError("cosine is undefined for a zero vector")
    return float(np.dot(u, v) / denominator)

orig_cos1 = cos_sim(query, doc1)
orig_cos2 = cos_sim(query, doc2)

documents = np.array([doc1, doc2])
mean = documents.mean(axis=0)
q_centered = query - mean
d1_centered = doc1 - mean
d2_centered = doc2 - mean

centered_cos1 = cos_sim(q_centered, d1_centered)
centered_cos2 = cos_sim(q_centered, d2_centered)

print(f"Original cos(query, doc1): {orig_cos1:.4f}")
print(f"Original cos(query, doc2): {orig_cos2:.4f}")
print(f"Centered cos(query, doc1): {centered_cos1:.4f}")
print(f"Centered cos(query, doc2): {centered_cos2:.4f}")

assert orig_cos1 > orig_cos2
assert centered_cos1 > 0
assert centered_cos2 < 0

Original cos(query, doc1): 0.9998
Original cos(query, doc2): 0.8000
Centered cos(query, doc1): 0.9802
Centered cos(query, doc2): -0.9802

After centering, doc2 sits more than 90 degrees from the query. The winner doesn't change in this example, but its angles do. A document equal to the fitted mean would become a zero vector, whose cosine is undefined. If you reduce dimensions for cosine retrieval, treat the compressed space as a new representation and validate its rankings, zero-vector handling, and 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:

import numpy as np

names = np.array(["billing_policy", "billing_help", "login_2fa", "deploy_policy", "runbook_note"])
documents = np.array([
    [-20.0, 0.0],
    [ 0.0, 0.0],
    [ 0.0, 3.0],
    [20.0, 0.0],
    [ 0.0, -3.0],
])
query = np.array([0.0, 3.2])

full_distance = np.linalg.norm(documents - query, axis=1)
full_winner = names[full_distance.argmin()]

mean = documents.mean(axis=0)
centered = documents - mean
singular_values, vt = np.linalg.svd(centered, full_matrices=False)[1:]
pc1 = vt[:1]
compressed_documents = centered @ pc1.T
compressed_query = (query - mean) @ pc1.T
compressed_distance = np.linalg.norm(compressed_documents - compressed_query, axis=1)
compressed_ties = names[np.isclose(compressed_distance, compressed_distance.min())].tolist()
kept = float((singular_values ** 2)[0] / (singular_values ** 2).sum())

print(f"PC1 variance kept: {kept:.3f}")
print("full-space nearest: ", full_winner)
print("one-PC nearest ties: ", compressed_ties)

PC1 variance kept: 0.978
full-space nearest:  login_2fa
one-PC nearest ties:  ['billing_help', 'login_2fa', 'runbook_note']

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. After repeated tuning, treat them as development data and use a separate untouched benchmark for evaluation.

import numpy as np

names = np.array(["billing_A", "billing_B", "login_A", "login_B", "deploy_A", "deploy_B"])
reviewed_topic = np.array(["billing", "billing", "login", "login", "deploy", "deploy"])
vectors = np.array([
    [1.0, 1.0], [1.2, 1.8],
    [4.0, 4.2], [5.0, 3.8],
    [8.0, 1.0], [9.0, 1.8],
])

squared = ((vectors[:, None, :] - vectors[[0, 2, 4]][None, :, :]) ** 2).sum(axis=2)
labels = squared.argmin(axis=1)
centers = np.vstack([vectors[labels == index].mean(axis=0) for index in range(3)])
labels = ((vectors[:, None, :] - centers[None, :, :]) ** 2).sum(axis=2).argmin(axis=1)

for cluster in range(3):
    topics, counts = np.unique(reviewed_topic[labels == cluster], return_counts=True)
    dominant = topics[counts.argmax()]
    purity = counts.max() / counts.sum()
    members = names[labels == cluster].tolist()
    print(f"cluster {cluster}: proposed={dominant:8} purity={purity:.2f} members={members}")

cluster 0: proposed=billing  purity=1.00 members=['billing_A', 'billing_B']
cluster 1: proposed=login    purity=1.00 members=['login_A', 'login_B']
cluster 2: proposed=deploy   purity=1.00 members=['deploy_A', 'deploy_B']

On this teaching set, purity is 1.0 because we designed clean coordinates and reviewed every row. Purity alone rewards fragmentation: making every message its own cluster also gives purity 1.0. Record cluster sizes, coverage of each reviewed topic, mixed and boundary examples, multilingual cases, 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. Reproduce the miss at $[2.8,2.71]$. The C0/C1 perpendicular bisector is $3.4q_x+2.6q_y=16.54$; verify the center-distance tie at $[2.8,2.7]$. Then move slightly to each side and compare center routing with exact document lookup. Count coarse and fine comparisons separately.
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++ changes a distribution over seeds. Compare inertia and iteration distributions across the same restart budget; don't expect every run to improve both. Greedy library seeding and the original single-candidate rule are different experiments.
• When unscaled metadata dominates distance, clustering reflects units rather than semantics. Verify feature normalization before interpreting clusters.
• One probe can miss an exact neighbor near a cell boundary. Sweep $nprobe$ against flat-search ANN recall and measured latency; then evaluate semantic relevance separately.
• PCA preserves global variance, not task labels. Evaluating downstream retrieval metrics on held-out pairs is mandatory before deploying compressed representations.