每天进步一点点。

前言

说好的要做个有逼格的技术博客，虽然这篇依然没什么水平，但总算走出了第一步，希望以后每天都能进步一点点吧！

接触机器学习也一年多了，也学了很多算法，而PCA是数据预处理中一个很重要的算法，当时学习的时候也在网上看了很多资料，没想到一个简单的PCA也有着严密的数学推导，终于知道当年研究生考试为什么数学要考这么难了。

这篇文章主要是对PCA算法的一个总结，数学原理部分主要来自PCA的数学原理，并对其进行了总结（其实人家已经总结的很棒了），然后通过Scikit-learn中的一个例子讲解PCA的使用，最后用Python编写一个简单的PCA算法。

正文

基本概念

方差：用来衡量随机变量与其数学期望（均值）之间的偏离程度。统计中的方差（样本方差）是各个数据分别与其平均数之差的平方的和的平均数。

协方差：度量两个随机变量关系的统计量，协方差为0的两个随机变量是不相关的。

协方差矩阵：在统计学与概率论中，协方差矩阵的每个元素是各个向量元素之间的协方差。特殊的，矩阵对角线上的元素分别是向量的方差。

—— 摘自百度百科

特别是协方差矩阵的概念非常重要。

PCA介绍

主成分分析（Principal Component Analysis），是一种用于探索高维数据的技术。PCA通常用于高维数据集的探索与可视化。还可以用于数据压缩，数据预处理等。PCA可以把可能具有线性相关性的高维变量合成为线性无关的低维变量，称为主成分（principal components），新的低维数据集会尽可能的保留原始数据的变量，可以将高维数据集映射到低维空间的同时，尽可能的保留更多变量。

注意：降维就意味着信息的丢失，这一点一定要明确，如果用原始数据在模型上没有效果，期望通过降维来进行改善这是不现实的，不过鉴于实际数据本身常常存在的相关性，我们可以想办法在降维的同时将信息的损失尽量降低。当你在原数据上跑了一个比较好的结果，又嫌它太慢模型太复杂时候才可以采取PCA降维。

PCA数学理论

基

二维空间默认(1,0)和(0,1)为一组基。其实任何两个线性无关的二维向量都可以成为一组基。因为正交基有较好的性质，所以一般使用的基都是正交的。

例如：只有一个(3,2)本身是不能够精确表示一个向量的。这里的3实际表示的是向量在x轴上的投影值是3，在y轴上的投影值是2。也就是说我们其实隐式引入了一个定义：以x轴和y轴上正方向长度为1的向量为标准，默认选择(1,0)和(0,1)为基。

内积

向量A和B的内积公式为：

\[A\cdot B=|A||B|cos(a)\]

我们希望基的模是1，因为从内积的意义可以看到，如果基的模是1，那么就可以方便的用向量点乘基而直接获得其在新基上的坐标了，即 向量在基上的投影=向量与基的内积=坐标。

基变换的矩阵表示

将一组向量的基变换表示为矩阵的相乘。一般的，如果我们有M个N维向量，想将其变换为由R个N维向量（R个基）表示的新空间中，那么首先将R个基按行组成矩阵P，然后将待变换向量按列组成矩阵X，那么两矩阵的乘积就是变换结果。R可以小于N，而R决定了变换后数据的维数。也就是说，我们可以将一N维数据变换到更低维度的空间中去，变换后的维度取决于基的数量。因此这种矩阵相乘可以表示降维变换：

\[Y_{R \times M}=P_{R \times N} \times X_{N \times M}\]

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

优化目标

如何选择基才是最优的。或者说，如果我们有一组N维向量，现在要将其降到R维，那么我们应该如何选择R个基才能最大程度保留原有的信息？

对于一个二维空间：要在二维平面中选择一个方向，将所有数据都投影到这个方向所在直线上，用投影值表示原始记录。这是一个实际的二维降到一维的问题。那么如何选择这个方向才能尽量保留最多的原始信息呢？一种直观的看法是：希望投影后的投影值尽可能分散，而这种分散程度，可以用数学上的方差来表述。

对于上面二维降成一维的问题来说，找到那个使得方差最大的方向就可以了。不过对于更高维，还有一个问题需要解决。考虑三维降到二维问题，与之前相同，首先我们希望找到一个方向使得投影后方差最大，这样就完成了第一个方向的选择，继而我们选择第二个投影方向。如果我们还是单纯只选择方差最大的方向，很明显，这个方向与第一个方向应该是“几乎重合在一起”，显然这样的维度是没有用的，因此，应该有其他约束条件。从直观上说，让两个不同维度尽可能表示更多的原始信息，我们是不希望它们之间存在（线性）相关性的，因为相关性意味着两个维度不是完全线性独立，必然存在重复表示的信息。

数学上用协方差表示两个维度的相关性，当协方差为0时，表示两个维度完全独立。为了让协方差为0，我们选择第二个基时只能在与第一个基正交的方向上选择，因此最终选择的两个方向一定是正交的。

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

协方差矩阵

上面推导出优化目标，那么具体该怎么实现呢，下面就用到了协方差矩阵。回顾一下，协方差矩阵的每个元素是各个向量元素之间的协方差，特殊的，矩阵对角线上的元素分别是各个向量的方差。

设原始矩阵为X(N×M)，表示M个N维向量，其协方差矩阵为C(N×N)；P(R×N)为变换矩阵；Y(R×M)为目标矩阵, 其协方差矩阵为D。我们要求降维后的矩阵Y的每一维包含的数据足够分散，也就是每一行（维）方差足够大，而且要求行之间的元素线性无关，也就是要求行之间的协方差全部为0，这就要求协方差矩阵D的对角线元素足够大，除对角线外元素都为0。相当于对C进行协方差矩阵对角化。

具体推导如下：

\[D=\frac{1}{M}YY'=\frac{1}{M}PXX'P'=PCP'\]

C是X的协方差矩阵，是实对称矩阵，整个PCA降维过程其实就是一个实对称矩阵对角化的过程。

PCA具体算法步骤

设有M个N维数据:

将原始数据按列组成N行M列矩阵X
将X的每一行进行零均值化，即减去每一行的均值
求出X的协方差矩阵C
求出协方差矩阵C的特征值及对应的特征向量，C的特征值就是Y的每维元素的方差，也是D的对角线元素，从大到小沿对角线排列构成D。
将特征向量按对应特征值大小从上到下按行排列成矩阵，根据实际业务场景，取前R行组成矩阵P
Y=PX即为降到R维后的目标矩阵

Scikit-learn PCA实例分析

Scikit-learn是Python下著名的机器学习库，关于它我在这里就不多做介绍了，反正很好很强大。

首先数据选用经典的手写字符数据。

from sklearn import datasets
digits = datasets.load_digits()
x = digits.data                                              #输入数据
y = digits.target                                            #输出数据

PCA的调用也很简单。

from sklearn import decomposition
pca = decomposition.PCA()
pca.fit(x)

可视化，matplotlib是Python下的绘图库，功能也是十分强大。

import matplotlib.pyplot as plt
plt.figure()
plt.plot(pca.explained_variance_, 'k', linewidth=2)
plt.xlabel('n_components', fontsize=16)
plt.ylabel('explained_variance_', fontsize=16)
plt.show()

pca.explained_variance_ 就是上面协方差矩阵D的对角线元素，如下图所示：

fig-01

至于到底降到多少维度，主要取决于方差，具体的方法可以采用交叉验证。

实现PCA算法

纸上得来终觉浅，绝知此事要躬行。

最后基于Nnmpy, Pandas, Matploylib实现PCA并可视化结果。整体代码很简单，按照上面总结的算法步骤一步一步地计算，Numpy和Pandas的功能很强大，你能想到的运算几乎都有。

整体的代码风格还是比较Pythonic的，主要为了养成一个良好的编程习惯。

首先定义异常类：

class DimensionValueError(ValueError):
    """定义异常类"""
    pass

定义PCA类：

class PCA(object):
    """定义PCA类"""
    def __init__(self, x, n_components=None):
        self.x = x
        self.dimension = x.shape[1]

        if n_components and n_components >= self.dimension:
            raise DimensionValueError("n_components error")

        self.n_components = n_components

接下来就是计算协方差矩阵，特征值，特征向量，为了方便下面的计算，我把特征值和特征向量整合在一个dataframe内，并按特征值的大小降序排列：

    def cov(self):
        """求x的协方差矩阵"""
        x_T = np.transpose(self.x)                           #矩阵转置
        x_cov = np.cov(x_T)                                  #协方差矩阵
        return x_cov

    def get_feature(self):
        """求协方差矩阵C的特征值和特征向量"""
        x_cov = self.cov()
        a, b = np.linalg.eig(x_cov)
        m = a.shape[0]
        c = np.hstack((a.reshape((m,1)), b))
        c_df = pd.DataFrame(c)
        c_df_sort = c_df.sort(columns=0, ascending=False)    
        return c_df_sort

最后就是降维，用了两种方式，指定维度降维和根据方差贡献率自动降维，默认方差贡献率为99%：

    def reduce_dimension(self):
        """指定维度降维和根据方差贡献率自动降维"""
        c_df_sort = self.get_feature()
        varience = self.explained_varience_()

        if self.n_components:                                #指定降维维度
            p = c_df_sort.values[0:self.n_components, 1:]
            y = np.dot(p, np.transpose(self.x))              
            return np.transpose(y)

        varience_sum = sum(varience)                         
        varience_radio = varience / varience_sum

        varience_contribution = 0
        for R in xrange(self.dimension):
            varience_contribution += varience_radio[R]       
            if varience_contribution >= 0.99:
                break

        p = c_df_sort.values[0:R+1, 1:]                      #取前R个特征向量
        y = np.dot(p, np.transpose(self.x))                  
        return np.transpose(y)

完整的代码已经push到了 Github，欢迎star和fork。

后记

希望读完本文对你有所帮助。

由于本人水平有限，文中难免有疏漏之处，还请大家批评指正。

本作品采用知识共享署名 3.0 中国大陆许可协议进行许可。