写在前面: 本来这篇应该是上周四更新, 但是上周四写了一篇深度学习的反向传播法的过程, 就推迟更新了本来想参考 PRML 来写, 但是发现里面涉及到比较多的数学知识, 写出来可能不好理解, 我决定还是用最通俗的方法解释 PCA, 并举一个实例一步步计算, 然后再进行数学推导, 最后再介绍一些变种以及相应的程序(数学推导及变种下次再写好了)
正文:
在数据处理中, 经常会遇到特征维度比样本数量多得多的情况, 如果拿到实际工程中去跑, 效果不一定好一是因为冗余的特征会带来一些噪音, 影响计算的结果; 二是因为无关的特征会加大计算量, 耗费时间和资源所以我们通常会对数据重新变换一下, 再跑模型数据变换的目的不仅仅是降维, 还可以消除特征之间的相关性, 并发现一些潜在的特征变量
一 PCA 的目的
PCA 是一种在尽可能减少信息损失的情况下找到某种方式降低数据的维度的方法通常来说, 我们期望得到的结果, 是把原始数据的特征空间 (n 个 d 维样本) 投影到一个小一点的子空间里去, 并尽可能表达的很好 (就是说损失信息最少) 常见的应用在于模式识别中, 我们可以通过减少特征空间的维度, 抽取子空间的数据来最好的表达我们的数据, 从而减少参数估计的误差注意, 主成分分析通常会得到协方差矩阵和相关矩阵这些矩阵可以通过原始数据计算出来协方差矩阵包含平方和与向量积的和相关矩阵与协方差矩阵类似, 但是第一个变量, 也就是第一列, 是标准化后的数据如果变量之间的方差很大, 或者变量的量纲不统一, 我们必须先标准化再进行主成分分析
二 PCA VS MDA
提到 PCA, 可能有些人会想到 MDA(Multiple Discriminate Analysis, 多元判别分析法), 这两者都是线性变换, 而且很相似只不过在 PCA 中, 我们是找到一个成分 (方向) 来把我们的数据最大化方差, 而在 MDA 中, 我们的目标是最大化不同类别之间的差异(比如说, 在模式识别问题中, 我们的数据包含多个类别, 与两个主成分的 PCA 相比, 这就忽略了类别标签)
换句话说, 通过 PCA, 我们把整个数据集 (不含类别标签) 投射到一个不同的子空间中, 在 MDA 中, 我们试图决定一个合适的子空间来区分不同类别再换种方式说, PCA 是找到数据传播最广的时候的最大方差的轴 axis,MDA 是最大化类别与类别之间的区别
上文我们提到了子空间, 那么怎么样去寻找好的子空间呢?
假设我们的目标是减少 d 维的数据集, 将其投影到 k 维的子空间上 (看 k<d) 所以, 我们如何来确定 k 呢? 如何知道我们选择的特征空间能够很好的表达原始数据呢?
下文中我们会计算数据中的特征向量 (主成分), 然后计算散布矩阵(scatter_matrix) 中(也可以从协方差矩阵中计算)每个特征向量与特征值相关, 即特征向量的长度或大小如果发现每个特征值都很小, 那就可以说明我们的原始数据就已经是一个好的空间了但是, 如果有些特征值比其他值要大得多, 我们只需要关注那些特别大的特征值, 因为这些值包含了数据分布情况的绝大部分信息反之, 那些接近于 0 的特征值包含的信息几乎没有, 在新的特征空间里我们可以忽略不计
三 PCA 的过程
通常来说有以下六步:
1. 去掉数据的类别特征(label), 将去掉后的 d 维数据作为样本
2. 计算 d 维的均值向量(即所有数据的每一维向量的均值)
3. 计算所有数据的散布矩阵(或者协方差矩阵)
4. 计算特征值 (e1,e2,...,ed) 以及相应的特征向量(lambda1,lambda2,...,lambda d)
5. 按照特征值的大小对特征向量降序排序, 选择前 k 个最大的特征向量, 组成 d*k 维的矩阵 W(其中每一列代表一个特征向量)
6. 运用 d*K 的特征向量矩阵 W 将样本数据变换成新的子空间(用数学式子表达就是
, 其中 x 是 d*1 维的向量, 代表一个样本, y 是 K*1 维的在新的子空间里的向量)
四具体步骤
1. 数据准备 ---- 生成三维样本向量
首先随机生成 40*3 维的数据, 符合多元高斯分布假设数据被分为两类, 其中一半类别为 w1, 另一半类别为 w2
- #coding:utf-8
- import numpy as np
- np.random.seed(4294967295)
- mu_vec1 = np.array([0,0,0])
- cov_mat1 = np.array([[1,0,0],[0,1,0],[0,0,1]])
- class1_sample = np.random.multivariate_normal(mu_vec1, cov_mat1, 20).T
- assert class1_sample.shape == (3,20)# 检验数据的维度是否为 3*20, 若不为 3*20, 则抛出异常
- mu_vec2 = np.array([1,1,1])
- cov_mat2 = np.array([[1,0,0],[0,1,0],[0,0,1]])
- class2_sample = np.random.multivariate_normal(mu_vec2, cov_mat2, 20).T
- assert class1_sample.shape == (3,20)# 检验数据的维度是否为 3*20, 若不为 3*20, 则抛出异常
运行这段代码后, 我们就生成了包含两个类别的样本数据, 其中每一列都是一个三维的向量, 所有数据是这样的矩阵:
结果:
2. 作图查看原始数据的分布
- from matplotlib import pyplot as plt
- from mpl_toolkits.mplot3d import Axes3D
- from mpl_toolkits.mplot3d import proj3d
- fig = plt.figure(figsize=(8,8))
- ax = fig.add_subplot(111, projection='3d')
- plt.rcParams['legend.fontsize'] = 10
- ax.plot(class1_sample[0,:], class1_sample[1,:], class1_sample[2,:], 'o', markersize=8, color='blue', alpha=0.5, label='class1')
- ax.plot(class2_sample[0,:], class2_sample[1,:], class2_sample[2,:], '^', markersize=8, alpha=0.5, color='red', label='class2')
- plt.title('Samples for class 1 and class 2')
- ax.legend(loc='upper right')
- plt.show()
结果:
3. 去掉数据的类别特征
- all_samples = np.concatenate((class1_sample, class2_sample), axis=1)
- assert all_samples.shape == (3,40)# 检验数据的维度是否为 3*20, 若不为 3*20, 则抛出异常
4. 计算 d 维向量均值
- mean_x = np.mean(all_samples[0,:])
- mean_y = np.mean(all_samples[1,:])
- mean_z = np.mean(all_samples[2,:])
- mean_vector = np.array([[mean_x],[mean_y],[mean_z]])
- print('Mean Vector:\n', mean_vector)
结果:
- print('Mean Vector:\n', mean_vector)
- Mean Vector:,
- array([[ 0.68047077],
- [ 0.52975093],
- [ 0.43787182]]))
5. 计算散步矩阵或者协方差矩阵
a. 计算散步矩阵
散布矩阵公式:
其中 m 是向量的均值:
(第 4 步已经算出来是 mean_vector)
- scatter_matrix = np.zeros((3,3))
- for i in range(all_samples.shape[1]):
- scatter_matrix += (all_samples[:,i].reshape(3,1) - mean_vector).dot((all_samples[:,i].reshape(3,1) - mean_vector).T)
- print('Scatter Matrix:\n', scatter_matrix)
结果:
- print('Scatter Matrix:\n', scatter_matrix)
- ('Scatter Matrix:,
- array([[ 46.81069724, 13.95578062, 27.08660175],
- [ 13.95578062, 48.28401947, 11.32856266],
- [ 27.08660175, 11.32856266, 50.51724488]]))
b. 计算协方差矩阵
如果不计算散布矩阵的话, 也可以用 python 里内置的 numpy.cov()函数直接计算协方差矩阵因为散步矩阵和协方差矩阵非常类似, 散布矩阵乘以 (1/N-1) 就是协方差, 所以他们的特征空间是完全等价的 (特征向量相同, 特征值用一个常数(1/N-1, 这里是 1/39) 等价缩放了)协方差矩阵如下所示:
- cov_mat = np.cov([all_samples[0,:],all_samples[1,:],all_samples[2,:]])
- print('Covariance Matrix:\n', cov_mat)
结果:
- >>> print('Covariance Matrix:\n', cov_mat)
- Covariance Matrix:,
- array([[ 1.20027429, 0.35784053, 0.69452825],
- [ 0.35784053, 1.23805178, 0.29047597],
- [ 0.69452825, 0.29047597, 1.29531397]]))
6. 计算相应的特征向量和特征值
- # 通过散布矩阵计算特征值和特征向量
- eig_val_sc, eig_vec_sc = np.linalg.eig(scatter_matrix)
- # 通过协方差矩阵计算特征值和特征向量
- eig_val_cov, eig_vec_cov = np.linalg.eig(cov_mat)
- for i in range(len(eig_val_sc)):
- eigvec_sc = eig_vec_sc[:,i].reshape(1,3).T
- eigvec_cov = eig_vec_cov[:,i].reshape(1,3).T
- assert eigvec_sc.all() == eigvec_cov.all()
- print('Eigenvector {}: \n{}'.format(i+1, eigvec_sc))
- print('Eigenvalue {} from scatter matrix: {}'.format(i+1, eig_val_sc[i]))
- print('Eigenvalue {} from covariance matrix: {}'.format(i+1, eig_val_cov[i]))
- print('Scaling factor:', eig_val_sc[i]/eig_val_cov[i])
- print(40 * '-')
结果:
- Eigenvector 1:
- [[-0.84190486]
- [-0.39978877]
- [-0.36244329]]
- Eigenvalue 1 from scatter matrix: 55.398855957302445
- Eigenvalue 1 from covariance matrix: 1.4204834860846791
- Scaling factor: 39.0
- ----------------------------------------
- Eigenvector 2:
- [[-0.44565232]
- [ 0.13637858]
- [ 0.88475697]]
- Eigenvalue 2 from scatter matrix: 32.42754801292286
- Eigenvalue 2 from covariance matrix: 0.8314755900749456
- Scaling factor: 39.0
- ----------------------------------------
- Eigenvector 3:
- [[ 0.30428639]
- [-0.90640489]
- [ 0.29298458]]
- Eigenvalue 3 from scatter matrix: 34.65493432806495
- Eigenvalue 3 from covariance matrix: 0.8885880596939733
- Scaling factor: 39.0
- ----------------------------------------
其实从上面的结果就可以发现, 通过散布矩阵和协方差矩阵计算的特征空间相同, 协方差矩阵的特征值 * 39 = 散布矩阵的特征值
当然, 我们也可以快速验证一下特征值 - 特征向量的计算是否正确, 是不是满足方程 (其中 为协方差矩阵, v 为特征向量, lambda 为特征值)
- for i in range(len(eig_val_sc)):
- eigv = eig_vec_sc[:,i].reshape(1,3).T
- np.testing.assert_array_almost_equal(scatter_matrix.dot(eigv), eig_val_sc[i] * eigv,decimal=6, err_msg='', verbose=True)
得出结果未返回异常, 证明计算正确
注: np.testing.assert_array_almost_equal 计算得出的结果不一样会返回一下结果:
- >>> np.testing.assert_array_almost_equal([1.0,2.33333,np.nan],
- ... [1.0,2.33339,np.nan], decimal=5)
- ...
- <type 'exceptions.AssertionError'>:
- AssertionError:
- Arrays are not almost equal
- (mismatch 50.0%)
- x: array([ 1. , 2.33333, NaN])
- y: array([ 1. , 2.33339, NaN])
可视化特征向量
- from matplotlib import pyplot as plt
- from mpl_toolkits.mplot3d import Axes3D
- from mpl_toolkits.mplot3d import proj3d
- from matplotlib.patches import FancyArrowPatch
- class Arrow3D(FancyArrowPatch):
- def __init__(self, xs, ys, zs, *args, **kwargs):
- FancyArrowPatch.__init__(self, (0,0), (0,0), *args, **kwargs)
- self._verts3d = xs, ys, zs
- def draw(self, renderer):
- xs3d, ys3d, zs3d = self._verts3d
- xs, ys, zs = proj3d.proj_transform(xs3d, ys3d, zs3d, renderer.M)
- self.set_positions((xs[0],ys[0]),(xs[1],ys[1]))
- FancyArrowPatch.draw(self, renderer)
- fig = plt.figure(figsize=(7,7))
- ax = fig.add_subplot(111, projection='3d')
- ax.plot(all_samples[0,:], all_samples[1,:], all_samples[2,:], 'o', markersize=8, color='green', alpha=0.2)
- ax.plot([mean_x], [mean_y], [mean_z], 'o', markersize=10, color='red', alpha=0.5)
- for v in eig_vec_sc.T:
- a = Arrow3D([mean_x, v[0]], [mean_y, v[1]], [mean_z, v[2]], mutation_scale=20, lw=3, arrowstyle="-|>", color="r")
- ax.add_artist(a)
- ax.set_xlabel('x_values')
- ax.set_ylabel('y_values')
- ax.set_zlabel('z_values')
- plt.title('Eigenvectors')
- plt.show()
结果:
7. 根据特征值对特征向量降序排列
我们的目标是减少特征空间的维度, 即通过 PCA 方法将特征空间投影到一个小一点的子空间里, 其中特征向量将会构成新的特征空间的轴然而, 特征向量只会决定轴的方向, 他们的单位长度都为 1, 可以用代码检验一下:
- for ev in eig_vec_sc:
- numpy.testing.assert_array_almost_equal(1.0, np.linalg.norm(ev))
因此, 对于低维的子空间来说, 决定丢掉哪个特征向量, 就必须参考特征向量相应的特征值通俗来说, 如果一个特征向量的特征值特别小, 那它所包含的数据分布的信息也很少, 那么这个特征向量就可以忽略不计了常用的方法是根据特征值对特征向量进行降序排列, 选出前 k 个特征向量
- # 生成 (特征向量, 特征值) 元祖
- eig_pairs = [(np.abs(eig_val_sc[i]), eig_vec_sc[:,i]) for i in range(len(eig_val_sc))]
- #对 (特征向量, 特征值) 元祖按照降序排列
- eig_pairs.sort(key=lambda x: x[0], reverse=True)
- #输出值
- for i in eig_pairs:
- print(i[0])
结果:
- 84.5729942896
- 39.811391232
- 21.2275760682
8. 选出前 k 个特征值最大的特征向量
本文的例子是想把三维的空间降维成二维空间, 现在我们把前两个最大特征值的特征向量组合起来, 生成 d*k 维的特征向量矩阵 W
- matrix_w = np.hstack((eig_pairs[0][1].reshape(3,1), eig_pairs[1][1].reshape(3,1)))
- print('Matrix W:\n', matrix_w)
结果:
- >>> print('Matrix W:\n', matrix_w)
- Matrix W:,
- array([[-0.62497663, 0.2126888 ],
- [-0.44135959, -0.88989795],
- [-0.643899 , 0.40354071]]))
9. 将样本转化为新的特征空间
最后一步, 把 2*3 维的特征向量矩阵 W 带到公式
中, 将样本数据转化为新的特征空间
- matrix_w = np.hstack((eig_pairs[0][1].reshape(3,1), eig_pairs[1][1].reshape(3,1)))
- print('Matrix W:\n', matrix_w)
- transformed = matrix_w.T.dot(all_samples)
- assert transformed.shape == (2,40), "The matrix is not 2x40 dimensional."
- plt.plot(transformed[0,0:20], transformed[1,0:20], 'o', markersize=7, color='blue', alpha=0.5, label='class1')
- plt.plot(transformed[0,20:40], transformed[1,20:40], '^', markersize=7, color='red', alpha=0.5, label='class2')
- plt.xlim([-4,4])
- plt.ylim([-4,4])
- plt.xlabel('x_values')
- plt.ylabel('y_values')
- plt.legend()
- plt.title('Transformed samples with class labels')
- plt.show()
结果:
到这一步, PCA 的过程就结束了其实 python 里有已经写好的模块, 可以直接拿来用, 但是我觉得不管什么模块, 都要懂得它的原理是什么 matplotlib 有 matplotlib.mlab.PCA(),sklearn 也有专门一个模块 Dimensionality reduction 专门讲 PCA, 包括传统的 PCA, 也就是我上文写的, 以及增量 PCA, 核 PCA 等等, 除了 PCA 以外, 还有 ZCA 白化等等, 在图像处理中也经常会用到, 内容太多, 下次再写
最后推荐一个博客, 动态展示了 PCA 的过程: http://setosa.io/ev/principal-component-analysis/ 写的也很清楚, 可以看一下; 再推荐一个维基百科的, 讲的真的是详细啊 https://en.wikipedia.org/wiki/Principal_component_analysis
来源: https://juejin.im/entry/5a7a45af6fb9a063592bb87e