```{r setup,include=F}
library(knitr)
library(rmdformats)
library(DT)
library(dplyr)
library(ggplot2)
library(gridExtra)
```

# Quick Recap 

- Unsupervised Learning
  - No knowledge of outcomes / labels
  - Useful in data reduction, clustering, visualization, exploratory analysis, etc.
- Principal Components Analysis (PCA)
  - Low-dimension representation of data that preserves the most information/variability from the full dataset
  - Principal components: unit vectors defining the lines/directions that maximize variance (a.k.a. directions, rotations, loadings)
  - Scores: (one-dimensional) linear combinations of original data with respect to principal components
- Covariance method for PCA 
  - Eigendecomposition of sample covariance: $\Sigma = \frac{1}{n-1} X^T X = V\Lambda V^T$
  - Eigenvalues of $\Sigma$ are contained in $\Lambda=\text{diag}(\lambda_1,\dots,\lambda_p)$ 
  - Eigenvectors are contained in orthogonal matrix $V$ 
- Interpretations 
  - The $m$th column of $V$ (denoted as $\Phi_m$) is the loading vector for PC $m$ 
  - The scores for PC $m$ can be calculated by $w_{\cdot m} = X\Phi_m$ 
  - The amount of variability explained by PC $m$ is measured by $\lambda_m$ (i.e. the variance of $w_{\cdot m}$)
- Singular value decomposition (SVD)
  - Directly decompose $X = U \Sigma V^T$
  - Eigenvectors / directions form the columns of $V$ 
  - Scores are computed by $U\Sigma$
  - Can be more efficient than covariance method for large $p$

# Problem 1: Eigenfaces

This exercise is adapted from https://github.com/k41m4n/eigenfaces.

Images can be represented as arrays of pixels. As we saw in lecture, PCA can be used for constructing lower dimensional "compressed" representations of high dimensional matrices, which can be useful for images. In this exercise we will use a classic facial recognition dataset known as the Olivetti data. These are pictures of faces (in greyscale) taken in the early 90s at AT&T Laboratories. The data contain ten different images of each of 40 individuals, taken with different expressions, angles, lighting, etc.

##### a. Use the code below to read in the data and visualize it.

```{r p2-1,eval=T,fig.width=7,cache=T}
## read in the faces data
## should be 400 x 4096
## images in rows
## pixels in columns
faces = read.csv('./olivetti_X_lab.csv',header=T)
dim(faces)
faces[1:10,1:10]

## create an individual subject identifier
## the images are ordered by subject, with 10 images for each
subject_id = rep(1:10,each=10)

## plot the first 40 images

## function for plotting
plot_faces = function(x){
  ## make each row of data into a matrix
  ## structured like the pixels in the image
  subi = matrix(x,nrow = 64, byrow = TRUE)
  subi = t(apply(subi,2,rev))
  ## plot the image
  image(subi,col=grey(seq(0, 1, length=256)), xaxt="n", yaxt="n")
}

par(mfrow=c(4, 10))
par(mar=c(0.05, 0.05, 0.05, 0.05))
apply(faces[1:40,],1,plot_faces) 

## look at the "mean" face
par(mfrow=c(1,1))
avFace = colMeans(faces)
plot_faces(avFace)

```

##### b.  Apply PCA to visualize some of the eigenfaces. Which one do you think is most creepy? Report the minimum number of PCs, $M$, that are needed to explain 75% of the variance in the original data.
```{r p2-2,eval=T,cache=T}
## apply PCA ##
pca_faces = prcomp(faces,scale=T)

## plot eigenfaces ##
par(mfrow=c(4, 4))
par(mar=c(0.05, 0.05, 0.05, 0.05))
## prcomp() uses "rotation" to refer to the matrix of PC loadings/directions
apply(pca_faces$rotation[,1:16],2,plot_faces)

## make a screeplot
par(mfrow=c(1,1))
screeplot(pca_faces,type='lines')

## find min number of components that explain 75% variance
lambda = pca_faces$sdev^2
M = min(which(cumsum(lambda)/sum(lambda)>0.75))
print(M)
```

##### c. Reconstruct the images based on the first $M$ PCs only. Are the subjects of the images recognizeable? What happens if we change $M$ to be the minimum number of PCs that explains 95% of the variance? How many total numbers are we using to represent the compressed images with 95% of variance retained? How does this compare to the original data?

```{r p2-3,eval=T,fig.width=7}
## reconstruct images
faces_r = matrix(avFace,
                 nrow=nrow(faces),
                 ncol=ncol(faces),
                 byrow=T) + pca_faces$x[,1:M]%*%t(pca_faces$rotation[,1:M])

## plot the first 40 reconstructed images
par(mfrow=c(4, 10))
par(mar=c(0.05, 0.05, 0.05, 0.05))
apply(faces_r[1:40,],1,plot_faces) 

## increase M and repeat
M_new = min(which(cumsum(lambda)/sum(lambda)>0.95))

faces_r = matrix(avFace,nrow=nrow(faces),ncol=ncol(faces),byrow=T)+pca_faces$x[,1:M_new]%*%t(pca_faces$rotation[,1:M_new])
par(mfrow=c(4, 10))
par(mar=c(0.05, 0.05, 0.05, 0.05))
apply(faces_r[1:40,],1,plot_faces) 

## numbers used to represent compressed images
prod(dim(pca_faces$x[,1:M_new])) + prod(dim(pca_faces$rotation[,1:M_new])) + length(avFace)

## numbers used to represent original images
prod(dim(faces))

```

# Problem 2: Generating data from eigenvalues/eigenvectors

In this problem, we will use code to understand how the eigenvalues and eigenvectors of a covariance matrix impact the data distribution. The eigenvectors are the directions along which the distribution is aligned. The eigenvalues are the variances along these directions.

##### a. Run the code below to visualize eigenvalues and eigenvectors for a bivariate standard normal random variable.
```{r p3-p1-1,eval=T}
set.seed(1)
library(MASS)  # for generating multivariate normals

# Function to draw a line segment from a to b.
draw_line = function(a,b,col="red",lwd=5,lty=1) { 
  lines(c(a[1],b[1]),c(a[2],b[2]),col=col,lwd=lwd,lty=lty)
}

# Function to plot a bivariate normal with the given mean
# and covariance as defined by the eigenvectors and eigen values
#   mu = mean vector
#   Cov = covariance matrix = ULU'
plot_bvn = function(n,mu,eigval1,eigval2,eigvec1,eigvec2,return_pts=FALSE) {
  U = cbind(eigvec1, eigvec2)
  L = diag(c(eigval1,eigval2))
  Cov = U %*% L %*% t(U)
  X = mvrnorm(n, mu, Cov)
  plot(X[,1],X[,2],xlim=c(-6,6),ylim=c(-6,6),cex=0.25)
  draw_line(c(0,0), sqrt(eigval1)*eigvec1)
  draw_line(c(0,0), sqrt(eigval2)*eigvec2)
  if (return_pts==T){
      return(X)
  }
}

n = 10000  # number of samples to draw
mu = c(0,0)  # mean of bivariate normals to plot

# Standard normal
plot_bvn(n,mu,eigval1=1,eigval2=1,eigvec1=c(1,0),eigvec2=c(0,1))

```

##### b. Using the code below, adjust the eigenvalues and observe what happens. What do the black dots in the plots represent? What do the red lines represent?

Black dots are the data points, red lines are the the scales (squareroot eigenvalues) times the corresponding eigenvectors of the covariance matrix.

```{r p3-p1-2,eval=T}
# Adjust eigenvalues (i.e., variance in direction of eigenvectors)
plot_bvn(n, mu, 2^2, 1, c(1,0), c(0,1))
plot_bvn(n, mu, 2^2, 0.5^2, c(1,0), c(0,1))
plot_bvn(n, mu, 1, 3^2, c(1,0), c(0,1))
```

##### c. Using the code below, adjust the eigenvectors and observe what happens. What choice of eigenvalues and eigenvectors would make $\mid X1-X2 \mid$ small with high probability?

To make $\mid X1-X2 \mid$ small, we must make the direction of the first eigenvector point toward a 1:1 relationship between the two variables and make its eigenvalue large and the other eigenvalue small so that there is little variability around the 1:1 line.

```{r p3-p1-3,eval=T}
# Adjust eigenvectors (i.e., directions)
plot_bvn(n, mu, 2^2, 1, c(1,0), c(0,1))
plot_bvn(n, mu, 2^2, 1, c(1,1)/sqrt(2), c(-1,1)/sqrt(2))
# (The eigenvectors must be orthogonal and have length 1.)
plot_bvn(n, mu, 2^2, 1, c(2,1)/sqrt(5), c(-1,2)/sqrt(5))
plot_bvn(n, mu, 2^2, 1, c(3,1)/sqrt(10), c(-1,3)/sqrt(10))

# make |X1-X2| small with high probability
plot_bvn(n, mu, 2^2, 0.001, c(1,1)/sqrt(2), c(-1,1)/sqrt(2))

# Adjust both the eigenvalues and eigenvectors
plot_bvn(n, mu, 2^2, 0.5^2, c(3,1)/sqrt(10), c(-1,3)/sqrt(10))
plot_bvn(n, mu, 0.8^2, 0.5^2, c(3,1)/sqrt(10), c(-1,3)/sqrt(10))
plot_bvn(n, mu, 0.8^2, 1.5^2, c(3,1)/sqrt(10), c(-1,3)/sqrt(10))

```


##### d. Use the code below to project each point in the data onto each eigenvector. Run the block of code a few times to get a visual sense of what the code is doing. What are the interpretations of the green and blue lines?

Blue line represent moving a data point onto the line represented by the eigenvector and the green line represents the projection value (i.e., the PC score) which is the distance from the origin to the blue point.

```{r p3-p1-4,eval=T}
# We can project each point onto each eigvec:
eigval1 = 2^2
eigval2 = 1
eigvec1 = c(3,1)/sqrt(10)
eigvec2 = c(-1,3)/sqrt(10)
X = plot_bvn(n, mu, eigval1, eigval2, eigvec1, eigvec2,return_pts = T)
score1 = X %*% eigvec1
score1[1:5]  # score1[i] = distance from (0,0) to projection of X[i,] along eigvec1
i = 1  # Let's look at point i
draw_line(X[i,], score1[i]*eigvec1, col="blue", lwd=5)
draw_line(c(0,0), score1[i]*eigvec1, col="green", lwd=3, lty=1)
```

# Problem 3: Covariance method of computing PCA

In the previous problem, we defined the scales and directions and generated a collection of data points with those scales along those directions. Conversely, we can take a collection of data points, and find the scales and directions (i.e., the sqrt(eigvals) and eigvecs). The directions can then be used to compute the scores. This is PCA!

##### a. The function below (pca_cov) performs PCA using the covariance method. Add a comment to each line of pca_cov to explain/justify it.
```{r p3-p2-1,eval=T}
# Do PCA using the covariance method
pca_cov = function(X) {
  ## find sample covariance matrix of X
  Cov_hat = cov(X)  
  ## do eigenvalue decomp
  eig = eigen(Cov_hat)
  # Note: eigen() returns the eigvecs/vals in sorted order, from largest to smallest eigval.
  ## extract directions/loadings from eigenvalue decomp
  directions = eig$vectors
  ## extract scales
  scales = sqrt(eig$values)
  ## compute PC scores
  scores = X %*% directions
  ## return output
  return(list(scores=scores, directions=directions, scales=scales))
}

```

##### b. Now run PCA on the data. Do the red lines here exactly match the red lines in the plot above based on the true scales and directions? How do they differ?

They are the same only up to a sign flip.

```{r p3-p2-2,eval=T}
# Run PCA and plot the PCA directions 
pca = pca_cov(X)
plot(X[,1], X[,2], xlim=c(-6,6),ylim=c(-6,6),cex=0.25)
draw_line(c(0,0),pca$scales[1]*pca$directions[,1])
draw_line(c(0,0),pca$scales[2]*pca$directions[,2])

# Plot PCA scores
plot(pca$scores[,1],pca$scores[,2], xlim=c(-6,6),ylim=c(-6,6),cex=0.25)

# Compare estimated eigvals/vecs to true eigvals/vecs
eigval1
pca$scales[1]^2
eigval2
pca$scales[2]^2
eigvec1
pca$directions[,1]
eigvec2
pca$directions[,2]

# Each direction is only determined up to a sign flip.
```
