PCA主成分分析

PCA是什么

PCA（Principal Component Analysis） ，即主成分分析法，是一种常用的数据分析方法。PCA通过线性变换将原始数据变换为一组各维度线性无关的表示，可用于提取数据的主要特征分量，常用于高维数据的降维。

为什么进行数据降维？

一般情况下，在数据挖掘和机器学习中，数据被表示为向量。很多机器学习算法的复杂度和数据的维数有着密切关系，甚至与维数呈指数级关联。机器学习中常会处理成千上万维的数据，这会带来较大的资源消耗，因此我们必须对数据进行降维。

而降维当然意味着信息的丢失，不过鉴于实际数据本身常常存在的相关性，我们可以想办法在降维的同时将信息的损失尽量降低。PCA 即为一种具有严格数学基础且已被广泛应用的降维方法。

PCA理论基础

一、协方差矩阵

1. 方差与协方差定义

统计学中，方差常用来度量某个随机变量的离散程度，协方差度量两个随机变量的相似程度。方差的计算公式为：

$Var(x)$ = $\frac{1}{n-1}$ ${\sum }_{i=1}^n(x_i-\overline{x}{)}^2$

协方差的计算公式为：

$Cov(x,y)$ = $\frac{1}{n-1}$ ${\sum }_{i=1}^n(x_i-\overline{x})(y_i-\overline{y})$

其中 $\overline{x}$ 和 $\overline{y}$ 分别代表两个随机变量对应的观测样本平均值。
又有：

$Cov(x,x)=Var(x)$

2. 协方差矩阵

给定 $d$ 个随机变量 $x_k,k=1,2......,d$ ， 其中 $x_{ki}$ 表示随机变量 $x_k$ 中的第 $i$ 个观测样本,每个随机变量对应的样本量为 $n$ 。
这些随机变量的方差为：

$Var(x_k)$ = $\frac{1}{n-1}$ ${\sum }_{i=1}^n(x_{ki}-{\overline{x}}_k{)}^2$

协方差为：

$Cov(x_m,x_k)$ = $\frac{1}{n-1}$ ${\sum }_{i=1}^n(x_{mi}-{\overline{x}}_m)(y_{ki}-{\overline{y}}_k)$

协方差矩阵为：
$$\Sigma =\left(\begin{array}{ccc}Cov(x_1,x_1)& \cdots & Cov(x_1,x_d)\\ ⋮& \ddots & ⋮\\ Cov(x_d,x_1)& \cdots & Cov(x_d,x_d)\end{array}\right)$$
其中，对角线上元素为各个随机变量的方差，而非对角线上的元素为两两随机变量之间的协方差。由于 $Cov(x_m,x_k)=Cov(x_k,x_m)$ ,$\Sigma$ 为 $d$ x $d$ 的对称矩阵，即：
$$\Sigma =\left(\begin{array}{ccc}Var(x_1)& \cdots & Cov(x_1,x_d)\\ ⋮& \ddots & ⋮\\ Cov(x_d,x_1)& \cdots & Var(x_d)\end{array}\right)$$
设有m个n维数据，将其按列展开组成矩阵X，即：

$$X=\left(\begin{array}{cccc}{a}_1& a_2& \cdots & a_n\\ b_1& b_2& \cdots & b_n\end{array}\right)$$
有：
$$\frac{1}{m}X{X}^T=\left(\begin{array}{cc}\frac{1}{m}{\sum }_{i=1}^m{a}_i^2& \frac{1}{m}{\sum }_{i=1}^m{a}_i{b}_i\\ & \\ \frac{1}{m}{\sum }_{i=1}^m{a}_i{b}_i& \frac{1}{m}{\sum }_{i=1}^m{b}_i^2\end{array}\right)$$

即X矩阵对应的协方差矩阵$\Sigma$为：
$$\Sigma =\frac{1}{m}X{X}^T$$

3. 协方差矩阵的特征值分解

对于任意对称矩阵 $\Sigma$ ，存在一特征值分解（EVD）:

$\Sigma =U\Lambda U^T$

其中 $U$ 的每一列都是相互正交的特征向量，且为单位向量，满足

$U^{T}U=I$

$\Lambda$ 对角线上的元素由大到小排列的特征值，非对角线上元素为0

4. 相关系数

对于二维随机变量 $(x,y)$ ，各自方差为：

$Var(x)={\sigma }_x^2$ , $Var(y)={\sigma }_y^2$

则随机变量 x 和 y 的相关系数为：

${\rho }_{xy}=\frac{Cov(x,y)}{{\sigma }_x{\sigma }_y}$

二、向量及基变换

1. 内积运算的意义

两个维数相同的向量的内积为：

$(a_1,a_2,...,a_n{)}^T\cdot (b_1,b_2,...,b_n{)}^T=a_1{b}_1+a_2{b}_2+....+a_n{b}_n$

以二维向量为例讨论内积运算的意义，假设二维向量$A=(x_1,y_1),B=(x_2,y_2)$,则它们的内积为：

$A\cdot B=|A||B|cos(\theta )$

若令 $|B|=1$ ，则内积变成：

$A\cdot B=|A|cos(\theta )$

即，设向量B的模为1，则A与B的内积值等于A在B所在直线上投影的矢量长度。

2. 基

我们已知，要准确描述向量，首先要确定一组基，给定向量在各个基所在直线上的投影值。我们总默认以 (1,0) 和 (0,1) 为基。

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

例如，(1,1)和(-1,1)也可以成为一组基。一般来说，我们希望基的模是1，因为从内积的意义可以看到，如果基的模是1，那么就可以方便的用向量点乘基而直接获得其在新基上的坐标了。同时，常选择正交基，因为正交基具有较好的性质

3. 基变换的矩阵表示

一般的，如果我们有M个N维向量，想将其变换为由R个N维向量表示的新空间中，那么首先将R个基按行组成矩阵A，然后将向量按列组成矩阵B，那么两矩阵的乘积AB就是变换结果，其中AB的第m列为A中第m列变换后的结果。

数学表示为：

$$\left(\begin{array}{c}{e}_1\\ e_2\\ ⋮\\ e_R\end{array}\right)\left(\begin{array}{cccc}{a}_1& a_2& \cdots & a_M\end{array}\right)=\left(\begin{array}{cccc}{e}_1{a}_1& e_1{a}_2& \cdots & e_1{a}_M\\ e_2{a}_1& e_2{a}_2& \cdots & e_2{a}_M\\ ⋮& ⋮& \ddots & ⋮\\ e_R{a}_1& e_R{a}_2& \cdots & e_R{a}_M\end{array}\right)$$

其中 $e_i$ 是一个行向量，表示第 i 个基， $a_j$ 是一个列向量，表示第 j 个原始数据记录。

特别要注意的是，这里R可以小于N，而R决定了变换后数据的维数。也就是说，我们可以将一N维数据变换到更低维度的空间中去，变换后的维度取决于基的数量。因此这种矩阵相乘的表示也可以表示降维变换。

最后，上述分析同时给矩阵相乘找到了一种物理解释：两个矩阵相乘的意义是将右边矩阵中的每一列列向量变换到左边矩阵中每一行行向量为基所表示的空间中去。更抽象的说，一个矩阵可以表示一种线性变换。

PCA步骤

如果基的数量少于向量本身的维数，则可以达到降维的效果。那么如何选择基？

1. 方差与协方差在PCA的应用

以二维数据降一维为例：

即，现有一组二位数据分布在二维平面内，先要在二维平面内找到一个基，将所有数据投影到基所在的直线上，用投影值表示原始记录。

那么如何选择基才能在降维后尽量保留更多的原始信息？我们希望投影后的投影值尽可能分散。

上文已提到，方差常用来表示数据的分散程度，所以问题就转化成了：寻找一个一维基，使得所有数据投影到这个基上后，投影值坐标的方差最大。

以三维降二维为例：
与二维降一维一样，我们还希望找到第一个基，使数据在该方向上投影的方差最大。但是 如何选择第二个基？

如果找第二个基时还选择使投影方差最大的，那么第二个基与第一个基几乎是重合的，这样的维度显然没有意义。

数据在两个基上的投影值组成两个不同的字段，想让这两个字段尽可能保留更多的原始信息，则不希望他们之间存在线性相关性，因为相关性表示两个字段不完全独立，有重复的信息。上文提到，协方差常用来表示两个数据之间的相关性。先用两个字段的协方差表示相关性，为使两个字段完全独立，则它们的协方差应为0，这时第二个基只能与第一个基正交。

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

2. 协方差矩阵对角化

现在已有优化目标，那么需要怎么做，才能得到我们想要的基?
假设P是一组基按行展开的矩阵，X是原始数据按列展开的矩阵，Y是X对P做基变换后的数据组成的矩阵，原始数据矩阵X对应的协方差矩阵为C，矩阵Y对应的协方差矩阵为D：

有：
$$P\cdot X=Y$$
即：

$$\left(\begin{array}{c}{e}_1\\ e_2\\ ⋮\\ e_R\end{array}\right)\left(\begin{array}{cccc}{x}_1& x_2& \cdots & x_m\end{array}\right)=\left(\begin{array}{cccc}{e}_1{x}_1& e_1{x}_2& \cdots & e_1{x}_m\\ e_2{x}_1& e_2{x}_2& \cdots & e_2{x}_m\\ ⋮& ⋮& \ddots & ⋮\\ e_R{x}_1& e_R{x}_2& \cdots & e_R{x}_m\end{array}\right)$$

其中 $e_i$ 是一个行向量，表示第 i 个基， $a_j$ 是一个列向量，表示第 j 个原始数据记录。

有:
$$C=\frac{1}{m}X{X}^T$$
$$D=\frac{1}{m}Y{Y}^T$$
$$D=\frac{1}{m}Y{Y}^T=\frac{1}{m}(PX)(PX{)}^T=\frac{1}{m}PX{X}^T{P}^T=PC{P}^T$$

至此，优化目标变成了寻找一个矩阵P，满足是一个对角矩阵，并且对角元素按从大到小依次排列，那么P的前K行就是要寻找的基，用P的前K行组成的矩阵乘以X就使得X从N维降到了K维并满足上述优化条件。

如何进行协方差矩阵对角化？

由上文知道，协方差矩阵C是一个是对称矩阵，有如下结论：

$U$ 的每一列都是相互正交的特征向量，且为单位向量

$$U=\left(\begin{array}{cccc}{e}_1& e_2& \cdots & e_n\end{array}\right)$$

对角矩阵$\Lambda$ 对角线上的元素由大到小排列的特征值，非对角线上元素为0

有：

$$C=U\Lambda U^T$$

所以：
$$\Lambda =U^{T}CU=\Sigma =\left(\begin{array}{cccc}{\lambda }_1& & & \\ & {\lambda }_2& & \\ & & \ddots & \\ & & & {\lambda }_n\end{array}\right)$$

因为有：

$$D=PC{P}^T$$
所以:
$$P=U^T$$

即P的每一行都为C的一个特征向量，P的前k行组成的矩阵乘原始数据矩阵X可得到降维后的数据矩阵Y。

PCA算法总结

设有m条n维数据。

1）将原始数据按列组成n行m列矩阵X
2）将X的每一行（代表一个属性字段）进行零均值化，即减去这一行的均值
3）求出协方差矩阵
4）求出协方差矩阵的特征值及对应的特征向量
5）将特征向量按对应特征值大小从上到下按行排列成矩阵，取前k行组成矩阵P
6）即为降维到k维后的数据

PCA代码实现

代码转自枫叶千言

#PCA.py
from numpy import *


def loadDataSet(fileName, delim='\t'):
    #dataMat = []
    fr = open(fileName)
    stringArr = [line.strip().split(delim) for line in fr.readlines()]
    datArr = [list(map(float, line)) for line in stringArr]
    return mat(datArr)


def pca(dataMat, topNfeat=9999999):  # topNfeat 降维后的维度
    meanVals = mean(dataMat, axis=0)  # 按列求均值，即每一列求一个均值，不同的列代表不同的特征
    # print meanVals
    meanRemoved = dataMat - meanVals  # remove mean     #去均值，将样本数据的中心点移到坐标原点
    print
    meanRemoved
    covMat = cov(meanRemoved, rowvar=0)  # 计算协方差矩阵
    # print covMat
    eigVals, eigVects = linalg.eig(mat(covMat))  # 计算协方差矩阵的特征值和特征向量
    # print eigVals
    # print eigVects
    eigValInd = argsort(eigVals)  # sort, sort goes smallest to largest  #排序，将特征值按从小到大排列
    # print eigValInd
    eigValInd = eigValInd[:-(topNfeat + 1):-1]  # cut off unwanted dimensions      #选择维度为topNfeat的特征值
    # print eigValInd
    redEigVects = eigVects[:, eigValInd]  # reorganize eig vects largest to smallest   #选择与特征值对应的特征向量
    print
    redEigVects
    lowDDataMat = meanRemoved * redEigVects  # transform data into new dimensions    #将数据映射到新的维度上，lowDDataMat为降维后的数据
    print
    lowDDataMat
    reconMat = (lowDDataMat * redEigVects.T) + meanVals  # 对原始数据重构，用于测试
    print
    reconMat
    return lowDDataMat, reconMat


def replaceNanWithMean():  # 均值代替那些样本中的缺失值
    datMat = loadDataSet('secom.data', ' ')
    numFeat = shape(datMat)[1]
    for i in range(numFeat):
        meanVal = mean(
            datMat[nonzero(~isnan(datMat[:, i].A))[0], i])  # values that are not NaN (a number) # .A表示把矩阵转化为数组array
        # nonzero(~isnan(datMat[:,i].A))[0] 返回非0元素所在行的索引；
        # >>> nonzero([True,False,True])
        #    (array([0, 2]),) 第0个和第3个元素非0
        # ~isnan()返回Ture or False
        datMat[nonzero(isnan(datMat[:, i].A))[0], i] = meanVal  # set NaN values to mean
    return datMat

#RunPca.py
import matplotlib.pyplot as plt
fig = plt.figure()
ax = fig.add_subplot(111)
import PCA
dataMat = PCA.loadDataSet('a.txt')
lowDMat, reconMat = PCA.pca(dataMat, 1)
ax.scatter(dataMat[:, 0].flatten().A[0], dataMat[:, 1].flatten().A[0], marker='^', s=90)
ax.scatter(reconMat[:, 0].flatten().A[0], reconMat[:, 1].flatten().A[0], marker='o', s=50, c='red')
plt.show()  # 由两维降为1维数据，降维后为一条红色直线，该方向是样本方差最大的方向，即样本离散程度最大的方向，该方向，将原来的2维数据融合为1维上

测试结果：

参考

PCA数学原理
PCA降维–代码实现