```python
%matplotlib inline
import matplotlib as mpl
import matplotlib.pyplot as plt
```


```python
import numpy as np
import numpy.linalg as la
```

# Dimension Reduction


```python
np.random.seed(123)
np.set_printoptions(3)
```

### PCA from scratch

Principal Components Analysis (PCA) basically means to find and rank all the eigenvalues and eigenvectors of a covariance matrix. This is useful because high-dimensional data (with $p$ features) may have nearly all their variation in a small number of dimensions $k$, i.e. in the subspace spanned by the eigenvectors of the covariance matrix that have the $k$ largest eigenvalues. If we project the original data into this subspace, we can have a dimension reduction (from $p$ to $k$) with hopefully little loss of information.

For zero-centered vectors,

\begin{align}
\text{Cov}(X, Y) &= \frac{\sum_{i=1}^n(X_i - \bar{X})(Y_i - \bar{Y})}{n-1} \\
  &= \frac{\sum_{i=1}^nX_iY_i}{n-1} \\
  &= \frac{XY^T}{n-1}
\end{align}

and so the covariance matrix for a data set X that has zero mean in each feature vector is just $XX^T/(n-1)$. 

In other words, we can also get the eigendecomposition of the covariance matrix from the positive semi-definite matrix $XX^T$.

We will take advantage of this when we cover the SVD later in the course.

Note: Here we are using a matrix of **row** vectors


```python
n = 100
x1, x2, x3 = np.random.normal(0, 10, (3, n))
x4 = x1 + np.random.normal(size=x1.shape)
x5 = (x1 + x2)/2 + np.random.normal(size=x1.shape)
x6 = (x1 + x2 + x3)/3 + np.random.normal(size=x1.shape)
```

#### For PCA calculations, each column is an observation


```python
xs = np.c_[x1, x2, x3, x4, x5, x6].T
```


```python
xs[:, :10]
```

#### Center each observation


```python
xc = xs - np.mean(xs, 1)[:, np.newaxis]
```


```python
xc[:, :10]
```

#### Covariance

Remember the formula for covariance

$$
\text{Cov}(X, Y) = \frac{\sum_{i=1}^n(X_i - \bar{X})(Y_i - \bar{Y})}{n-1}
$$

where $\text{Cov}(X, X)$ is the sample variance of $X$.


```python
cov = (xc @ xc.T)/(n-1)
```


```python
cov
```

#### Check


```python
np.cov(xs)
```

#### Eigendecomposition


```python
e, v = la.eigh(cov)
```


```python
idx = np.argsort(e)[::-1]
```


```python
e = e[idx]
v = v[:, idx]
```

#### Explain the magnitude of the eigenvalues

Note that $x4, x5, x6$ are linear combinations of $x1, x2, x3$ with some added noise, and hence the last 3 eigenvalues are small. 


```python
plt.stem(e)
pass
```

#### The eigenvalues and eigenvectors give a factorization of the covariance matrix


```python
v @ np.diag(e) @ v.T
```

### Geometry of PCA



### Algebra of PCA



Note that $Q^T X$ results in a new data set that is uncorrelated.


```python
m = np.zeros(2)
s = np.array([[1, 0.8], [0.8, 1]])
x = np.random.multivariate_normal(m, s, n).T
```


```python
x.shape
```

#### Calculate covariance matrix from centered observations


```python
xc = (x - x.mean(1)[:, np.newaxis])
```


```python
cov = (xc @ xc.T)/(n-1)
```


```python
cov
```

#### Find eigendecoposition


```python
e, v = la.eigh(cov)
idx = np.argsort(e)[::-1]
e = e[idx]
v = v[:, idx]
```

#### In original coordinates


```python
plt.scatter(x[0], x[1], alpha=0.5)
for e_, v_ in zip(e, v.T):
    plt.plot([0, e_*v_[0]], [0, e_*v_[1]], 'r-', lw=2)
plt.xlabel('x', fontsize=14)
plt.ylabel('y', fontsize=14)
plt.axis('square')
plt.axis([-3,3,-3,3])
pass
```

#### After change of basis


```python
yc = v.T @ xc
```


```python
plt.scatter(yc[0,:], yc[1,:], alpha=0.5)
for e_, v_ in zip(e, np.eye(2)):
    plt.plot([0, e_*v_[0]], [0, e_*v_[1]], 'r-', lw=2)
plt.xlabel('PC1', fontsize=14)
plt.ylabel('PC2', fontsize=14)
plt.axis('square')
plt.axis([-3,3,-3,3])
pass
```

#### Check


```python
from sklearn.decomposition import PCA
```


```python
pca = PCA()
```

#### Note that the PCA from scikit-learn works with feature vectors, not observation vectors


```python
z = pca.fit_transform(x.T)
```

#### Eigenvectors from PCA


```python
pca.components_
```

#### Eigenvalues from PCA


```python
pca.explained_variance_
```

#### Explained variance

This is just a consequence of the invariance of the trace under change of basis. Since the original diagnonal entries in the covariance matrix are the variances of the featrues, the sum of the eigenvalues must also be the sum of the orignal varainces. In other words, the cumulateive proportion of the top $n$ eigenvaluee is the "explained variance" of the first $n$ principal components. 


```python
e/e.sum()
```


```python
pca.explained_variance_ratio_
```

#### The principal components are identical to our home-brew version, up to a flip in direction of eigenvectors


```python
plt.scatter(z[:, 0], z[:, 1], alpha=0.5)
for e_, v_ in zip(e, np.eye(2)):
    plt.plot([0, e_*v_[0]], [0, e_*v_[1]], 'r-', lw=2)
plt.xlabel('PC1', fontsize=14)
plt.ylabel('PC2', fontsize=14)
plt.axis('square')
plt.axis([-3,3,-3,3])
pass
```


```python
plt.subplot(121)
plt.scatter(-z[:, 0], -z[:, 1], alpha=0.5)
for e_, v_ in zip(e, np.eye(2)):
    plt.plot([0, e_*v_[0]], [0, e_*v_[1]], 'r-', lw=2)
plt.xlabel('PC1', fontsize=14)
plt.ylabel('PC2', fontsize=14)
plt.axis('square')
plt.axis([-3,3,-3,3])
plt.title('Scikit-learn PCA (flipped)')

plt.subplot(122)
plt.scatter(yc[0,:], yc[1,:], alpha=0.5)
for e_, v_ in zip(e, np.eye(2)):
    plt.plot([0, e_*v_[0]], [0, e_*v_[1]], 'r-', lw=2)
plt.xlabel('PC1', fontsize=14)
plt.ylabel('PC2', fontsize=14)
plt.axis('square')
plt.axis([-3,3,-3,3])
plt.title('Homebrew PCA')
plt.tight_layout()
pass
```

bf53d179cb8d8a4e9d1565dfb8c743537eb85278
ipynb
Jupyter Notebook
notebooks/B03A_PCA.ipynb
cliburn/bios-823-2021
2a8e8f38b60c947adaeafdcee15f5396a0d01b52
MIT
2021-07-15T02:20:30.000Z
2022-01-17T20:55:29.000Z
notebooks/B03A_PCA.ipynb
cliburn/bios-823-2021
2a8e8f38b60c947adaeafdcee15f5396a0d01b52
MIT
notebooks/B03A_PCA.ipynb
cliburn/bios-823-2021
2a8e8f38b60c947adaeafdcee15f5396a0d01b52
MIT
2021-08-23T16:29:44.000Z
2022-02-10T17:28:29.000Z
Qwen/Qwen-72B
1. YES
2. YES