# Principal Component Analysis

## 1. The goal - dimensionality reduction

Principal Component Analysis (sometimes called Proper Orthogonal Decomposition) is a technique for reducing the dimensions of data. 

If we think of our data as points in a coordinate system (easy with 2-D data, but harder to visualise in n-D), what PCA does is project the data onto a new set of axes, where the basis vectors for the new axes are a **linear combination** of the original coordinate system basis vectors. 

The key to PCA is that the new basis vectors are chosen in a defined way: the basis vector of the first principle component is aligned so as to capture the maximum possible variation in the data; the second principle component captures the second largest variation in the data; and so on.






Above is a cartoon illustration. The data points have a large variation in both the x and y basis vectors; note, also, that 2 data points that differ in x will also differ in y - these changes are correlated. The principal components (right picture) are chosen so that the maximum variation in the data is obtained along the PC1 basis vector; also, there is no correlation between a change in PC2 and in PC1.

## 2. Link to the Singular Value Decomposition
Principal Component Analyis and the Singular Value Decomposition are closely related. The easiest way to illustrate this is by assuming the answer: for mean-centred data (each dimension has zero mean), the basis vectors required for PCA are the columns of $\mathbf{U}$ in the SVD.

Our data is stored in matrix $\mathbf{X}$,

\begin{equation}
\mathbf{X} = \begin{bmatrix} 
\mid & \mid &  & \mid  \\ 
\mathbf{x_1} & \mathbf{x_2} & \cdots  & \mathbf{x_m} \\
\mid & \mid &  & \mid \end{bmatrix}\quad\text{.}
\end{equation}

Each column $\mathbf{x_k}$ is a new data point, and each row is a different dimension. For the 2-D example above, $\mathbf{X}$ would have 2 rows and as many columns as we have points.

We perform the SVD on $\mathbf{X}$, so that,

\begin{equation}
\mathbf{X} = \mathbf{U}\,\mathbf{\Sigma}\,\mathbf{V^T}\quad\text{.}
\end{equation}

Now, we recall that $\mathbf{U}$ (and also $\mathbf{V}$) is orthonormal: each column is a vector of unit magnitude, and all the column vectors are orthogonal. If we want to project our original data $\mathbf{X}$ on to the basis vectors formed by the columns of $\mathbf{U}$ (to obtain a new view of our data which we will call $\mathbf{Y}$) then all we have to do is,

\begin{equation}
\mathbf{Y}=\mathbf{U^T}\mathbf{X}\quad\text{.}
\end{equation}

To make further progress, we introduce *variance* and *covariance*. If we have two samples of data, each with zero mean, stored in vector $\mathbf{a}$ and a vector $\mathbf{b}$, then the variances are:

\begin{equation}
\sigma^2_\mathbf{a} = \frac{1}{n-1}\mathbf{a}\mathbf{a^T} \\
\sigma^2_\mathbf{b} = \frac{1}{n-1}\mathbf{b}\mathbf{b^T}
\end{equation}

(the n-1 is because we treat the vectors as a sample of data, and we seek an 'unbiased estimator').

The covariance of $\mathbf{a}$ and $\mathbf{b}$ is

\begin{equation}
\sigma^2_\mathbf{ab} = \frac{1}{n-1}\mathbf{a}\mathbf{b^T}\quad\text{.}
\end{equation}

If the data in $\mathbf{a}$ and $\mathbf{b}$ are correlated, the covariance will be high. If $\mathbf{a}$ and $\mathbf{b}$ are not correlated, the covariance will be low.

We can obtain the covariance of our transformed data matrix $\mathbf{Y}$ as follows,

\begin{equation}
\mathbf{C_Y}=\frac{1}{n-1}\mathbf{Y}\mathbf{Y^T}\quad\text{.}
\end{equation}

What do we want the covariance matrix, $\mathbf{C_Y}$, to look like? The entries on the diagonal will be the variances in each of our (transformed) measurement dimensions. If we have chosen these dimensions correctly, i.e. we are using the principal components of the data, then all of the off-diagonal entries will be zero (uncorrelated).

Since $\mathbf{Y} = \mathbf{U}^T\mathbf{X}$,

\begin{align}
\mathbf{C_Y} &= \frac{1}{n-1}\mathbf{Y}\mathbf{Y^T} \\
 &= \frac{1}{n-1} (\mathbf{U}^T\mathbf{X})(\mathbf{U}^T\mathbf{X})^T \\
 &= \frac{1}{n-1} \mathbf{U}^T\:\mathbf{X}\mathbf{X}^T\:\mathbf{U} \\
 &= \frac{1}{n-1} \mathbf{U}^T\:\mathbf{U}\mathbf{\Sigma}\mathbf{V^T}(\mathbf{U}\mathbf{\Sigma}\mathbf{V^T})^T\:\mathbf{U} \\
 &= \frac{1}{n-1} \mathbf{U}^T\:\mathbf{U}\mathbf{\Sigma}\mathbf{V^T}\mathbf{V}\mathbf{\Sigma}\mathbf{U^T}\:\mathbf{U} \\
 & = \frac{1}{n-1} \mathbf{\Sigma}^2
\end{align}

This is exactly what we want. $\mathbf{C_Y}$ is a matrix which only has entries on the diagonal. This means that the basis vectors on to which we have projected $\mathbf{X}$ must be the principal components. We also know that $\mathbf{\Sigma}$ has the singular values in order largest to smallest, so the columns of $\mathbf{U}$ are our principal component basis vectors arranged such that the first column gives the basis with the largest variance in $\mathbf{Y}$ and so on.

## 3. Simple demonstration of PCA

We can use Python to demonstrate PCA.


```python
import numpy as np 
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D
%matplotlib inline
```

First, we construct an artifical set of data like the cartoon above.


```python
npts = 500                                        # number of points
angle = 20.                                       # this is the angle of our first principal component direction
angle_rad = angle/180.*np.pi                 
s = np.random.normal(size=npts)                   # the s and n coordinates are along, and perpendicular, to PC1
n = 0.1* np.random.normal(size=npts)
x = s*np.cos(angle_rad) - n*np.sin(angle_rad)     # obtain the x and y of each point
y = s*np.sin(angle_rad) + n*np.cos(angle_rad)
x = x - np.mean(x)                                # make sure x and y have zero mean
y = y - np.mean(y)
plt.scatter(x,y,alpha=0.4)
plt.axis("equal")
plt.xlabel("x")
plt.ylabel("y");
```

Construct our data matrix $\mathbf{X}$. The matrix will have two rows (corresponding to x and y) and the number of columns will be the number of data points.


```python
X = np.stack((x,y))                                # numpy has a handy function for combinine two row vectors into a matrix
```


```python
plt.scatter(X[0,:],X[1,:],alpha=0.4)
plt.axis("equal")
plt.xlabel("x")
plt.ylabel("y");
```

Obtain the SVD of $\mathbf{X}$


```python
U, S, VT = np.linalg.svd(X)                         
```

We can immediately check to see if the first principle component is at the expected angle.


```python
np.arctan(U[0,1]/U[0,0])/np.pi*180.             # U[0,1] is the y component of the first PC, U[0,0] is the x component
```

We now obtain $\mathbf{Y}$ by projecting our data matrix on to the column vectors of $\mathbf{U}$.


```python
Y = U.T @ X
```

And we can plot the data in this transformed basis


```python
plt.scatter(Y[0,:],Y[1,:],alpha=0.4)
plt.axis("equal")
plt.xlabel("PC 1")
plt.ylabel("PC 2");
```

## 4. 'Eigenwakes'

Demonstrations of PCA in the literature sometimes use an example with many photographs of faces. The principal component vectors (each of which is an image) are called [Eigenfaces](https://en.wikipedia.org/wiki/Eigenface). In the spirit of this example, we look at Eigenwakes!

First load a set of wakes 1188 wakes from CFD simulations of compressor stators.


```python
wakes = np.load("data/cfd_images.npy")
wakes = wakes.T                               # take the transpose so that each column of wakes is an individual wake
imsize = [56,56]                              # each wake is a 56 x 56 "image"
n_wakes = wakes.shape[1]
wakes.shape
```

The first column of wakes is a clean wake:


```python
nwake = 0
wake = np.reshape(wakes[:,nwake],imsize)
plt.imshow(wake, vmin=0, vmax=1.)
plt.colorbar()
plt.title("wake "+str(nwake));
```

Wake number 700, for example, has a large corner separation:


```python
nwake = 700
wake = np.reshape(wakes[:,nwake],imsize)
plt.imshow(wake, vmin=0, vmax=1.)
plt.colorbar()
plt.title("wake "+str(nwake));
```

To perform PCA, we need the mean of each row of our data matrix to be zero.
Obtain the mean wake:


```python
wake_mean = np.mean(wakes, axis = 1)
plt.imshow(np.reshape(wake_mean,imsize), vmin=0, vmax=1.)
plt.colorbar()
plt.title("mean wake");
```

Subtract the mean wake from each individual wake to obtain our data matrix, $\mathbf{X}$. Then perform the SVD.


```python
X = wakes - np.outer(wake_mean, np.ones([1,n_wakes]))
U, S, VT = np.linalg.svd(X)
```

We can plot the columns of $\mathbf{U}$ - our principal component vectors - as images.


```python
plt.figure(figsize=(15,8))
for i in range(4):
    plt.subplot(1,4,i+1)
    plt.imshow(np.reshape(U[:,i],imsize))
    plt.title('PC mode '+str(i+1))
```

Project $\mathbf{X}$ on to the columns of $\mathbf{U}$ to obtain $\mathbf{Y}$:


```python
Y = U.T @ X
```

What do the wakes look like if we reduce the number of dimensions to `nPC = 10` ? Remember that the number of dimensions of our original data was 56 x 56 = 3136.


```python
nPC = 10
```

Reconstruct wake 0 using only the first nPC Principal components:


```python
nwake = 0
wake_recon = wake_mean + U[:,:nPC] @ Y[:nPC, nwake]
plt.imshow(np.reshape(wake_recon,imsize), vmin=0, vmax=1.)
plt.title("Wake "+str(nwake)+", reconstructed with "+str(nPC)+" PC modes");
```

and wake 700:


```python
nwake = 700
wake_recon = wake_mean + U[:,:nPC] @ Y[:nPC, nwake]
plt.imshow(np.reshape(wake_recon,imsize), vmin=0, vmax=1.)
plt.title("Wwake "+str(nwake)+", reconstructed with "+str(nPC)+" PC modes");
```

We can plot the principal components (the values in $\mathbf{Y}$) that are needed to reconstruct the wakes.


```python
xx = np.arange(1, nPC+1)
plt.scatter(xx, Y[:nPC, 0],c='b')
plt.scatter(xx, Y[:nPC, 700],c='r')
plt.xlabel('PC number')
plt.ylabel('Value of PC');
```

In the above figure, the blue points are from wake 0 - a clean wake. The red points are from wake 700 - a corner separation.

Using a different technique (convolutional neural network image processing) all of the wakes in our data set have been labeled, and the results are in the cornsep.txt file. A value of 0 means the wake is clean, a value of 1 means a corner separation has been detected.


```python
cornsep = np.loadtxt("data/cornsep.txt", dtype=int)
```

Using these labels, we can plot the values of the principal components for all the wakes in $\mathbf{Y}$.


```python
for i in range(n_wakes):
    plt.scatter(xx,Y[:nPC,i],c=np.array(["b", "r"])[cornsep[i]],alpha=0.01)
plt.xlabel('PC number')
plt.ylabel('Value of PC');
```

You can see that, for some of the principal components, the blue and red dots do not overlap much - this is the case for PC2 for example. If a principal compoment, or a small set of principal components, can be used to differentiate between clean wakes and corner separated wakes, this could be used in an automated classification process.

To illustrate this, we plot each wake as a point in a 2D plot where the axes are PC1 and PC2:


```python
plt.scatter(Y[0,:],Y[1,:],c=np.array(["b", "r"])[cornsep],alpha=0.2)
plt.xlabel('PC 1')
plt.ylabel('PC 2');
```

or in a 3D plot where we use the first 3 principal components as coodinates: 


```python
fig = plt.figure()
ax = fig.add_subplot(111, projection='3d')
ax.scatter(Y[0,:],Y[1,:],Y[2,:],c=np.array(["b", "r"])[cornsep],alpha=0.2)
ax.set_xlabel('PC 1')
ax.set_ylabel('PC 2')
ax.set_zlabel('PC 3');

```

The PCA has helped us to see the low dimensional - reduced order - structure that is important in our original, high dimensional data.


```python

```

83d307c398ef6d1bdfca9f47c097b7571782cd42
ipynb
Jupyter Notebook
02 - Principal Component Analysis.ipynb
grahampullan/datascienceBB
87c1127d9deddd3e4e11f1a785b4f9e93a74fb19
MIT
02 - Principal Component Analysis.ipynb
grahampullan/datascienceBB
87c1127d9deddd3e4e11f1a785b4f9e93a74fb19
MIT
02 - Principal Component Analysis.ipynb
grahampullan/datascienceBB
87c1127d9deddd3e4e11f1a785b4f9e93a74fb19
MIT
Qwen/Qwen-72B
1. YES
2. YES