PCA基本原理

数据降维

当数据维度或特征变量的数量过大时，极有可能造成维数灾难，导致计算成本上升。而在实际应用中提取出的有用信息并不需要那么高的维度。

因此本文将介绍一种极其常见的降维方法：主成分分析（Principal Component Analysis， PCA）。

主成分分析

在主成分分析中，数据从原来的坐标系转换到新的坐标系中，尽可能减少信息冗余，新坐标系的选择是由数据本身决定的。

通常来说，第一个新坐标轴选择的是原始数据中方差最大的方向，对应的特征向量作为一个主成分。第二个新坐标轴的选择和第一个坐标轴正交且具有最大方差的方向，对应的特征向量作为另一个主成分。该过程一直重复，重复次数为原始数据中特征的数目。

注意，得到的不同轴之间是不相关的（即协方差为0）。此时新空间的维度小于原来的空间，所以把数据投影到新的空间后，尽可能地保留了原始数据集的信息，又大大降低数据的复杂度。

正如下面这张图，把二维数据降维到了一维空间中的一条线上，此时方差最大。

那么，这些主成分应该如何求出呢？一个矩阵的主成分是它协方差矩阵的特征向量，特征值的大小则能表明该成分的重要性，因此关键还是协方差矩阵。

协方差矩阵

协方差可以用来衡量两个变量之间的相关程度。当变量不相关时，协方差为0，正负号表示两变量之间呈正相关或负相关。其计算公式如下：

cov(X,Y)=∑in(Xi−Xˉ)(Yi−Yˉ)n−1cov(X,Y) = \sum^n _i\frac {(X_i - \bar{X})(Y_i - \bar{Y})}{n-1}

其中Xˉ\bar{X}与Yˉ\bar{Y}为矩阵均值，每个变量都减去它的均值。

若将所有变量两两求协方差，可得到一个协方差矩阵。一个包含3个特征的协方差矩阵如下：

C=[cov(x1,x1)cov(x1,x2)cov(x1,x3)cov(x2,x1)cov(x2,x2)cov(x2,x3)cov(x3,x1)cov(x3,x2)cov(x3,x3)]C = \begin{bmatrix} cov(x_1, x_1) & cov(x_1, x_2) & cov(x_1, x_3) \\ cov(x_2, x_1) & cov(x_2, x_2) & cov(x_2, x_3) \\ cov(x_3, x_1) & cov(x_3, x_2) & cov(x_3, x_3) \end{bmatrix}

若矩阵大小为M×NM \times N，其中MM为样本数，NN为特征数量，对于下面这个数据集而言，M=4,N=3M=4, N=3。

长 宽 高
5 3 2
6 2 1
4 3 4
8 4 3

简单计算得长宽高的均值分别为6.5, 2.5, 3.5，接着我们可以计算得到3×33\times3的协方差矩阵：

1
2
3
4
5
6
7
import numpy as np
a = np.array([[5, 3, 2], [6, 2, 1], [4, 3, 4], [8, 4, 3]])
print(np.cov(a.T))   # np.cov()计算协方差矩阵
# 输出
[[ 2.91666667  0.66666667 -0.5       ]
 [ 0.66666667  0.66666667  0.66666667]
 [-0.5         0.66666667  1.66666667]]

特征值与特征向量

前面提到，主成分是协方差矩阵的特征向量，而特征值越大，表明该成分的方差越大。

线性代数中，设AA为矩阵，vv为特征向量，λ\lambda为特征值，此公式体现出了三者之间的关系：

Av→=λv→A \overrightarrow{v} = \lambda \overrightarrow{v}

我们的原始数据是4×34 \times 3的矩阵，假设要降到一维4×14 \times 1， 则需要将数据投影到第一主成分上，[4×3]×[3×1]=[4×1][4 \times 3] \times [3 \times 1] = [4 \times 1]，因此第一主成分特征向量维度为3×13 \times 1。

但注意，使用的是原始数据减去该特征的平均值后再乘以第一主成分：

[5−6.53−2.52−3.56−6.52−2.51−3.54−6.53−2.54−3.58−6.54−2.53−3.5]\begin{bmatrix} 5-6.5 & 3-2.5 & 2-3.5 \\ 6-6.5 & 2-2.5 & 1-3.5 \\ 4-6.5 & 3-2.5 & 4-3.5 \\ 8-6.5 & 4-2.5 & 3-3.5 \end{bmatrix}

接下来通过np.linalg.eig()计算出协方差矩阵的特征值和特征向量：

1
2
3
4
5
6
7
8
9
10
import numpy as np
a = np.array([[5, 3, 2], [6, 2, 1], [4, 3, 4], [8, 4, 3]])
eig_values, eig_vectors = np.linalg.eig(np.cov(a.T))
print(eig_values)
# [0.09417879 3.1723972  1.98342402]
print(eig_vectors, eig_vectors.shape)
# 输出 
[[-0.2799539   0.95365046  0.11034772]
   [ 0.84864579  0.19210102  0.49284636]
   [-0.44880524 -0.23162038  0.86309087]]  (3, 3)

由特征值大小可以发现第二个特征(3.172)为第一主成分，第三个特征(1.983)为第二主成分，第一个特征(0.094)为第三主成分。

通过计算每个组成部分所解释的方差百分比，可知前两个维度已经能代表绝大多数情况了。因此在对特征值排序之后，选择前两个维度，并获得相应的特征向量构成矩阵。

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
# 对特征值与特征向量排序
idx = np.argsort(-eig_values)
eig_values = eig_values[idx]
eig_vectors = eig_vectors[:, idx]

# 计算方差百分比
total_eig_values = np.sum(eig_values, axis=0)
eig_values_ratio = eig_values / total_eig_values
print(eig_values_ratio)  # [0.60426613 0.37779505 0.01793882]

# 取前2大的特征值与特征向量
eig_values = eig_values[:2]
eig_vectors = eig_vectors[:, :2]
print(eig_values, eig_vectors)

# 输出
[[ 0.95365046  0.11034772]
 [ 0.19210102  0.49284636]
 [-0.23162038  0.86309087]] 

[3.1723972  1.98342402]

最后将矩阵与特征向量相乘，实现三维到二维的转换：

1
2
3
trans = np.matmul(a, eig_vectors)
fig = plt.scatter(trans[:, 0], trans[:, 1], marker='o')
plt.show()

可以看到，二维坐标原点发生了移动，新坐标系的选择是由数据本身决定的。

因此，PCA的大致思路如下：

  1. 样本中心化：算出数据在每一个维度上的平均值，让该维数值减去这个平均值，中心化不会改变求得的新空间，但会减少计算量。
  2. 对中心化后的数据，算出这些数据的协方差矩阵。
  3. 对协方差矩阵进行对角化，即算出协方差矩阵的特征值与特征向量。
  4. 取特征值大的一些特征向量构成一个矩阵 P，这个矩阵是一个投影矩阵，把原始空间的数据投影到新的空间。