### Principal Component Analysis


```python
import scipy.io as sio
import numpy as np
import matplotlib.pyplot as plt
mat_contents = sio.loadmat('mnistAll.mat')
```

#### Load data


```python
train_img = mat_contents['mnist'][0][0][0]
train_label = mat_contents['mnist'][0][0][2]
idx3 = np.where(mat_contents['mnist'][0][0][2]==3)
ndx = idx3[0][:1000]; d = 28*28;
n = len(ndx)
X = np.double(train_img[:,:,ndx].reshape(d,n))
```


```python
plt.subplot(1,2,1)
plt.imshow(X[:,0].reshape(28,28))
plt.subplot(1,2,2)
mu = np.mean(X,axis=1)
plt.imshow(mu.reshape(28,28))
```

#### processing


```python
def pcaPmtk(x):
    n,d=x.shape
    K = np.linalg.matrix_rank(x)x
    mu = np.mean(x,axis=1)
    xc = x-np.transpose(np.tile(mu,(d,1)))
    w,v = np.linalg.eigh(np.cov(xc))
    wsort = np.sort(w)[::-1]
    ind = np.argsort(w)[::-1]
    B = v[:,ind[:K]]
    z = np.matmul(np.transpose(xc),B)
    xrecon = np.matmul(z,np.transpose(B))+np.tile(mu,(d,1))
    return B,z,w
```


```python
v,z,evals = pcaPmtk(X)
```

#### Principal Axis


```python
for i in range(4):
    plt.subplot(2,2,i+1)
    plt.imshow(v[:,i].reshape(28,28))
```

#### Reconstructed Image


```python
# number of principal axes
Ks = [5,10,20,np.linalg.matrix_rank(X)]
for j in range(4):
    xre = np.matmul(z[124,:Ks[j]],np.transpose(v[:,:Ks[j]]))+mu
    plt.subplot(2,2,j+1)
    plt.imshow(xre.reshape(28,28))
```

### Proof


#### Cost function
$$J(w,z_{i1})=\frac{1}{N}\sum^N_{i=1}(x_i-z_{i1}w_1)^\top(x_i-z_{i1}w_1)=\frac{1}{N}\sum^N_{i=1}(x_i^\top x_i-2x_i^\top z_{i1}w_1+z_{i1}^2w^\top_1 w_1)$$

#### Minimizing the Cost function, 
Since $w_1$ is orthnormal, $w_1^\top w_1 = 1$

\begin{align*}
\frac{\partial J}{\partial z_{j1}} & = \frac{1}{N}\sum^N_{i=1}(-2x_i^\top \frac{\partial z_{i1}}{\partial z_{j1}}w_1+2z_i\frac{\partial z_{i1}}{\partial z_{j1}}w_1^\top w_1) = \frac{1}{N}\sum^N_{i=1} \delta_{ij}(-2)(x_i^\top w_1-z_{i1}w_1^\top w_1) 
\\
&= \frac{1}{N}(x_j^\top w_1-z_{j1}w_1^\top w_1) = 0 
\qquad\text{  Therefore,   } \qquad z_{j1}=w_1^\top x_j
\end{align*}

#### Substitute back to the Cost function,
\begin{align*}
J(z_{i1}) &= \frac{1}{N}\sum^N_{i=1}(x^\top_i x_i - 2z_{i1}w_1^\top x_i + z_{i1}^2) = \frac{1}{N}\sum^N_{i=1}(x_i^\top x_i - 2z_{i1}^2+z_{i1}^2) = \frac{1}{N}\sum^N_{i=1}(x_i^\top x_i-z_{i1}^2) = constant - \frac{1}{N}\sum^N_{i=1}z^2_{i1}
\end{align*}

#### Covariance matrix
Consider $z_{i1} = w_1^\top x_i$,

$$
\frac{1}{N}\sum^N_{i=1}z_{i1}^2 = \frac{1}{N}\sum^N_{i=1}w_1^\top x_ix_i^\top w_1=w_1^\top \bigg(\frac{1}{N}\sum^N_{i=1}x_ix_i^\top \bigg)w_1 = w_1^\top \hat{\Sigma}w_1
$$

To minimize the cost function $J(z_{i1})$ is equivalent to maximize $\frac{1}{N}\sum^N_{i=1}z_{i1}^2$, so we can rewrite the cost function as $\tilde{J}(w_1)$ with a constraint of $w^\top_1w_1=1$ using Lagrange multiplier.
\begin{equation}
\tilde{J}(w_1) = w^\top_1\hat{\Sigma}w_1+\lambda_1(w^\top_1w_1-1)
\end{equation}
After minimizing the new cost function, obtain the eigenvector 
$$\frac{\partial \tilde{J}}{\partial w_1} = 2\hat{\Sigma}w_1-2\lambda_1w_1 = 0 \qquad \longrightarrow \qquad  \hat{\Sigma}w_1=\lambda_1w_1$$


```python

```


```python

```


```python

```


```python

```


```python

```

e0ebb4c1e4bdc26c29634e37625d45704821e8ef
ipynb
Jupyter Notebook
Machine Learning A Probabilistic Perspective/12LatentLinearModels/F12.5/.ipynb_checkpoints/12.5pcaImageDemo-checkpoint.ipynb
zcemycl/ProbabilisticPerspectiveMachineLearning
8291bc6cb935c5b5f9a88f7b436e6e42716c21ae
MIT
2019-11-20T10:20:29.000Z
2021-11-09T11:15:23.000Z
Machine Learning A Probabilistic Perspective/12LatentLinearModels/F12.5/.ipynb_checkpoints/12.5pcaImageDemo-checkpoint.ipynb
zcemycl/ProbabilisticPerspectiveMachineLearning
8291bc6cb935c5b5f9a88f7b436e6e42716c21ae
MIT
Machine Learning A Probabilistic Perspective/12LatentLinearModels/F12.5/.ipynb_checkpoints/12.5pcaImageDemo-checkpoint.ipynb
zcemycl/ProbabilisticPerspectiveMachineLearning
8291bc6cb935c5b5f9a88f7b436e6e42716c21ae
MIT
2020-05-27T03:56:38.000Z
2021-05-02T13:15:42.000Z
Qwen/Qwen-72B
1. YES
2. YES