机器学习-PCA降维【手撕】

降维算法

降维算法中的”降维“，指的是降低特征矩阵中特征的数量，其目的是为了让算法运算更快，效果更好，同时可以方便数据可视化。过高的维度特征维度的特征矩阵无法通过可视化，数据的性质也就比较难理解。其中主要用到的降维方法为PCA和SVD

PCA

在降维中，我们会减少特征的数量，这意味着删除数据，数据量变少则表示模型可以获取的信息会变少，模型的表现可能会因此受影响。同时，在高维数据中，必然有一些特征是不带有效的信息的（比如噪音），或者有一些特征带有的信息和其他特征的信息是重复的。我们希望能够找到一种方法来帮助我们衡量特征上所带的信息量，让我们在降维过程中，能偶既减少特征的数量，又保留大部分有效信息——将那些带有重复信息的特征合并，并删除那些带无效信息的特征等等——逐渐创造出能够代表原特征矩阵大部分信息的，特征更少的，新特征矩阵。
为了让降维后的特征能够包含最大差异性的主成分方向，我们通过计算数据矩阵的协方差矩阵，然后得到协方差矩阵的特征值特征向量，选择特征值最大的k个特征所对应的特征向量组成的矩阵。这样就可以讲数据矩阵转换到新的空间中，实现数据特征的降维。

计算得到协方差矩阵的特征值特征向量有两种方法：特征值分解协方差矩阵、奇异值分解协方差矩阵，所以PCA主要有两种方法：
基于特征值分解协方差矩阵实现PCA算法、基于SVD分解协方差矩阵实现PCA算法。

特征值分解协方差矩阵

样本方差：

\bar{x} = \frac{1}{n} \sum^N_{i=1}x_i 

样本方差：

S^2 = \frac{1}{n-1}\sum^N_{i=1}(x_i - \bar{x})^2

样本X和Y的协方差：

Cov(X,Y) = E[(X-E(X))(Y-E(Y))] = \frac{1}{n-1} \sum^n_(i=1)(x_i - \bar{x})(y_i - \bar{y})

import numpy as np

data_npy = np.array(data)

# 使用 numpy 按步骤实现 PCA
data_cov = np.cov(data_npy, rowvar=False)  
# 将列也就是样本的特征看作变量，计算它们的协方差，我们得到的应该是一个 6 * 6 的矩阵
# 对协方差矩阵进行特征值分解
eigenvalue, eigenvector = np.linalg.eig(data_cov)  
# eigenvalue: 特征值，eigenvector: 特征向量print(eigenvalue,'\n', eigenvector)'''

# 将特征值排序
sorted_id = sorted(range(len(eigenvalue)), key=lambda k: eigenvalue[k], reverse=True)  
# 返回降序排列好的特征值对应的索引
# 假设我们降到 2 维，即取最前面的两个特征向量
w = np.array([eigenvector[sorted_id[0]], eigenvector[sorted_id[1]]])

SVD分解协方差矩阵

相较于特征值分解只适用于方阵，奇异值分解适合所有实数矩阵。设$m?n$实矩阵$A$，可将其分解为

A = U \Sigma V^T

其中，$U$是满足$U^{T}U=I$的m阶矩阵，$U$是满足$V^{T}V=I$的n阶矩阵。$\Sigma$是$m?n$的矩阵，其中$(\Sigma {)}_{ii}={\sigma }_i$，其位置元素均为0,${\sigma }_i$均为非负实数且满足${\sigma }_1\ge {\sigma }_2\ge {\sigma }_3\ge ...\ge 0$

其中$U$的列向量陈伟矩阵$A$的左奇异向量，$V$的列向量矩阵$A$的右奇异向量。${\sigma }_i$就是矩阵$A$的奇异值，其个数等于$A$的秩。

接下来我们看如何计算$U,\Sigma ,V$。
首先将要分解的$m\times n$的矩阵$A$乘上它的转置，我们就可以得到一个$m\times m$的方阵$(A{A}^T)$。那么是方阵的话我们就可以对其进行特征值分解，得到 $(A{A}^T)U={\lambda }_{i}U$。这里的$U$就是我们需要的那个$U$。
同上一步，接下来将$A$的转置乘上它本身，得到一个$n\times n$的方阵$(A^{T}A)$。同样进行特征值分解，得到$(A^{T}A)V={\lambda }_{i}V$ 就可计算得到$V$由：

A = U\Sigma V^T  \Longrightarrow AV = U\Sigma V ^TV \Longrightarrow AV = U \Sigma \Longrightarrow  
\begin{pmatrix}  
  a_{11} & \cdots & a_{1n} \\  
  \vdots & \ddots & \vdots \\  
  a_{m1} & \cdots & a_{mn}  
\end{pmatrix}  (v_1,v_2,...v_n) = (u_1, u_2,...u_n) \Sigma

通过最后得出矩阵方程我们就可以解出${\sigma }_i$，至此我们对$A$的奇异值分解就结束了。

A = np.array([[1, 2, 8], [-1, 5, 10], [0, 2, 9], [2, 1, 7], [0, 3, 10]])  # shape: (5, 3)
s = np.linalg.svd(A)  # 返回包含三个元素的tuple, 依次是U、Sigma和V^T
'''
(array([[-0.39725362, -0.28034428,  0.20685827, -0.83712199, -0.14153824],
        [-0.52944196,  0.68764079,  0.41136745,  0.16104935, -0.22733515],
        [-0.44192138, -0.11186842, -0.71914378,  0.15421894, -0.50123816],
        [-0.33830158, -0.65916459,  0.4473559 ,  0.49908567, -0.04289846],
        [-0.50213375,  0.039303  , -0.26587853,  0.02049123,  0.821709  ]]),
 array([20.78369486,  3.16909449,  0.99743079]),
 array([[-0.02619432, -0.29687997, -0.95455547],
        [-0.72144086,  0.66660144, -0.18752494],
        [ 0.69198045,  0.68374323, -0.23164251]]))
'''
# 根据公式验证结果
Sigma = np.array([[s[1][0], 0, 0], [0, s[1][1], 0], [0, 0, s[1][2]], [0, 0, 0], [0, 0, 0]])  # Sigma形状是(m, n), 即(5, 3)
A_ = np.linalg.multi_dot([s[0], Sigma, s[2]])
print(np.allclose(A, A_))
# True

参考

https://zhuanlan.zhihu.com/p/620573209