What is an intuitive explanation for how PCA turns from a geometric problem (with distances) to a linear algebra problem (with eigenvectors)?

Question

I've read a lot about PCA, including various tutorials and questions (such as this one, this one, this one, and this one).

The geometric problem that PCA is trying to optimize is clear to me: PCA tries to find the first principal component by minimizing the reconstruction (projection) error, which simultaneously maximizes the variance of the projected data.

When I first read that, I immediately thought of something like linear regression; maybe you can solve it using gradient descent if needed.

However, then my mind was blown when I read that the optimization problem is solved by using linear algebra and finding eigenvectors and eigenvalues. I simply do not understand how this use of linear algebra comes into play.

So my question is: How can PCA turn from a geometric optimization problem to a linear algebra problem? Can someone provide an intuitive explanation?

I am not looking for an answer like this one that says "When you solve the mathematical problem of PCA, it ends up being equivalent to finding the eigenvalues and eigenvectors of the covariance matrix." Please explain why eigenvectors come out to be the principal components and why the eigenvalues come out to be variance of the data projected onto them

I am a software engineer and not a mathematician, by the way.

Note: the figure above was taken and modified from this PCA tutorial.

Answer 1

Problem statement

The geometric problem that PCA is trying to optimize is clear to me: PCA tries to find the first principal component by minimizing the reconstruction (projection) error, which simultaneously maximizes the variance of the projected data.

That's right. I explain the connection between these two formulations in my answer here (without math) or here (with math).

Let's take the second formulation: PCA is trying the find the direction such that the projection of the data on it has the highest possible variance. This direction is, by definition, called the first principal direction. We can formalize it as follows: given the covariance matrix $\mathbf C$, we are looking for a vector $\mathbf w$ having unit length, $\|\mathbf w\|=1$, such that $\mathbf w^\top \mathbf{Cw}$ is maximal.

(Just in case this is not clear: if $\mathbf X$ is the centered data matrix, then the projection is given by $\mathbf{Xw}$ and its variance is $\frac{1}{n-1}(\mathbf{Xw})^\top \cdot \mathbf{Xw} = \mathbf w^\top\cdot (\frac{1}{n-1}\mathbf X^\top\mathbf X)\cdot \mathbf w = \mathbf w^\top \mathbf{Cw}$.)

On the other hand, an eigenvector of $\mathbf C$ is, by definition, any vector $\mathbf v$ such that $\mathbf{Cv}=\lambda \mathbf v$.

It turns out that the first principal direction is given by the eigenvector with the largest eigenvalue. This is a nontrivial and surprising statement.

---

Proofs

If one opens any book or tutorial on PCA, one can find there the following almost one-line proof of the statement above. We want to maximize $\mathbf w^\top \mathbf{Cw}$ under the constraint that $\|\mathbf w\|=\mathbf w^\top \mathbf w=1$; this can be done introducing a Lagrange multiplier and maximizing $\mathbf w^\top \mathbf{Cw}-\lambda(\mathbf w^\top \mathbf w-1)$; differentiating, we obtain $\mathbf{Cw}-\lambda\mathbf w=0$, which is the eigenvector equation. We see that $\lambda$ has in fact to be the largest eigenvalue by substituting this solution into the objective function, which gives $\mathbf w^\top \mathbf{Cw}-\lambda(\mathbf w^\top \mathbf w-1) = \mathbf w^\top \mathbf{Cw} = \lambda\mathbf w^\top \mathbf{w} = \lambda$. By virtue of the fact that this objective function should be maximized, $\lambda$ must be the largest eigenvalue, QED.

This tends to be not very intuitive for most people.

A better proof (see e.g. this neat answer by @cardinal) says that because $\mathbf C$ is symmetric matrix, it is diagonal in its eigenvector basis. (This is actually called spectral theorem.) So we can choose an orthogonal basis, namely the one given by the eigenvectors, where $\mathbf C$ is diagonal and has eigenvalues $\lambda_i$ on the diagonal. In that basis, $\mathbf w^\top \mathbf{C w}$ simplifies to $\sum \lambda_i w_i^2$, or in other words the variance is given by the weighted sum of the eigenvalues. It is almost immediate that to maximize this expression one should simply take $\mathbf w = (1,0,0,\ldots, 0)$, i.e. the first eigenvector, yielding variance $\lambda_1$ (indeed, deviating from this solution and "trading" parts of the largest eigenvalue for the parts of smaller ones will only lead to smaller overall variance). Note that the value of $\mathbf w^\top \mathbf{C w}$ does not depend on the basis! Changing to the eigenvector basis amounts to a rotation, so in 2D one can imagine simply rotating a piece of paper with the scatterplot; obviously this cannot change any variances.

I think this is a very intuitive and a very useful argument, but it relies on the spectral theorem. So the real issue here I think is: what is the intuition behind the spectral theorem?

---

Spectral theorem

Take a symmetric matrix $\mathbf C$. Take its eigenvector $\mathbf w_1$ with the largest eigenvalue $\lambda_1$. Make this eigenvector the first basis vector and choose other basis vectors randomly (such that all of them are orthonormal). How will $\mathbf C$ look in this basis?

It will have $\lambda_1$ in the top-left corner, because $\mathbf w_1=(1,0,0\ldots 0)$ in this basis and $\mathbf {Cw}_1=(C_{11}, C_{21}, \ldots C_{p1})$ has to be equal to $\lambda_1\mathbf w_1 = (\lambda_1,0,0 \ldots 0)$.

By the same argument it will have zeros in the first column under the $\lambda_1$.

But because it is symmetric, it will have zeros in the first row after $\lambda_1$ as well. So it will look like that:

$$\mathbf C=\begin{pmatrix}\lambda_1 & 0 & \ldots & 0 \\ 0 & & & \\ \vdots & & & \\ 0 & & & \end{pmatrix},$$

where empty space means that there is a block of some elements there. Because the matrix is symmetric, this block will be symmetric too. So we can apply exactly the same argument to it, effectively using the second eigenvector as the second basis vector, and getting $\lambda_1$ and $\lambda_2$ on the diagonal. This can continue until $\mathbf C$ is diagonal. That is essentially the spectral theorem. (Note how it works only because $\mathbf C$ is symmetric.)

---

Here is a more abstract reformulation of exactly the same argument.

We know that $\mathbf{Cw}_1 = \lambda_1 \mathbf w_1$, so the first eigenvector defines a 1-dimensional subspace where $\mathbf C$ acts as a scalar multiplication. Let us now take any vector $\mathbf v$ orthogonal to $\mathbf w_1$. Then it is almost immediate that $\mathbf {Cv}$ is also orthogonal to $\mathbf w_1$. Indeed:

$$ \mathbf w_1^\top \mathbf{Cv} = (\mathbf w_1^\top \mathbf{Cv})^\top = \mathbf v^\top \mathbf C^\top \mathbf w_1 = \mathbf v^\top \mathbf {Cw}_1=\lambda_1 \mathbf v^\top \mathbf w_1 = \lambda_1\cdot 0 = 0.$$

This means that $\mathbf C$ acts on the whole remaining subspace orthogonal to $\mathbf w_1$ such that it stays separate from $\mathbf w_1$. This is the crucial property of symmetric matrices. So we can find the largest eigenvector there, $\mathbf w_2$, and proceed in the same manner, eventually constructing an orthonormal basis of eigenvectors.

Answer 2

This is my take on the linear algebra behind PCA. In linear algebra, one of the key theorems is the Spectral Theorem. It states if S is any symmetric n by n matrix with real coefficients, then S has n eigenvectors with all the eigenvalues being real. That means we can write $S = ADA^{-1} $ with D a diagonal matrix with positive entries. That is $ D = \mbox{diag} (\lambda_1, \lambda_2, \ldots, \lambda_n)$ and there is no harm in assuming $\lambda_1 \geq \lambda_2 \geq \ldots \geq \lambda_n$ . A is the change of basis matrix. That is, if our original basis was $x_1,x_2, \ldots, x_n$, then with respect to the basis given by $A(x_1), A(x_2), \ldots A(x_n)$ , the action of S is diagonal. This also means that the $A(x_i)$ can be considered as a orthogonal basis with $||A(x_i)|| = \lambda_i$ If our covariance matrix was for n observations of n variables, we would be done. The basis provided by the $A(x_i)$ is the PCA basis . This follows from the linear algebra facts. In essence it is true because a PCA basis is a basis of eigenvectors and there are atmost n eigenvectors of a square matrix of size n.
Of course most data matrices are not square. If X is a data matrix with n observations of p variables, then X is of size n by p. I will assume that $ n>p$ (more observations than variables) and that $rk(X) = p $ (all the variables are linearly independent). Neither assumption is necessary, but it will help with the intuition. Linear algebra has a generalization from the Spectral theorem called the singular value decomposition. For such an X it states that $ X = U \Sigma V^{t} $ with U,V orthonormal (square) matrices of size n and p and $\Sigma = (s_{ij}) $ a real diagonal matrix with only non-negative entries on the diagonal. Again we may rearrange the basis of V so that $s_{11} \geq s_{22} \geq \ldots s_{pp}> 0 $ In matrix terms, this means that $ X(v_i) = s_{ii} u_i$ if $ i \leq p$ and $ s_{ii} = 0 $ if $ i> n$ . The $ v_i$ give the PCA decomposition. More precisely $ \Sigma V^{t} $ is the PCA decomposition. Why ?Again, linear algebra says that there can only be p eigenvectors. The SVD gives new variables (given by the columns of V) that are orthogonal and have decreasing norm.

Answer 3

There is a result from 1936 by Eckart and Young (https://ccrma.stanford.edu/~dattorro/eckart%26young.1936.pdf), which states the following

$\sum_1^r d_k u_k v_k^T = arg min_{\hat{X} \epsilon M(r)} ||X-\hat{X}||_F^2$

where M(r) is the set of rank-r matrices, which basically means first r components of SVD of X gives the best low-rank matrix approximation of X and best is defined in terms of the squared Frobenius norm - the sum of squared elements of a matrix.

This is a general result for matrices and at first sight has nothing to do with data sets or dimensionality reduction.

However if you don't think of $X$ as a matrix but rather think of the columns of the matrix $X$ representing vectors of data points then $\hat{X}$ is the approximation with the minimum representation error in terms of squared error differences.

Answer 4

" which simultaneously maximizes the variance of the projected data." Have you hear of Rayleigh quotient? Maybe that's one way of seeing this. Namely the rayleigh quotient of the covariance matrix gives you the variance of the projected data. (and the wiki page explains why eigenvectors maximise the Rayleigh quotient)

Answer 5

Lagrange multipliers are fine but you don't actually need that to get a decent intuitive picture of why eigenvectors maximize the variance (the projected lengths).

So we want to find the unit length $w$ such that $\|Aw\|$ is maximal, where $A$ is the centered data matrix and $\frac{A^TA}{n} = C$ is our covariance matrix.

Since squaring is monotonically increasing over non-negative real numbers, maximizing $\|Aw\|$ is equivalent to maximizing $\|Aw\|^2 = (Aw)^TAw = w^TA^TAw = n (w^TCw)$. And we can also ignore that $n$ since we're choosing the $w$ that maximizes that and $n$ is constant, so it won't affect which $w$ maximizes the expression.

But we don't actually need to enforce the unit length constraint with a Lagrange multiplier because we can turn any non-zero vector into a unit vector by dividing by its length. So, for any $w$ of non-zero length, the vector $\frac{w}{\|w\|}$ is always unit length.

So now we just need to maximize

$$ \frac{w^T}{\|w\|}C\frac{w}{\|w\|} = \frac{w^TCw}{\|w\|^2} = \left(\frac{1}{n}\right)\frac{\|Aw\|^2}{\|w\|^2} $$

That last expression shows that this is equivalent to maximizing the ratio of the squared length of $Aw$ to the squared length of $w$, where we let $w$ be of any length. Instead of forcing $w$ to be unit-length and maximizing the numerator of that ratio (the denomitator will be 1 if $w$ is forced to be unit length), we can let $w$ be whatever length it wants and then maximize that ratio. As someone else pointed out, this ratio is called the Rayleigh Quotient.

As with lots of maximization problems, we need to find where the gradient vanishes (where the derivative is equal to zero). Before we do that with our particular multivariate case, let's derive something general about where derivatives equal zero for quotients in one dimension.

Consider the quotient $\frac{f(x)}{g(x)}$. The derivative with respect to x of this, using the product rule and chain rule (or "quotient" rule) from basic calc, we get:

$$ \frac{f'(x)}{g(x)} - \frac{f(x)g'(x)}{g(x)^2} $$

If we set this equal to zero (to find maxima and minima) and then rearrange a bit, we get

$$ \frac{f'(x)}{g'(x)} = \frac{f(x)}{g(x)} $$

So when the ratio of the rates of change equals the ratio of the current values, the derivative is zero and you're at a minimum or maximum.

Which actually makes a lot of sense when you think about it. Think informally about small changes in $f$ and $g$ that happen when you take a small step in $x$, then you'll go

$$ \frac{f(x)}{g(x)} \xrightarrow{\text{small step in x}} \frac{f(x) + \Delta f}{g(x) + \Delta g} $$

Since we're interested in the case where there's no net change, we want to know when

$$ \frac{f(x)}{g(x)} \approx \frac{f(x) + \Delta f}{g(x) + \Delta g} $$

$\approx$ because this is all informal with finite small changes instead of limits. The above is satisfied when

$$ \frac{\Delta f}{\Delta g} \approx \frac{f(x)}{g(x)} $$

If you currently have 100 oranges and 20 apples, you have 5 oranges per apple. Now you're going to add some oranges and apples. In what case will the ratio (quotient) of oranges to apples be preserved? It would be preserved when, say, you added 5 oranges and 1 apple because $\frac{100}{20} = \frac{105}{21}$. When you went from (100, 20) to (105, 21), the ratio didn't change because the ratio of the changes in quantity was equal to the ratio of the current quantities.

What we'll use is (after one more rearrangement), now using formal symbols again, the following condition:

$$ f'(x) = \frac{f(x)}{g(x)}g'(x) $$

"The instantaneous rate of change in the numerator must be equal to the rate of change in the denominator scaled by the ratio of the current values".

In our multivariate case, we want the whole gradient to be zero. That is, we want every partial derivative to be zero. Let's give a name to our numerator:

$$ f(w) = \|Aw\|^2 $$

$f$ is a multivariate function. It's a function from a vector $w$ to a scalar, $\|Aw\|^2$.

Let's make $A$ and $w$ explicit to illustrate.

$$ A = \begin{bmatrix} a & e & i \\ b & f & j \\ c & g & k \\ d & h & l \\ \end{bmatrix} $$

and

$$ w = \begin{bmatrix} x \\ y \\ z \\ \end{bmatrix} $$

If you write out $\|Aw\|^2$ explicitly and take the partial derivative with respect to $y$ for instance (notated as $f_y$), you will get

$$ \begin{align} f_y & = \frac{d}{dy}(\|Aw\|^2) \\ & = \frac{d}{dy}((ax + ey + iz)^2 + (bx + fy + jz)^2 + \dots) \\ & = 2e(ax + ey + iz) + 2f(bx + fy + jz) + \dots \\ & = 2\left<\begin{bmatrix}e & f & g & h\end{bmatrix}, Aw\right> \end{align} $$

So that's 2 times the inner product of the 2nd column of $A$ (corresponding to $y$ being in the 2nd row of $w$) with the vector $Aw$. This makes sense because, e.g., if the 2nd column is pointing in the same direction as $Aw$'s current position, you'll increase its squared length the most. If it's orthogonal, your rate will be 0 because you'll be (instantaneously) rotating $Aw$ instead of moving forward.

And let's give a name to the denominator in our quotient: $g(w) = \|w\|^2$. It's easier to get

$$ g_y = 2y $$

And we know what condition we want on each of our partial derivatives simulatenously to have the gradient vector equal to the zero vector. In the case of the partial w.r.t. $y$, that will become

$$ f_y = \frac{f(w)}{g(w)}g_y $$

Keep in mind every term there is a scalar. Plugging in $f_y$ and $g_y$, we get the condition:

$$ 2\left<\begin{bmatrix}e & f & g & h\end{bmatrix}, Aw\right> = \frac{\|Aw\|^2}{\|w\|^2} 2y $$

If we go ahead and derive partial derivatives $f_x$ and $f_z$ too, and arrange them into a column vector, the gradient, we get

$$ \nabla f = \begin{bmatrix} f_x \\ f_y \\ f_z \end{bmatrix} = \begin{bmatrix} 2\left<\begin{bmatrix}a & b & c & d\end{bmatrix}, Aw\right> \\ 2\left<\begin{bmatrix}e & f & g & h\end{bmatrix}, Aw\right> \\ 2\left<\begin{bmatrix}i & j & k & l\end{bmatrix}, Aw\right> \end{bmatrix} = 2A^TAw $$

The three partial derivatives of $f$ turn out to be equal to something we can write as a matrix product, $2A^TAw$.

Doing the same for $g$, we get

$$ \nabla g = 2w $$

Now we just need to simultaneously plug in our quotient derivative condition for all three partial derivatives, producting three simultaneous equations:

$$ 2A^TAw = \frac{\|Aw\|^2}{\|w\|^2} 2w $$

Cancelling the 2's, subbing in $C$ for $A^TA$ and letting the $n$'s cancel, we get

$$ Cw = \left(\frac{w^TCw}{w^Tw}\right)w $$

So the 3 simultaneous conditions we got from our derivative of ratios thing, one for each of the 3 partial derivatives of the expression (one for each component of $w$), produces a condition on the whole of $w$, namely that it's an eigenvector of $C$. We have a fixed ratio (the eigenvalue) scaling each partial derivative of $g$ (each component of an eigenvector) by the same amount, producing the partials of $f$ (the components of the output of the linear transformation done by $C$).

Answer 6

@amoeba gives neat formalization and proof of:

We can formalize it as follows: given the covariance matrix C, we are looking for a vector w having unit length, ‖w‖=1, such that w^TCw is maximal.

But I think there is one intuitive proof to:

It turns out that the first principal direction is given by the eigenvector with the largest eigenvalue. This is a nontrivial and surprising statement.

We can interpret w^TCw as a dot product between vector w and Cw, which is obtain by w going through transformation C:

w^TCw = ‖w‖ * ‖Cw‖ * cos(w, Cw)

Since w has fix length, to maximize w^TCw, we need:

1. maximize ‖Cw‖
2. maximize cos(w, Cw)

It turn out if we take w to be eigenvector of C with the largest eigenvalue, we can archive both simultaneously:

1. ‖Cw‖ is max, (if w deviate from this eigenvector, decomposite it along orthogonal eigenvectors, you should see ‖Cw‖ decrease.)
2. w and Cw in same direction, cos(w, Cw) = 1, max

Since eigenvectors are orthogonal, together with the other eigenvectors of C they forms a set of principal components to X.

---

proof of 1

decomposite w into orthogonal primary and secondary eigenvector v1 and v2, suppose their length is v1 and v2 respectively. we want to proof

(λ₁w)² > ((λ₁v1)² + (λ₂v2)²)

since λ₁ > λ₂, we have

((λ₁v1)² + (λ₂v2)²)

< ((λ₁v1)² + (λ₁v2)²)

= (λ₁)² * (v1² + v2²)

= (λ₁)² * w²