机器学习之 PCA(主成分分析)

时间 2019-11-08

标签机器学习 pca 成分分析繁體版

原文原文链接

以前相关文章目录：html

下面介绍一种降维算法，即PCA(主成分分析)。python

在机器学习中，有一种问题被称为维数灾难，在实际机器学习项目中，咱们要处理的的样本数据的维数多是成千上万，甚至几十万或者更多的，这种状况下，直接对原始样本数据进行训练建模会耗费大量时间，对应的资源消耗是不可接受，这个时候，咱们就须要对数据进行降维，降维固然意味着信息的丢失，不过鉴于实际数据自己经常存在的相关性，咱们能够想办法在降维的同时将信息的损失尽可能下降.git

PCA是一种具备严格数学基础而且已被普遍采用的降维方法算法

在具体给出PCA算法以前，先回顾下线性代数中的相关知识，首先，须要强调一点，下面所说的向量，在没有特殊说明下，都指的是列向量数组

1.向量的内积和投影

两个维数相同的向量的内积被定义为：dom

(a_1,a_2,\cdots,a_n)\cdot (b_1,b_2,\cdots,b_n)^\mathsf{T}=a_1b_1+a_2b_2+\cdots+a_nb_n

接下来，看下内积的几何意义，假设A，B是两个n维向量，为了简单起见咱们假设A和B均为二维向量，则。则在二维平面上A和B能够用两条发自原点的有向线段表示,以下图：
如今咱们从A点向B所在直线引一条垂线。咱们知道垂线与B的交点叫作A在B上的投影，再设A与B的夹角是a，则投影的矢量长度为，其中 $|A|=\sqrt{x_1^2+y_1^2}$ 是向量A的模，也就是A线段的标量长度机器学习

咱们知道，内积还可表示为 $A\cdot B=|A||B|cos(a)$ ，也就是说，A与B的内积等于A到B的投影长度乘以B的模
因此，设向量B的模为1，则A与B的内积值等于A向B所在直线投影的矢量长度！

基与基变换

咱们知道，一个二维向量能够对应二维笛卡尔直角坐标系中从原点出发的一个有向线段。例以下面这个向量：函数

图中的那个向量，咱们能够表示为（3，2），其中，3，表示的是向量在X轴的投影值为3，2，表示的是在Y轴的投影值为2。post

也就是说咱们其实隐式引入了一个定义：以x轴和y轴上正方向长度为1的向量为标准。那么一个向量(3,2)实际是说在x轴投影为3而y轴的投影为2。注意投影是一个矢量，因此能够为负。学习

因此，向量(x,y)实际上表示线性组合： $$x(1,0)^\mathsf{T}+y(0,1)^\mathsf{T}$$ 而此处的(1,0),(0,1)叫作二维空间中的一组基。也是一组正交基。

咱们之因此默认选择(1,0)和(0,1)为基，固然是比较方便，由于它们分别是x和y轴正方向上的单位向量，所以就使得二维平面上点坐标和向量一一对应，很是方便。但实际上任何两个线性无关的二维向量均可以成为一组基，所谓线性无关在二维平面内能够直观认为是两个不在一条直线上的向量。

例如，(1,1)和(-1,1)也能够成为一组基。通常来讲，咱们但愿基的模是1，由于从内积的意义能够看到，若是基的模是1，那么就能够方便的用向量点乘基而直接得到其在新基上的坐标了！实际上，对应任何一个向量咱们总能够找到其同方向上模为1的向量，只要让两个份量分别除以模就行了。例如，上面的基能够变为$$(\frac{1}{\sqrt{2}},\frac{1}{\sqrt{2}})和(-\frac{1}{\sqrt{2}},\frac{1}{\sqrt{2}})。$$

如今，咱们想得到(3,2)在新基上的坐标，即在两个方向上的投影矢量值，那么根据内积的几何意义，咱们只要分别计算(3,2)和两个基的内积，不可贵到新的坐标为 $(\frac{5}{\sqrt{2}},-\frac{1}{\sqrt{2}})$ 。下图给出了新的基以及(3,2)在新基上坐标值的示意图：

另外这里要注意的是，咱们列举的例子中基是正交的（即内积为0，或直观说相互垂直），但能够成为一组基的惟一要求就是线性无关，非正交的基也是能够的。不过由于正交基有较好的性质，因此通常使用的基都是正交的。

基变换的矩阵表示

经过前面的描述，咱们知道，将(3,2)变换为新基上的坐标，就是用(3,2)去和新的一组基中的每个去作内积运算。咱们能够用矩阵相乘的形式简洁的表示这个变换：

\begin{pmatrix} 1/\sqrt{2} & 1/\sqrt{2} \\ -1/\sqrt{2} & 1/\sqrt{2} \end{pmatrix} \begin{pmatrix} 3 \\ 2 \end{pmatrix} = \begin{pmatrix} 5/\sqrt{2} \\ -1/\sqrt{2} \end{pmatrix}

推广到多个向量，假设有m个向量，只要将二维向量按列排成一个两行m列矩阵，而后用“基矩阵”乘以这个矩阵，就获得了全部这些向量在新基下的值。例如(1,1)，(2,2)，(3,3)，想变换到刚才那组基上，则能够这样表示：

\begin{pmatrix} 1/\sqrt{2} & 1/\sqrt{2} \\ -1/\sqrt{2} & 1/\sqrt{2} \end{pmatrix} \begin{pmatrix} 1 & 2 & 3 \\ 1 & 2 & 3 \end{pmatrix} = \begin{pmatrix} 2/\sqrt{2} & 4/\sqrt{2} & 6/\sqrt{2} \\ 0 & 0 & 0 \end{pmatrix}

由此咱们能够能够看到，基变换能够表示为矩阵相乘

通常的，若是咱们有M个N维向量，想将其变换为由R个N维向量表示的新空间中，那么首先将R个基按行组成矩阵A，而后将向量按列组成矩阵B，那么两矩阵的乘积AB就是变换结果，其中AB的第m列为A中第m列变换后的结果。数学表示为：

\begin{pmatrix} p_1 & p_2 & \cdots & p_M \end{pmatrix}^\mathsf{T} \begin{pmatrix} a_1 & a_2 & \cdots & a_M \end{pmatrix} = \begin{pmatrix} p_1a_1 & p_1a_2 & \cdots & p_1a_M \\ p_2a_1 & p_2a_2 & \cdots & p_2a_M \\ \vdots & \vdots & \ddots & \vdots \\ p_Ra_1 & p_Ra_2 & \cdots & p_Ra_M \end{pmatrix}

若是R小于N时，这就是将一N维数据变换到更低维度的空间去，所以这种矩阵相乘的表示也能够表示降维变换。

两个矩阵相乘的意义是将右边矩阵中的每一列列向量变换到左边矩阵中每一行行向量为基所表示的空间中去

优化目标

经过上面咱们知道，经过选择不一样的一组基，能够将同一组数据给出不一样的表示，并且若是一组基中基的数量少于向量自己的维数，则能够达到降维的效果。

那么，接下来咱们的问题就是如何找到一组最优的基，也就是说，若是咱们有一组N维向量，如今要将它降到K维（K小于N），那么，咱们该如何选择K个基才能最大程度的保留原始N维向量的信息。

下面，以一个具体例子来展开，假设咱们的数据有5条记录，以下：

\begin{pmatrix} 1 & 1 & 2 & 4 & 2 \\ 1 & 3 & 3 & 4 & 4 \end{pmatrix}

每一列为一条数据，一行为一个字段，也就是一个特征。

首先，为了后续处理方便，咱们将数据进行均值归零化，也就是将每一个字段内（也就是每一行）全部值都减去字段的均值，处理之后，新的数据每一个字段的均值都为0。变换后的数据以下：

\begin{pmatrix} -1 & -1 & 0 & 2 & 0 \\ -2 & 0 & 0 & 1 & 1 \end{pmatrix}

在坐标系中的分布以下：

如今，咱们要将这些数据降到一维，可是有要尽量的保留原始的信息。如何去处理呢？

经过上面的讨论咱们知道，这个问题其实是要在二维平面中选择一个方向，将全部的数据否投影到这个方向所在的直线上，那么，如何选择这个方向（或者说基）才能更多的保留的原始信息呢？一种直观的见解是：咱们但愿使投影后的数值尽量分散

以上图为例，能够看出若是向x轴投影，那么最左边的两个点会重叠在一块儿，中间的两个点也会重叠在一块儿，因而自己四个各不相同的二维点投影后只剩下两个不一样的值了，这是一种严重的信息丢失，同理，若是向y轴投影最上面的两个点和分布在x轴上的两个点也会重叠。因此看来x和y轴都不是最好的投影选择。咱们直观目测，若是向经过第一象限和第三象限的斜线投影，则五个点在投影后仍是能够区分的。

下面，咱们经过具体的数学方法来表述下这个问题

方差

上面咱们讨论出，但愿投影后的数值尽量分散，这种分散程度，在数学上能够用方差来表述。字段a的方差表示为：

Var = \frac{1}{m} \sum_{i=1}^m{(a_i-\mu)^2}

因为上面咱们已经对每一个字段进行均值归零化，因此：

因此，上述咱们的问题被表述为：寻找一个一维基，使得全部数据变换为这个基上的坐标表示后，方差值最大

协方差

对于上面二维降成一维的问题来讲，找到那个使得方差最大的方向就能够了。不过对于更高维，还有一个问题须要解决。考虑三维降到二维问题。与以前相同，首先咱们但愿找到一个方向使得投影后方差最大，这样就完成了第一个方向的选择，继而咱们选择第二个投影方向。

若是咱们仍是单纯只选择方差最大的方向，很明显，这个方向与第一个方向应该是“几乎重合在一块儿”，显然这样的维度是没有用的，所以，应该有其余约束条件。从直观上说，让两个字段尽量表示更多的原始信息，咱们是不但愿它们之间存在（线性）相关性的，由于相关性意味着两个字段不是彻底独立，必然存在重复表示的信息。

数学上能够用两个字段的协方差表示其相关性：

因为已经让每一个字段均值为0，则：

能够看到，在字段均值为0的状况下，两个字段的协方差简洁的表示为其内积除以元素数m。

当协方差为0时，表示两个字段彻底独立。为了让协方差为0，咱们选择第二个基时只能在与第一个基正交的方向上选择(正交的化，内积为0，正好和协方差为0相对应)。所以最终选择的两个方向必定是正交的。

至此，咱们获得了降维问题的优化目标：将一组N维向量降为K维（K大于0，小于N），其目标是选择K个单位（模为1）正交基，使得原始数据变换到这组基上后，各字段两两间协方差为0，而字段的方差则尽量大（在正交的约束下，取最大的K个方差）。

协方差矩阵

咱们看到，最终要达到的目的与字段内方差及字段间协方差有密切关系。所以咱们但愿能将二者统一表示，仔细观察发现，二者都可以表示为内积的形式，而内积又与矩阵相乘密切相关。因而咱们来了灵感：

假设咱们只有a和b两个字段，那么咱们将它们按行组成矩阵X：

X=\begin{pmatrix} a_1 & a_2 & \cdots & a_m \\ b_1 & b_2 & \cdots & b_m \end{pmatrix}

而后咱们用X乘以X的转置，并乘上系数1/m：

这个矩阵对角线上的两个元素分别是两个字段的方差，而其它元素是a和b的协方差。二者被统一到了一个矩阵的。这就是两个字段间的协方差矩阵

根据矩阵相乘的运算法则，这个结论很容易被推广到通常状况：设咱们有m个n维数据记录，将其按列排成n乘m的矩阵X，设 $C=\frac{1}{m}XX^\mathsf{T}$ ，则C是一个对称矩阵，其对角线分别个各个字段的方差，而第i行j列和j行i列元素相同，表示i和j两个字段的协方差

协方差矩阵对角化

根据上述推导，咱们发现要达到优化目前，等价于将协方差矩阵对角化：即除对角线外的其它元素化为0，而且在对角线上将元素按大小从上到下排列，这样咱们就达到了优化目的，这样说可能还不是很明晰，咱们进一步看下原矩阵与基变换后矩阵协方差矩阵的关系：

设原始数据矩阵X对应的协方差矩阵为C，P是一组基（每一个基都是一个列向量）组成的矩阵，设,则Y为X对P所包含的那组基作基变换后的数据。设Y的协方差矩阵为D，则有：

\begin{array}{l l l} D & = & \frac{1}{m}YY^\mathsf{T} \\ & = & \frac{1}{m}(P^TX)(P^TX)^\mathsf{T} \\ & = & \frac{1}{m}P^TXX^\mathsf{T}P \\ & = & P^T(\frac{1}{m}XX^\mathsf{T})P \\ & = & P^TCP \end{array}

如今事情很明白了！咱们要找的P不是别的，而是能让原始协方差矩阵对角化的P。换句话说，优化目标变成了寻找一个矩阵P，知足 $P^\mathsf{T}CP$ 是一个对角矩阵，而且对角元素按从大到小依次排列，那么P的前K列就是要寻找的基，用P的前K列组成的矩阵的转置乘以X就使得X从N维降到了K维并知足上述优化条件。

上文知道，协方差矩阵C是一个是实对称矩阵，在线性代数上，实对称矩阵有一系列很是好的性质：
1）实对称矩阵的特征值都为实数
2）实对称矩阵的全部特征向量正交。
3）设特征值 $\lambda$ 重数为r，则必然存在r个线性无关的特征向量对应于 $\lambda$ ，所以能够将这r个特征向量单位正交化。

由上面两条可知，一个n行n列的实对称矩阵必定能够找到n个单位正交特征向量，设这n个特征向量为 $e_1,e_2,\cdots,e_n$ ，咱们将其按列组成矩阵：

E=\begin{pmatrix} e_1 & e_2 & \cdots & e_n \end{pmatrix}

则对协方差矩阵C有以下结论：

E^\mathsf{T}CE=\Lambda=\begin{pmatrix} \lambda_1 & & & \\ & \lambda_2 & & \\ & & \ddots & \\ & & & \lambda_n \end{pmatrix}

其中 $\Lambda$ 为对角矩阵，其对角元素为各特征向量对应的特征值（可能有重复）。

到此，咱们就已经找到了咱们须要的矩阵P = E

P是协方差矩阵的特征向量单位化后按列排列出的矩阵，其中每一列都是C的一个特征向量。若是设P按照 $\Lambda$ 中特征值的从大到小，将特征向量从左到右排列，则用P的前K行组成的矩阵的装置乘以原始数据矩阵X，就获得了咱们须要的降维后的数据矩阵Y。

上述，就是整个PCA的数学原理讨论

PCA算法

总结一下PCA的算法步骤：

设有m条n维数据。

1）将原始数据按列组成n行m列矩阵X

2）将X的每一行（表明一个属性字段）进行零均值化，即减去这一行的均值

3）求出协方差矩阵C=\frac{1}{m}XX^\mathsf{T}

4）求出协方差矩阵的特征值及对应的特征向量

5）将特征向量按对应特征值大小从左到右按列排列成矩阵，取前k列组成矩阵P

6）Y=P^TX即为降维到k维后的数据

实例

这里以上文提到的

为例，咱们用PCA方法将这组二维数据其降到一维。

由于这个矩阵的每行已是零均值，这里咱们直接求协方差矩阵：

C=\frac{1}{5}\begin{pmatrix} -1 & -1 & 0 & 2 & 0 \\ -2 & 0 & 0 & 1 & 1 \end{pmatrix}\begin{pmatrix} -1 & -2 \\ -1 & 0 \\ 0 & 0 \\ 2 & 1 \\ 0 & 1 \end{pmatrix}=\begin{pmatrix} \frac{6}{5} & \frac{4}{5} \\ \frac{4}{5} & \frac{6}{5} \end{pmatrix}

而后求其特征值和特征向量，具体求解方法再也不详述，能够参考相关资料。求解后特征值为：

其对应的特征向量分别是：

c_1\begin{pmatrix} 1 \\ 1 \end{pmatrix},c_2\begin{pmatrix} -1 \\ 1 \end{pmatrix}

其中对应的特征向量分别是一个通解，c_1和c_2可取任意实数。那么标准化后的特征向量为：

\begin{pmatrix} 1/\sqrt{2} \\ 1/\sqrt{2} \end{pmatrix},\begin{pmatrix} -1/\sqrt{2} \\ 1/\sqrt{2} \end{pmatrix}

所以咱们的矩阵P是：

P=\begin{pmatrix} 1/\sqrt{2} & -1/\sqrt{2} \\ 1/\sqrt{2} & 1/\sqrt{2} \end{pmatrix}

能够验证协方差矩阵C的对角化：

P^\mathsf{T}CP=\begin{pmatrix} 1/\sqrt{2} & 1/\sqrt{2} \\ -1/\sqrt{2} & 1/\sqrt{2} \end{pmatrix}\begin{pmatrix} 6/5 & 4/5 \\ 4/5 & 6/5 \end{pmatrix}\begin{pmatrix} 1/\sqrt{2} & -1/\sqrt{2} \\ 1/\sqrt{2} & 1/\sqrt{2} \end{pmatrix}=\begin{pmatrix} 2 & 0 \\ 0 & 2/5 \end{pmatrix}

最后咱们用P的第一行乘以数据矩阵，就获得了降维后的表示：

Y=\begin{pmatrix} 1/\sqrt{2} \\ 1/\sqrt{2} \end{pmatrix}^T\begin{pmatrix} -1 & -1 & 0 & 2 & 0 \\ -2 & 0 & 0 & 1 & 1 \end{pmatrix}=\begin{pmatrix} -3/\sqrt{2} & -1/\sqrt{2} & 0 & 3/\sqrt{2} & -1/\sqrt{2} \end{pmatrix}

代码实现

import numpy as np
import pandas as pd
复制代码

## X: 须要降维的原始数据 topNfeat：须要降到的维数
def PCA(X, topNfeat=9999999):
    # 1.原始数据默认都是每一行为一个样本数据，每一列为一个特征，
    # 因此进行转置，让每一列表明一个样本数据
    X = X.T
    # 2.将数据的每一行（表明一个属性字段）进行零均值化
    meanValues = np.mean(X, axis=1) #计算每一行的均值
    meanValues = np.mat(meanValues).T #将一个向量转换成n*1的矩阵
    meanRemoved = X - meanValues   #均值归零
    
    #3.求出协方差矩阵
    covMat = np.cov(meanRemoved) #cov计算协方差时除的是 （样本个数-1）,也就是自由度
# covMat = meanRemoved @ meanRemoved.T / (meanRemoved.shape[1])
    print("协方差矩阵：\n",covMat)
    
    #4求出协方差矩阵的特征值及对应的特征向量
    eigVals, eigVects = np.linalg.eig(covMat)  #eigVals:特征值 eigVects：特征向量
    print("特征值\n",eigVals)
    print("特征向量\n",eigVects)
    
    #5将特征向量按对应特征值大小从左到右按列排列成矩阵，取前k列组成矩阵
    # argsort函数返回的是数组值从小到大的索引值,参数中加个-号，变为从大到小
    eigValInd = np.argsort(-eigVals)
    eigValInd = eigValInd[0:topNfeat]  #取出前topNfeat个最大的特征值所对应的索引
    redEigVects = eigVects[:,eigValInd]  #redEigVects 即为须要的变换矩阵，即P
    print("变换矩阵：\n",redEigVects)
    
    #6 Y=P^TX即为降维到k维后的数据
    X_PCA = redEigVects.T @ X
    return X_PCA
    
    
复制代码

PCA(X, topNfeat = 1)
复制代码

协方差矩阵：
 [[1.5 1. ]
 [1.  1.5]]
特征值
 [2.5 0.5]
特征向量
 [[ 0.70710678 -0.70710678]
 [ 0.70710678  0.70710678]]
变换矩阵：
 [[0.70710678]
 [0.70710678]]





array([[-2.12132034, -0.70710678,  0.        ,  2.12132034,  0.70710678]])
复制代码

sklearn

sklearn中为咱们已经封装好了对应的PCA接口，下面咱们使用PCA对sklearn中自带的一个手写数字数据集进行降维。

from sklearn import datasets
复制代码

digits = datasets.load_digits()
X = digits.data
y = digits.target
X.shape,y.shape
复制代码

((1797, 64), (1797,))
复制代码

首先划分数据集为训练集和测试集

from sklearn.model_selection import train_test_split
X_train, X_test, y_train, y_test = train_test_split(X, y, random_state=666)
X_train.shape
复制代码

(1347, 64)
复制代码

首先，在不进行pca降维的状况下，使用knn算法进行模型的创建和预测

%%time

from sklearn.neighbors import KNeighborsClassifier

knn_clf = KNeighborsClassifier()
knn_clf.fit(X_train, y_train)
knn_clf.score(X_test,y_test)  #准确率为
复制代码

Wall time: 69.8 ms
复制代码

knn_clf.score(X_test,y_test)  #准确率为
复制代码

0.9866666666666667
复制代码

接下来看下 PCA降维

#decomposition(分解)
from sklearn.decomposition import PCA
复制代码

# n_components 要降到的维数，第一次咱们直接降到2维
pca = PCA(n_components =2)
pca.fit(X_train)
X_train_reducation = pca.transform(X_train)
X_test_reducation = pca.transform(X_test) #注意测试集也须要作降维处理
X_train_reducation.shape
复制代码

(1347, 2)
复制代码

pca.explained_variance_ratio_ 
#表明降维后的各主成分的方差值占总方差值的比例，这个比例越大，则越是重要的主成分。
复制代码

array([0.14566817, 0.13735469])
复制代码

%%time 

knn_clf = KNeighborsClassifier()
knn_clf.fit(X_train_reducation, y_train) # #降维后，所用时间明显减小
复制代码

Wall time: 2 ms
复制代码

knn_clf.score(X_test_reducation,y_test) 
#因为咱们将64维的数据直接降到了两维，因此信息损失严重，准确率只有0.6
复制代码

0.6066666666666667
复制代码

#若是咱们直接传入一个小数值，则表示要保留的信息占原始信息的比例
#好比下面，保留95%的主要信息
pca = PCA(0.95)
pca.fit(X_train)
X_train_reducation = pca.transform(X_train)
X_test_reducation = pca.transform(X_test)
复制代码

pca.n_components_  #n_components_能够查看pca后降到的具体维数
复制代码

28
复制代码

knn_clf = KNeighborsClassifier()
knn_clf.fit(X_train_reducation, y_train)
knn_clf.score(X_test_reducation,y_test)   
复制代码

0.98
复制代码

能够看到，准确度基本没变化，可是维度却从以前的64维降到28维，从而节约了大量资源。这就是PCA主成分分析法的原理和代码实现

参考文章： www.cnblogs.com/mikewolf200…

欢迎关注个人我的公众号 AI计算机视觉工坊，本公众号不按期推送机器学习，深度学习，计算机视觉等相关文章，欢迎你们和我一块儿学习，交流。