Principal Components Analysis
36-740/620
Lecture 3 (1 September 2026)
\[
\newcommand{\X}{\mathbf{x}}
\newcommand{\w}{\mathbf{w}}
\newcommand{\V}{\mathbf{v}}
\newcommand{\S}{\mathbf{s}}
\newcommand{\Expect}[1]{\mathbb{E}\left[ #1 \right]}
\newcommand{\Var}[1]{\mathrm{Var}\left[ #1 \right]}
\newcommand{\SampleVar}[1]{\widehat{\mathrm{Var}}\left[ #1 \right]}
\newcommand{\Cov}[1]{\mathrm{Cov}\left[ #1 \right]}
\newcommand{\TrueRegFunc}{\mu}
\newcommand{\EstRegFunc}{\widehat{\TrueRegFunc}}
\DeclareMathOperator{\tr}{tr}
\DeclareMathOperator*{\argmin}{argmin}
\DeclareMathOperator{\dof}{DoF}
\DeclareMathOperator{\det}{det}
\newcommand{\TrueNoise}{\epsilon}
\newcommand{\EstNoise}{\widehat{\TrueNoise}}
\]
In our last episode…
- Data is a vector
- Break the vector up into additive components
- Each component is a vector, which is a pattern
- Which components?
- Before: components were eigenvectors of the smoother matrix
- Made sense to work with them because we were smoothing
- Can we get the data to tell us what the “right” components are?
Finding a “principal” component
Projections
- \(\vec{x}_i \cdot \vec{w} =\) length of \(\vec{x}_i\)’s projection on to the direction of \(\vec{w}\)
- \((\vec{x}_i \cdot \vec{w})\vec{w} =\) the actual projected vector
Projections
Projections
- \(\vec{x}_i \cdot \vec{w} =\) length of \(\vec{x}_i\)’s projection on to the direction of \(\vec{w}\)
- \((\vec{x}_i \cdot \vec{w})\vec{w} =\) the actual projected vector
- \(\vec{x}_i - (\vec{x}_i \cdot \vec{w})\vec{w} =\) vector difference between data and projected vectors
- \(\|\vec{x}_i - (\vec{x}_i \cdot \vec{w})\vec{w}\|^2=\) squared error from replacing data vector with projected vector
How well does the projection approximate the original?
Do it for one vector first:
\[\begin{eqnarray}
{\|\vec{x_i} - (\vec{w}\cdot\vec{x_i})\vec{w}\|}^2
& =& \left(\vec{x_i} - (\vec{w}\cdot\vec{x_i})\vec{w}\right)\cdot\left(\vec{x_i} - (\vec{w}\cdot\vec{x_i})\vec{w}\right)\\
& = & \vec{x_i}\cdot\vec{x_i} -\vec{x_i}\cdot (\vec{w}\cdot\vec{x_i})\vec{w}\\
\nonumber & & - (\vec{w}\cdot\vec{x_i})\vec{w}\cdot\vec{x_i} + (\vec{w}\cdot\vec{x_i})\vec{w}\cdot(\vec{w}\cdot\vec{x_i})\vec{w}\\
& = & {\|\vec{x_i}\|}^2 -2(\vec{w}\cdot\vec{x_i})^2 + (\vec{w}\cdot\vec{x_i})^2\vec{w}\cdot\vec{w}\\
& = & \|\vec{x_i}\|^2 - (\vec{w}\cdot\vec{x_i})^2
\end{eqnarray}\]
(This is the Pythagorean theorem)
How well does the projection approximate the original?
Average across all the data vectors:
\[\begin{eqnarray}
MSE(\vec{w}) & = & \frac{1}{n}\sum_{i=1}^{n}{ {\|\vec{x_i} - (\vec{w}\cdot\vec{x_i})\vec{w}\|}^2}\\
& = & \frac{1}{n}\sum_{i=1}^{n}{\left( \|\vec{x_i}\|^2 -{(\vec{w}\cdot\vec{x_i})}^2 \right)}\\
& = & \frac{1}{n}\sum_{i=1}^{n}{\|\vec{x_i}\|^2} -\frac{1}{n}\sum_{i=1}^{n}{(\vec{w}\cdot\vec{x_i})^2}
\end{eqnarray}\]
- First bit doesn’t depend on \(\vec{w}\), so doesn’t matter for minimizing
- So we want to maximize \[
L(\vec{w}) = \frac{1}{n}\sum_{i=1}^{n}{{(\vec{w}\cdot\vec{x_i})}^2}
\]
Minimizing MSE is maximizing variance
\[\begin{eqnarray}
L(w) & = & \frac{1}{n}\sum_{i=1}^{n}{{(\vec{w}\cdot\vec{x_i})}^2}\\
& = & {\left(\frac{1}{n}\sum_{i=1}^{n}{\vec{x_i}\cdot\vec{w}}\right)}^2 +
\SampleVar{\vec{w}\cdot\vec{x_i}}
\end{eqnarray}\]
(Because: \(\Expect{Z^2} = (\Expect{Z})^2 + \Var{Z}\))
Centering simplifies our book-keeping (as promised): \[
\frac{1}{n}\sum_{i=1}^{n}{\vec{x_i} \cdot \vec{w}} = 0
\]
\(\therefore\) \[
L(\vec{w}) = \SampleVar{\vec{w}\cdot\vec{x_i}}
\]
Minimizing MSE is maximizing variance
The direction which gives us the best approximation of the data is the direction with the greatest variance
OK, how do we find this magic direction?
Matrix form: all the lengths of projections is \(\mathbf{x}\mathbf{w}\) \([n\times 1]\)
\[\begin{eqnarray}
\SampleVar{\vec{w}\cdot\vec{x_i}} & = & \frac{1}{n}\sum_{i}{{\left(\vec{x_i} \cdot \vec{w}\right)}^2}\\
& = & \frac{1}{n}{\left(\X \w\right)}^{T} \left(\X \w\right)\\
& = & \frac{1}{n} \w^T \X^T \X \w\\
& = & \w^T \frac{\X^T \X}{n} \w\\
\end{eqnarray}\]
- Fact: \(\V \equiv \frac{\X^T \X}{n} =\) sample covariance matrix of the vectors
- Because we’ve centered; \(v_{jk} = n^{-1}\sum_{i=1}^{n}{x_{ij} x_{ik}} =\) sample covariance of coordinate \(x^{(j)}\) with coordinate \(x^{(k)}\)
- Same reason \(\X^T \X\) shows up in formulas for linear regression
OK, how do we find this magic direction?
- We need to maximize \[\begin{equation}
\SampleVar{\vec{w}\cdot\vec{x_i}} = \w^T \V \w
\end{equation}\]
- Constraint: \(\vec{w}\) has length 1 \(\Leftrightarrow\) \(\w^T \w = 1\)
- Introduce a Lagrange multiplier \(\lambda\) \[\begin{eqnarray}
\mathcal{L}(\w,\lambda) & \equiv & \w^T\V\w - \lambda(\w^T \w -1)\\
\frac{\partial \mathcal{L}}{\partial \lambda} & = & \w^T \w -1\\
\frac{\partial \mathcal{L}}{\partial \w} & = & 2\V\w - 2\lambda\w
\end{eqnarray}\]
- Set derivatives to zero: \[\begin{eqnarray}
\w^T \w & = & 1\\
\V \w & = & \lambda \w
\end{eqnarray}\]
The magic direction is an eigenvector
\[\begin{eqnarray}
\w^T \w & = & 1\\
\V \w & = & \lambda \w
\end{eqnarray}\]
THIS IS AN EIGENVALUE/EIGENVECTOR EQUATION!
The value of the the solution is \[
\SampleVar{\vec{w}\cdot\vec{x_i}} = \w^T \V \w = \w^T \lambda \w = \lambda
\] so the maximum is the leading eigenvector of \(\V\)
About the sample covariance matrix
\(\V\) is a special matrix: symmetric and non-negative definite
- Eigenvalues are all real (b/c symmetric)
- If \(\lambda_i \neq \lambda_j\), then \(v_i \perp v_j\) (b/c symmetric)
- Eigenvectors form a basis, and are orthonormal
- Or can be chosen to be orthonormal
- All \(\lambda_i \geq 0\) (b/c non-negative definite)
\[\begin{eqnarray}
\text{Lead eigenvector of}\ \V & = & 1^{\mathrm{st}}\ \text{principal component}\\
& = & \text{Direction of maximum variance}\\
& = & \text{Best 1D approximation to the data}
\end{eqnarray}\]
Multiple principle components
- What about approximating by a plane, hyper-plane, hyper-hyper-plane, etc.?
- Intuition: take the direction of maximum variance \(\perp\) the first principal component
- Then direction of maximum variance \(\perp\) the first two principal components
- These are the eigenvectors of \(\V\), in order of decreasing \(\lambda\)
- The gory details can be found at the end of these slides
Some properties of the PCs
- The principal components are orthonormal
- \(\vec{w}_i \cdot \vec{w}_i = 1\)
- \(\vec{w}_i \cdot \vec{w}_j = 0\) (unless \(i=j\))
- Or in matrix form: \(\w^T\w = \mathbf{I}\)
- PC1 is the direction of maximum variance through the data
- That variance is \(\lambda_1\), biggest eigenvalue of \(\V\)
- PC \(i+1\) is the direction of maximum variance \(\perp\) PC1, PC2, \(\ldots\) PC \(i\)
- That variance is \(\lambda_{i+1}\)
Some properties of the eigenvalues
- All eigenvalues \(\geq 0\)
- If \(n \geq p\), then generally \(p\) non-zero eigenvalues
- If the data are exactly in a \(q\)-dimensional subspace, then exactly \(q\) non-zero eigenvalues
- If \(n < p\), at most \(n\) non-zero eigenvalues
- Two points define a line, three define a plane, …
Some properties of PCA as a whole
- If we use all \(p\) principal components, we have the eigendecomposition of \(\V\): \[
\V = \w \mathbf{\Lambda} \mathbf{w}^T
\] \(\mathbf{\Lambda}=\) diagonal matrix of eigenvalues \(\lambda_1, \ldots \lambda_p\)
- If we use all \(p\) principal components, \[\begin{eqnarray*}
\S & = & \X \w\\
\X & = & \S\w^T
\end{eqnarray*}\]
- If we use only the top \(q\) PCs, we get:
- the best rank-\(q\) approximation to \(\V\)
- the best dimension-\(q\) approximation to \(\X\)
Some properties of PC scores
- Average score on each PC \(=0\) (b/c we centered the data)
- Variance of score on PC \(i\) \(=\lambda_i\) (by construction)
- Covariance of score on PC \(i\) with score on PC \(j\) \(=0\)
\[\begin{eqnarray}
\Var{\text{scores}} & = & \frac{1}{n} \S^T \S\\
& = & \frac{1}{n} (\X\w)^T(\X\w)\\
& = & \frac{1}{n}\w^T \X^T \X \w\\
& = & \w^T \V\w ~\text{ by definition of} ~ \V\\
& = & \w^T ( \w \mathbf{\Lambda} \mathbf{w}^T) \w ~\text{by eigendecomposition}\\
& = & (\w^T \w) \mathbf{\Lambda} (\w^T\w)\\
& = & \mathbf{\Lambda}
\end{eqnarray}\]
Another way to think about PCA
- The original coordinates are correlated
- There is always another coordinate system with uncorrelated coordinates
- We’re rotating to that coordinate system
- Rotating to new coordinates \(\Rightarrow\) multiplying by an orthogonal matrix
- That matrix is \(\mathbf{w}\)
- The new coordinates are the scores
PCA can be used for any multivariate data
- Nothing in the math of PCA cares about where the data came from
- \(n\) measurements on \(p\) variables is all that matters
- Areas of application:
- Reducing multiple measurements (original idea)
- Dealing with collinearity or high-dimensional covariates in regression
- Recommendation engines (e.g. Netflix (Feuerverger, He, and Khatri 2012))
- \(n\) users ratings or engagement with \(p\) items of content, used to predict which users will like/engage with which items
- Your feed on Facebook, Twitter, etc. is a recommendation engine that’s trying to maximize your engagement
- Gene expression levels in molecular biology (Wall, Rechtsteiner, and Rocha 2003)
- Fashion (and cf. “eigendresses”)
- … and, of course, spatio-temporal data
PCA with spatial data
- \(n\) locations for \(p\) variables
- Each PC is a \(p\)-dimensional vector (in the space of variables)
- Scores are distributed over physical space
- That is, each location has a score (on each PC)
- vs. \(n\) variables at \(p\) locations
- Each PC is a spatial pattern
- One score for each original variable
- See backup for examples of doing this with the states
- These are almost interchangeable
- Again, see backups for the math
In R
prcomp is the usual built-in PCA command- Slightly complicated object returned
- Understand it with a (spatial) data example
USA, \(\approx 1977\)
Dataset pre-loaded in R:
head(state.x77)
## Population Income Illiteracy Life Exp Murder HS Grad Frost Area
## Alabama 3615 3624 2.1 69.05 15.1 41.3 20 50708
## Alaska 365 6315 1.5 69.31 11.3 66.7 152 566432
## Arizona 2212 4530 1.8 70.55 7.8 58.1 15 113417
## Arkansas 2110 3378 1.9 70.66 10.1 39.9 65 51945
## California 21198 5114 1.1 71.71 10.3 62.6 20 156361
## Colorado 2541 4884 0.7 72.06 6.8 63.9 166 103766
Principal components of the USA, \(\approx 1977\)
state.pca <- prcomp(state.x77, scale. = TRUE)
str(state.pca)
## List of 5
## $ sdev : num [1:8] 1.897 1.277 1.054 0.841 0.62 ...
## $ rotation: num [1:8, 1:8] 0.126 -0.299 0.468 -0.412 0.444 ...
## ..- attr(*, "dimnames")=List of 2
## .. ..$ : chr [1:8] "Population" "Income" "Illiteracy" "Life Exp" ...
## .. ..$ : chr [1:8] "PC1" "PC2" "PC3" "PC4" ...
## $ center : Named num [1:8] 4246.42 4435.8 1.17 70.88 7.38 ...
## ..- attr(*, "names")= chr [1:8] "Population" "Income" "Illiteracy" "Life Exp" ...
## $ scale : Named num [1:8] 4464.49 614.47 0.61 1.34 3.69 ...
## ..- attr(*, "names")= chr [1:8] "Population" "Income" "Illiteracy" "Life Exp" ...
## $ x : num [1:50, 1:8] 3.79 -1.053 0.867 2.382 0.241 ...
## ..- attr(*, "dimnames")=List of 2
## .. ..$ : chr [1:50] "Alabama" "Alaska" "Arizona" "Arkansas" ...
## .. ..$ : chr [1:8] "PC1" "PC2" "PC3" "PC4" ...
## - attr(*, "class")= chr "prcomp"
Principal components of the USA, \(\approx 1977\)
The weight/loading matrix \(\w\) gets called $rotation (why?):
signif(state.pca$rotation[, 1:2], 2)
## PC1 PC2
## Population 0.130 0.410
## Income -0.300 0.520
## Illiteracy 0.470 0.053
## Life Exp -0.410 -0.082
## Murder 0.440 0.310
## HS Grad -0.420 0.300
## Frost -0.360 -0.150
## Area -0.033 0.590
Each column is an eigenvector of \(\V\)
Principal components of the USA, \(\approx 1977\)
signif(state.pca$sdev, 2)
## [1] 1.90 1.30 1.10 0.84 0.62 0.55 0.38 0.34
Standard deviations along each principal component \(=\sqrt{\lambda_i}\)
If we keep \(k\) components, \[
R^2 = \frac{\sum_{i=1}^{k}{\lambda_i}}{\sum_{j=1}^{p}{\lambda_j}}
\]
(Denominator \(=\tr{\V}\) [why?])
Principal components of the USA, \(\approx 1977\)
signif(state.pca$x[1:10, 1:2], 2)
## PC1 PC2
## Alabama 3.80 -0.23
## Alaska -1.10 5.50
## Arizona 0.87 0.75
## Arkansas 2.40 -1.30
## California 0.24 3.50
## Colorado -2.10 0.51
## Connecticut -1.90 -0.24
## Delaware -0.42 -0.51
## Florida 1.20 1.10
## Georgia 3.30 0.11
Columns are \(\vec{x}_i \cdot \vec{w}_1\) and \(\vec{x}_i \cdot \vec{w}_2\)
PC1 is kinda southern
- size of state abbreviation \(\propto\) projection on to PC1
- coordinates = state capitols, except for AK and HI
PC1 is kinda the legacy of slavery
- Correlation of PC1 with having been a slave state in 1861 is 0.78
- Correlation of PC1 with having been in the Confederacy is 0.8
A famous spatial example
- For each of \(n\) variant genes (“alleles”):
- Measure prevalence of allele at \(p\) different locations
- Originally by blood tests, now gene sequencing
- Each observation is a fraction \(\in [0,1]\)
- Accumulated over many thousands of genes (\(n\)) and of locations (\(p\))
- Outstanding reference: Cavalli-Sforza, Menozzi, and Piazza (1994)
- or the good-parts version, Cavalli-Sforza (2000)
Some maps from Cavalli-Sforza, Menozzi, and Piazza (1993)
World PC1
(\(\approx 35\%\) of between-population variance)
Some maps from Cavalli-Sforza, Menozzi, and Piazza (1993)
World PC2
(\(\approx 18\%\) of between-population variance)
Some maps from Cavalli-Sforza, Menozzi, and Piazza (1993)
World PC3
(\(\approx 12\%\) of between-population variance)
Some maps from Cavalli-Sforza, Menozzi, and Piazza (1993)
- green = PC1, blue = PC2, red = PC3
- No sharp breaks…
- This is \(\approx 65\%\) of between-population variance
- (\(\approx 90\%\) of variance across people is within populations)
- Similar analyses on smaller geographic scales say more about history
PCA with multiple time series
- One approach: \(n\) time points for \(p\) variables
- Each PC is a \(p\)-dimensional vector \(\Rightarrow\) a pattern across variables
- E.g., across space
- Flip perspective: \(n\) variables for \(p\) time-points
- Each PC is a \(n\)-dimensional vector \(\Rightarrow\) a pattern across time
Irish wind data
Irish wind data
## year month day RPT VAL ROS KIL SHA BIR DUB CLA MUL CLO
## 1 61 1 1 15.04 14.96 13.17 9.29 13.96 9.87 13.67 10.25 10.83 12.58
## 2 61 1 2 14.71 16.88 10.83 6.50 12.62 7.67 11.50 10.04 9.79 9.67
## 3 61 1 3 18.50 16.88 12.33 10.13 11.17 6.17 11.25 8.04 8.50 7.67
## 4 61 1 4 10.58 6.63 11.75 4.58 4.54 2.88 8.63 1.79 5.83 5.88
## 5 61 1 5 13.33 13.25 11.42 6.17 10.71 8.21 11.92 6.54 10.92 10.34
## 6 61 1 6 13.21 8.12 9.96 6.67 5.37 4.50 10.67 4.42 7.17 7.50
## BEL MAL time
## 1 18.50 15.04 1961-01-01 12:00:00
## 2 17.54 13.83 1961-01-02 12:00:00
## 3 12.75 12.71 1961-01-03 12:00:00
## 4 5.46 10.88 1961-01-04 12:00:00
## 5 12.92 11.83 1961-01-05 12:00:00
## 6 8.12 13.17 1961-01-06 12:00:00
Irish wind data
Irish wind data — one time series
Irish wind data — all the time series
PCA: \(n = 6574\), \(p=12\)
wind.pca.1 <- prcomp(wind[, 4:15])
wind.pca.1$sdev
## [1] 15.149749 4.806761 3.848214 2.840283 2.796445 1.932717 1.809999
## [8] 1.559231 1.408849 1.355770 1.164033 1.079990
PC1: The eigenvector
plot(-wind.pca.1$rotation[, 1], ylim = c(0, 1))
text(1:12, -wind.pca.1$rotation[, 1], pos = 3, labels = colnames(wind)[4:15])
A pattern over space
PC1: The eigenvector
A function of space
PC1: The scores
A function of time
Try to describe the first component here
PCA with spatio-temporal data
- … we just did this!
- \(n\) time points for \(p\) locations
- PCs are spatial patterns
- We get a time series of scores for each component
- \(n\) spatial locations, each with a time series of length \(p\)
- PCs are time series / temporal patterns
- Trends, or components of trends?
- We get a score for each location for each component = how much that location participated in that component (over all time)
Interpreting PCA results
- PCs are linear combinations of the original coordinates
- \(\Rightarrow\) PCs change as you add or remove coordinates
- Put in 1000 measures of education and PC1 is education…
- Sometimes this is even what you want (Zeller and Carmines 1980)
- Very tempting to reify the PCs
- i.e., to “make them a thing”
- sometimes totally appropriate…
- sometimes not at all appropriate…
- Be very careful when the only evidence is the PCA (Glymour 1998)
- Smoothing artifacts can be deadly (Novembre and Stephens 2008)
- In essence: Yule and Slutsky strike again
PCA is exploratory analysis, not statistical inference
- We assumed no model
- Pro: That’s the best linear approximation to these data, no matter what
- Con: doesn’t tell us where the data came from, what other data will look like, or how much our results are driven by noise
- Prediction: PCA predicts nothing
- Inference: If \(\V \rightarrow \Var{X}\) then PCs \(\rightarrow\) eigenvectors of \(\Var{X}\) and \(\lambda\)s \(\rightarrow\) eigenvalues of \(\Var{X}\)
- But PCA doesn’t need this assumption
- Doesn’t tell us about uncertainty
- We will see ways to tackle this by simulation
Some alternatives to PCA
- Independent component analysis (ICA)
- PCA analyzes data into uncorrelated components
- But uncorrelated \(\neq\) independent (unless everything’s Gaussian)
- ICA tries to break data into statistically-independent additive components
- Measure dependence between components, minimize
- Good overview: Stone (2004)
- Nonlinear approximation
- PCA finds low-dimensional linear approximation
- What if the real structure has curves?
- Locally linear embedding, spectral component analysis, …
- Could-do-worse overview: Shalizi (n.d.), appendix
- Yet more alternatives…
Summing up
- PCA rotates to new, uncorrelated coordinates
- Using the first \(q\) PCs gives the best \(q\)-dimensional approximation to the data
- These are the \(q\) directions of largest variance
- We can make either the basis vectors or the scores into spatio-temporal patterns
- Interpretation needs domain knowledge, and some caution
- PCA does no inference or prediction
Backup: Historical asides
- PCA was invented by Pearson (1901)
- In the minimize-distance-to-a-plane formulation
- PCA was (independently) re-invented by Hotelling (1933a), Hotelling (1933b)
- In the maximize-variance formulation; also provided the name
- Also implied (incorrectly!) that \(X\) needs to be Gaussian
- I am not sure of the subsequent history, and in particular of
- The continuous-coordinate version of PCA is the “Karhunen-Loeve transform” (KL transform, KLT), worked out by Karhunen and Loeve in the 1940s, see Loève (1955), ch. X, sec. 34.5
- Hand-wavy version below
- Loeve didn’t know about either Pearson or Hotelling; I don’t know about Karhunen
- In many branches of the physical sciences, PCA/KLT is also called “empirical orthogonal functions”
- On Cavalli-Sforza’s life and work, see Stone and Lurquin (2005)
Backup: Karhunen-Loeve in one (hand-wavy) slide
- Given: \(X(t)\) defined on a finite domain (e.g., interval) \(D\), \(\Expect{X(t)}=0\) (centering)
- Desired \(X(t) = \sum_{i=1}^{\infty}{S_i \phi_i(t)}\) with random coefficients \(S_i\)
- Implies \(\Expect{S_i}=0\)
- Require orthonormality of the functions: \(\int_{D}{\phi_i(t) \phi_j(t) dt} = \delta_{ij}\)
- Implies \(\Cov{X(t), X(s)} = \sum_{i=1}^{\infty}{\sum_{j=1}^{\infty}}{\phi_i(t) \phi_j(s) \Cov{S_i, S_j}}\)
- Require minimal mean squared error, so we pick the \(\phi_i\) to minimize \(\Expect{\int_{D}{\left(X(t) - \sum_{i=1}^{q}{S_i \phi_i(t)}\right)^2 dt}}\) at each \(q\), subject to the orthonormality constraint
- Implies \(\int_{D}{\Cov{X(t), X(s)} \phi_i(s) ds} = \lambda_i \phi_i(t)\)
- An eigenproblem again, but now solved by an eigenfunction, not an eigenvector
- Eigenproblems for integral operators are classic early-20th-century math (Courant and Hilbert 1953, ch. III)
- Implies \(\Cov{S_i, S_j} = \lambda_i \delta_{ij}\)
- \(\therefore\) \(\Cov{X(t), X(s)} = \sum_{i}^{\infty}{\phi_i(t) \phi_i(s) \lambda_i}\)
- We need some sort of estimate of the covariance function…
Backup: Even more abstract K-L
- Given: function \(X\) in some space with an inner product \(\langle X, Y \rangle\), with \(\Expect{X} = 0\) (centering)
- Associated norm \(\|X\|^2 = \langle X, X \rangle\)
- Desired: \(X = \sum_{i=1}^{\infty}{S_i \phi_i}\) with random coefficients \(S_i\)
- Require orthonormality, \(\langle \phi_i, \phi_j \rangle = \delta_{ij}\)
- Implies \(\Cov{X(t), X(s)} = \sum_{i}{\sum_{j}{\phi_i(t) \phi_j(s) \Cov{S_i, S_j}}}\)
- Require minimal error, minimize \(\Expect{\| X - \sum_{i=1}^{q}{S_i \phi_i}\|^2}\) over choice of \(\phi_i\)
- Define \(\gamma_t(s) \equiv \Cov{X(t), X(s)}\)
- Doing the minimization, with the Lagrange multiplier to enforce orthonormality, gives us \(\langle \gamma_t, \phi_i \rangle = \lambda_i \phi_i(t)\) which is an eigenvalue equation
- \(\langle \gamma_t, af + bg \rangle = a \langle \gamma_t, f \rangle + b \langle \gamma_t, g \rangle\) so this is a linear operator
- And we get the implication that \(\Cov{S_i, S_j} = \lambda_i \delta_{ij}\)
Backup: \(\V\) is symmetric and non-negative definite
- Symmetry: \(V_{ij} = \Cov{X^{(i)}, X^{(j)}} = \Cov{X^{(j)}, X^{(i)}} = \V_{ji}\)
- Non-negative-definiteness: Pick any \(p\)-dimensional vector \(a\), so \(a \cdot X\) is a scalar, therefore \(\Var{a \cdot X} \geq 0\) \[\begin{eqnarray}
\Var{a \cdot X} & = & a \cdot \Var{X} a\\
& = & a \cdot \V a
\end{eqnarray}\] But “\(\V\) is non-negative definite” means “\(a \cdot \V a \geq 0\) for all vectors \(a\)”
To see the matrix trick: \[\begin{eqnarray}
\Var{a \cdot X} & = & \Cov{a \cdot X, a \cdot X}\\
& = & \Cov{\sum_{i=1}^{p}{a_i X_i}, \sum_{j=1}^{p}{a_j X_j}}\\
& = & \sum_{i=1}^{p}{a_i \sum_{j=1}^{p}{\Cov{X_i, X_j} a_j}}\\
& = & \sum_{i=1}^{p}{a_i \sum_{j=1}^{p}{\V_{ij} a_j}}\\
& = & \sum_{i=1}^{p}{a_i (\V a)_i}\\
& = & a \cdot \V a
\end{eqnarray}\]
(Backup) The gory details for multiple PCs
Use \(k\) directions in a \(p\times k\) matrix \(\w\)
Require: \(\mathbf{w}^T\mathbf{w} = \mathbf{I}\), the basis vectors are orthonormal
\(\X \w =\) matrix of projection lengths \([n\times k]\)
\(\X \w \w^T =\) matrix of projected vectors \([n\times p]\)
\(\X - \X \w \w^T =\) matrix of vector residuals \([n\times p]\)
\((\X-\X\w\w^T)(\X-\X\w\w^T)^T =\) matrix of inner products of vector residuals \([n\times n]\)
\(\tr{((\X-\X\w\w^T)(\X-\X\w\w^T)^T)} =\) sum of squared errors \([1\times 1]\)
Backup: The gory details (cont’d.)
\[\begin{eqnarray}
MSE(\w) & = & \frac{1}{n} \tr{((\X-\X\w\w^T)(\X^T - \w\w^T \X^T))}\\
& = & \frac{1}{n} \tr{(\X \X^T - \X\w\w^T\X^T - \X\w\w^T\X^T + \X\w\w^T\w\w^T\X^T)}\\
& = & \frac{1}{n}\left(\tr{(\X\X^T)} - 2\tr{(\X\w\w^T\X^T)} + \tr{(\X\w\w^T\X^T)}\right)\\
& = & \frac{1}{n}\tr{(\X\X^T)} - \frac{1}{n}\tr{(\X\w\w^T\X^T)}
\end{eqnarray}\]
so maximize \(\frac{1}{n}\tr{(\X\w\w^T\X^T)}\)
Backup: The gory details (cont’d.)
“trace is cyclic” so \[
\tr{(\X\w\w^T\X^T)} = \tr{(\X^T\X\w\w^T)} = \tr{(\w^T\X^T\X\w)}
\] so we want to maximize \[
\tr{\left(\w^T \frac{\X^T \X}{n}\w\right)}
\] under the constraint \[
\w^T \w = \mathbf{I}
\]
This is the same form we saw before, so it has the same sort of solution: each column of \(\w\) must be an eigenvector of \(\V\).
Backup: The Lagrange multiplier trick
- Want to solve \[
\max_{w}{L(w)}
\] with constraint \(f(w) = c\)
- Option I: Use \(f(w)=c\) to eliminate one coordinate of \(w\)
- Then \(w=g(v,c)\) for some function \(g\) and unconstrained \(v\)
- Do \(\max_{v}{L(g(v,c))}\)
- An unconstrained problem with one less variable
- Drawback: algebra is hard!
- Drawback (more general): who wants to study operations research?
- Option II: Solve an unconstrained problem with more variables
- Drawback: Only a French mathematician could think this made sense
Backup: The Lagrange multiplier trick (cont’d.)
- Define a Lagrangian \[
\mathcal{L}(w,\lambda) \equiv L(w) - \lambda(f(w) - c)
\]
- \(=L(w)\) when the constraint holds
- \(\lambda\) is the Lagrange multiplier which enforces the constraint \(f(w)=c\)
- Now do an unconstrained optimization over \(w\) and \(\lambda\): \[
\max_{w, \lambda}{\mathcal{L}(w,\lambda)}
\]
Backup: The Lagrange multiplier trick (cont’d.)
\[
\max_{w, \lambda}{\mathcal{L}(w,\lambda)}
\]
- Take derivatives: \[\begin{eqnarray}
\frac{\partial \mathcal{L}}{\partial \lambda} & = & -(f(w)-c)\\
\frac{\partial \mathcal{L}}{\partial w} & = & \frac{\partial L}{\partial w} - \lambda\frac{\partial f}{\partial w}
\end{eqnarray}\]
- Set to 0 at the maximum: \[\begin{eqnarray}
f(w) & = & c\\
\frac{\partial L}{\partial w} & = & \lambda\frac{\partial f}{\partial w}
\end{eqnarray}\]
- We’ve automatically recovered the constraint!
- One equation per unknown \(\Rightarrow\) Solve for \(\lambda\), \(w\) together
Backup: The Lagrange multiplier trick (cont’d.)
- More equality constraints \(\Rightarrow\) more Lagrange multipliers
- Inequality constraints, \(g(w) \leq d\), are trickier
- Is the unconstrained maximum inside the feasible set?
- Yes: problem solved
- No: constrained maximum is on the boundary
- Boundary is usually \(g(w) = d\) \(\Rightarrow\) treat like an equality
- There are subtleties; sometimes we need to learn some operations research
Backup: The Lagrange multiplier trick (cont’d.)
- To think through: Why not to do this instead? \[
\max_{\lambda, w}{L(w) + \lambda(f(w)-c))^2}
\]
Backup: Some more maps from Cavalli-Sforza, Menozzi, and Piazza (1993)
- PC1 for Europe and Southwest Asia
- Agriculture starts in the Fertile Crescent
- Farmers (not just farming) spread out
Backup: Some more maps from Cavalli-Sforza, Menozzi, and Piazza (1993)
- PC3 for Europe and Southwest Asia
- Something centered on the plains (“steppes”) north of the Caucasus mountains and the Black Sea
- Archaeology: Domestication of the horse, chariots
- Linguistics: origins of the Indo-European languages
- Most likely explanation: Indo-European-speaking barbarians using horses and chariots to expand out from the steppes (Anthony 2007)
- Similar patterns for Central and South Asia
Backup: Principal Components Regression
- Say we want to linearly regress \(Y\) on \(X\), with \(X\) being \(p\)-dimensional
- If \(p\) is too big, for whatever reason, we could do PCA on \(X\) and get a \(q\)-dimensional \(S\)
- Then we regress \(Y\) on \(S\) and get our predictions that way
- This is principal components regression
- Variant: \(X=(Z,U)\) where \(Z\) is the variable we really care about and \(U\) is a big mess of controls
- Do PCA on \(U\), get \(q\)-dimensional scores \(S\), and then regress \(Y\) on \(Z\) and \(S\)
- This is less common
- Major application of both forms of PC regression: \(n < p\), we’ve measured too many variables on each observation to get a fit by ordinary least squares
- Also called high-dimensional regression
- Another major application: dealing with collinearity or multi-collinearity
- Multicollinearity means that the \(p\) predictor variables really live in a lower, \(q\)-dimensional subspace
- Could drop predictor variables until they’re no longer collinear
- Or you could just find the scores on the \(q\) non-zero principal components
- No guarantee that the directions of maximum variance through the predictors are the directions along which \(Y\) varies \(\Rightarrow\) PC regression might not be very good
- There are other ways of dealing with high-dimensional or collinear predictors…
- Ridge regression, the lasso, the elastic net, etc.
- But some recent work suggest PC regression is surprisingly good compared to those modern methods (Dhillon et al. 2013)
Backup: Orthogonal matrices
- Matrix \(\mathbf{o}\) is orthogonal iff \(\mathbf{o}^T = \mathbf{o}^{-1}\)
- \(\Leftrightarrow\) the columns of \(\mathbf{o}\) are orthonormal vectors
- Why we say “orthogonal matrix” rather than “orthonormal matrix” is lost in the mists of 19th century German mathematics
- Every rotation (around the origin) corresponds to an orthogonal matrix
- E.g. to rotate by an angle \(\theta\) in two dimensions, the matrix is \(\mathbf{o} = \left[ \begin{array}{cc} \cos{\theta} & \sin{\theta} \\ -\sin{\theta} & \cos{\theta} \end{array}\right]\)
- If that looks funny compared to what you saw in an earlier class, remember we’re writing vectors as \([1 \times p]\) matrices, so the rotation works as \(\X\mathbf{o}\)
- You can check that \(\mathbf{o}^T \mathbf{o} = \left[\begin{array}{cc} 1 & 0 \\ 0 & 1\end{array}\right] = \mathbf{I}\) (because \(\cos^2{\theta} + \sin^2{\theta} = 1\))
- Are there orthogonal matrices which aren’t rotations?
Backup: Exchanging locations and variables
Recall the states:
state.pca <- prcomp(state.x77, scale. = TRUE)
signif(state.pca$rotation[, 1:2], 2)
## PC1 PC2
## Population 0.130 0.410
## Income -0.300 0.520
## Illiteracy 0.470 0.053
## Life Exp -0.410 -0.082
## Murder 0.440 0.310
## HS Grad -0.420 0.300
## Frost -0.360 -0.150
## Area -0.033 0.590
states are locations, PCs are patterns of variables
Backup: Exchanging locations and variables
Each score is spatially distributed
Backup: Exchanging locations and variables
Turn the data on its side
state.vars.pca <- prcomp(t(scale(state.x77))) # What's t()?
length(state.vars.pca$sdev) # Why 8?
## [1] 8
head(signif(state.vars.pca$rotation[, 1:2]), 4)
## PC1 PC2
## Alabama -0.2801370 0.03161830
## Alaska 0.0147876 0.56532600
## Arizona -0.0700666 0.00872764
## Arkansas -0.1653660 0.03283480
signif(state.vars.pca$x[, 1], 2)
## Population Income Illiteracy Life Exp Murder HS Grad Frost
## -2.60 2.90 -6.80 4.90 -6.70 4.80 4.30
## Area
## -0.69
Backup: Exchanging locations and variables
The states turned on their sides…
- This is the same as the other map (up to rounding errors)
- This is no coincidence (see below)
Backup: Exchanging locations and variables = PCA of \(\X\) vs. PCA of \(\X^T\)
- Starting from \(\X\) (an \(n\times p\) matrix)
- Eigenvectors \(\vec{w}_i\), eigenvalues \(\lambda_i\)
- \(\X = \S \w^T\)
- Eigendecomposition of the variance matrix is \(n^{-1} \X^T \X = \w \mathbf{\Lambda} \w^T\) (a \(p\times p\) matrix)
- Starting from \(\X^T\) (a \(p\times n\) matrix)
- Eigenvectors \(\vec{u}_i\), eigenvalues \(\psi_i\). diagonal matrix \(\mathbf{\Psi}\)
- Eigendecomposition of this variance is \(p^{-1} \X \X^T = \mathbf{u}\mathbf{\Psi}\mathbf{u}^T\) (an \(n\times n\) matrix)
\[\begin{eqnarray}
\mathbf{u}\mathbf{\Psi}\mathbf{u}^T & = & p^{-1} \X \X^T\\
\mathbf{u}\mathbf{\Psi}\mathbf{u}^T & = & p^{-1} \S \w^T (\S \w^T)^T\\
\mathbf{u}\mathbf{\Psi}\mathbf{u}^T & = & p^{-1} \S \w^T \w \S^T\\
\mathbf{u}\mathbf{\Psi}^{1/2} \mathbf{\Psi}^{1/2}\mathbf{u}^T & = & p^{-1/2} \S \w^T \w \S^T p^{-1/2}\\
(\mathbf{u}\mathbf{\Psi}^{1/2}) (\mathbf{u}\mathbf{\Psi}^{1/2})^T & = & p^{-1/2} \S \S^T p^{-1/2}\\
\mathbf{u} & = & p^{-1/2} \mathbf{\Psi}^{-1/2} \S
\end{eqnarray}\]
New PC1 vector \(\propto\) old scores on PC1, etc.
Backup: No, really, PCA doesn’t do statistical inference
- PCA resembles a statistical model called factor analysis
- The factor model is \(\vec{X} = \mathbf{\omega} \vec{F} + \vec{\epsilon}\)
- \(\vec{X}\) is the \(p\)-dimensional random vector we observe
- \(\vec{F}\) is the \(q\)-dimensional random vector of hidden (“latent”) factors or factor scores
- Usually assume \(q < p\)
- Usually assume \(\Var{\vec{F}} = \mathbf{I}\)
- \(\mathbf{\omega}\) is a \([p\times q]\) matrix of loadings
- \(\vec{\epsilon}\) is random noise uncorrelated with \(\vec{F}\)
- This is a statistical model which can generate new data
- It’s can also make predictions:
- Distribution of new data points
- If \(\vec{X}\) has some missing values, can still estimate \(\vec{F}\) and then use that to predict the unobserved entries in \(\vec{X}\)
- This looks similar, but PCs \(\neq\) factors
- \(\mathbf{\omega} \neq \mathbf{w}\)
- Scores on PCs aren’t even estimates of factor values
- We’ll come back to factor models later, but you could do worse than to read the factor-analysis chapter of Shalizi (n.d.)
Other alternatives to PCA
- Slow feature analysis
- PCs of multiple time series are trend-ish
- Trends ought to change slowly
- SFA finds components with high correlation over time
- Forecastable component analysis (Goerg 2013)
- Finds highly predictable components
References
Anthony, David W. 2007. The Horse, the Wheel and Language: How Bronze-Age Riders from the Eurasian Steppes Shaped the Modern World. Princeton: Princeton University Press.
Cavalli-Sforza, Luigi L. 2000. Genes, Peoples, and Languages. New York: North Point Press.
———. 1994. The History and Geography of Human Genes. Princeton: Princeton University Press.
Courant, Richard, and David Hilbert. 1953. Methods of Mathematical Physics. New York: Wiley.
Dhillon, Paramveer S., Dean P. Foster, Sham M. Kakade, and Lyle H. Ungar. 2013. “A Risk Comparison of Ordinary Least Squares Vs Ridge Regression.” Journal of Machine Lerning Research 14:1505–11. http://jmlr.org/papers/v14/dhillon13a.html.
Feuerverger, Andrey, Yu He, and Shashi Khatri. 2012. “Statistical Significance of the Netflix Challenge.” Statistical Science 27:202–31. https://doi.org/10.1214/11-STS368.
Glymour, Clark. 1998. “What Went Wrong? Reflections on Science by Observation and The Bell Curve.” Philosophy of Science 65:1–32. https://doi.org/10.1086/392624.
Goerg, Georg M. 2013. “Forecastable Component Analysis (Foreca).” In Proceedings of the 30th International Conference on Machine Learning [Icml 2013], edited by Sanjoy Dasgupta and David McAllester, 28:64–72. 2. http://proceedings.mlr.press/v28/goerg13.html.
Hotelling, Harold. 1933a. “Analysis of a Complex of Statistical Variables into Principal Components [Part 1 of 2].” Journal of Educational Psychology 24:417–41. https://doi.org/10.1037/h0071325.
———. 1933b. “Analysis of a Complex of Statistical Variables into Principal Components [Part 2 of 2].” Journal of Educational Psychology 24:498–520. https://doi.org/10.1037/h0070888.
Loève, Michel. 1955. Probability Theory. 1st ed. New York: D. Van Nostrand Company.
Novembre, John, and Matthew Stephens. 2008. “Interpreting Principal Component Analyses of Spatial Population Genetic Variation.” Nature Genetics 40:646–49. https://doi.org/10.1038/ng.139.
Stone, James V. 2004. Independent Component Analysis: A Tutorial Introduction. Cambridge, Massachusetts: MIT Press.
Stone, Linda, and Paul F. Lurquin. 2005. A Genetic and Cultural Odyssey: The Life and Work of L. Luca Cavalli-Sforza. New York: Columbia University Press. https://doi.org/10.7312/ston13396.
Wall, Michael E., Andreas Rechtsteiner, and Luis M. Rocha. 2003. “Singular Value Decomposition and Principal Component Analysis.” In A Practical Approach to Microarray Data Analysis, edited by D. P. Berrar, W. Dubitsky, and M. Granzow, 91–109. Norwell, Massachusetts: Kluwer. https://arxiv.org/abs/physics/0208101.
Zeller, Richard A., and Edward G. Carmines. 1980. Measurement in the Social Sciences: The Link Between Theory and Data. Cambridge, England: Cambridge University Press.
