主成分分析(PCA)

假设$x$为$m$ 维随机变量，其均值为$\mu$，协方差矩阵为$\Sigma$

考虑由$m$维随机变量$x$到$m$维随机变量$y$的线性变换

$$y_i={\alpha }_i^{T}x=\sum _{k=1}^m{\alpha }_{ki}{x}_k,i=1,2,\cdots ,m$$

其中${\alpha }_i^T=({\alpha }_{1i},{\alpha }_{2i},\cdots ,{\alpha }_{mi})$

如果该线性变换满足以下条件，则称之为总体主成分：

${\alpha }_i^T{\alpha }_i=1,i=1,2,\cdots ,m$；

$\mathrm{cov}(y_i,y_j)=0(i\ne j)$;

变量$y_1$是$x$的所有线性变换中方差最大的

$y_2$是与$y_1$不相关的$x$的所有线性变换中方差最大的

一般地，$y_i$是与$y_1,y_2,\cdots ,y_{i-1},(i=1,2,\cdots ,m)$都不相关的$x$的所有线性变换中方差最大的

这时分别称$y_1,y_2,\cdots ,y_m$为$x$的第一主成分、第二主成分、…、第$m$主成分

---

假设$x$是$m$维随机变量

其协方差矩阵是$\Sigma$

$\Sigma$的特征值分别是${\lambda }_1\ge {\lambda }_2\ge \cdots \ge {\lambda }_m\ge 0$

特征值对应的单位特征向量分别是${\alpha }_1,{\alpha }_2,\cdots ,{\alpha }_m$，则$x$的第2主成分可以写作

$$y_i={\alpha }_i^{T}x=\sum _{k=1}^m{\alpha }_{ki}{x}_k,i=1,2,\cdots ,m$$

并且，$x$的第$i$主成分的方差是协方差矩阵$\Sigma$的第$i$个特征值，即$$\mathrm{var}(y_i)={\alpha }_i^T\Sigma {\alpha }_i={\lambda }_i$$

---

主成分有以下性质：

主成分$y$的协方差矩阵是对角矩阵$$\mathrm{cov}(y)=\Lambda =\mathrm{diag}({\lambda }_1,{\lambda }_2,\cdots ,{\lambda }_m)$$

主成分$y$的方差之和等于随机变量$x$的方差之和

$$\sum _{i=1}^m{\lambda }_i=\sum _{i=1}^m{\sigma }_{ii}$$

其中${\sigma }_{ii}$是$x_2$的方差，即协方差矩阵$\Sigma$的对角线元素

主成分$y_k$与变量$x_2$的相关系数$\rho (y_k,x_i)$称为因子负荷量（factor loading）

它表示第$k$个主成分$y_k$与变量$x$的相关关系，即$y_k$对$x$的贡献程度

$$\rho (y_k,x_i)=\frac{\sqrt{{\lambda }_k}{\alpha }_{ik}}{\sqrt{{\sigma }_{ii}}},k,i=1,2,\cdots ,m$$

---

样本主成分分析就是基于样本协方差矩阵的主成分分析

给定样本矩阵

$$X=\left[\begin{array}{llll}{x}_1& x_2& \cdots & x_n\end{array}\right]=\left[\begin{array}{cccc}{x}_{11}& x_{12}& \cdots & x_{1n}\\ & & & \\ x_{21}& x_{22}& \cdots & x_{2n}\\ & & & \\ ⋮& ⋮& & ⋮\\ & & & \\ x_{m1}& x_{m2}& \cdots & x_{mn}\end{array}\right]$$

其中$x_j=(x_{1j},x_{2j},\cdots ,x_{mj}{)}^T$是$x$的第$j$个独立观测样本，$j=1,2，\dots ,n$。

$X$的样本协方差矩阵

$$\begin{array}{c}S=[s_{ij}{]}_{m\times m},s_{ij}=\frac{1}{n-1}{\sum }_{k=1}^n(x_{ik}-{\overline{x}}_i)(x_{jk}-{\overline{x}}_j)\\ \\ i=1,2,\cdots ,m,j=1,2,\cdots ,m\end{array}$$

给定样本数据矩阵$X$，考虑向量$x$到$y$的线性变换$$y=A^{T}x$$

这里

$$A=\left[\begin{array}{llll}{a}_1& a_2& \cdots & a_m\end{array}\right]=\left[\begin{array}{cccc}{a}_{11}& a_{12}& \cdots & a_{1m}\\ & & & \\ a_{21}& a_{22}& \cdots & a_{2m}\\ & & & \\ ⋮& ⋮& & ⋮\\ & & & \\ a_{m1}& a_{m2}& \cdots & a_{mm}\end{array}\right]$$

如果该线性变换满足以下条件，则称之为样本主成分。样本第一主成分$y_1=a_1^{T}x$是在$a_1^T{a}_1=1$条件下

使得$a_1^T{x}_j(j=1,2,\cdots ,n)$的样本方差$a_1^{T}S{a}_1$最大的$x$的线性变换

样本第二主成分$y_2=a_2^{T}x$

是在$a_2^T{a}_2=1$和$a_2^T{x}_j$与$a_1^T{x}_j(j=1,2,\cdots ,n)$的样本协方差

$a_1^{T}S{a}_2=0$条件下

使得$a_2^T{x}_j(j=1,2,\cdots ,n)$的样本方差$a_2^{T}S{a}_2$最大的$x$的线性变换

一般地，样本第$i$主成分$y_i=a_i^{T}x$是在$a_i^T{a}_i=1$和$a_i^T{x}_j$与$a_k^T{x}_j(k<i,j=1,2,\cdots ,n)$的样本协方差

$a_k^{T}S{a}_i=0$条件下

使得$a_i^T{x}_j(j=1,2,\cdots ,n)$的样本方差$a_k^{T}S{a}_i$最大的$x$的线性变换

---

主成分分析方法主要有两种，可以通过相关矩阵的特征值分解或样本矩阵的奇异值分解进行

相关矩阵的特征值分解算法

针对$m\times n$样本矩阵$X$，求样本相关矩阵

$$R=\frac{1}{n-1}X{X}^T$$

再求样本相关矩阵的$k$个特征值和对应的单位特征向量，构造正交矩阵

$$V=(v_1,v_2,\cdots ,v_k)$$

$V$的每一列对应一个主成分，得到$k\times n$样本主成分矩阵

$$Y=V^{T}X$$

矩阵$X$的奇异值分解算法

针对$m\times n$样本矩阵$X$

$$X^{\prime }=\frac{1}{\sqrt{n-1}}{X}^T$$

对矩阵$X^{\prime }$进行截断奇异值分解，保留$k$个奇异值、奇异向量，得到

$$X^{\prime }=US{V}^T$$

$V$的每一列对应一个主成分，得到$k\times n$样本主成分矩阵$Y$

$$Y=V^{T}X$$

---

PCA（principal components analysis）即主成分分析技术旨在利用降维的思想，把多指标转化为少数几个综合指标。

PCA的算法相当简单。 在确保数据被归一化之后，输出仅仅是原始数据的协方差矩阵的奇异值分解。

现在我们有主成分（矩阵U），我们可以用这些来将原始数据投影到一个较低维的空间中

对于这个任务，我们将实现一个计算投影并且仅选择顶部K个分量的函数，有效地减少了维数

我们也可以通过反向转换步骤来恢复原始数据。

第一主成分的投影轴基本上是数据集中的对角线

当我们将数据减少到一个维度时，我们失去了该对角线周围的变化，所以在我们的再现中，一切都沿着该对角线

---

代码实现

import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
import seaborn as sb
from scipy.io import loadmat

X = data['X']
def pca(X):
    # normalize the features
    X = (X - X.mean()) / X.std()
    
    # compute the covariance matrix
    X = np.matrix(X)
    cov = (X.T * X) / X.shape[0]
    
    # perform SVD
    U, S, V = np.linalg.svd(cov)
    
    return U, S, V
def project_data(X, U, k):
    U_reduced = U[:,:k]
    return np.dot(X, U_reduced)
def recover_data(Z, U, k):
    U_reduced = U[:,:k]
    return np.dot(Z, U_reduced.T)
U, S, V = pca(X)
Z = project_data(X, U, 1)
X_recovered = recover_data(Z, U, 1)
fig, ax = plt.subplots(figsize=(12,8))
ax.scatter(list(X_recovered[:, 0]), list(X_recovered[:, 1]))
plt.show()