---
title: Linear Dimension Reduction (mostly Principal Components)
author: 36-462/662, Spring 2022
date: 24 March 2022 (Lecture 18)
output: slidy_presentation
bibliography: locusts.bib
---


```{r, include=FALSE}
# General set-up options
library(knitr)
opts_chunk$set(size="small",background="white", highlight=FALSE, cache=TRUE,
               autodep=TRUE, warning=FALSE, message=FALSE,
               tidy.opts=list(comment=FALSE))
```



\[
\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]}
\DeclareMathOperator{\tr}{tr}
\DeclareMathOperator*{\argmin}{argmin}
\]

## Previously...

- Represent items in our data base as vectors of numerical variables
- Want to look for patterns in those features
- First, or as an example, low-dimensional structure in the vectors

## What the data looks like

- Data is $\mathbf{x} = n\times p$ matrix
    + Could be multivariate ($p$-dimensional) data on $n$ items in our data base
	  * Or multivariate data from $n$ moments in time
      * Or multivariate data from $n$ spatial locations
	  * Or scalar data at $n$ time points and $p$ locations
	  * Or multivariate data from $n_S$ locations at $n_T$ moments ($n=n_S n_T$)
	  * Etc., etc.
    + Write $\vec{x}_i$ for row $i$ ($1\times p$ matrix)
- First, center the data, just to reduce book-keeping
    + i.e., subtract the mean of each column
	+ Optional: scale each column to equal variance


## Finding _the_ "principal" component

- We don't want to keep track of $p$ dimensions
- We want _one_ dimension
- We also don't want to distort the data too much

> - Pick a direction in the $p$-dimensional space, and a length-1 vector $\vec{w}$
> - What's the best $\vec{w}$?


## Projections

> - $\vec{x}_i \cdot \vec{w} =$ length of $\vec{x}_i$'s projection on to the direction of $\vec{w}$
>    + (Strictly the _signed_ length)
>    + (Remember we're requiring $\|\vec{w}\| = 1$)
> - $(\vec{x}_i \cdot \vec{w})\vec{w} =$ the actual projected vector


## 

```{r demo.sim, include=FALSE}
# Make up date for the running example
# Make an arc
demo.theta <- runif(10,min=0,max=pi/2)
demo.x <- cbind(cos(demo.theta),sin(demo.theta))
# Center the coordinates
demo.x <- scale(demo.x, center=TRUE, scale=FALSE)
```

```{r demo.plot, echo=FALSE}
# Square plot!
par(pty="s")
# Plot the points
plot(demo.x,xlab=expression(x^1),ylab=expression(x^2), xlim=c(-1,1),
     ylim=c(-1,1),
     main="Projections of data")
# Pick a direction (not a very good one), by a unit vector
demo.w <- c(cos(-3.5*pi/6), sin(-3.5*pi/6))
# Draw an arrow for the unit vector
arrows(0,0,demo.w[1],demo.w[2],col="blue")
text(demo.w[1],demo.w[2],pos=4,labels=expression(w))
# Draw a dashed line for the whole line
abline(0,b=demo.w[2]/demo.w[1],col="blue",lty="dashed")
# Get the length of the projection of each data point on to the line
projection.lengths <- demo.x %*% demo.w
# Get the actual projections
projections <- projection.lengths %*% demo.w
# Draw those as points along the line
points(projections,pch=16,col="blue")
# Draw lines from each data point to its projection
segments(x0=demo.x[,1], y0=demo.x[,2],x1=projections[,1],y1=projections[,2], col="grey")
legend("topright", legend=c("Data vectors", "Projections onto arbitrary w",
                            "Errors with arbitrary w"),
       pch=c(1, 16, NA),
       lty=c("blank", "blank", "solid"),
       col=c("black", "blue", "lightblue"), cex=0.75)
```


## How well does the projection approximate the original? (1)

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}


## How well does the projection approximate the original? (2)

- Average across all the data vectors:

\begin{eqnarray}
MSE(\vec{w}) & = & \frac{1}{n}\sum_{i=1}^{n}{\|\vec{x_i}\|^2 -{(\vec{w}\cdot\vec{x_i})}^2}\\
& = & \frac{1}{n}\left(\sum_{i=1}^{n}{\|\vec{x_i}\|^2} -\sum_{i=1}^{n}{(\vec{w}\cdot\vec{x_i})^2}\right)
\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 (1)

\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}
($\Expect{Z^2} = (\Expect{Z})^2 + \Var{Z}$)

But
\[
\frac{1}{n}\sum_{i=1}^{n}{\vec{x_i} \cdot \vec{w}} = 0
\]
(Why?)


## Minimizing MSE is maximizing variance (2)

so
\[
L(\vec{w}) = \SampleVar{\vec{w}\cdot\vec{x_i}}
\]


## Minimizing MSE is maximizing variance (3)

> The direction which gives us the best approximation of the data is the direction along which the projections have the greatest variance


## OK, how do we find this magic direction? (1)

Matrix form: all _lengths_ of projections  $=\mathbf{x}\mathbf{w}$
$[n\times 1]$

\begin{eqnarray}
\SampleVar{\vec{w}\cdot\vec{x_i}} & = & \frac{1}{n}\sum_{i=1}^{n}{{\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


## OK, how do we find this magic direction? (2)

- 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$

- Add a [**Lagrange multiplier**](http://www.stat.cmu.edu/~cshaizi/dm/22/lectures/07/lecture-07.html) $\lambda$ times the constraint equation, and take derivatives


\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!

At the solution,
\[
\SampleVar{\vec{w}\cdot\vec{x_i}} = \w^T \V \w = \w^T \lambda \w = \lambda
\]
so the maximum comes from the _leading_ eigenvector of $\V$


## The first principal component (PC1) and scores on the first principal component

\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}

- Write the leading eigenvector as $\vec{w}_1$, or, as a $p\times 1$ matrix, $\mathbf{w}_1$
- The **score** of the $i^{\mathrm{th}}$ data point on PC1 is $\vec{x}_i \cdot \vec{w}$, the (signed) length of $\vec{x}_i$'s projection on to $\vec{w}$
- The $[n\times 1]$ matrix of all scores on PC1 is $\mathbf{x} \mathbf{w}_1$
- The scores are the new feature values we'd use after dimension reduction
- The $[n\times p]$ matrix of all _approximations_ using just PC1 is $\mathbf{x} \mathbf{w}_1 \mathbf{w}_1^T$


##

```{r demo.plot.plus.pc1, echo=FALSE}
<<demo.plot>>
# Extract PC1 "by hand"; we'll see more efficient commands to use instead later
# Find the variance-covariance matrix of the data
v <- var(demo.x)
# First principal component is the leading eigenvector
pc.1 <- eigen(v)$vectors[,1]
# Draw the arrow and the dashed line, but in a different color
arrows(0,0,pc.1[1],pc.1[2],col="red")
abline(0,b=pc.1[2]/pc.1[1],col="red",lty="dashed")
text(pc.1[1], pc.1[2], pos=4, labels="PC1")
# Get the projection lengths and actual projections
projection.lengths <- demo.x %*% pc.1
projections <- projection.lengths %*% pc.1
# Draw the points
points(projections, pch=16, col="red")
segments(x0=demo.x[,1], y0=demo.x[,2],
         x1=projections[,1], y1=projections[,2],
         col="pink")
legend("bottomleft", legend=c("Projection on PC1", "Error with PC1"),
       pch=c(16, NA),
       lty=c("blank", "solid"),
       col=c("red", "pink"), cex=0.75)
```




## Multiple principle components

- What about approximating by a plane, hyper-plane, hyper-hyper-plane, etc.?
- If we want to use $q$ principal components, we'll need to make $\w$ a $[p\times q]$ matrix
  + Each column of $\w$ will be the vector in a different direction of the linear subspace we're projecting down to
  + $\X\w$ will be the $n\times q$ matrix of new feature values
- 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$
  + Proof: See the "gory details" back-up slides
- But clearly we need to know something about the eigenvalues and eigenvectors


## About the sample covariance matrix (1)

- $\V$ is a $p\times p$ matrix, so it has $p$ eigenvalues $\lambda_1, \ldots \lambda_p$, and $p$ eigenvectors, $\vec{w}_1, \ldots \vec{w}_p$
  + Some eigenvalues might be repeated
  + We can _always_ chose eigenvectors to have length 1 (normalize them)
  + We can stack the eigenvectors in to a $p\times p$ matrix $\w$ where each column is a different $\vec{w}_i$
- $\V$ is a covariance matrix, so it has two special properties;
  1. It's **symmetric**
     + Meaning: $\V = \V^T$
     + Since $v_{ij} = \Cov{X_i, X_j} = \Cov{X_j, X_i} = v_{ji}$
  2. It's **non-negative definite**
     + Meaning: $\vec{a} \cdot (\V \vec{a}) \geq 0$,  for _any_ vector $\vec{a}$
     + You're showing that all covariance matrices are non-negative definite in HW 9
     + Compatible with some entries in $\V$ being negative numbers
	 
	 
## About the sample covariance matrix (2)

- Because $\V$ is symmetric:
  + All its eigenvalues are real (not complex)
    * Somewhat nontrivial math
  + Eigenvectors with distinct eigenvalues are orthogonal: _If_ $\lambda_i \neq \lambda_j$, then $\vec{w}_i \perp \vec{w}_j$
    * See back-up for proof
  + Eigenvectors with equal eigenvalues can be _chosen_ to be orthogonal
    * See backup for proof
  + So there are $p$ orthogonal, normalized ("orthonormal") eigenvectors
  + In $p$ dimensions, any set of $p$ orthogonal vectors is a basis
  + $\Rightarrow$ any vector $\vec{a}$ is a linear combination of the eigenvectors, $\vec{a} = \sum_{j=1}^{p}{(\vec{a} \cdot \vec{w}_j) \vec{w}_j}$
- Because $\V$ is non-negative definite:
  + All $\lambda_i \geq 0$
     * For any eigenvector $\vec{w}_i$, $\vec{w}_i \cdot (\V \vec{w}_i) = \lambda_i (\vec{w}_i \cdot \vec{w}_i)$ but this must be $\geq 0$ because $\V$ is non-negative definite


## Some properties of the eigenvalues

- All eigenvalues $\geq 0$
  + Again, this is because $\V$ is a variance matrix
- In general, $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 the PC vectors

- The principal components are orthonormal
    + $\vec{w}_i \cdot \vec{w}_i = 1$
    + $\vec{w}_i \cdot \vec{w}_j = 0$ (unless $i=j$)
    + In matrix form, $\w^T\w = \mathbf{I}$
	  * Remember, $\w =$ the $p\times p$ matrix whose columns are the PC vectors
	+ Again, this is because $\V$ is symmetric
- PC1 is the direction of maximum variance through the data
    + That variance is $\lambda_1$, biggest eigenvalue of $\V$
	  * As we saw in the proof above
- PC $(i+1)$ is the direction of maximum variance $\perp$ PC1, PC2, $\ldots$ PC $i$
    + That variance is $\lambda_{i+1}$
	  * See the "gory details" backup slides



## The **eigendecomposition** of $\V$ {.smaller}

- We saw that $\V$'s eigenvectors are an orthonormal basis, so for any $\vec{a}$,
\[
\vec{a} = \sum_{j=1}^{p}{\vec{w}_j (\vec{w}_j \cdot \vec{a})}
\]
- $\therefore$ 
\[
\V \vec{a} = \sum_{j=1}^{p}{(\V \vec{w}_j)  (\vec{w}_j \cdot \vec{a})} = \sum_{j=1}^{p}{\vec{w}_j \lambda_j (\vec{w}_j \cdot \vec{a})}
\]
- Again, define $\w$ as the $p\times p$ matrix whose columns are the eigenvectors
  + So $\w^T$ is the $p\times p$ matrix whose rows are the eigenvectors
- Define $\mathbf{\Lambda}$ as the $p\times p$ diagonal matrix whose diagonal is $\lambda_1, \lambda_2, \ldots \lambda_p$
- What we've just show is that for any vector $\vec{a}$,
\[
\V \vec{a} = \w \mathbf{\Lambda} \w^T \vec{a}
\]
so
\[
\V = \w \mathbf{\Lambda} \w^T
\]
- This is the **eigendecomposition** of $\V$



## Some properties of the PC scores (1)

Define the $[n\times p]$ matrix of the **scores** of all data points on all PCs:
\[
\S \equiv \X \w
\]

- $s_{ij} =$ score of the $i^{\mathrm{th}}$ data point on PC $j$
- Average score on each PC $=0$ (b/c we centered the data)
    + Here, "average" is over the $n$ data points, not over the $p$ PCs
	+ That is, for each $j \in 1:p$, $\frac{1}{n}\sum_{i=1}^{n}{s_{ij}} = 0$
- This implies that
\begin{eqnarray}
\Var{\text{scores}} & = & \frac{1}{n} \S^T \S\\
\end{eqnarray}


## Some properties of the PC scores (2)

\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 = \w^T (\w \mathbf{\Lambda} \mathbf{w}^T) \w
& = & (\w^T \w) \mathbf{\Lambda} (\w^T\w)\\
& = & \mathbf{\Lambda}
\end{eqnarray}

- Variance of score on PC $j$ $=\lambda_j$
  + $=$ variance over the $n$ data points; each has its own score on PC $i$
  + That is, $\frac{1}{n}\sum_{i=1}^{n}{s_{ij}^2} = \lambda_j$
- Covariance of score on PC $j$ with score on PC $k$ $=0$
  + Again, covariance over the $n$ data points, each point gets its own score on both PCs
  + That is, $\frac{1}{n}\sum_{i=1}^{n}{s_{ij} s_{ik}} = 0$



## 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
\]
   + Again, $\mathbf{\Lambda}=$ diagonal matrix of eigenvalues $\lambda_1, \ldots \lambda_p$ and $\w =$ matrix whose columns are the PCs
- If we use all $p$ principal components,
\[
\X = \S\w^T
\]
- If we use only the top $q$ PCs, we get (see backup):
    + the best rank-$q$ approximation to $\V$
    + the best dimension-$q$ approximation to $\X$
	+ the best rank-$q$ approximation to $\X$




## 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
	
	
##  Interpreting PCA results

- PCs are linear combinations of the original coordinates
   + Usually a good idea to scale variables with different units of measurement first
- $\Rightarrow$ PCs change as you add or remove coordinates
    + Put in 1000 measures of education and PC1 is education...
- 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
	
	

## USA, $\approx 1977$

Dataset pre-loaded in R:

```{r}
head(state.x77)
```

- What're the units here?
- [Inflation](https://data.bls.gov/cgi-bin/cpicalc.pl?cost1=1.00&year1=202202&year2=197402): $1 in 1974 $\approx$ $ $`r signif(1/0.17,3)`$ in 2022 
    + Or $1 in 2022 $\approx$ $ $0.17$ in 1974
- Reference: US income per person in 2020 $\approx$ $ 64k in 2020 [(source)](https://data.worldbank.org/indicator/NY.GDP.PCAP.CD?locations=US&most_recent_value_desc=true)	

## Principal components of the USA, $\approx 1977$

```{r}
state.pca <- prcomp(state.x77,scale.=TRUE)
str(state.pca)
```

## The principal component vectors

The weight/loading matrix $\w$ gets called `$rotation` (why?):

```{r}
signif(state.pca$rotation[,1:2], 2)
```
Each column is an eigenvector of $\V$

> -Break for in-class exercise

> 1. What kind of state would get a large positive score on the 1st PC, and what kind of state would get a large negative score?
> 2. What kind of state would get a large positive score on the 2nd PC, and what kind of state would get a large negative score?


## Scores on the principal components

```{r}
signif(state.pca$x[1:10, 1:2], 2)
```

Columns here are $\vec{x}_i \cdot \vec{w}_1$ and $\vec{x}_i \cdot \vec{w}_2$

So, for instance,
\[
s_{\text{Alabama}, 1} = `r signif(state.pca$rot[1,1],2)` x_{\text{Alabama}, \text{Population}} + `r signif(state.pca$rot[2,1],2)` x_{\text{Alabama}, \text{Income}} + \ldots `r signif(state.pca$rot[8,1],2)` x_{\text{Alabama}, \text{Area}}
\]

(after centering and scaling the features)


## PC1 is kinda southern (1)

```{r}
signif(head(sort(state.pca$x[,1])),2)
signif(tail(sort(state.pca$x[,1])),2)
```

## PC1 is kinda southern (2)

```{r, echo=FALSE}
plot.states_scaled <- function(sizes,min.size=0.4,max.size=2,text.color="black",
                               ...) {
  plot(state.center,type="n",...)
  out.range = max.size - min.size
  in.range = max(sizes)-min(sizes)
  scaled.sizes = out.range*((sizes-min(sizes))/in.range)
  text(state.center,state.abb,cex=scaled.sizes + min.size, col=text.color)
  invisible(scaled.sizes)
}

plot.states_scaled(state.pca$x[,1],min.size=0.3,max.size=1.5,
                   xlab="longitude",ylab="latitude")
```

- size of state abbreviation $\propto$ projection on to PC1

- coordinates = state capitols, except for AK and HI

## PC1 is kinda southern (3)

```{r, echo=FALSE}
state.confed <- rep(0,times=50)
names(state.confed) <- rownames(state.x77)
state.confed[c("South Carolina", "Mississippi", "Florida", "Alabama",
               "Georgia", "Louisiana", "Texas", "Virginia", "Arkansas",
               "Tennessee", "North Carolina")] <- 1
state.slave <- state.confed
state.slave[c("Kentucky","Missouri","Maryland","Delaware","West Virginia")] <- 1

state.colors <- rep("blue", 50)
state.colors[state.confed==1] <- "grey"
state.colors[(state.slave==1) & (state.confed==0)] <- "lightblue"

plot.states_scaled(state.pca$x[,1],min.size=0.3,max.size=1.5,
                   xlab="longitude",ylab="latitude",
                   text.color=state.colors)
```

> - (Correlation of PC1 with having been in the Confederacy is `r signif(cor(state.confed, state.pca$x[,1]), 2)`)


## The eigenvalues / variances along each PC

```{r}
signif(state.pca$sdev, 2)
```

Standard deviations along each principal component $=\sqrt{\lambda_i}$

If we keep $q$ components,
\[
R^2 = \frac{\sum_{i=1}^{q}{\lambda_i}}{\sum_{j=1}^{p}{\lambda_j}}
\]

(Denominator $=\tr{\V}$ --- why?)

## Scree plot (1)

```{r}
plot(state.pca$sdev^2, xlab="PC", ylab="Variance", type="b")
```

## Scree plot (2)

```{r}
plot(cumsum(state.pca$sdev^2), xlab="Number of PCs", ylab="Cumulative variance",
     type="b", ylim=c(0, sum(state.pca$sdev^2)))
```




## PCA as exploratory analysis, not statistical inference

- We assumed no model
    + Pro: That's the best linear approximation, no matter what
    + Con: doesn't tell us 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}$
    + But PCA doesn't **need** this assumption
    + Doesn't tell us about uncertainty
- There are PCA-_like_ models (like factor models), which we'll come back to later, but those models $\neq$ PCA


	



## Some  statistical-learning / data-mining-y applications

- Latent semantic indexing
- Multidimensional scaling
- Principal components regression
- Netflix / recommender systems









## Latent semantic indexing [@Indexing-by-latent-semantic-analysis;@solution-to-platos-problem] {.smaller}

- Build bag-of-words representation of our documents, $p \approx 10^4$ or even $10^5$
    + That is, represent documents by vectors saying how often each word appears
- Do PCA
- Keep $q \ll p$ components
    + These are the "latent semantic dimensions" of the documents
- Use the _scores_ on those components as the new features for
    + Classification
    + Prediction
    + Similarity-based search
    + Analogies
    + Many, many other things which people now think you need to use complicated neural networks to do
- Worked example in the backup slides




## Multidimensional scaling (MDS) {.smaller}

- Given: High-dimensional vectors with known distances
- Desired: Low-dimensional vectors with nearly the same distances
- One approach: PCA
    + Because: triangle inequality; gory details in @Hand-Mannila-Smyth, sec. 3.7

```{r, echo=FALSE}
par(pty="s")
plot(state.pca$x[,1:2],
     xlab="PC1", ylab="PC2", type="n")
text(state.pca$x[,1:2], labels=state.abb)
```
	
	
	
## Principal components regression

- Run PCA on the features in $\mathbf{x}$
- Take the top $q$ principal components
- Then regress $y$ on the scores on PC1, PC2, $\ldots$ PC$q$
- Can lose information, but:
  + New features are uncorrelated
  + Can be a good way to deal with multicollinearity
  + Can be a good way to deal with high-dimensional, $n < p$ problems
  + Surprisingly competitive with ridge regression [@Dhillion-et-al-ridge-regression-vs-pca-regression]



## Netflix and recommender systems in general

- The data: user ratings of movies (1--5), say $n$ users by $p$ movies
    + (Really, _accounts_ rather than _users_...)
- Goal: predict how a user will rate a movie they haven't seen before
    + So Netflix can recommend movies users will like, and keep them paying
- Originally, Netflix released data $10^8$ ratings from $4.8\times 10^5$ users of $1.8 \times 10^4$ movies
  + Challenge: improve the RMSE of ratings on a testing set by 10\%
  + Released in 2006; prize of $ ${10}^6$ for first team to hit the threshold by 2011 (with smaller prizes for incremental progress)
  + The first person to use PCA got an immediate improvement of $`r signif(100*(0.9525-0.9167)/(0.9525), 3)`$% (!)
- Questions to decide before using PCA:
  + Are users features for movies, or are movies features for users?
  + What's the best way to use the new features (= scores on PCs) for prediction?
  + How do we handle missing values?
- We'll come back to this when we look at recommender systems in more detail







## Summing up

- We have multivariate data $\X$ (dimension = $[n\times p]$)
- We want the best $q$-dimensional _linear_ approximation
- Solution: **Principal components analysis**
    + Take the $q$ leading eigenvectors of $\V \equiv \frac{1}{n}\X^T \X =$ sample/empirical covariance matrix
	+ These eigenvectors = the principal components = directions of largest variance = biggest contrasts within the data
	+ Project the data on to the principal components
	+ Equivalently: rotate to new, uncorrelated coordinates
- Assemble the eigenvectors into $\w_q$ ($[p\times q]$)
    + **Scores** for the data are $\X\w_q \equiv \S_q$ ($[n\times q]$)
	+ **Approximations** are $(\X\w_q) \w_q^T = \S_q \w_q^T$ ($[n\times p]$)
- No inference or prediction, just linear approximation
- If we just want to (approximately) preserve distances, random linear projections can work really well too (see backup)







## Backup: Gory details of PCA in matrix form

- Use $q$ directions in a $p\times q$ matrix $\w$
  + i.e., each column of $\w$ will be a unit vector in $p$-dimensional space
- _Require:_ $\mathbf{w}^T\mathbf{w} = \mathbf{I}$, the basis vectors are orthonormal
- By analogy with one principal component:
  + $\X \w =$ matrix of projection lengths  $[n\times q]$
  + $\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

\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

"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: Distinct eigenvalues of a symmetric matrix have orthogonal eigenvectors

> For any symmetric matrix with eigenvalues $\lambda_1, \ldots \lambda_p$ and eigenvectors $\vec{w}_1, \ldots \vec{w}_p$, if $\lambda_i \neq \lambda_j$, then $\vec{w}_i \perp \vec{w}_j$

Proof: Remember that $\vec{w}_i \perp \vec{w}_j$ iff $\vec{w}_i \cdot \vec{w}_j = 0$, so let's get at that inner product.  Now, for any two vectors $\vec{a}$, $\vec{b}$ and square matrix $\mathbf{c}$ (not necessarily symmetric),
\[
\vec{a} \cdot (\mathbf{c} \vec{b}) = (\mathbf{c}^T \vec{a}) \cdot \vec{b}
\]
(To see this, write $\vec{a}$ and $\vec{b}$ as $p\times 1$ matrices, so we've got $\vec{a}^T \mathbf{c} \vec{b} = \vec{a}^T (\mathbf{c}^T)^T \vec{b} = (\mathbf{c}^T \vec{a})^T \vec{b}$.)
Now let's apply this to our eigenvectors:
\begin{eqnarray}
\vec{w}_i \cdot (\V \vec{w}_j) & = & (\V^T \vec{w}_i) \cdot \vec{w}_j\\
\vec{w}_i \cdot (\V \vec{w}_j) & = & (\V \vec{w}_i) \cdot \vec{w}_j
\vec{w}_i \cdot (\lambda_j \vec{w}_j) & = & (\lambda_i \vec{w}_i) \cdot \vec{w}_j\\
\lambda_j (\vec{w}_i \cdot \vec{w}_j) & = & \lambda_i (\vec{w}_i \cdot \vec{w}_j)
\end{eqnarray}
(The second line uses the assumption that $\V$ is symmetric, the third line uses the fact that $\vec{w}_i$ and $\vec{w}_j$ are both eigenvectors of $\V$, the last line uses linearity of inner products.)

Since, by assumption, $\lambda_i \neq \lambda_j$, the only way the two sides of the equation can balance is if $\vec{w}_i \cdot \vec{w}_j = 0$, as was to be shown.

## Backup: Distinct eigenvectors with equal eigenvalues can be made or chosen to be orthogonal

- If $\vec{w}_i$ and $\vec{w}_j$ are both eigenvectors of some matrix $\mathbf{c}$, and both have eigenvalue $\lambda$, then $a\vec{w}_i + b\vec{w}_j$ is also an eigenvector of $\mathbf{c}$, with eigenvalue $\lambda$:
\[
\mathbf{c}\left(a\vec{w}_i + b\vec{w}_j\right) = a(\mathbf{c}\vec{w}_i) + b(\mathbf{c}\vec{w}_j) = \lambda(a\vec{w}_i + b\vec{w}_j)
\]
- So suppose $\vec{w}_i$ and $\vec{w}_j$ are both eigenvectors of $\mathbf{c}$ with eigenvalue $\lambda$, and that $\vec{w}_i \not\perp \vec{w}_j$; we'll construct a new eigenvector, also with eigenvalue $\lambda$, which is orthogonal to $\vec{w}_i$
  + Without loss of generality ("wolog"), assume $\|\vec{w}_i\|=\|\vec{w}_j\|=1$
- Define $\vec{u} \equiv \vec{w}_j - (\vec{w}_i \cdot \vec{w}_j) \vec{w}_i$
  + Claim 1: $\vec{u}$ is an eigenvector with eigenvalue $\lambda$ (because it's a linear combination of two eigenvectors with eigenvalue $\lambda$)
  + Claim 2: $\vec{u} \perp \vec{w}_i$:
\begin{eqnarray}
\vec{u} \cdot \vec{w}_i & = & \left(\vec{w}_j - (\vec{w}_i \cdot \vec{w}_j) \vec{w}_i\right) \cdot \vec{w}_i\\
& = & \vec{w}_j \cdot \vec{w}_i - (\vec{w}_i \cdot \vec{w}_j) (\vec{w}_i \cdot \vec{w}_i)\\
& = & \vec{w}_j \cdot \vec{w}_i - \vec{w}_i \cdot \vec{w}_j\\
& = & 0
\end{eqnarray}
(using $\|\vec{w}_i\|=1$ and symmetry of inner products)
   + Claim 3: $\vec{w}_j^{\prime} = \frac{\vec{u}}{\|\vec{u}\|}$ is also an eigenvector with eigenvalue $\lambda$, and also orthogonal to $\vec{w}_i$
     * Because it's proportional to $\vec{u}$
- So we can replace $\vec{w}_j$ with $\vec{w}_j^{\prime}$ in our list of eigenvectors
- By induction, this works if we've got more than 2 non-orthogonal eigenvectors with the same eigenvalue


## Backup: Low-rank approximation

- If we use _all_ the scores, we just re-express the data:
\[
\X = \S\w^T = (\X\w) \w^T
\]
    + Remember $\w^T = \w^{-1}$
- Suppose we just keep the top $q$ PCs, so $\w_{q}$ is $p\times q$
\[
\S_{q} = \X\w_{q} ~(\text{dimension} ~ [n\times q])
\]
- Approximation to the data:
\[
\X_{q} = \S_{q} \w_{q}^T ~ ~(\text{dimension} ~ [n\times q][q\times p])
\]
- $\X_q$ is a rank-$q$ matrix, and it's the closest rank-$q$ matrix to $\X$





## Backup Example: The axis of historical complexity (1)

@Turchin-et-al-single-axis-of-complexity

```{r, include=FALSE}
soccomp <- read.csv("http://www.stat.cmu.edu/~cshalizi/dst/18/hw/03/soccomp.irep1.csv")
complexities <- 4:12 # Column numbers
```

`r length(complexities)` variable measuring social complexity, over time, for dozens of locations across the world

Dates from `r min(soccomp$Time)` to `r max(soccomp$Time)`

```{r}
summary(soccomp[,complexities], digits=2)
```

## Backup Example: The axis of historical complexity (2)

```{r}
soccomp.pca <- prcomp(soccomp[,complexities], scale=TRUE)
```

What are the parts of the return value?

```{r}
str(soccomp.pca)
```

- `sdev` = standard deviations along PCs $=\sqrt{\lambda_i}$
- `rotation` = matrix of eigenvectors $=\w$
- `x` = scores on all PCs $=\S$

## Backup Example: The axis of historical complexity (3)

```{r}
# R^2 for k components = sum(first k eigenvalues)/sum(all eigenvalues)
# cumsum(x) returns the cumulative sums along the vector x
plot(cumsum(soccomp.pca$sdev^2)/sum(soccomp.pca$sdev^2),
     xlab="Number of components used", ylab=expression(R^2),
     ylim=c(0,1))
# Add the lines indicating the desired levels
abline(h=0.75, lty="dashed")
abline(h=0.90, lty="dotted")
```

_One_ principal component already keeps more than 75% of the variance (3 components keep just under 90%)



## Backup Example: The axis of historical complexity (4)



```{r}
# Plot the coordinates/weights/loadings of PC1
  # But suppress the usual labels on the horizontal axis
plot(soccomp.pca$rotation[,1], xlab="", main="PCs of the complexity measures",
     xaxt="n", ylab="PC weight", ylim=c(-1,1), type="l")
# Label the horizontal axis by the names of the variables
  # The las option turns the axis labels on their side (so they don't overlap
  # each other)
axis(side=1, at=1:9, labels=colnames(soccomp)[complexities], las=2,
     cex.axis=0.5)
# Now add PC2 and PC3, in distinct colors
lines(soccomp.pca$rotation[,2], col="red")
lines(soccomp.pca$rotation[,3], col="blue")
# Legend
legend("bottomright", legend=c("PC1", "PC2", "PC3"), lty="solid",
       col=c("black", "red", "blue"))
# Guide-to-the-eye horizontal line at 0
abline(h=0, col="grey")
```

- PC1: Nearly equal positive weight on all variables
    + High scores: large population, large areas, large capital cities, much hierarchy, very differentiated, etc.
	+ Low scores: small _and_ unsophisticated
- PC2: negative weight on population, area, levels of hierarchy, positive weight on gov't, infrastructure, information, money
    + High scores: small but sophisticated
	+ Low scores: big but dumb

## Backup Example: The axis of historical complexity (5)

```{r}
# Add scores to data frame
soccomp$PC1 <- soccomp.pca$x[,1]
# Plot vs. time (all areas)
with(soccomp, plot(PC1 ~ Time))
```

## Backup Example: The axis of historical complexity (6)


```{r, echo=FALSE}
# Set up a grid of smaller plots
par(mfrow=c(2,3))
# We're doing the same thing to each area, so list the areas we want...
desired.NGAs <- c("Cahokia", "Kachi Plain", "Latium",
                    "Middle Yellow River Valley", "Niger Inland Delta",
                    "Susiana")
# ... and then loop over them, making a plot of PC1 against time
for (the.NGA in desired.NGAs) {
    with(soccomp[soccomp$NGA==the.NGA,],
         plot(PC1 ~ Time, ylim=c(-6, 4), main=paste(the.NGA)))
    # ylim here enforces a common vertical scale for all plots --- optional,
    # but makes them easier to compare
}
# Turn off the grid of smaller plots for succeeding code chunks
par(mfrow=c(1,1))
```

- Illinois, North India, Italy around Rome, central China, West Africa, Persia/Iran
- Spot the collapses of civilization

## Backup Example: The axis of historical complexity (7)


```{r}
plot(soccomp.pca$x[,1], soccomp.pca$x[,2], xlab="Score on PC1",
     ylab="Score on PC2")
```

- Correlation of scores between PC1 and PC2 is `r signif(cor(soccomp.pca$x[,1], soccomp.pca$x[,2]))`
- Uncorrelated but _dependent_


## Backup Example: Latent Semantic Indexing (1)

A random selection of stories from the _Times_, about either art or music

```{r nytac.setup, include=FALSE}
source("http://www.stat.cmu.edu/~cshalizi/dm/19/hw/01/nytac-and-bow.R")
art.stories <- read.directory("nyt_corpus/art")
music.stories <- read.directory("nyt_corpus/music")
art.BoW.list <- lapply(art.stories, table)
music.BoW.list <- lapply(music.stories, table)
nyt.BoW.frame <- make.BoW.frame(c(art.BoW.list, music.BoW.list),
                                row.names=c(paste("art", 1:length(art.BoW.list),
                                                  sep="."),
                                            paste("music", 1:length(music.BoW.list),
                                                  sep=".")))
nyt.nice.frame <- div.by.euc.length(idf.weight(nyt.BoW.frame))
```

- `nyt.nice.frame` holds all the stories, with some standard text-processing tweaks applied
  + (Up-weight counts of words which appear in few documents, "inverse document frequency" weighting)
  + (Normalize all feature vectors to length 1)

```{r}
dim(nyt.nice.frame)
```

- Again, The basic command for PCA in `R` is `prcomp()`:
```{r}
nyt.pca <- prcomp(nyt.nice.frame)
```




## Backup Example: Latent Semantic Indexing (2)

First component:

```{r}
nyt.latent.sem <- nyt.pca$rotation
signif(sort(nyt.latent.sem[,1],decreasing=TRUE)[1:30],2)
signif(sort(nyt.latent.sem[,1],decreasing=FALSE)[1:30],2)
```

- Music vs. art
- What's up with "p", "m", "y", personal pronouns?


## Backup Example: Latent Semantic Indexing (3)

Second component:

```{r}
signif(sort(nyt.latent.sem[,2],decreasing=TRUE)[1:30],2)
signif(sort(nyt.latent.sem[,2],decreasing=FALSE)[1:30],2)
```

## Backup Example: Latent Semantic Indexing (4)

```{r}
story.classes <- c(rep("art", times=length(art.stories)),
                   rep("music", times=length(music.stories)))
plot(nyt.pca$x[,1:2],
     pch=ifelse(story.classes=="music","m","a"),
     col=ifelse(story.classes=="music","blue","red"))
```

- Notice: each document (= data point, row of $\X$) gets its own score on every PC (I'm just using the first two for visualization)
  + What part of `prcomp()`'s return value contains the scores?
- Notice: really easy to separate using just 2 principal components!


## Bonus: PCA and distance preservation

- In PCA, we project from $\vec{x}_i$ to $\mathbf{u}\vec{x}_i$
    + $\mathbf{u} =$ projection operator for the principal components
	+ What's $\mathbf{u}$ in terms of the eigenvectors?
- We reduce from $p$ dimensions down to $q$
- PCA tries to make the projected points close to the original points, so
\[
\frac{1}{n}\sum_{i=1}^{n}{\| \vec{x}_i - \mathbf{u}\vec{x}_i\|^2}
\]
is minimized
- So distances between projected points track distances between original points
  + (By the triangle inequality)


## Bonus: Distance-preserving projections

- Suppose we just want to preserve distances
- A projection matrix $\mathbf{u}$ would be **distance preserving** or an **isometry** if
\[
\|\mathbf{u}\vec{x}_i - \mathbf{u}\vec{x}_j\|^2 = \|\vec{x}_i - \vec{x}_j\|^2
\]
- _Exact_ distance-preservation is a lot to ask
    + Unless $p > n$
- Approximate distance preservation: we say that $\mathbf{u}$ is "$\epsilon$-close to being an isometry", or an **$\epsilon$-isometry**, when
\[
(1-\epsilon) \|\vec{x}_i - \vec{x}_j\|^2 \leq \|\mathbf{u}\vec{x}_i - \mathbf{u}\vec{x}_j\|^2 \leq (1+\epsilon) \|\vec{x}_i - \vec{x}_j\|^2
\]
- A remarkable math fact: If $q = O(\frac{\log{n}}{\epsilon^2})$ then there is _always_ an $\epsilon$-isometry

## Bonus: The random projection trick

- Make $Z_{ij} \sim \mathcal{N}(0,1)$ independent standard Gaussians
- Set
\[
\mathbf{U} = \frac{1}{\sqrt{q}}\left[\begin{array}{ccc} Z_{11} & Z_{12} & \ldots & Z_{1p}\\
Z_{21} & Z_{22} & \ldots & Z_{2p}\\
\vdots & \vdots & \ldots & \vdots\\
Z_{q1} & Z_{q2} & \ldots & Z_{qp}\end{array}\right]
\]
- $\mathbf{U}\vec{x}_i$ is like projecting on to $q$ _random_ $p$-dimensional vectors
    + What is the expected squared length of those vectors?
- Now you can check that for _any_ vector $\vec{x}$,
\[
\Expect{\|\mathbf{U}\vec{x}\|^2} = \|\vec{x}\|^2
\]
- Hint: $\Expect{\|\mathbf{U}\vec{x}\|^2} = \Expect{\sum_{i=1}^{q}{\left(\frac{1}{\sqrt{q}}\sum_{j=1}^{p}{Z_{ij} x_j}\right)^2}}$ and $\|\vec{x}\|^2 = \sum_{j=1}^{p}{x_j^2}$

## Bonus: The random projection trick

- Expectation is linear and the $Z_{ij}$ are uncorrelated so
\[
\Expect{\|\mathbf{U}\vec{x}\|^2} = \frac{1}{q}\sum_{i=1}^{q}{\sum_{j=1}^{p}{x_j^2}} = \|\vec{x}\|^2
\]

- **Random projections are distance-preserving on average**

## Bonus: Random projections are nearly distance-preserving with high probability

- Pick your favorite $\epsilon > 0$ (approximation level) and $\delta > 0$ (error probability)
- **Johnson-Lindenstrauss theorem**: With probability at least $1-\delta$, $\mathbf{U}$ is an $\epsilon$-isometry for our $n$ data points if
\[
q \geq \frac{8}{\epsilon^2}\log{\frac{n}{\sqrt{\delta}}}
\]
    + E.g., for 95% confidence, need $q \geq \frac{8}{\epsilon^2}\left(\log{(n)} + `r signif(log(1/sqrt(0.05)),4)`\ldots \right)$
	+ For 95% confidence with $n=205$ (like in HW 10) and $\epsilon=0.2$, we need $q \geq `r ceiling((8/((0.2)^2))*log(205/sqrt(0.05)))`$
	+ If we had $n=2\times 10^6$, we'd need $q \geq `r ceiling((8/((0.2)^2))*log(2e6/sqrt(0.05)))`$
- The remarkable bits: $p$ doesn't matter and neither do the actual data vectors
    + For the HW 10 data, $p=`r 350*450*3`$
- You can also turn this around to see what $\epsilon$ you can achieve with a fixed budget of $q$
    + With $q = 100$, $n=205$, $\delta=0.05$, we get $\epsilon \approx `r signif(sqrt((8/100)*log(205/sqrt(0.05))),2)`$ which isn't great but is _some_ signal
    + Again, independent of $p$






## MOAR READING

- PCA goes back to @Pearson-invents-PCA
    + This makes it clear that PCA just needs the sample variance matrix for the features (what we've been writing $\V$), and that it's _like_ doing linear regression without picking out one dimension as the thing to be predicted
       * In addition to PCA, Pearson invented the correlation coefficient, the $\chi^2$ test for dependence in contingency tables (and for goodness of fit), and a lot more of modern statistics.  He did all this as part of a very racist and classist (but explicitly anti-sexist!) project to breed better human beings [@Porter-on-KPearson].
  + In the subsequent 120 years, PCA got re-invented multiple times, most influentially by @Hotelling-re-invents-PCA-1;@Hotelling-re-invents-PCA-2
    * Lots of presentations inspired by Hotelling leap right to maximizing variance, but I think it makes more sense to focus on trying to come close to the data (as in Pearson's original paper)
	* The latest re-invention is reported (credulously) in [https://www.wired.com/story/the-style-maven-astrophysicists-of-silicon-valley/]
  + One re-invention was actually an important extension, the **Karhunen-Loeve transformation** or **Karhunen-Loeve expansion** [@Loeve-probability]:
    * We used PCA to write $p$-dimensional vectors as combinations of $p$ orthogonal basis vectors, so $\vec{x}_i = \sum_{j=1}^{p}{s_{ij} \vec{w}_j}$ with $\vec{w}_j \cdot \vec{w}_k = \delta_{jk}$ and $n^{-1}\sum_{i=1}^{n}{s_{ij} s_{ik}} = \lambda_j \delta_{jk}$
	* Suppose instead our data was $n$ continuous _functions_, $x_i(t)$, say with $t \in [0,T]$.  We'd like to find a set of **basis functions** $\phi_j(t)$ where we can write $x_i(t) = \sum_{j}{s_{ij} \phi_j(t)}$ (notice the score $s_{ij}$ doesn't change over time), with $\int_{0}^{T}{\phi_j(t) \phi_k(t) dt} = \delta_{jk}$ (basis functions are orthogonal) and $n^{-1}\sum_{i=1}^{n}{s_{ij} s_{ik}} = \lambda_j \delta_{jk}$ (scores on different basis functions are uncorrelated)
	* K&L (independently) showed how to do this, and also showed how to make it work if we imagine all of the data $x_i(t)$ as being randomly drawn from some distribution over functions
	* The same trick works for functions defined over 2, 3, many, ... dimensions
  + In many fields, PCA is called "empirical orthogonal functions" (see e.g. @Eshel-spatiotemporal) because of the connection to the KL transformation
- The random projection trick relies on the Johnson-Lindenstrauss (JL) theorem, which goes back to @Johnson-Lindenstrauss
    + The theorem is often called the JL _lemma_, because it was just a lemma (intermediate result) in their original work [@Johnson-Lindenstrauss, Lemma 1, pp. 190--191]
        * Johnson and Lindenstrauss were actually solving a kind of interpolation problem, about how to extend functions defined on a limited set of points to functions defined everywhere on a continuous space, and the JL theorem was just a tool for them
    + The slickest proof of the Johnson-Lindenstrauss theorem I've seen is @Dasgupta-Gupta-on-Johnson-Lindenstrauss
        * See also [http://blog.geomblog.org/2011/11/intuitive-argument-for-jl-lemma.html] by Suresh Venkatasubramanian
    + The version of the JL theorem I've quoted (for the constants) comes from @Boucheron-Lugosi-Massart-concentration, ch. 2
  




## References

