PCA（主成分分析法）的Python代码实现（numpy，sklearn）

• 语言描述
• 算法描述
• 示例
• • 1 使用numpy降维
  • 2 直接使用sklearn中的PCA进行降维

语言描述

PCA设法将原来众多具有一定相关性的属性（比如p个属性），重新组合成一组相互无关的综合属性来代替原属性。通常数学上的处理就是将原来p个属性做线性组合，作为新的综合属性。

PCA 中的线性变换等价于坐标变换，变换的目的是使 $n$ 个样本点在新坐标轴 $y_1$ 上的离散程度（方差）最大，这样变量 $y_1$ 就代表了原始数据的绝大部分信息，即使忽略 $y_2$ 也无损大局，从而把两个指标压缩成一个指标。从几何上看，找主成分的问题就是找出 $N$ 维空间中椭球体的主轴问题。从数学上也可以证明，它们分别是相关矩阵的 $k$ 个较大的特征值所对应的特征向量。

算法描述

输入：$n$ 个 $m$ 维 $(n\times m)$ 样本集 $X=(x_1,x_2,\cdots ,x_m)$，低维空间维数 $k$ 。

主成分的计算步骤如下：

1. 对所有样本进行中心化：$x_i\leftarrow x_i-\frac{1}{n}{\sum }_{i=1}^n{x}_i$

2. 计算样本的协方差矩阵$X^{T}X$ $(m\times m)$

3. 计算特征值与特征向量

   • 解特征方程$|\lambda E-X^{T}X|=0$，常用雅可比 ( Jacobi ) 法求出特征值，并使其按大小顺序排列，即${\lambda }_1\ge {\lambda }_2\ge \cdots \ge {\lambda }_m\ge 0$ 。
   • 分别求出对应于特征值 ${\lambda }_i$ 的特征向量 $e_i(i=1,2,\cdots ,m)$，要求 $||e_i||=1$，即${\sum }_{j=1}^m{e}_{ij}^2=1$，其中 $e_{ij}$ 表示向量 $e_i$ 的第 $j$ 个分量。
   • 计算主成分贡献率及累计贡献率。
     • 贡献率的公式：$f_i=\frac{{\lambda }_i}{{\sum }_{i=1}^m{\lambda }_i}$
     • 累计贡献率：${\alpha }_k={\sum }_{i=1}^k{f}_i$
     • 一般取累计贡献率达 85% ~ 95% 的特征值 ${\lambda }_1,{\lambda }_2,\cdots ,{\lambda }_l$ 所对应的第一、第二、第$l(l\le m)$个主成分。

4. 计算主成分值

   前 $k$ 个主成分值 $z=(X{e}_1,X{e}_2,\cdots ,X{e}_k)=(z_1,z_2,\cdots ,z_k)$ $(n\times k即n个k维)$

与通过保留原属性集的一个子集来减少属性集的大小不同，PCA通过创建一个能替换的、较小的变量集“组合”属性的基本要素。原数据可以投影到较小的集合中。PCA常常能够揭示先前未被察觉的联系，并允许解释不寻常的结果。

示例

学号	语文	数学	物理	化学	英语	历史
1	84	65	61	72	79	81
2	64	77	77	76	55	70
3	65	67	63	49	57	67
4	74	80	69	75	63	74
5	84	74	70	80	74	82

1 使用numpy降维

>>> import numpy as np

# 输入待降维数据 (5 * 6) 矩阵，6个维度，5个样本值
>>> A = np.array([[84,65,61,72,79,81],[64,77,77,76,55,70],[65,67,63,49,57,67],[74,80,69,75,63,74],[84,74,70,80,74,82]])
>>> print(A)
[[84 65 61 72 79 81]
 [64 77 77 76 55 70]
 [65 67 63 49 57 67]
 [74 80 69 75 63 74]
 [84 74 70 80 74 82]]

# 对每一个属性的样本求均值
>>> MEAN = np.mean(A, axis=0) # 沿轴0调用mean函数
>>> print(MEAN)
[74.2 72.6 68.  70.4 65.6 74.8]

# 去中心化
>>> X = np.subtract(A, MEAN)
>>> print(X)
[[  9.8  -7.6  -7.    1.6  13.4   6.2]
 [-10.2   4.4   9.    5.6 -10.6  -4.8]
 [ -9.2  -5.6  -5.  -21.4  -8.6  -7.8]
 [ -0.2   7.4   1.    4.6  -2.6  -0.8]
 [  9.8   1.4   2.    9.6   8.4   7.2]]
>>> print(X.T) #矩阵的转置
[[  9.8 -10.2  -9.2  -0.2   9.8]
 [ -7.6   4.4  -5.6   7.4   1.4]
 [ -7.    9.   -5.    1.    2. ]
 [  1.6   5.6 -21.4   4.6   9.6]
 [ 13.4 -10.6  -8.6  -2.6   8.4]
 [  6.2  -4.8  -7.8  -0.8   7.2]]
 
# 计算协方差矩阵
>>> COV = np.dot(X.T, X)
>>> print(COV)
[[ 380.8  -55.6  -95.   248.6  401.4  252.2]
 [ -55.6  165.2  131.   179.8 -107.8  -20.4]
 [ -95.   131.   160.   170.  -132.   -34. ]
 [ 248.6  179.8  170.   605.2  214.8  215.4]
 [ 401.4 -107.8 -132.   214.8  443.2  263.6]
 [ 252.2  -20.4  -34.   215.4  263.6  174.8]]
 
# 计算特征值和特征向量
>>> W, V = np.linalg.eig(COV)
>>> print(W) # 特征值
[1.22517276e+03 6.54041238e+02 3.95721181e+01 1.04138814e+01
 1.50877843e-14 5.51899893e-14]
>>> print(V) # 特征向量
[[-0.53264253  0.20279107 -0.34433806  0.39437042 -0.61869481 -0.55543331]
 [ 0.00876193 -0.46059524 -0.81597078  0.02185232  0.25842516  0.34848844]
 [ 0.04593605 -0.47328385  0.37877077  0.70892582 -0.03144886  0.21014772]
 [-0.51955599 -0.64238594  0.24891406 -0.45230979 -0.15412561 -0.22434743]
 [-0.55131936  0.32775478  0.09651389 -0.13044526  0.29446728  0.67491022]
 [-0.37445103  0.05145202  0.0297077   0.34614812  0.66255449  0.14160509]]
 
# 计算主成分贡献率以及累计贡献率
>>> sum_lambda = np.sum(W) # 特征值的和
>>> print(sum_lambda)
1929.1999999999994
>>>f = np.divide(W, sum_lambda) # 每个特征值的贡献率（特征值 / 总和）
>>> print(f)
[6.35067780e-01 3.39021998e-01 2.05121906e-02 5.39803100e-03
 7.82074656e-18 2.86077075e-17]
>>> f[0]+f[1] # 前两大的特征值的累计贡献率
0.974089778403108
>>> f[0]+f[1]+f[2] # 前三大的特征值的累计贡献率
0.9946019690025047
# 0.97 > 0.85，只需要选取前两大特征值即可，以从6维降到2维
# 前两大特征值对应的特征向量为：
>>> e1 = V.T[0]
>>> print(e1)
[-0.53264253  0.00876193  0.04593605 -0.51955599 -0.55131936 -0.37445103]
>>> e2 = V.T[1]
>>> print(e2)
[ 0.20279107 -0.46059524 -0.47328385 -0.64238594  0.32775478  0.05145202]

# 计算主成分值（已去中心化）
>>> z1 = np.dot(X, e1)
>>> print(z1)
[-16.14860528  10.61676743  23.40212697  -0.43966353 -17.43062559]
>>> z2 = np.dot(X, e2)
>>> print(z2)
[ 12.48396235 -15.67317428  13.607117    -7.77054621  -2.64735885]

# 输出降维后的结果（已去中心化）
>>> RES = np.array([z1,z2])
>>> print(RES)
[[-16.14860528  10.61676743  23.40212697  -0.43966353 -17.43062559]
 [ 12.48396235 -15.67317428  13.607117    -7.77054621  -2.64735885]]
>>> print(RES.T)
[[-16.14860528  12.48396235]
 [ 10.61676743 -15.67317428]
 [ 23.40212697  13.607117  ]
 [ -0.43966353  -7.77054621]
 [-17.43062559  -2.64735885]]

2 直接使用sklearn中的PCA进行降维

>>> import numpy as np
>>> from sklearn.decomposition import PCA

# 输入待降维数据 (5 * 6) 矩阵，6个维度，5个样本值
>>> A = np.array([[84,65,61,72,79,81],[64,77,77,76,55,70],[65,67,63,49,57,67],[74,80,69,75,63,74],[84,74,70,80,74,82]])
>>> print(A)
[[84 65 61 72 79 81]
 [64 77 77 76 55 70]
 [65 67 63 49 57 67]
 [74 80 69 75 63 74]
 [84 74 70 80 74 82]]
 
# 直接使用PCA进行降维
>>> pca = PCA(n_components=2) #降到 2 维
>>> pca.fit(A)
PCA(n_components=2)
>>> pca.transform(A) # 降维后的结果
array([[-16.14860528, -12.48396235],
       [ 10.61676743,  15.67317428],
       [ 23.40212697, -13.607117  ],
       [ -0.43966353,   7.77054621],
       [-17.43062559,   2.64735885]])
>>> pca.explained_variance_ratio_ # 降维后的各主成分的方差值占总方差值的比例，即方差贡献率
array([0.63506778, 0.339022  ])
>>> pca.explained_variance_ # 降维后的各主成分的方差值
array([306.29319053, 163.51030959])