Principal component analysis (PCA)

I want to discuss another important topic in machine learning, which is Principal Component Analysis (PCA). PCA is a statistical technique that is used for dimensionality reduction while preserving as much variance as possible in the data. It transforms the original features into a new set of uncorrelated features called principal components.

You may ask, why do we do this ?

The reason is that transforming features into uncorrelated features allows us to break down the problem into smaller pieces and reason about them independently.

If you remember from my previous post on Logistic Regression, we trained a model to predict whether a patient has breast cancer or not based on the 30 features of the Wisconsin Breast Cancer dataset. Holding a mental model of all these features interacting with each other to give a result of Malign or Benign is intractable if you are not a domain expert in the field. Our logistic regression model somehow captures those nuances but we don't know how to reason about it.

So imagine you share this trained model with your team in the organisation and you get bombarded with following questions:

What is really going on with this model? Is it really predicting breast cancer, or just memorizing the data? If we bring in a new patient, how can we be sure that the model will predict correctly?

These are all valid questions, and to answer them with some confidence, we need to understand the data better. PCA is one of the many tools to analyse the data. What it will do is show us the structure of the data the model is learning from: which directions carry the most information, and whether the two tumor classes even separate in a lower-dimensional view.

Now, I used 2 phrases - "directions which carry most information" and "whether the two classes even separate in lower dimensional view". These are not wishy-washy statements. Bear with me and I will demonstrate what I mean.

Digression: Co-Variance and Variance

Imagine we have 2 random variables $X, Y \in \R^n$ with $\mathbb{E}[X], \mathbb{E}[Y]$ as their mean values. Then their co-variance is given by

$\Sigma_{XY} = \text{Cov}(X, Y) = \mathbb{E}[X - \mathbb{E}[X]] \cdot \mathbb{E}[Y - \mathbb{E}[Y]]$

which evaluates to $\Sigma_{XY} = \mathbb{E}[XY] - \mathbb{E}[X]\mathbb{E}[Y]$

To get an intuition of the expression $\mathbb{E}[X - \mathbb{E}[X]] \cdot \mathbb{E}[Y - \mathbb{E}[Y]]$, think of it in terms of vectors, we are essentially computing an dot product. The higher the dot product, the more closely aligned they are. Hence, more closely correlated the random variables $X, Y$ are.

Now, what happens if we take $Y = X$ ?

$$ \Sigma_{XX} = \mathbb{E}[(X - \mathbb{E}[X]) \cdot (X - \mathbb{E}[X])^T]\\ = \mathbb{E}[XX^T] - \mathbb{E}[X]\mathbb{E}[X^T] $$

There are some nice properties of $\Sigma_{XX}$.

1. It is symmetric, i.e. $\Sigma_{XX} = \Sigma_{XX}^T$
2. Diagonal elements are simply $\text{Var}(X_i)$ terms.
3. Off-diagonal elements describe how each pair $(X_i, X_j)$ vary wrt each other.
4. And, $tr(\Sigma_{XX}) = \sum_{i} Var(X_i)$

Directions carrying most information

Now, imagine that we want to compress random variable $X \in \R^n$ into $X' \in \R^k$ where $k << n$ and still want to retain most of the variance. How can we do this ?

We first center $X$ by computing $Z = X - \mu$ which gives us $\mathbb{E}[Z] = 0$. (We could also scale each feature by its standard deviation, which is a good idea when features are on wildly different scales - but it doesn't change the structure of the derivation, so let's keep it simple.)

and it can be shown that, $$ \Sigma = \mathbb{E}[ZZ^T]\\ = \frac{1}{n} \sum_{i=1}^{n} Z^{(i)} (Z^{(i)})^T $$

where $Z^{(i)} \in \R^n$ is the i-th (centered) datapoint written as a column vector.

And so if we pick a unit vector direction $u \in \R^n$ and project $Z$ on it we get the projection $(Z^{(i)})^T u$ of i-th centered datapoint.

therefore, taking an average of squared projections over all datapoints gives us $$ \frac{1}{n} \sum_{i=1}^{n} ((Z^{(i)})^T u)^2 = \frac{1}{n} \sum_{i=1}^{n} ((Z^{(i)})^T u)^T ((Z^{(i)})^T u)\\ = \frac{1}{n} \sum_{i=1}^{n} u^T (Z^{(i)} (Z^{(i)})^T) u\\ = u^T \left( \frac{1}{n} \sum_{i=1}^{n} Z^{(i)} (Z^{(i)})^T \right) u\\ = u^T \Sigma u $$

So the variance captured along direction $u$ is exactly $u^T \Sigma u$. Now, given the constraint $||u||_2 = 1$, our problem becomes:

$$ u_1 = \argmax_{||u||_2 = 1} \; u^T \Sigma u $$

This is a constrained optimization problem, and the standard tool for it is the method of Lagrange multipliers. We encode the constraint $u^T u = 1$ into the objective:

$$ \mathcal{L}(u, \lambda) = u^T \Sigma u - \lambda (u^T u - 1) $$

Taking the gradient with respect to $u$ and setting it to zero:

$$ \frac{\partial \mathcal{L}}{\partial u} = 2\Sigma u - 2\lambda u = 0 \implies \Sigma u = \lambda u $$

And there it is - the direction that maximizes variance is an eigenvector of the covariance matrix, and the variance along it is the corresponding eigenvalue:

$$ u_1^T \Sigma u_1 = u_1^T (\lambda_1 u_1) = \lambda_1 $$

Since $\Sigma$ is symmetric (property 1 from before), the spectral theorem guarantees it has $n$ real, orthogonal eigenvectors $u_1, \dots, u_n$ with real eigenvalues $\lambda_1 \geq \lambda_2 \geq \dots \geq \lambda_n \geq 0$. If we sort them by eigenvalue, the first eigenvector $u_1$ is our answer: it is the direction of maximum variance, i.e. the first principal component.

Extending to $k$ dimensions

The single-direction result generalizes beautifully: instead of one direction, we want the best $k$-dimensional subspace to project onto. The intuition is that we simply repeat the same game - find the direction of maximum variance, then find the next best direction that is orthogonal to the first (so it captures variance the first one missed, rather than re-explaining it), and so on $k$ times. Because each new direction must be orthogonal to the previous ones, the variances they capture simply add up, and the optimal subspace turns out to be the one spanned by the top-$k$ eigenvectors $u_1, \dots, u_k$, capturing $\sum_{j=1}^{k} \lambda_j$ of the total $tr(\Sigma) = \sum_{j=1}^{n} \lambda_j$ variance.

This gives us a principled way to choose $k$: pick the smallest $k$ such that the explained variance ratio $\sum_{j=1}^{k} \lambda_j / \sum_{j=1}^{n} \lambda_j$ is large enough (say, 95%).

Empirical evidence

We don't need to go far to see PCA in action - we already used it, quietly, in the logistic regression post. There, we trained a model on the 30 features of the Wisconsin Breast Cancer dataset and reached ~98% validation accuracy. But we hit a wall when trying to show the decision boundary: the boundary is a hyperplane in 30-dimensional space, which no one can draw.

The way out was PCA. We projected the (normalized) data onto its top-2 principal components - the two orthogonal directions carrying the most variance - and plotted the model's predicted probability across that plane:

Three things to notice in this plot, each confirming a claim from the math above:

1. Two directions carry almost everything. The 30 features collapse into 2 coordinates and the data is still spread out and structured - the classes occupy visibly different regions. This is the "directions carrying most information" claim made concrete: the top-2 eigenvectors of the covariance matrix capture the bulk of the variance, so almost nothing important is lost in the projection.
2. The classes separate even in 2-D. Malignant and benign tumors form largely distinct clusters, with only a handful of points in the transition band. This is the second promise from the intro, kept: a linear model achieving ~98% accuracy is consistent with what we see - if the classes separate along the top-variance directions, a hyperplane has an easy job.
3. The boundary is still a straight line. Since the model is linear, its decision boundary remains a line in any linear projection - the $p = 0.5$ contour in the plot. The PCA view didn't distort the geometry; it just made it visible.

There is also a subtler lesson hiding in the data. Many of the 30 features are derived from one another (area is roughly radius$^2$, compactness is perimeter$^2$/area $- 1$), so the covariance matrix has strongly correlated off-diagonal entries. This is exactly the redundancy PCA is built to exploit: correlated features mean the variance concentrates in a few directions, which is why 30 features compress so cleanly into 2. It also explains why individual logistic regression coefficients were hard to interpret - collinear features trade off against each other - while the PCA view gives a picture that doesn't depend on that arbitrary split.

One honest caveat: PCA is unsupervised, so the top-variance directions are not guaranteed to be the class-separating directions. Here they happened to align, which is why the 2-D view is so clean - but that is a property of this dataset, not a theorem.

Running PCA ourselves

To back these observations with numbers, we ran the full PCA pipeline from the blog post on the same data - center, SVD, sort eigenvalues - and looked at what the components actually capture.

First, the explained variance. The scree plot shows how much variance each component carries individually (left) and cumulatively (right):

The numbers behind the plot:

Two lessons here. First, the variance drops off fast - the first component alone carries more than the next five combined, which is the correlated-features story from above made quantitative. Second, and just as important: reaching the 95% threshold takes 10 components, not 2. This is the honest caveat in numbers - the top-2 view is a visualization tool, not a full compression. It shows us the structure of the data, but a third of the variance still lives in the discarded directions.

Next, the projection itself - the same top-2 view as the decision boundary plot, but showing the data alone:

Finally, the loadings - which original features drive each component. Each principal component is a linear combination of all 30 features, and the eigenvector entries tell us the weights:

PC 1 is dominated by the size/shape cluster - concave points, concavity, compactness, perimeter, radius, area - which is exactly the set of mutually correlated features we flagged earlier. PC 2 picks up a genuinely different direction: fractal_dimension and smoothness, features that measure texture rather than size. This is orthogonality in action - the second component cannot re-explain what the first already captured, so it is forced to find independent structure.

As a final sanity check, we also verified the SVD shortcut numerically: eigendecomposing the covariance matrix $\Sigma = \frac{1}{n} Z^T Z$ directly gives eigenvalues matching $\lambda_j = d_j^2/n$ to within $10^{-15}$. The full notebook with all the code is available here.

The algorithm

Putting it all together, PCA is:

1. Center the data: $Z = X - \mu$ (optionally scale by per-feature standard deviation).
2. Compute the covariance matrix $\Sigma = \frac{1}{n} Z^T Z$.
3. Eigendecompose $\Sigma$ and sort eigenvectors by decreasing eigenvalue.
4. Take the top $k$ eigenvectors as columns of $U$.
5. Project: $X' = ZU \in \R^{n \times k}$.