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

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.

• In the long thread behind your first link, there is @amoeba's answer with animation, which explains the core thing. PCA is the rotaion of the data axes (columns) till they become uncorrelated as data vectors (variables). Such rotation matrix is found via eigendecomposition or singular value decomposition and is called eigenvector matrix. – ttnphns Jun 9 '16 at 3:50
• Besides, even if you are not a mathematician (I'm not too) you've probably heard about that linear algebra and Euclidean geometry are very intimately tied fields of maths; they are even studied together as a discipline called analytical geometry. – ttnphns Jun 9 '16 at 3:54
• optimization problem Yes PCA problem could be solved via (iterative, convergent) optimization approaches, I believe. But since it has closed form solution via maths why not use that simpler, efficient solution? – ttnphns Jun 9 '16 at 3:59
• You ask to provide an intuitive explanation. I wonder why intuitive and clear answer by amoeba, where I've linked to, won't suit you. You ask _why_ eigenvectors come out to be the principal components... Why? By definition! Eigenvectors are the principal directions of a data cloud. – ttnphns Jun 14 '16 at 1:24
• @ttnphns: I actually think that the question is reasonable. Here is how I understand it. PCA wants to find the direction of maximal variance of the projection. This direction is called (by definition) the first principal direction. On the other hand, an eigenvector of the covariance matrix $C$ is (by definition) such vector $w$ that $Cw=\lambda w$. So why is the first principal direction given by the eigenvector with the largest eigenvalue? What's the intuition here? It's certainly not by definition. I've been thinking about it, and I know how to prove it, but it's hard to explain intuitively. – amoeba Jun 14 '16 at 22:19

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.

• "Lagrange multiplier" is really clear for me. However, could you tell me why we need a unit length constraint? Thanks – Haitao Du Jun 17 '16 at 20:21
• @hxd1011 There is exactly this question here already, but briefly: that's because otherwise you can multiply $w$ by any number and $w^\top Cw$ will increase by the square of this number. So the problem becomes ill-defined: the maximum of this expression is infinite. In fact, the variance of the projection on the direction of $w$ is $w^\top Cw$ only if $w$ is unit length. – amoeba Jun 17 '16 at 20:26
• I guess $n-1$ might be a bit more familiar to most readers; I replaced it here. Thanks. – amoeba Jun 17 '16 at 20:44
• @amoeba: Thank you for the answer. I'm confused by some of your notation. You use w to indicate the unit-length vector that turns out to be the first eigenvector (principal component). When I run PCA in R (e.g. prcomp(iris[,1:4], center=T, scale=T)), I see unit-length eigenvectors with a bunch of floats like (0.521, -0.269, 0.580, 0.564). However, in your answer under "Proofs", you write It is almost immediate that to maximize this expression one should simply take w=(1,0,0,…,0), i.e. the first eigenvector. Why does the eigenvector in your proof look so well-formed like that? – stackoverflowuser2010 Jun 20 '16 at 2:06
• Hi @user58865, thanks for the nudge: I simply forgot to answer the first time. The thin is, $w^\top_1 C v$ is a scalar - it's just a number. Any number is "symmetric" :) and is equal to its transpose. Does it make sense? – amoeba Oct 8 '16 at 8:49

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.

This is my take on the linear algebra behind PCA. In linear algebra, one of the key theorems is the $\textit{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.

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

@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 wTCw 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 wTCw as a dot product between vector w and Cw, which is obtain by w going through transformation C:

wTCw = ‖w‖ * ‖Cw‖ * cos(w, Cw)

Since w has fix length, to maximize wTCw, 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

1w)2 > ((λ1v1)2 + (λ2v2)2)

since λ1 > λ2, we have

((λ1v1)2 + (λ2v2)2)

< ((λ1v1)2 + (λ1v2)2)

= (λ1)2 * (v12 + v22)

= (λ1)2 * w2