《机器学习实战》之十三——利用PCA来简化数据

PCA目录

一、前言
二、降维技术
三、PCA
- 1、PCA的数学原理
- - - （1）向量内积与投影
  - （2）基
  - （3）基变换
  - （4）方差与协方差
- 2、PCA算法步骤
- 3、在numpy中实现PCA
- - （1）零均值化
  - （2）求协方差矩阵
  - （3）求特征值、特征矩阵
  - （4）保留主要的成分[即保留值比较大的前n个特征]
- 4、选择主成分个数
四、示例：利用PCA对半导体制造数据降维
五、总结
参考资料

一、前言

我们通过电视观看体育比赛，在电视的纯平显示器上有一个球。显示器大概包含了100万像素，而球则可能是由较少的像素组成的，比如1000像素。在大部分体育比赛中，我们关注的是给定时刻球的位置。人的大脑要想了解比赛的进展，就需要了解球在运动场中的位置，在这个场景中，人们实时地将显示器上的百万像素转换成了一个三维的图像，该图像就给出了运动场上球的位置。在这个过程中，人们已经将数据从一百万维降至了三维。

在上述的体育比赛中，人们面对的原本是百万像素的数据，但是只有球的三维位置才最重要，这就被称为 降维（dimensionality reduction）

二、降维技术

降维是对数据高维度特征的一种预处理方法。将高维度的数据保留下最重要的一些特征，去除噪声和不重要的特征，从而实现提升数据处理速度的目的。在实际的生产和应用中，降维在一定的信息损失范围内，可以为我们节省大量的时间和成本。降维也成为了应用非常广泛的数据预处理方法。

1、降维作用：

使得数据集更易使用
降低很多算法的计算开销
.去除噪声
使得结果易懂

2、常用的降维方法：

主成分分析PCA（Principal Component Analysis）。
在PCA中，数据从原来的坐标系转换到了新的坐标系，新坐标系的选择是由数据本身决定的。第一个新坐标轴选择的是原始数据中方差最大的方向，第二个坐标轴的选择和第一个坐标轴正交且有最大方差的方向。该过程一直重复，重复次数为原始数据中特征的数目。我们会发现，大部分方差都包含在最前面的几个新坐标轴中。因此，我们可以忽略余下的坐标轴，即对数据进行降维处理。
因子分析（Factor Analysis）
我们假设在观察数据的生成中有一些观察不到的隐变量。假设观察数据是这些隐变量和某些噪声的线性组合。那么隐变量的数据可能比观察数据的数目少，也就是说通过找隐变量就可以实现数据的降维。
独立成分分析ICA（Independent Component Analysis）
ICA假设数据是从N个数据源生成的，这一点和因子分析有些类似。假设数据为多个数据源的混合观察结果，这些数据源之间在统计上是相互独立的，而在PCA中只假设数据是不相关的。同因子分析一样，如果数据源的数目少于观察数据的数目，则可以实现降维过程。

在上述3种降维技术中，PCA的应用目前最为广泛。所以，我们主要来介绍PCA降维技术。

三、PCA

1、PCA的数学原理

一般情况下，在数据挖掘和机器学习中，数据被表示为向量。例如某个淘宝店2012年全年的流量及交易情况可以看成一组记录的集合，其中每一天的数据是一条记录，格式如下：

( 日期 , 浏览量 , 访客数 , 下单数 , 成交数 , 成交金额 ) \mathit{\left ( 日期, 浏览量, 访客数, 下单数, 成交数, 成交金额 \right )} (日期,浏览量,访客数,下单数,成交数,成交金额)

其中“日期”是一个记录标志而非度量值，而数据挖掘关心的大多是度量值，因此如果我们忽略日期这个字段后，我们得到一组记录，每条记录可以被表示为一个五维向量，其中一条看起来大约是这个样子：

( 500 , 240 , 25 , 13 , 2312.15 ) T \left ( 500,240,25,13,2312.15 \right )^{T} (500,240,25,13,2312.15)T

注意这里我用了转置，因为习惯上使用列向量表示一条记录（后面会看到原因），本文后面也会遵循这个准则。不过为了方便有时我会省略转置符号，但我们说到向量默认都是指列向量。

（1）向量内积与投影

下面先来看一个高中就学过的向量运算：内积。两个维数相同的向量的内积被定义为：

( a 1 , a 2 , . . . , a n ) T ⋅ ( b 1 , b 2 , . . . , b n ) = a 1 b 1 + a 2 b 2 + . . . + a n b n {\color{Red} \mathit{\left ( a_{1},a_{2},...,a_{n}\right )^{T}\cdot \left ( b_{1},b_{2},...,b_{n}\right )=a_{1}b_{1}+a_{2}b_{2}+...+a_{n}b_{n}}} (a1,a2,...,an)T⋅(b1,b2,...,bn)=a1b1+a2b2+...+anbn

内积运算将两个向量映射为一个实数。其计算方式非常容易理解，但是其意义并不明显。下面我们分析内积的几何意义。假设A和B是两个n维向量，我们知道n维向量可以等价表示为n维空间中的一条从原点发射的有向线段，为了简单起见我们假设A和B均为二维向量，则A=(x1,y1)，B=(x2,y2)。则在二维平面上A和B可以用两条发自原点的有向线段表示，见下图：

好，现在我们从A点向B所在直线引一条垂线。我们知道垂线与B的交点叫做A在B上的投影，再设A与B的夹角是a，则投影的矢量长度为 ∣ A ∣ c o s ( a ) {\color{Red} \left | A \right |cos(a)} ∣A∣cos(a)，其中是向量A的模 ∣ A ∣ = x 1 2 + y 1 2 {\color{Red} \left | A \right |=\sqrt{x_{1}^{2}+y_{1}^{2}}} ∣A∣=x12+y12 ，也就是A线段的标量长度。

注意这里我们专门区分了矢量长度和标量长度，标量长度总是大于等于0，值就是线段的长度；而矢量长度可能为负，其绝对值是线段长度，而符号取决于其方向与标准方向相同或相反。

到这里还是看不出内积和这东西有什么关系，不过如果我们将内积表示为另一种我们熟悉的形式：

A ⋅ B = ∣ A ∣ ∣ B ∣ c o s ( a ) {\color{Red} A\cdot B=\left | A \right |\left | B \right |cos(a)} A⋅B=∣A∣∣B∣cos(a)

现在事情似乎是有点眉目了：A与B的内积等于A到B的投影长度乘以B的模。再进一步，如果我们假设B的模为1，即让|B|=1|B|=1，那么就变成了：

A ⋅ B = ∣ A ∣ c o s ( a ) {\color{Red} A\cdot B=\left | A \right |cos(a)} A⋅B=∣A∣cos(a)

也就是说，设向量B的模为1，则A与B的内积值等于A向B所在直线投影的矢量长度！这就是内积的一种几何解释，也是我们得到的第一个重要结论。在后面的推导中，将反复使用这个结论。

（2）基

下面我们继续在二维空间内讨论向量。上文说过，一个二维向量可以对应二维笛卡尔直角坐标系中从原点出发的一个有向线段。例如下面这个向量：

在代数表示方面，我们经常用线段终点的点坐标表示向量，例如上面的向量可以表示为(3,2)，这是我们再熟悉不过的向量表示。

不过我们常常忽略，只有一个(3,2)本身是不能够精确表示一个向量的。我们仔细看一下，这里的3实际表示的是向量在x轴上的投影值是3，在y轴上的投影值是2。也就是说我们其实隐式引入了一个定义：以x轴和y轴上正方向长度为1的向量为标准。那么一个向量(3,2)实际是说在x轴投影为3而y轴的投影为2。注意投影是一个矢量，所以可以为负。

更正式的说，向量(x,y)实际上表示线性组合：

x ( 1 , 0 ) T + y ( 0 , 1 ) T {\color{Red} x\left ( 1,0 \right )^{T}+y\left ( 0,1 \right )^{T}} x(1,0)T+y(0,1)T

不难证明所有二维向量都可以表示为这样的线性组合。此处(1,0)和(0,1)叫做二维空间中的一组基。

所以，要准确描述向量，首先要确定一组基，然后给出在基所在的各个直线上的投影值，就可以了。只不过我们经常省略第一步，而默认以(1,0)和(0,1)为基。

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

例如，(1,1)和(-1,1)也可以成为一组基。一般来说，我们希望基的模是1，因为从内积的意义可以看到，如果基的模是1，那么就可以方便的用向量点乘基而直接获得其在新基上的坐标了！实际上，对应任何一个向量我们总可以找到其同方向上模为1的向量，只要让两个分量分别除以模就好了。例如，上面的基可以变为 ( 1 2 , 1 2 ) {\color{Red} (\frac{1}{\sqrt{2}},\frac{1}{\sqrt{2}})} (2 1,2 1)和 ( − 1 2 , 1 2 ) {\color{Red} (-\frac{1}{\sqrt{2}},\frac{1}{\sqrt{2}})} (−2 1,2 1)

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

另外这里要注意的是，我们列举的例子中基是正交的（即内积为0，或直观说相互垂直），但可以成为一组基的唯一要求就是线性无关，非正交的基也是可以的。不过因为正交基有较好的性质，所以一般使用的基都是正交的。

（3）基变换

下面我们找一种简便的方式来表示基变换。还是拿上面的例子，想一下，将(3,2)变换为新基上的坐标，就是用(3,2)与第一个基做内积运算，作为第一个新的坐标分量，然后用(3,2)与第二个基做内积运算，作为第二个新坐标的分量。实际上，我们可以用矩阵相乘的形式简洁的表示这个变换：

( 1 2 1 2 − 1 2 1 2 ) ( 3 2 ) = ( 5 2 − 1 2 ) {\color{Red} \begin{pmatrix} \frac{1}{\sqrt{2}} &\frac{1}{\sqrt{2}} \\ \frac{-1}{\sqrt{2}} & \frac{1}{\sqrt{2}} \end{pmatrix}\begin{pmatrix} 3\\ 2 \end{pmatrix}=\begin{pmatrix} \frac{5}{\sqrt{2}}\\ \frac{-1}{\sqrt{2}} \end{pmatrix}} (2 12 −12 12 1)(32)=(2 52 −1)

太漂亮了！其中矩阵的两行分别为两个基，乘以原向量，其结果刚好为新基的坐标。可以稍微推广一下，如果我们有m个二维向量，只要将二维向量按列排成一个两行m列矩阵，然后用“基矩阵”乘以这个矩阵，就得到了所有这些向量在新基下的值。例如(1,1)，(2,2)，(3,3)，想变换到刚才那组基上，则可以这样表示：

( 1 2 1 2 − 1 2 1 2 ) ( 1 2 3 1 2 3 ) = ( 2 2 4 2 6 2 0 0 0 ) {\color{Red} \begin{pmatrix} \frac{1}{\sqrt{2}} &\frac{1}{\sqrt{2}} \\ \frac{-1}{\sqrt{2}} & \frac{1}{\sqrt{2}} \end{pmatrix}\begin{pmatrix} 1 & 2 &3 \\ 1 & 2 &3 \end{pmatrix}=\begin{pmatrix} \frac{2}{\sqrt{2}} &\frac{4}{\sqrt{2}} &\frac{6}{\sqrt{2}} \\ 0 & 0 &0 \end{pmatrix}} (2 12 −12 12 1)(112233)=(2 202 402 60)
于是一组向量的基变换被干净的表示为矩阵的相乘。

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

数学表示为：
( p 1 p 2 . . . p R ) ( a 1 a 2 . . . a M ) = ( p 1 a 1 p 1 a 2 . . . p 1 a M p 2 a 1 p 2 a 2 . . . p 2 a M . . . . . . . . . . . . p R a 1 p R a 2 . . . p R a M ) {\color{Red}\begin{pmatrix} p_{1}\\ p_{2}\\ .\\ .\\ .\\ p_{R} \end{pmatrix}\begin{pmatrix} a_{1} &a_{2} &. &. &. &a_{M} \end{pmatrix}=\begin{pmatrix} p_{1}a_{1} &p_{1}a_{2} &. &. &. &p_{1}a_{M} \\ p_{2}a_{1}&p_{2}a_{2} &. &. &. &p_{2}a_{M} \\ .& .& .& & &. \\ .& .& & .& &. \\ .& .& & & .&. \\ p_{R}a_{1}&p_{R}a_{2} & .& .& .&p_{R}a_{M} \end{pmatrix} } ⎝⎜⎜⎜⎜⎜⎜⎛p1p2...pR⎠⎟⎟⎟⎟⎟⎟⎞(a1a2...aM)=⎝⎜⎜⎜⎜⎜⎜⎛p1a1p2a1...pRa1p1a2p2a2...pRa2............p1aMp2aM...pRaM⎠⎟⎟⎟⎟⎟⎟⎞
其中 p i p_{i} pi 是一个行向量，表示第i个基， a j a_{j} aj 是一个列向量，表示第 j个原始数据记录。

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

最后，上述分析同时给矩阵相乘找到了一种物理解释：两个矩阵相乘的意义是将右边矩阵中的每一列列向量变换到左边矩阵中每一行行向量为基所表示的空间中去。更抽象的说，一个矩阵可以表示一种线性变换。很多同学在学线性代数时对矩阵相乘的方法感到奇怪，但是如果明白了矩阵相乘的物理意义，其合理性就一目了然了。

（4）方差与协方差

上面我们讨论了选择不同的基可以对同样一组数据给出不同的表示，而且如果基的数量少于向量本身的维数，则可以达到降维的效果。但是我们还没有回答一个最最关键的问题：如何选择基才是最优的。或者说，如果我们有一组N维向量，现在要将其降到K维（K小于N），那么我们应该如何选择K个基才能最大程度保留原有的信息？

要完全数学化这个问题非常繁杂，这里我们用一种非形式化的直观方法来看这个问题。

为了避免过于抽象的讨论，我们仍以一个具体的例子展开。假设我们的数据由五条记录组成，将它们表示成矩阵形式：

( 1 1 2 4 2 1 3 3 4 4 ) {\color{Red} \begin{pmatrix} 1 & 1 &2 &4 &2 \\ 1 & 3&3 &4 &4 \end{pmatrix} } (1113234424)
其中每一列为一条数据记录，而一行为一个字段。为了后续处理方便，我们首先将每个字段内所有值都减去字段均值，其结果是将每个字段都变为均值为0（这样做的道理和好处后面会看到）。

我们看上面的数据，第一个字段均值为2，第二个字段均值为3，所以变换后：

( − 1 − 1 0 2 0 − 2 0 0 1 1 ) {\color{Red} \begin{pmatrix} -1 & -1 &0 &2 &0 \\ -2 & 0&0 &1 &1 \end{pmatrix} } (−1−2−10002101)
我们可以看下五条数据在平面直角坐标系内的样子：

现在问题来了：如果我们必须使用一维来表示这些数据，又希望尽量保留原始的信息，你要如何选择？

通过上面对基变换的讨论我们知道，这个问题实际上是要在二维平面中选择一个方向，将所有数据都投影到这个方向所在直线上，用投影值表示原始记录。这是一个实际的二维降到一维的问题。

那么如何选择这个方向（或者说基）才能尽量保留最多的原始信息呢？一种直观的看法是：希望投影后的投影值尽可能分散。

以上图为例，可以看出如果向x轴投影，那么最左边的两个点会重叠在一起，中间的两个点也会重叠在一起，于是本身四个各不相同的二维点投影后只剩下两个不同的值了，这是一种严重的信息丢失，同理，如果向y轴投影最上面的两个点和分布在x轴上的两个点也会重叠。所以看来x和y轴都不是最好的投影选择。我们直观目测，如果向通过第一象限和第三象限的斜线投影，则五个点在投影后还是可以区分的。

下面，我们用数学方法表述这个问题。

- 方差

上文说到，==我们希望投影后投影值尽可能分散，而这种分散程度，可以用数学上的方差来表述。==此处，一个字段的方差可以看做是每个元素与字段均值的差的平方和的均值，即：

V a r ( a ) = 1 m ∑ i = 1 m ( a i − μ ) 2 {\color{Red} Var(a)=\frac{1}{m}\sum_{i=1}^{m}(a_{i}-\mu )^2} Var(a)=m1i=1∑m(ai−μ)2

由于上面我们已经将每个字段的均值都化为0了，因此方差可以直接用每个元素的平方和除以元素个数表示：

V a r ( a ) = 1 m ∑ i = 1 m a i 2 {\color{Red} Var(a)=\frac{1}{m}\sum_{i=1}^{m}a_{i}^2} Var(a)=m1i=1∑mai2

于是上面的问题被形式化表述为：寻找一个一维基，使得所有数据变换为这个基上的坐标表示后，方差值最大。

- 协方差

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

如果我们还是单纯只选择方差最大的方向，很明显，这个方向与第一个方向应该是“几乎重合在一起”，显然这样的维度是没有用的，因此，应该有其他约束条件。从直观上说，让两个字段尽可能表示更多的原始信息，我们是不希望它们之间存在（线性）相关性的，因为相关性意味着两个字段不是完全独立，必然存在重复表示的信息。

数学上可以用两个字段的协方差表示其相关性，由于已经让每个字段均值为0，则：

C o v ( a , b ) = 1 m ∑ i = 1 m a i b i {\color{Red}Cov(a,b)=\frac{1}{m}\sum_{i=1}^{m}a_{i}b_{i}} Cov(a,b)=m1i=1∑maibi

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

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

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

- 协方差矩阵

上面我们导出了优化目标，但是这个目标似乎不能直接作为操作指南（或者说算法），因为它只说要什么，但根本没有说怎么做。所以我们要继续在数学上研究计算方案。

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

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

X = ( a 1 a 2 . . . a m b 1 b 2 . . . b m ) {\color{Red} X=\begin{pmatrix} a_{1} &a_{2} &. &. &. &a_{m} \\ b_{1}&b_{2} &. &. &. &b_{m} \end{pmatrix}} X=(a1b1a2b2......ambm)

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

1 m X X T = ( 1 m ∑ i = 1 m a i 2 1 m ∑ i = 1 m a i b i 1 m ∑ i = 1 m a i b i 1 m ∑ i = 1 m b i 2 ) {\color{Red} \frac{1}{m}XX^{T}=\begin{pmatrix} \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{pmatrix}} m1XXT=(m1∑i=1mai2m1∑i=1maibim1∑i=1maibim1∑i=1mbi2)

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

根据矩阵相乘的运算法则，这个结论很容易被推广到一般情况：

设我们有m个n维数据记录，将其按列排成n乘m的矩阵X，设 C = 1 m X X T {\color{Red} C=\frac{1}{m}XX^{T}} C=m1XXT，则C是一个对称矩阵，其对角线分别是各个字段的方差，而第 i i i 行 j j j 列和 j j j 行 i i i 列元素相同，表示 i i i 和 j j j 两个字段的协方差。

- 协方差矩阵对角化

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

设原始数据矩阵X对应的协方差矩阵为C，而P是一组基按行组成的矩阵，设Y=PX，则Y为X对P做基变换后的数据。设Y的协方差矩阵为D，我们推导一下D与C的关系：

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

至此，我们离“发明”PCA还有仅一步之遥！

现在所有焦点都聚焦在了协方差矩阵对角化问题上，有时，我们真应该感谢数学家的先行，因为矩阵对角化在线性代数领域已经属于被玩烂了的东西，所以这在数学上根本不是问题。

由上文知道，协方差矩阵C是一个是对称矩阵，在线性代数上，实对称矩阵有一系列非常好的性质：

(1) 实对称矩阵不同特征值对应的特征向量必然正交。
(2) 设特征向量 λ \lambda λ 重数为 r r r，则必然存在 r r r 个线性无关的特征向量对应于 λ \lambda λ ，因此可以将这 r r r 个特征向量单位正交化。

由上面两条可知，一个n行n列的实对称矩阵一定可以找到n个单位正交特征向量，设这n个特征向量为 e 1 , e 2 , ⋯ , e n e_{1},e_{2},⋯,e_{n} e1,e2,⋯,en，我们将其按列组成矩阵：

E = ( e 1 , e 2 , ⋯ , e n ) {\color {Red} E=\left ( e_{1},e_{2},⋯,e_{n} \right ) } E=(e1,e2,⋯,en)

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

E T C E = Λ = ( λ 1 λ 2 . . . λ n ) {\color{Red} E^{T}CE =\Lambda =\begin{pmatrix} &\lambda _{1} & & & & \\ & &\lambda _{2} & & & \\ & & &. & & \\ & & & &. & \\ & & & & &. \\ & & & & & \lambda _{n} \end{pmatrix}} ETCE=Λ=⎝⎜⎜⎜⎜⎜⎜⎛λ1λ2...λn⎠⎟⎟⎟⎟⎟⎟⎞

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

以上结论不再给出严格的数学证明，对证明感兴趣的朋友可以参考线性代数书籍关于“实对称矩阵对角化”的内容。

到这里，我们发现我们已经找到了需要的矩阵P：

P = E T {\color {Red} P=E^{T}} P=ET

P是协方差矩阵的特征向量单位化后按行排列出的矩阵，其中每一行都是 C 的一个特征向量。如果设 P 按照 Λ \Lambda Λ 中特征值的从大到小，将特征向量从上到下排列，则用P的前K行组成的矩阵乘以原始数据矩阵X，就得到了我们需要的降维后的数据矩阵Y。

至此我们完成了整个PCA的数学原理讨论

2、PCA算法步骤

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

3、在numpy中实现PCA

还是老规矩，我们先打开“testSet.txt”文件，看一下数据集是什么样的。该数据集中有两列，都是数值类型的。

新建一个名为pca.py的文件，并在其中写入如下代码：

# -*- coding: utf-8 -*-import numpy as np
import matplotlib.pyplot as pltdef loadDataSet(fileName, delim='\t'):fr = open(fileName)stringArr = [line.strip().split(delim) for line in fr.readlines()]datArr = [list(map(float, line)) for line in stringArr]return np.mat(datArr)"""
函数说明：数据可视化
Parameters：fileName：数据文件名
Returns：无
"""
def ShowData(fileName):dataMat = loadDataSet(fileName)fig = plt.figure()ax = fig.add_subplot(111)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()if __name__ == '__main__':ShowData('testSet.txt')

运行结果如下所示，这图中展示的是数据的分布情况：

（1）零均值化

假如原始数据集为矩阵dataMat，dataMat中每一行代表一个样本，每一列代表同一个特征。零均值化就是求每一列的平均值，然后该列上的所有数都减去这个均值。也就是说，这里零均值化是对每一个特征而言的，零均值化后，每个特征的均值变成0。实现代码如下：

"""
函数说明：零均值化
Parameters：dataMat：数据集
Returns：newData：去除均值后的数据meanVal：平均值
"""
def zeroMean(dataMat):meanVal = np.mean(dataMat, axis=0)  #按列求平均值newData = dataMat - meanVal  #减去平均值return newData, meanVal

（2）求协方差矩阵

    newData, meanVal = zeroMean(dataMat)covMat = np.cov(newData, rowvar=0)

numpy中的cov函数用于求协方差矩阵，参数rowvar很重要！若rowvar=0，说明传入的数据一行代表一个样本，若非0，说明传入的数据一列代表一个样本。因为newData每一行代表一个样本，所以将rowvar设置为0。covMat即所求的协方差矩阵。

（3）求特征值、特征矩阵

    eigVals, eigVects = np.linalg.eig(np.mat(covMat))

调用numpy中的线性代数模块linalg中的eig函数，可以直接由covMat求得特征值和特征向量：eigVals存放特征值，行向量。eigVects存放特征向量，每一列代表一个特征向量。特征值和特征向量是一一对应的。

（4）保留主要的成分[即保留值比较大的前n个特征]

第三步得到了特征值向量eigVals，假设里面有m个特征值，我们可以对其排序，排在前面的n个特征值所对应的特征向量就是我们要保留的，它们组成了新的特征空间的一组基n_eigVect。将零均值化后的数据乘以n_eigVect就可以得到降维后的数据。代码如下：

    eigValIndice=np.argsort(eigVals)            #对特征值从小到大排序n_eigValIndice=eigValIndice[-1:-(n+1):-1]   #最大的n个特征值的下标n_eigVect=eigVects[:,n_eigValIndice]        #最大的n个特征值对应的特征向量lowDDataMat=newData*n_eigVect               #低维特征空间的数据reconMat=(lowDDataMat*n_eigVect.T)+meanVal  #重构数据return lowDDataMat,reconMat

代码中有几点要说明一下，首先argsort对特征值是从小到大排序的，那么最大的n个特征值就排在后面，所以 eigValIndice[-1:-(n+1):-1] 就是根据特征值排序结果的逆序，取出这个n个最大的特征向量。【python里面，list[a：b：c] 代表从下标a开始到b，步长为c。】

OK，这四步下来就可以从高维的数据dataMat得到低维的数据lowDDataMat，另外，程序也返回了重构数据reconMat，有些时候reconMat更便于数据分析。

总的代码如下：

# -*- coding: utf-8 -*-import numpy as np
import matplotlib.pyplot as pltdef loadDataSet(fileName, delim='\t'):fr = open(fileName)stringArr = [line.strip().split(delim) for line in fr.readlines()]datArr = [list(map(float, line)) for line in stringArr]return np.mat(datArr)"""
函数说明：数据可视化
Parameters：fileName：数据文件名
Returns：无
"""
def ShowData(fileName):dataMat = loadDataSet(fileName)lowDMat, reconMat = pca(dataMat, 1)fig = plt.figure()ax = fig.add_subplot(111)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()
"""
函数说明：零均值化
Parameters：dataMat：数据集
Returns：newData：去除均值后的数据meanVal：每个特征的平均值
"""
def zeroMean(dataMat):meanVal = np.mean(dataMat, axis=0)  #按列求平均值newData = dataMat - meanVal  #减去平均值return newData, meanVal
"""
函数说明：
Parameters：dataMat：数据集topNfeat：最大的特征数
Returns：lowDDataMat：低维特征空间的数据reconMat：重构数据
"""
def pca(dataMat, topNfeat=9999999):newData, meanVal = zeroMean(dataMat) #零均值化covMat = np.cov(newData, rowvar=0) #求协方差矩阵eigVals, eigVects = np.linalg.eig(np.mat(covMat))  #求特征值、特征矩阵eigValIndice = np.argsort(eigVals) #对特征值从小到大排序n_eigValIndice = eigValIndice[-1:-(topNfeat+1):-1] #最大的topNfeat个特征值的下标n_eigVect = eigVects[:, n_eigValIndice]  #最大的topNfeat个特征值对应的特征向量lowDDataMat = newData * n_eigVect #低维特征空间数据reconMat = (lowDDataMat * n_eigVect.T) + meanVal #重构数据return lowDDataMat, reconMatif __name__ == '__main__':ShowData('testSet.txt')

这里对testSet.txt文件中的1000个数据点组成的数据，进行PCA操作，topNfeat取值为1，得到的结果如下图所示：

如果topNfeat取值为2，那么就是不剔除任何特征，那么重构后的数据会和原始的数据重合。

4、选择主成分个数

文章写到这里还没有完，应用PCA的时候，对于一个1000维的数据，我们怎么知道要降到几维的数据才是合理的？即 topNfeat 要取多少，才能保留最多信息同时去除最多的噪声？一般，我们是通过方差百分比来确定 topNfeat 的，这一点在Ufldl教程中说得很清楚，并且有一条简单的公式，下面是该公式：

根据这条公式，可以写个函数，函数传入的参数是百分比percentage和特征值向量，然后根据percentage确定n，代码如下：

"""
函数说明：按照百分比来获取特征ges
Parameters：eigVals：特征值percentage：百分比
Returns：num：特征个数
"""
def percentage2n(eigVals, percentage):sortArray = np.sort(eigVals)  #升序sortArray = sortArray[-1::-1]  #降序arraySum = sum(sortArray)  #求特征值的和tmpSum = 0num = 0for i in sortArray:tmpSum += i   #最大k个特征值的总和num += 1if tmpSum >= arraySum*percentage:return num

那么pca函数也可以重写成百分比版本，默认百分比99%。

"""
函数说明：
Parameters：dataMat：数据集percentage：百分比
Returns：lowDDataMat：低维特征空间的数据reconMat：重构数据
"""
def pca(dataMat, percentage=0.99):newData, meanVal = zeroMean(dataMat) #零均值化covMat = np.cov(newData, rowvar=0) #求协方差矩阵eigVals, eigVects = np.linalg.eig(np.mat(covMat))  #求特征值、特征矩阵n = percentage2n(eigVals, percentage) #要达到percent的方差百分比，需要前n个特征向量eigValIndice = np.argsort(eigVals) #对特征值从小到大排序n_eigValIndice = eigValIndice[-1:-(n+1):-1] #最大的topNfeat个特征值的下标n_eigVect = eigVects[:, n_eigValIndice]  #最大的topNfeat个特征值对应的特征向量lowDDataMat = newData * n_eigVect #低维特征空间数据reconMat = (lowDDataMat * n_eigVect.T) + meanVal #重构数据return lowDDataMat, reconMat

四、示例：利用PCA对半导体制造数据降维

半导体是在一些极为先进的工厂中制造出来的。工厂或制造设备不仅需要花费上亿美元，而且还需要大量的工人。制造设备仅能在几年内保持其先进性，随后就必须更换了。单个集成电路的加工时间会超过一个月。在设备生命期有限，花费又极其巨大的情况下，制造过程中的每一秒都价值巨大。如果制造过程中存在瑕疵，我们就必须尽早发现，从而确保宝贵的时间不会花费在缺陷产品的生产上。

一些工程上的通用解决方案是通过早期测试和频繁测试来发现有缺陷的产品，但仍然有一些存在瑕疵的产品通过了测试。如果机器学习技术能够用于进一步减少错误，那么它就会为制造商节省大量的资金。

打开数据集，我们看到它包含了很多的特征，具体地讲，它拥有590个特征。我们的任务就是对这些特征进行降维处理。

该数据包含了很多的缺失值。这些缺失值是以NaN（Not a Number的缩写）标识的。对于这些缺失值的处理办法，这里，我们用平均值来代替缺失值，平均值根据那些非NaN得到。
打开pac.py文件，写入如下的代码：

"""
函数说明：将数据集中的NaN替换为平均值
Parameters：无
Returns：datMat：处理之后的数据集
"""
def replaceNanWithMean():datMat = loadDataSet('secom.data',' ')numFeat = np.shape(datMat)[1]for i in range(numFeat):meanVal = np.mean(datMat[np.nonzero(~np.isnan(datMat[:,i].A))[0], i])datMat[np.nonzero(np.isnan(datMat[:,i].A))[0], i] = meanVal      return datMatif __name__ == '__main__':
#    ShowData('testSet.txt')dataMat = replaceNanWithMean()meanRemoved, meanVal = zeroMean(dataMat) #零均值化covMat = np.cov(meanRemoved, rowvar=0)  #计算协方差矩阵eigVals, eigVects = np.linalg.eig(np.mat(covMat)) #特征值分析print("90%的主成分方差总和：",sum(eigVals)*0.9)     # 计算90%的主成分方差总和print("前6个主成分的方差和：",sum(eigVals[:6]))     # 计算前6个主成分所占的方差和percentage = []for i in range(21):p1 = eigVals[i]/sum(eigVals)p2 = sum(eigVals[:i])/sum(eigVals)percentage.append(p1*100)print(i,"主成分的方差百分比：{:.2%}".format(float(p1)))print(i,"主成分的累积方差百分比：{:.2%}".format(float(p2)))plt.rcParams['font.sans-serif']=['SimHei'] #用来正常显示中文标签plt.rcParams['axes.unicode_minus']=False #用来正常显示负号plt.plot(percentage)    plt.xticks(np.linspace(0,20,5))  #设置x坐标轴的刻度，范围0-20，生成共5个点plt.title('方差百分比')plt.xlabel('主成分数目')plt.ylabel('方差的百分比')plt.show()

运行结果如下：

绘图输出了前20个主成分的总方差的百分比，从图中可以看出，大部分的方差都包含在前面的几个主成分中，舍弃后面的主成分并不会损失太多的信息。我们还可以打印出累计方差的百分比。跟书上是一样的。

上述分析能够得到所有用到的主成分数目，然后我们可以将该数目输入到PCA算法中，最后得到约简后的数据就可以在分类器中使用了。

五、总结

降维技术使得数据变得更容易使用，并且它们往往能去除数据中的噪声
PCA可以从数据中识别出主要特征，他是通过沿着数据最大方差方向旋转坐标轴来实现的，选择方差最大的方向作为第一条坐标轴，后续坐标轴则与前面的坐标轴正交。

参考资料

【1】PCA的数学原理 http://blog.codinglabs.org/articles/pca-tutorial.html
【2】《机器学习实战》第十三章利用PCA来简化数据
【3】UFLDL教程： http://ufldl.stanford.edu/wiki/index.php/主成分分析
【4】博客： https://blog.csdn.net/u012162613/article/details/42177327
【5】Markdown 文本居中、字体颜色以及数学公式： https://blog.csdn.net/silver1225/article/details/89036250#_140