The Mathematics Behind Unsupervised Learning: PCA and K-Means

The mathematics under the toolkit: what k-means minimises, why Lloyd’s algorithm converges, why scaling decides the clustering, how maximising variance produces eigenvectors, and where PCA and k-means turn out to meet.

Unsupervised learning removes the one thing supervised learning leans on. There are no labels, no correct answers, and no loss function comparing a prediction against a truth. What remains is geometry. Every method in the field, from k-means to PCA to the recommender at the end of the pipeline, is an argument about distance, variance, and direction in a space where the data points are the only evidence available.

This article works through the mathematics of the two foundational techniques. K-means asks where the groups are, and answers it by minimising a sum of squared distances. PCA asks which directions matter, and answers it by maximising variance. Both reduce to optimisation problems with closed forms or provably convergent algorithms, and both explain, once you see the mathematics, the practical rules that otherwise look like folklore: why scaling changes everything, why inertia can never choose k for you, and why the components come out uncorrelated whether you wanted that or not.

What k-means is actually minimising

Start with a dataset of n points in d dimensions, and suppose we want k clusters. A clustering is two things at once: an assignment of every point to a cluster, and a centroid for each cluster. K-means treats the quality of that pair as a single number, the within-cluster sum of squares, usually called inertia:

J=j=1kxiCjxiμj2J = \sum_{j=1}^{k} \sum_{x_i \in C_j} \lVert x_i – \mu_j \rVert^2

Here CjC_j is the set of points assigned to cluster j, and μj\mu_j is that cluster’s centroid. Every point contributes the squared distance to the centre of its own cluster, and J adds them all up.

Read the formula as a definition of what “good” means. A clustering is good when points sit close to their own centre, and J is small exactly when the clusters are tight. Nothing about labels or categories appears anywhere in the expression, which is precisely why the method is unsupervised: it is optimising a geometric quantity that any unlabelled dataset possesses.

The objective also explains why the centroid is a mean rather than any other summary. Hold the assignments fixed and ask which point μj\mu_j minimises the inner sum. Differentiating with respect to μj\mu_j and setting the result to zero gives:

μjxiCjxiμj2=2xiCj(xiμj)=0\frac{\partial}{\partial \mu_j} \sum_{x_i \in C_j} \lVert x_i – \mu_j \rVert^2 = -2 \sum_{x_i \in C_j} (x_i – \mu_j) = 0

Solving gives the arithmetic mean of the cluster’s members:

μj=1|Cj|xiCjxi\mu_j = \frac{1}{|C_j|} \sum_{x_i \in C_j} x_i

The mean is not a design choice; it is the unique minimiser of squared distance. Had the objective used absolute distances instead of squared ones, the same derivation would have produced the median, which is exactly what k-medians does. The loss function chooses the statistic.

Lloyd’s algorithm and why it converges

Minimising J over all possible assignments is NP-hard, since the number of ways to partition n points into k groups explodes combinatorially. Lloyd’s algorithm, which is what KMeans runs, sidesteps the search with alternating minimisation. It repeats two steps:

assign: ci=argminjxiμj2\text{assign: } c_i = \arg\min_j \lVert x_i – \mu_j \rVert^2
update: μj=1|Cj|xiCjxi \text{update: } \mu_j = \frac{1}{|C_j|} \sum_{x_i \in C_j} x_i

The convergence argument is short and worth following, because it explains both the algorithm’s reliability and its weakness. Each step minimises J with respect to one set of variables while holding the other fixed. The assignment step sends every point to its nearest centroid, which cannot increase any point’s contribution to J. The update step moves each centroid to the mean, which we just proved is the minimiser for fixed assignments, so it cannot increase J either. Therefore J never rises. Since J is bounded below by zero and the number of possible assignments is finite, the algorithm must terminate.

Notice what the argument does not say. It proves J decreases monotonically to a local minimum, not that it reaches the global one. The algorithm’s final answer depends on where the centroids started, which is why n_init=10 exists: it runs the whole procedure ten times from different initialisations and keeps the run with the lowest J. The mathematics guarantees you will reach a valley; only repetition improves the odds it is a deep one.

Why scaling decides the outcome

The objective is built from the Euclidean norm, and expanding it reveals the problem:

xiμj2=m=1d(ximμjm)2\lVert x_i – \mu_j \rVert^2 = \sum_{m=1}^{d} (x_{im} – \mu_{jm})^2

Every feature contributes its own squared difference, and the sum treats those contributions as equals. But features are not equals in raw units. In the wine data, proline ranges over roughly a thousand units while hue ranges over about one. A difference of 200 in proline contributes 40,000 to the distance, and a difference of 0.5 in hue contributes 0.25. The proline term is five orders of magnitude larger, so the sum is effectively proline alone, and the clustering describes one column rather than the data.

Standardisation removes the unit dependence by expressing every feature in standard deviations from its own mean:

zim=ximxmσmz_{im} = \frac{x_{im} – \bar{x}_m}{\sigma_m}

Now each feature has variance one and contributes to the distance on the same terms. This is the mathematics behind the empirical result in the code-along workbook, where the identical k-means jumped from roughly 70 percent cluster purity to near-perfect after a single StandardScaler. The algorithm never changed; the geometry it was measuring did.

The same reasoning explains the alternatives. When only the direction of a vector carries meaning and its length does not, as with word-frequency profiles, Normalizer rescales each row to unit length so that distance measures angle. When the matrix is sparse and centring it would destroy that sparsity by filling every zero, MaxAbsScaler divides by the largest absolute value and leaves the zeros untouched. Each scaler is a different answer to the same question: what should distance mean here?

Why inertia can never choose k

Practitioners learn to distrust inertia as a criterion for k without always learning why. The reason is a theorem, not a heuristic. Adding a cluster can never increase J, because any clustering with k clusters can be reproduced with k+1 by splitting one cluster in two, and splitting cannot raise the sum of squared distances:

Jk+1Jkfor all kJ_{k+1} \leq J_k \quad \text{for all } k

The inequality is strict in practice, so J falls with every additional cluster until it reaches zero at k = n, where every point is its own centroid and every distance is zero. A perfect score, and a completely useless model. Minimising inertia over k therefore always recommends the maximum possible k.

This is why the elbow method reads the shape of the curve instead of its minimum. The quantity of interest is the marginal gain:

Δk=Jk1Jk\Delta_k = J_{k-1} – J_k

While k is below the true number of groups, each new cluster splits a genuinely heterogeneous blob and Δk\Delta_k is large. Once k passes the true number, additional clusters can only subdivide groups that are already coherent, and Δk\Delta_k collapses to the small, steady reductions you would get from splitting noise. The elbow is where the sequence of Δk\Delta_k drops sharply, and that bend is a signal about the data’s structure that the value of J alone cannot express.

PCA: the variance of a direction

PCA starts from a different question. Rather than grouping points, it asks which directions in the feature space carry information. The mathematics defines information as variance, and the entire method follows from taking that definition seriously.

First, centre the data so each feature has mean zero. The covariance matrix is then:

Σ=1n1XX\Sigma = \frac{1}{n-1} X^{\top} X

Σ\Sigma is a d by d symmetric matrix whose diagonal holds each feature’s variance and whose off-diagonal entries hold the covariances between pairs.

Now take any unit vector w representing a direction, and project every point onto it. The projection of point xix_i is the scalar wxiw^{\top} x_i, and the variance of those projected values is:

Var(Xw)=1n1Xw2=wΣw\operatorname{Var}(Xw) = \frac{1}{n-1} \lVert Xw \rVert^2 = w^{\top} \Sigma w

This quadratic form is the heart of PCA. It takes a direction and returns how spread out the data is along it. A direction through the data’s long axis returns a large number; a direction through its thin axis returns a small one. Finding the most informative direction is now a well-posed optimisation problem.

The constraint, the Lagrangian, and the eigenvector

We want the direction of maximum variance, but the problem needs a constraint, since scaling w up scales wΣww^{\top} \Sigma w up with it and the maximum would be unbounded. Restricting w to unit length fixes this:

maxw;wΣwsubject toww=1\max_{w} ; w^{\top} \Sigma w \quad \text{subject to} \quad w^{\top} w = 1

Introduce a Lagrange multiplier $\lambda$ and form the Lagrangian:

(w,λ)=wΣwλ(ww1)\mathcal{L}(w, \lambda) = w^{\top} \Sigma w – \lambda (w^{\top} w – 1)

Differentiating with respect to w and setting the gradient to zero gives:

2Σw2λw=02 \Sigma w – 2 \lambda w = 0

which rearranges into the eigenvector equation:

Σw=λw\Sigma w = \lambda w

This is the moment PCA stops being a statistical procedure and becomes linear algebra. The directions of stationary variance are exactly the eigenvectors of the covariance matrix, and nothing was assumed to make that happen; it fell out of maximising a quadratic form under a norm constraint.

The eigenvalue tells us what we gained. Left-multiplying the eigenvector equation by ww^{\top} and using ww=1w^{\top} w = 1:

wΣw=λww=λw^{\top} \Sigma w = \lambda w^{\top} w = \lambda

The variance captured by a principal component is its eigenvalue. Ranking components by variance and ranking eigenvalues are the same operation, which is why explained_variance_ratio_ is simply the eigenvalues normalised to sum to one:

explained variance ratiom=λml=1dλl \text{explained variance ratio}m = \frac{\lambda_m}{\sum{l=1}^{d} \lambda_l}

Intrinsic dimension and the trace

The denominator above deserves attention. For a symmetric matrix the sum of the eigenvalues equals the trace, and the trace of the covariance matrix is the sum of the individual feature variances:

m=1dλm=tr(Σ)=m=1dVar(xm)\sum_{m=1}^{d} \lambda_m = \operatorname{tr}(\Sigma) = \sum_{m=1}^{d} \operatorname{Var}(x_m)

Total variance is conserved. PCA does not create or destroy information; it redistributes the same total across a new set of axes, packing as much as possible into the first few. When the wine data’s first three components hold two-thirds of the variance, the statement is precise: two-thirds of the total spread of all thirteen original features lives in a three-dimensional subspace, and the remaining ten directions share the last third between them.

That is what intrinsic dimension means mathematically. The data nominally occupies d dimensions but effectively occupies far fewer, and the eigenvalue spectrum measures the gap. Choosing how many components to keep is choosing where on the cumulative curve to stop, and the curve, not a rule of thumb, is the evidence.

Decorrelation is not a bonus, it is the theorem

PCA’s second gift looks like luck until you write it down. Because Σ\Sigma is real and symmetric, the spectral theorem guarantees it can be diagonalised by an orthogonal matrix of its eigenvectors:

Σ=WΛW,WW=I \Sigma = W \Lambda W^{\top}, \qquad W^{\top} W = I

Transform the data with $Z = XW$ and compute the covariance of the result:

ΣZ=1n1ZZ=1n1WXXW=WΣW=Λ\Sigma_Z = \frac{1}{n-1} Z^{\top} Z = \frac{1}{n-1} W^{\top} X^{\top} X W = W^{\top} \Sigma W = \Lambda

The covariance matrix of the transformed data is Λ\Lambda, which is diagonal. Every off-diagonal entry, meaning every covariance between two distinct components, is exactly zero. This is why petal length and petal width can correlate at +0.96 in the raw iris data while the first two principal components correlate at 0.00. The decorrelation is not an empirical tendency to be checked; it is an algebraic identity, guaranteed by the orthogonality of the eigenvectors of a symmetric matrix.

It also explains why PCA is a standard pre-step for models that suffer under correlated inputs. Multicollinearity is a property of the covariance matrix’s off-diagonal structure, and PCA transforms that structure away by construction.

Loadings: what a component means

An eigenvector is a vector of d numbers, one per original feature, and those numbers are the loadings. Since the component is the linear combination:

zi1=m=1dw1m,ximz_{i1} = \sum_{m=1}^{d} w_{1m} , x_{im}

each loading w1mw_{1m} is the weight of feature m in the recipe. When several related features carry heavy weights with the same sign, the component is measuring whatever they have in common, which is how PC1 on the iris data becomes interpretable as overall petal size. The interpretation is not mystical; it is reading the coefficients of a weighted sum.

SVD: the same mathematics, computed properly

In practice scikit-learn does not build Σ\Sigma and extract its eigenvectors. It uses the singular value decomposition of the centred data matrix directly:

X=USVX = U S V^{\top}

Substituting into the covariance definition shows the two routes lead to the same place:

Σ=1n1XX=1n1VSUUSV=V(S2n1)V \Sigma = \frac{1}{n-1} X^{\top} X = \frac{1}{n-1} V S^{\top} U^{\top} U S V^{\top} = V \left( \frac{S^2}{n-1} \right) V^{\top}

Comparing this with Σ=WΛW\Sigma = W \Lambda W^{\top} identifies the correspondence exactly: the right singular vectors V are the principal components, and the eigenvalues relate to the singular values by:

λm=sm2n1\lambda_m = \frac{s_m^2}{n-1}

The preference for SVD is numerical. Forming XXX^{\top} X squares the condition number of the problem and can lose precision that the decomposition of X itself retains.

This also resolves the sparse-data question with a clean mathematical answer. PCA requires centring, and centring a sparse matrix means subtracting a nonzero mean from every entry, including the zeros, which destroys sparsity and can turn a manageable matrix into one that will not fit in memory. TruncatedSVD factorises X without centring it, so the zeros survive and the decomposition stays cheap. The trade is real: without centring, the first component often captures the mean direction rather than a direction of variation. On TF-IDF matrices, where the data is non-negative and sparse by nature, that trade is worth making, which is the mathematics behind the rule that dense data takes PCA and sparse text takes TruncatedSVD.

Where the two methods meet

K-means and PCA look unrelated, one partitioning points and the other rotating axes, but they optimise closely connected quantities. Total variance decomposes into a within-cluster part and a between-cluster part:

ixix2total=jxiCjxiμj2within: J+j|Cj|,μjx2between\underbrace{\sum_i \lVert x_i – \bar{x} \rVert^2}{\text{total}} = \underbrace{\sum_j \sum{x_i \in C_j} \lVert x_i – \mu_j \rVert^2}{\text{within: } J} + \underbrace{\sum_j |C_j| , \lVert \mu_j – \bar{x} \rVert^2}{\text{between}}

The left side is fixed by the data. Therefore minimising J, which is what k-means does, is identical to maximising the between-cluster term, which is a variance quantity of exactly the kind PCA maximises. The formal result is that the cluster indicators of k-means, once relaxed from binary values to continuous ones, are spanned by the top principal components. Running PCA before k-means is not merely a speed trick, then. It is projecting onto the subspace where the separating variance already lives, which is why the TF-IDF pipeline of TruncatedSVD followed by k-means works as well as it does.

What the mathematics tells you to do

Five practical rules in the workbook are consequences of the derivations above rather than accumulated habit. Scale before clustering, because the objective sums squared differences and unscaled features enter that sum with wildly unequal weight. Use n_init, because Lloyd’s algorithm provably reaches a local minimum and nothing more. Never pick k by minimising inertia, because Jk+1JkJ_{k+1} \leq J_k makes the minimum useless and only the shape of the marginal gains carries information. Read the cumulative explained variance to choose components, because the eigenvalues sum to the conserved total variance and the curve shows exactly what each additional axis buys. And reach for TruncatedSVD on text, because centring is what PCA needs and what sparsity cannot survive.

The deeper point is the one the field is built on. With no labels available, variance and distance are all the information the data contains about its own structure. K-means minimises distance within groups, PCA maximises variance along directions, and the variance decomposition shows those are two views of a single quantity. Once you see that, the toolkit stops being a list of algorithms and becomes a small set of consequences flowing from a single choice: measuring structure with squared Euclidean geometry.

See you soon.

View Comments (2)

Leave a Reply

Subscribe to My Newsletter

Subscribe to my email newsletter to get the latest posts delivered right to your email. Pure inspiration, zero spam.

Discover more from Discuss Data Science, Machine Learning and Analytics

Subscribe now to keep reading and get access to the full archive.

Continue reading