简介
降维是由一些问题带来的：
- 可以缓解由维度诅咒（高维）带来的问题；
- 可以用来压缩数据，将损失数据最小化；
- 可以将高维数据降到低维进行可视化。
主成分分析（Principal components analysis，简称PCA）是最重要的降维方法之一。一般我们提到降维最容易想到的算法就是PCA，下面我们就对PCA的原理做一个总结。
核心思想
PCA降维的核心思想是：一个矩阵的主成分是它的协方差矩阵的特征向量，及其对应的特征值排序。
PCA将一系列可能相关联的高维变量减少为一系列被称为主成分的低维线性不相关变量。这些低维数据会尽可能地保留原始数据的方差。比如我们有如下的数据分布：
图中数据从原点到右上角呈现散点分布，我们可以通过x轴和y轴来描述整个数据分布的构型。但是如果我们只想使用一个描述符来描述整个数据怎么办？可以预想这个描述符能够最大化解释整个数据集。比如下图：
在图中，我们发现沿着红色虚线的方向能够最大程度将数据分开，或者说在该方向上数据集的方差达到最大；而如果考虑黑色虚线的方向，此时就很难将数据分开。实际上，在这里红色虚线是第一个主成分，黑色虚线是第二个主成分，它们俩必须正交。如果要选择一个主成分那必须是红色虚线的方向。所以主成分分析事实上就是重新依次寻找方差最大的特征向量作为一个主成分，每个后面的主成分会保留剩余方差的最大值，唯一的限制是它必须和其他的主成分正交。
需要注意的是，当一个数据集中的方差围绕每个维度不均匀分布时，PCA降维才有效。考虑一个具有球形凸包的三维数据集，由于每个维度的方差都相等，没有一个维度可以在不损失信息的情况下被丢弃，此时，PCA将无法有效地运用到该数据集中。
那么这些主成分该如何求出来呢？记住：一个矩阵的主成分是它的协方差矩阵的特征向量，及其对应的特征值排序。关键还是协方差矩阵。
协方差矩阵
协方差是一种衡量两个变量之间相关性程度的方法。如果协方差为0，则变量不相关，正负号表示变量之间是正相关或负相关。它的计算公式如下：
每个变量都减去它的均值。如果将所有变量两两求协方差，就得到一个协方差矩阵。一个包含三个特征的协方差矩阵如下：
假如我们有一个这样的数据集：
特征v1和v2的均值分别是2和5，接着我们可以计算每一对变量的协方差并得到如下的协方差矩阵：
|
1 2 3 4 5 |
>>>import numpy as np >>>a = np.array([[1, 4], [2, 5], [3, 6]]) >>>np.cov(a.T) # 协方差矩阵 array([[1., 1.], [1., 1.]]) |
np.cov()方法提供了非常方便的计算，最后得到2*2的协方差矩阵。
特征值和特征向量
回顾一下，主成分是数据协方差矩阵的特征向量，由对应的特征值进行排序。一个矩阵的特征向量是一个非零的向量。它的定义如下：
其中，v是一个特征向量，A是一个方阵，λ是一个特征值。也就是说，一个矩阵乘以它的一个特征向量（线性变换），等于对这个特征向量进行缩放。通过np.linalg.eig()我们很快就能计算出特征值和特征向量，比如计算前面协方差矩阵的特征值和特征向量：
|
1 2 3 4 5 6 7 8 |
>>>import numpy as np >>>a = np.array([[1, 4], [2, 5], [3, 6]]) >>>eigval, eigvect = np.linalg.eig(np.cov(a.T)) >>>eigval # 特征值 array([2., 0.]) >>>eigvect # 特征向量 array([[ 0.70710678, -0.70710678], [ 0.70710678, 0.70710678]]) |
原始的数据是3*2的矩阵，假设我们要降到一维，也就是3*1，这个过程就是[3*2]*[2*1]=[3*1]，其中[2*1]就是第一主成分的特征向量。下面展示的就是将数据投影到第一主分上的结果：
需要注意的是，这里使用原始数据减去该特征的平均值后再乘以第一主成分。
事实上，主成分除了上面这种技巧来求解特征值和特征向量外，还有使用数据矩阵的奇异值分解(SVD)的方法来寻找协方差的特征向量和特征值的平方根。这种方法即被scikit-learn库的PCA所使用。SVD公式如下：
其中，U是协方差矩阵的左-奇异向量，表示特征向量；V是右-奇异向量；Σ的对角元素是协方差矩阵特征值的平方根。
PCA的Numpy实现
下面我们用一个实例来学习下Numpy中的PCA类使用。
首先随机构建一组包含三个特征的数据集，该数据集可以分成四类:
|
1 2 3 4 5 6 7 8 9 10 11 12 13 14 |
import numpy as np import matplotlib.pyplot as plt from mpl_toolkits.mplot3d import Axes3D from sklearn.datasets.samples_generator import make_blobs # 样本10000个，特征维度3 X, y = make_blobs(n_samples=10000, n_features=3, centers=[[3, 3, 3], [0, 0, 0], [1, 1, 1], [2, 2, 2]], cluster_std=[0.2, 0.1, 0.2, 0.2], random_state=9) # 画图 fig = plt.figure() ax = Axes3D(fig, rect=[0, 0, 1, 1], elev=30, azim=20) plt.scatter(X[:, 0], X[:, 1], X[:, 2], marker='o') plt.show() |
我们生成随机数据并可视化，三维数据的分布图如下：
接下来我们使用前面的方法计算特征值和特征向量：
|
1 2 3 4 |
# 计算特征值和特征向量 eigvalue, eigvector = np.linalg.eig(np.cov(X.T)) print(eigvalue) print(eigvector) |
得到如下结果：
|
1 2 3 4 5 6 |
# 特征值 [3.78521638 0.03202212 0.03272613] # 特征向量 [[-0.57601622 -0.74821695 -0.32920617] [-0.57846191 0.65765428 -0.48256882] [-0.57757002 0.08753423 0.81163454]] |
从特征值大小可以发现第一个特征（3.785）为第一主成分，第三个特征（0.0327）为第二主成分，第二个特征（0.0320）为第三主成分。假设我们最后将数据降维到二维，就是取第一个特征和第三个特征进行计算：
|
1 2 3 4 5 6 |
# 获取第一个和第三个特征 a = np.hstack((eigvector[:, 0].reshape(3, -1), eigvector[:, -1].reshape(3, -1))) print(a) X = X - X.mean(axis=0) X_new1 = X.dot(a) print(X_new1) |
结果为：
|
1 2 3 4 5 6 7 8 9 10 11 12 |
# 第一个和第三个特征向量构成的矩阵 [[-0.57601622 -0.32920617] [-0.57846191 -0.48256882] [-0.57757002 0.81163454]] # 降维后的结果 [[-1.29049617 0.01162118] [ 2.5902227 -0.04141849] [-2.81225258 -0.05286925] ... [-2.52492314 -0.0935418 ] [-2.98206456 0.03861322] [-2.28089246 -0.13618233]] |
这样我们就实现了三维到二维的转换，可以把结果可视化：
|
1 2 |
fig = plt.scatter(X_new1[:, 0], X_new1[:, 1], marker='o') plt.show() |
得到下图：
从图中可以清晰地看出数据点分成了四类。
PCA的scikit-learn实现
scikit-learn集成了PCA方法，调用起来也更加方便。需要注意的是，在scikit-learn中使用了奇异值分解来计算特征向量和特征值。PCA的调用方法非常简单，这里列出全部代码。
|
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 |
import numpy as np import matplotlib.pyplot as plt from sklearn.decomposition import PCA from sklearn.datasets.samples_generator import make_blobs X, y = make_blobs(n_samples=10000, n_features=3, centers=[[3, 3, 3], [0, 0, 0], [1, 1, 1], [2, 2, 2]], cluster_std=[0.2, 0.1, 0.2, 0.2], random_state=9) # 降维到二维 pca = PCA(n_components=2) pca.fit(X) # 输出特征值 print(pca.explained_variance_) 输出特征向量 print(pca.components_) # 降维后的数据 X_new = pca.transform(X) print(X_new) fig = plt.figure() plt.scatter(X_new[:, 0], X_new[:, 1], marker='o') plt.show() |
结果如下：
|
1 2 3 4 5 6 7 8 9 10 11 12 13 |
# 特征值 [3.78521638 0.03272613] # 特征向量 [[ 0.57601622 0.57846191 0.57757002] [-0.32920617 -0.48256882 0.81163454]] # 降维后的结果 [[ 1.29049617 0.01162118] [-2.5902227 -0.04141849] [ 2.81225258 -0.05286925] ... [ 2.52492314 -0.0935418 ] [ 2.98206456 0.03861322] [ 2.28089246 -0.13618233]] |
作图结果如下：
需要注意的是PCA()方法中的参数设置：
|
1 |
pca = PCA(n_components=0.95) |
当n_components=0.95表示指定了主成分累加起来至少占95%的那些成分。
|
1 |
pca = PCA(n_components='mle') |
当n_components='mle'表示让MLE算法自己选择降维维度的效果。
总结
在本文中，我们介绍了一种使用非常多的降维方法—PCA方法。该方法的核心思想是一个矩阵的主成分是它的协方差矩阵的特征向量，及其对应的特征值排序。除此之外，我们介绍了PCA方法的Numpy 实现过程和更简单的scikit-learn库的调用方式。