24.10  Principal Components

Chapter 13 §10 proved that the $k$-dimensional subspace through the origin closest to a cloud of points is spanned by the leading left singular vectors, and then said plainly what it had not proved: the centering, the covariance matrix, and the reading of ${\sigma }_i^2/\sum _l{\sigma }_l^2$ as a fraction of variance. This section supplies all three — and the variance fraction turns out to be a theorem, not a convention — and then marks exactly where linear algebra stops and statistics begins. The honest answer is: more is provable than one might expect, and less than the phrase “explained variance” usually suggests.

Throughout, the field is $\mathbb{R}$. The data are $N\ge 1$ points ${\mathbf{x}}_1,\dots ,{\mathbf{x}}_N\in {\mathbb{R}}^n$, written as the columns of a matrix, and $\Vert \cdot \Vert$ is the Euclidean norm.

24.10.1. Centering

An affine subspace of ${\mathbb{R}}^n$ is a coset $A=\mathbf{p}+U=\{\mathbf{p}+\mathbf{z}:\mathbf{z}\in U\}$ of a subspace $U$, and its dimension is $\dim U$; Chapter 18 §01 calls $U$ the direction space and shows in Proposition 18.1.4 that every affine hull has this form. The subspaces of ${\mathbb{R}}^n$ are the affine subspaces containing $\mathbf{0}$.

Chapter 13 §10 fitted subspaces, which all pass through $\mathbf{0}$. That is the wrong family for data: a cloud of points far from the origin is badly described by any line through the origin, and the quantity minimized is then dominated by the distance from $\mathbf{0}$ to the cloud rather than by its shape. Widening the competition to all affine subspaces fixes this, and it costs nothing: the extra freedom is used up in one place, and we can say in advance where.

Subtract the mean once, and the affine problem becomes the linear problem Chapter 13 solved.

Two clauses deserve comment. Columns are the data points, so $n$ is the number of measured quantities and $N$ the number of records; the opposite convention is also in use, and every formula below would be transposed under it. And $\overline{\mathbf{x}}{\mathbf{1}}^{\mathrm{\top }}$ is the $n\times N$ matrix all of whose columns equal $\overline{\mathbf{x}}$, since $\mathbf{1}\in {\mathbb{R}}^N$ is the all-ones vector.

The defining property of the centered matrix is that its own mean is zero: $${\mathbf{X}}_c\mathbf{1}=\mathbf{X}\mathbf{1}-\overline{\mathbf{x}}({\mathbf{1}}^{\mathrm{\top }}\mathbf{1})=N\overline{\mathbf{x}}-N\overline{\mathbf{x}}=\mathbf{0},$$ using ${\mathbf{1}}^{\mathrm{\top }}\mathbf{1}=N$. In words, each row of ${\mathbf{X}}_c$ sums to $0$.

Change one clause and the object changes. Suppose one subtracts the mean of each column instead of forming the mean of the columns, replacing ${\mathbf{x}}_j$ by ${\mathbf{x}}_j-(\frac{1}{n}{\mathbf{1}}_n^{\mathrm{\top }}{\mathbf{x}}_j){\mathbf{1}}_n$ with ${\mathbf{1}}_n\in {\mathbb{R}}^n$. That removes the average of the measured quantities inside each record, which is a meaningful operation for some data and a meaningless one for most, and in any case is not the one below. The clause that matters is which index is averaged over.

24.10.2. The best-fitting affine subspace

Distances to an affine subspace reduce to distances to its direction space, and the reduction is Chapter 11’s Best Approximation Theorem.

Part (b) is the whole of “the best affine subspace passes through the mean”. It holds for every direction space $U$, before any optimization over $U$ has happened, which is why the affine problem costs nothing extra.

24.10.3. The sample covariance matrix

Fix a direction $\mathbf{w}\in {\mathbb{R}}^n$ with $\mathbf{w}\ne \mathbf{0}$ and look at the data through it: each point ${\mathbf{x}}_j$ has the single coordinate $t_j={\mathbf{w}}^{\mathrm{\top }}{\mathbf{x}}_j$ along $\mathbf{w}$. The mean of these numbers is $$\frac{1}{N}\sum _j{\mathbf{w}}^{\mathrm{\top }}{\mathbf{x}}_j={\mathbf{w}}^{\mathrm{\top }}\overline{\mathbf{x}},$$ and their spread about that mean is the number $$\begin{array}{rl}\frac{1}{N}\sum _{j=1}^N({\mathbf{w}}^{\mathrm{\top }}{\mathbf{x}}_j-{\mathbf{w}}^{\mathrm{\top }}\overline{\mathbf{x}}{)}^2& =\frac{1}{N}\sum _{j=1}^N({\mathbf{w}}^{\mathrm{\top }}{\mathbf{c}}_j{)}^2\\ & =\frac{1}{N}{\mathbf{w}}^{\mathrm{\top }}{\mathbf{X}}_c{\mathbf{X}}_c^{\mathrm{\top }}\mathbf{w},\end{array}$$ the last step because $\sum _j({\mathbf{w}}^{\mathrm{\top }}{\mathbf{c}}_j{)}^2={\Vert {\mathbf{X}}_c^{\mathrm{\top }}\mathbf{w}\Vert }^2={\mathbf{w}}^{\mathrm{\top }}{\mathbf{X}}_c{\mathbf{X}}_c^{\mathrm{\top }}\mathbf{w}$. The direction $\mathbf{w}$ enters only through a quadratic form, and the matrix of that form does not depend on $\mathbf{w}$. That is a recurring expression, so it gets a name.

The sample covariance matrix is the single matrix that answers “how spread out are the data in direction $\mathbf{w}$?” for every $\mathbf{w}$ at once.

Clause by clause. The factor $1/N$ makes $\mathbf{S}$ an average, so it does not grow when data are duplicated. The entry $s_{aa}$ is the ordinary variance of the $a$-th measured quantity, and $s_{ab}$ for $a\ne b$ says whether quantities $a$ and $b$ tend to deviate from their means in the same direction (positive) or in opposite directions (negative). The centering is inside the definition, not optional: without it the same formula would measure spread about $\mathbf{0}$ rather than about the data.

Note also what $\mathbf{S}$ does not depend on: the order of the columns, since the defining sum runs over all of them. It does depend on the units of each coordinate, which is the subject of the warning below.

Change the definition in one place and it stops being a covariance. Replacing ${\mathbf{X}}_c$ by $\mathbf{X}$ gives $\frac{1}{N}\mathbf{X}{\mathbf{X}}^{\mathrm{\top }}$, whose $(a,a)$ entry is the mean of $x_{aj}^2$ rather than the variance; for the data of (a) it is $\left(\begin{array}{cc}46/4& 30/4\\ 30/4& 26/4\end{array}\right)$, which is not $\mathbf{S}$, and its top eigenvector is a multiple of $(1,0.72\dots )$ rather than of $(1,1)$. The clause that fails is the subtraction of the mean, and with it every statement below.

Part (c) connects this section to Chapter 17. The question “which direction has the most spread?” is the question “where is $R_{\mathbf{S}}$ largest?”, and Proposition 17.1.1 answers it: the maximum is ${\lambda }_1(\mathbf{S})={\sigma }_1^2/N$, attained at ${\mathbf{u}}_1$. “Which direction has the most spread among those orthogonal to ${\mathbf{u}}_1,\dots ,{\mathbf{u}}_{k-1}$?” is answered the same way by Proposition 17.1.3, with maximum ${\sigma }_k^2/N$ at ${\mathbf{u}}_k$. So the leading left singular vectors are not only the best-fitting directions of Theorem 24.10.2; they are also, one at a time, the directions of greatest remaining spread. The two descriptions agree because they are the two halves of one Pythagorean identity, which is the next theorem.

24.10.4. How much variance a subspace explains

So “the first $k$ principal directions explain $(\sum _{i\le k}{\sigma }_i^2)/(\sum _i{\sigma }_i^2)$ of the variance” is a theorem about the $N$ given points, not a convention: the numerator is the variance of the projected cloud, the denominator is the variance of the original cloud, and part (b) is what makes the two add up with nothing left over.

24.10.5. Computing it, and what is not proved

Two closing matters, one computational and one about scope.

Never form ${\mathbf{X}}_c{\mathbf{X}}_c^{\mathrm{\top }}$. Everything above is stated in terms of $\mathbf{S}$, and $\mathbf{S}$ is the natural object to state it with; but forming $\mathbf{S}$ and then finding its eigenvalues is the worst available way to compute the principal directions. The reason is Theorem 16.8.4: passing from $\mathbf{A}$ to ${\mathbf{A}}^{\ast }\mathbf{A}$ squares the condition number, ${\kappa }_2({\mathbf{A}}^{\ast }\mathbf{A})={\kappa }_2(\mathbf{A}{)}^2$, and Section 3 of this chapter turns that into a count of lost digits — Theorem 24.3.2 says the product route loses about $2{\log }_{10}{\kappa }_2$ decimal digits where a backward stable factorization of $\mathbf{A}$ itself loses about ${\log }_{10}{\kappa }_2$. The same arithmetic applies with $\mathbf{A}={\mathbf{X}}_c^{\mathrm{\top }}$ — but only when that matrix has full column rank, which is the hypothesis of Theorem 16.8.4. That asks for $\mathrm{rank}{\mathbf{X}}_c=n$, and ${\mathbf{X}}_c\mathbf{1}=\mathbf{0}$ already forces $\mathrm{rank}{\mathbf{X}}_c\le N-1$, so it fails whenever $N\le n$: for three records of ten measured quantities there is nothing to apply. The count that follows needs no hypothesis on the rank, and it is the one to remember.

Here is the count in this setting. Each entry of ${\mathbf{X}}_c{\mathbf{X}}_c^{\mathrm{\top }}$ is an inner product of two rows, each of length $N$, so by Theorem 24.1.2 the computed product is ${\mathbf{X}}_c{\mathbf{X}}_c^{\mathrm{\top }}+\mathbf{\Delta }$ with ${\Vert \mathbf{\Delta }\Vert }_2$ of order $c(n,N)u{\sigma }_1^2$, where $u$ is the unit roundoff of Section 1 and $c$ is a modest dimensional factor. It is a factor and not a constant because the entrywise bound has to be aggregated over the whole $n\times n$ product, by the same count as in Section 3, where the corresponding step reads $\Vert \mathbf{C}-{\mathbf{A}}^{\ast }\mathbf{A}\Vert \le {\gamma }_m{\Vert \mathbf{A}\Vert }_F^2\le n{\gamma }_m{\Vert \mathbf{A}\Vert }^2$: the length of the inner products sits inside $\gamma$, and passing from the Frobenius norm to the spectral norm puts their number in front. Here that makes $c(n,N)$ of size about $nN$.

If the implementation forms only the upper triangle and reflects it — an assumption about the computation, not something proved here — then $\mathbf{\Delta }$ is symmetric, and Corollary 17.3.6 says a symmetric perturbation of that size moves every eigenvalue by at most ${\Vert \mathbf{\Delta }\Vert }_2$. So every ${\sigma }_i^2$ below $c(n,N)u{\sigma }_1^2$ — that is, up to the factor $\sqrt{c}$, every ${\sigma }_i$ below $\sqrt{u}{\sigma }_1$ — is indistinguishable from $0$ in the computed product, while a backward stable factorization of ${\mathbf{X}}_c$ itself keeps it. Those are exactly the singular values a low-rank fit has to decide about. The right computation therefore works on ${\mathbf{X}}_c$ directly and never builds the product.

What this chapter does not do is assemble that computation. The ingredients are here — the Householder reductions of Section 3, the shifted $\mathbf{Q}\mathbf{R}$ iteration of Section 5 — but no algorithm for the singular values of a given matrix is written down or analyzed anywhere in this book. So the paragraph above pays the parenthesis of the last clause of Chapter 13 §10’s promise, that a singular value decomposition is “never” computed by forming ${\mathbf{A}}^{\ast }\mathbf{A}$ and finding its eigenvalues; the clause itself, how one is actually computed for a large matrix, is left unpaid.

The finite population, and the boundary. Chapter 13 §10 said that the statistics needs “a probabilistic model, a reason to prefer squared error, and an argument that the sample answer says something about the population”. Here is exactly what has been delivered.

A probabilistic model, for the finite population of the data themselves. Regard the $N$ columns as a population of $N$ equally likely outcomes, each of weight $1/N$. Then “expectation” means the average $\frac{1}{N}\sum _j$, which needs no measure theory and no limit, and the mean of the population is $\overline{\mathbf{x}}$. For two real-valued quantities $f,g$ on it, with averages $\overline{f}=\frac{1}{N}\sum _{j}f({\mathbf{x}}_j)$ and $\overline{g}=\frac{1}{N}\sum _{j}g({\mathbf{x}}_j)$, the covariance of the population is $$\mathrm{Cov}(f,g):=\frac{1}{N}\sum _{j=1}^N(f({\mathbf{x}}_j)-\overline{f})(g({\mathbf{x}}_j)-\overline{g}),$$ and the variance of $f$ is $\mathrm{Cov}(f,f)$. With that reading, the variance of $\mathbf{x}\mapsto {\mathbf{w}}^{\mathrm{\top }}\mathbf{x}$ is ${\mathbf{w}}^{\mathrm{\top }}\mathbf{S}\mathbf{w}$, which for a unit $\mathbf{w}$ is the variance in the direction $\mathbf{w}$ of Definition 24.10.2, hence the Rayleigh quotient of Proposition 24.10.3 (c); and the covariance of $\mathbf{x}\mapsto x_a$ with $\mathbf{x}\mapsto x_b$ is $s_{ab}$, which is the displayed formula of Definition 24.10.2 read backwards. So $\mathbf{S}$ is the covariance matrix of that population, and Theorem 24.10.4 is a statement about it. This much is a theorem, and it is what was proved above.

A reason to prefer squared error. Three, none of them a proof that squared error is the right criterion, which is a modeling question. (i) It does not depend on the coordinate frame: an orthogonal change of coordinates is an isometry, so it leaves every distance $d({\mathbf{x}}_j,A)$, hence $D(A)$, unchanged. The largest entry of the error has no such property — the matrix $\mathrm{diag}(1,0)$ has largest entry $1$, while the same matrix read in the frame $\frac{1}{\sqrt{2}}(1,\pm 1)$ is $\frac{1}{2}\left(\begin{array}{cc}1& 1\\ 1& 1\end{array}\right)$, with largest entry $\frac{1}{2}$ — and Chapter 13 §10’s warning shows that it really does prefer a different answer. (ii) It splits: Theorem 24.10.4 (b) says captured plus lost equals total, exactly, which is what makes a fraction of variance meaningful at all. (iii) It has a closed-form minimizer, computable by a single factorization (Theorem 13.10.3); for the sum of unsquared distances $\sum _{j}d({\mathbf{x}}_j,A)$, nothing proved in this book produces a minimizer at all.

An argument that the sample says something about a larger population. This is not proved here, and it cannot be with the tools of this book. Such an argument requires a probability model in which the $N$ points are drawn from a distribution on ${\mathbb{R}}^n$, an expectation defined by an integral rather than a finite sum, and theorems saying that $\mathbf{S}$ converges to the distribution’s covariance and that its eigenvectors converge to the distribution’s principal directions. Every one of those belongs to probability and statistics, which this book does not develop. Chapter 13 §10’s sentence promised it; the honest correction is that linear algebra supplies the finite-population model and the reasons for squared error, and hands the inference to a subject this book does not enter.

24.10.6. Exercises

A. Check your understanding

B. Practice