失眠网 > 全面理解主成分分析（PCA）和MNIST数据集的Python降维实现

全面理解主成分分析（PCA）和MNIST数据集的Python降维实现

时间：2020-07-06 05:47:44

注：本博文为原创博文，如需转载请注明原创链接！！！

这篇博文主要讲述主成分分析的原理并用该方法来实现MNIST数据集的降维。

一、引言

主成分分析是一种降维和主成分解释的方法。举一个比较容易理解的例子，如果将三维世界的可乐罐子踩一脚变成二维的，踩的过程就是降维。可以有很多种方法，比如将可乐罐子立起来从上向下踩，或者是将罐子平躺后再踩，甚至我们可以斜着踩或是选择不同的角度。那么如何踩这个可乐罐子可以保存更多的信息呢？显然不是竖着踩，而是平躺着踩下去才会保留可乐罐子的商标、颜色等一些重要信息。

正如上面这个例子，在实际的问题研究中，数据往往会涉及到很多的变量，将多变量的大量数据之间的规律研究得透彻这无疑是一件好事，但不可否认在一定程度上也加重了对数据进行分析和解释的任务量。一般情况，变量与变量之间是存在一定相关性的，提供的信息在一定情况下有所重叠。如果单独拿出个别的变量（指标）进行分析往往是孤立的，没有完全利用数据中的信息，会产生错误的结论。

由于各变量之间存在着一定的相关性，因此我们希望用较少的互不相关的新变量来反映原变量提供的绝大部分信息。主成分分析（Principal Component Analysis，简称PCA）是一种通过降维技术把多变量化为少数几个主成分的统计分析方法。

所谓降维就是对高维的数据特性进行预处理的一种方法，将高维数据中的一些重要的特征保留下来并去除一些不重要的噪声，从而提升了数据处理的速度。

二、理解主成分分析的误区

在详细介绍主成分分析的原理之前需要知道的误区：

主成分分析不是抽取数据中的某些特征而是在原有的特征中通过某种数学思想重建新的特征。新特征的特点就是数量比原有特征要少但大概率涵盖了原有的大部分特征的信息。主成分分析剔除了原有数据中的噪声和冗余信息，并不是舍弃原有数据中的特征信息。主成分分析的目的就是最小化噪声，最大化提取有用的信息。

三、主成分分析的核心思想

所谓主成分分析就是将样本的nnn维特征映射到kkk维上，这kkk维全新的正交特征称为数据的kkk个主成分。主成分分析把原有的数据变换到一个新的坐标系中，使得任何数据投影的最大化方差在第一个坐标（第一主成分）上，第二大方差在第二个坐标（第二主成分）上，以此类推，次数为原始数据中特征的数目。对于条目2，第二个新坐标系（第二主成分）是与第一个坐标系（第一主成分）正交的空间中的方差最大化。对于二维数据，第二个坐标系则垂直于第一坐标系。对于条目3，第三个坐标系（第三主成分）则分别与第一个和第二个坐标系正交且方差最大化。这就等价于线性代数中的施密特正交化的过程。在整个过程中通常用方差占比来衡量每个主成分能够解释信息的比例，如果前kkk个主成分几乎涵盖绝大部分的信息而后面的坐标系所含的方差几乎为000，那么就可以保留数据的前kkk个主成分，实现了降维的过程。在降维中的方差最大化就是投影误差最小化，这两个过程是等价的。找到方差最大化的过程就是主成分分析的本质。最大化数据的方差则数据降维后就越分散，能够表征的信息量则越多，特征的区分度就越明显。找到的新维度就是变异（方差）最大的新维度。

如下面两幅图所示：

左右两幅图是同为100010001000个二维数据样本的分布，显然左图中AAA直线代表的方向是覆盖数据的最大方差位置，右图是将该数据集的投影到红线的方向，将会是方差最大化且投影误差最小化的坐标系方向，能够最大程度的涵盖数据的信息，是数据的第一主成分特征，下面的章节中会以代码的形式展示。

四、算法流程

将数据进行“去中心化”。计算数据的协方差矩阵。计算协方差矩阵的特征值和所对应的特征向量。将特征值从大到小排序，保留前kkk个特征值，即前kkk个主成分。将数据转换到上述的kkk个特征向量构建的新空间中。

五、所需的数学公式

要理解主成分分析的整个数学原理则要先了解下述数学公式：

样本均值：

xˉ=1n∑i=1nxi\bar{x}=\frac{1}{n}\sum_{i=1}^nx_ixˉ=n1i=1∑nxi

样本方差：

S2=1n−1∑i=1n(xi−xˉ)2S^2=\frac{1}{n-1}\sum_{i=1}^n(x_i-\bar{x})^2S2=n−11i=1∑n(xi−xˉ)2

样本XXX和样本YYY的协方差：

Cov(X,Y)=E[(X−E[X])(Y−E[Y])]=1n−1∑i=1n(xi−xˉ)(yi−yˉ)\begin{aligned} \text{Cov}(X,Y)&=E[(X-E[X])(Y-E[Y])]\\ &=\frac{1}{n-1}\sum_{i=1}^n(x_i-\bar{x})(y_i-\bar{y}) \end{aligned} Cov(X,Y)=E[(X−E[X])(Y−E[Y])]=n−11i=1∑n(xi−xˉ)(yi−yˉ)

方差是针对一维特征的不同样本的取值的运算，协方差是针对二维以上的特征的不同取值之间的运算。方差是协方差的特殊情况。

上述公式中的方差和协方差的除数是n−1n-1n−1代表的是无偏估计。

矩阵的特征值和特征向量：

A⋅v=λ⋅v\mathbf{A}\cdot{\mathbf{v}}=\lambda\cdot{\mathbf{v}}A⋅v=λ⋅v

其中，vvv是矩阵AAA的特征向量，λ\lambdaλ是矩阵AAA的特征值，每个特征值所对应的特征向量相互正交。

对于矩阵特征值和特征向量公式的理解：一个矩阵点乘一个向量就相当于对该向量进行旋转或者拉伸等一系列线性变换。在我们求解矩阵的特征值和特征向量的时候就是要求解矩阵AAA能够使得这些特征向量只发生伸缩变换，而变换的效果等价于特征向量与某个常量λ\lambdaλ相乘。

矩阵的特征值分解：

A=QΣQ−1A=Q\Sigma Q^{-1}A=QΣQ−1

QQQ是矩阵AAA所对应的特征向量，Σ\SigmaΣ是矩阵AAA的特征值组成的对角矩阵。

六、主成分分析的应用实例和原理

下面利用一个生活中的常见问题来讲解主成分分析的原理。

有A，B，C，D，E五款汽车，其价格如下述表格所示：

放在一维坐标系上的表示如下：

555款车的均值为：

xˉ=10+3+1+7+25=4.6\bar{x}=\frac{10+3+1+7+2}{5}=4.6xˉ=510+3+1+7+2=4.6

为方便后续方差的计算，首先将数据进行“去中心化”，即xi−xˉx_i-\bar{x}xi−xˉ，“去中心化”后的数据显示为：

“去中心化”后的点在坐标系的表示：

显然上述图片可以看做是以000为中心的数据分布，这样做则更具有可观测性。通过上图可以清楚地看到AAA和DDD的售价比较高，CCC，EEE和BBB的售价稍微低一些。

555款车的方差为：

S2=14(5,42+(−1.6)2+(−3.6)2+2.42+(−2.6)2)=14.3S^2=\frac{1}{4}\bigg (5,4^2+(-1.6)^2+(-3.6)^2+2.4^2+(-2.6)^2 \bigg)=14.3S2=41(5,42+(−1.6)2+(−3.6)2+2.42+(−2.6)2)=14.3

现在又添加了这555款车的使用年限属性，如下表所示：

相应的坐标系表示为：

555款汽车的两个特征“去中心化”后的数据显示为：

相应的平面直角坐标系的表示如下：

现在需要从这555款车的两个特征中解析出一个新的特征来实现对数据的降维。

下面穿插一些数学知识点：

高中时候都学习过向量的知识。在二维世界里可以用两个正交基来表示任一点，即任何一个点都是两个正交基e1,e2\mathbf{e_1,e_2}e1,e2的线性组合：

A=(xy)=xe1+ye2A=\binom{x}{y}=x\mathbf{e_1}+y\mathbf{e_2}A=(yx)=xe1+ye2

下图所示的向量e1,e2\mathbf{e_1,e_2}e1,e2是平面直角坐标系的两个标准正交基：

坐标系以OOO为中心进行旋转则xxx和yyy的值也会发生变化。但是，无论坐标系旋转到何种角度，∣OA∣2=x2+y2\mathbf{|OA|}^2=x^2+y^2∣OA∣2=x2+y2是不变的。当xxx轴旋转到与向量OA\mathbf{OA}OA共线时，点AAA则无需e2\mathbf{e_2}e2向量的表示，降成一维的点，如下图所示：

此时，上图AAA点可以表示成：

A=3.82e1+0e2A=3.82\mathbf{e_1}+0\mathbf{e_2}A=3.82e1+0e2

直接去掉向量e1\mathbf{e_1}e1对点AAA并不影响：

上图AAA点可以表示成：

A=3.82e2A=3.82\mathbf{e_2}A=3.82e2

如果是两个点，情况就稍微复杂了点。有如下图所示的两个点A=(x1y1)A=\binom{x_1}{y_1}A=(y1x1)，C=(x2y2)C=\binom{x_2}{y_2}C=(y2x2)：

则点AAA和点CCC的向量表示为：

w=x1e1+y1e2,v=x2e1+y2e2\mathbf{w}=x_1\mathbf{e_1}+y_1\mathbf{e_2},\mathbf{v}=x_2\mathbf{e_1}+y_2\mathbf{e_2}w=x1e1+y1e2,v=x2e1+y2e2

由上图可知，点AAA和点CCC在xxx轴上投影的长度要多一写。为了降维，应该多分配给x1,x2x_1,x_2x1,x2，少分配给y1,y2y_1,y_2y1,y2。

对于上面的两个点A=(X1Y1)A=\binom{X_1}{Y_1}A=(Y1X1)，B=(X2Y2)B=\binom{X_2}{Y_2}B=(Y2X2)，在xxx轴的投影如下图所示：

随着以OOO为原点的坐标轴的不断旋转，X1X_1X1和X2X_2X2的值将不断发生变化，当旋转到如下图所示的情况，X1X_1X1和X2X_2X2的值最大。

尽量多的分配给X1X_1X1和X2X_2X2，则要借鉴最小二乘法的思想：

max⁡(X12+X22)=max⁡∑i=12Xi2\max{(X_1^2+X_2^2)}=\max{\sum_{i=1}^2{X_i^2}}max(X12+X22)=maxi=1∑2Xi2

假设坐标系的正交基向量为e1=(e11e12)\mathbf{e_1}=\binom{e_{11}}{e_{12}}e1=(e12e11)，e2=(e21e22)\mathbf{e_2}=\binom{e_{21}}{e_{22}}e2=(e22e21)。设向量OA\mathbf{OA}OA为a\mathbf{a}a，向量OB\mathbf{OB}OB为b\mathbf{b}b。根据向量点积的公式有：

X1=a⋅e1=(x1y1)⋅(e11e12)=x1e11+y1e12X_1=\mathbf{a}\cdot\mathbf{e_1}=\binom{x_1}{y_1}\cdot\binom{e_{11}}{e_{12}}=x_1e_{11}+y_1e_{12}X1=a⋅e1=(y1x1)⋅(e12e11)=x1e11+y1e12

X2=b⋅e1=(x2y2)⋅(e11e12)=x2e11+y2e12X_2=\mathbf{b}\cdot\mathbf{e_1}=\binom{x_2}{y_2}\cdot\binom{e_{11}}{e_{12}}=x_2e_{11}+y_2e_{12} X2=b⋅e1=(y2x2)⋅(e12e11)=x2e11+y2e12

则有：

X12+X22=(x1e11+y1e12)2+(x2e11+y2e12)2=x12e112+2x1y1e11e12+y12e122+x22e112+2x2y2e11e12+y22e122=(x12+x22)2e112+2(x1y1+x2y2)e11e12+(y12+y22)2e122\begin{aligned} X_1^2+X_2^2 & = (x_1e_{11}+y_1e_{12})^2+(x_2e_{11}+y_2e_{12})^2 \\ & = x_1^2e_{11}^2+2x_1y_1e_{11}e_{12}+y_1^2e_{12}^2+x_2^2e_{11}^2+2x_2y_2e_{11}e_{12}+y_2^2e_{12}^2 \\ &= (x_1^2+x_2^2)^2e_{11}^2+2(x_1y_1+x_2y_2)e_{11}e_{12}+(y_1^2+y2^2)^2e_{12}^2 \end{aligned} X12+X22=(x1e11+y1e12)2+(x2e11+y2e12)2=x12e112+2x1y1e11e12+y12e122+x22e112+2x2y2e11e12+y22e122=(x12+x22)2e112+2(x1y1+x2y2)e11e12+(y12+y22)2e122

上式转变为二次型形式：

X12+X22=e1⊤((x12+x22)2x1y1+x2y2x1y1+x2y2(y12+y22)2)⏟Pe1=e1⊤Pe1X_1^2+X_2^2=\mathbf{e_1^\top}\begin{matrix} \underbrace{ \begin{pmatrix} (x_1^2+x_2^2)^2 & x_1y_1+x_2y_2 \\ x_1y_1+x_2y_2 & (y_1^2+y2^2)^2 \end{pmatrix} } \\ \mathbf{P} \end{matrix} \mathbf{e_1}=\mathbf{e_1^\top} \mathbf{P}\mathbf{e_1} X12+X22=e1⊤((x12+x22)2x1y1+x2y2x1y1+x2y2(y12+y22)2)Pe1=e1⊤Pe1

这里的矩阵P\mathbf{P}P是一个二次型可以进行特征值（奇异值）分解：

P=AΣA⊤\mathbf{P}=\mathbf{A\Sigma{A^\top}}P=AΣA⊤

其中，A\mathbf{A}A为正交矩阵，即A−1=A⊤\mathbf{A^{-1}}=\mathbf{A^\top}A−1=A⊤，∣A∣=1\mathbf{|A|}=1∣A∣=1或者−1-1−1并且AA⊤=AA−1=I\mathbf{AA^{\top}}=\mathbf{AA^{-1}}=\mathbf{I}AA⊤=AA−1=I。

Σ\mathbf{\Sigma}Σ是对角矩阵，

Σ=(σ100σ2)\mathbf{\Sigma}=\begin{pmatrix} \sigma_1 & 0 \\ 0 & \sigma_2 \end{pmatrix} Σ=(σ100σ2)

其中，σ1\sigma_1σ1和σ2\sigma_2σ2是矩阵P\mathbf{P}P的特征值（奇异值），σ1>σ2\sigma_1>\sigma_2σ1>σ2。

将矩阵P\mathbf{P}P回代：

X12+X22=e1⊤Pe2=e1⊤AΣA⊤e1=(A⊤e1)⊤Σ(A⊤e1)\begin{aligned} X_1^2+X_2^2 &=\mathbf{e_1^\top{P}e_2} \\ &= \mathbf{e_1^\top{A}\Sigma{A^\top}e_1} \\ &= \mathbf{(A^\top{e_1})^\top\Sigma{(A^\top{e_1})}} \end{aligned} X12+X22=e1⊤Pe2=e1⊤AΣA⊤e1=(A⊤e1)⊤Σ(A⊤e1)

令n=A⊤e1\mathbf{n}=\mathbf{A^\top{e_1}}n=A⊤e1，由于A\mathbf{A}A是正交矩阵，所得到的n\mathbf{n}n也是单位向量，即：

n=(n1n2),n12+n22=1\mathbf{n}=\binom{n_1}{n_2},n_1^2+n_2^2=1n=(n2n1),n12+n22=1

则有

X12+X22=(A⊤e1)⊤Σ(A⊤e1)=n⊤Σn=(n1n2)(σ100σ2)(n1n2)=(n1σ1n2σ2)(n1n2)=σ1n12+σ2n22\begin{aligned} X_1^2+X_2^2 &= \mathbf{(A^\top{e_1})^\top\Sigma{(A^\top{e_1})}} \\ &=\mathbf{n^\top{\Sigma}n} \\ &= \begin{pmatrix}n_1 & n_2 \end{pmatrix} \begin{pmatrix} \sigma_1 & 0\\ 0 &\sigma_2 \end{pmatrix} \begin{pmatrix} n_1 \\ n_2 \end{pmatrix}\\ &=\begin{pmatrix} n_1\sigma_1 & n_2\sigma_2 \end{pmatrix} \begin{pmatrix} n_1 \\ n_2 \end{pmatrix}\\ &=\sigma_1n_1^2+\sigma_2n_2^2 \end{aligned} X12+X22=(A⊤e1)⊤Σ(A⊤e1)=n⊤Σn=(n1n2)(σ100σ2)(n1n2)=(n1σ1n2σ2)(n1n2)=σ1n12+σ2n22

此时，求X12+X22X_1^2+X_2^2X12+X22的最大值问题就转化成了如下的最优化问题。

max⁡(X12+X22)⟺{max⁡(σ1n12+σ2n22)n12+n22=1σ1>σ2\max(X_1^2+X_2^2)\Longleftrightarrow \left \{\begin{matrix} \max(\sigma_1n_1^2+\sigma_2n_2^2) \\ n_1^2+n_2^2=1 \\ \sigma_1>\sigma_2\\ \end{matrix} \right. max(X12+X22)⟺⎩⎨⎧max(σ1n12+σ2n22)n12+n22=1σ1>σ2

上述问题可以使用拉格朗日乘数法来解答。其实当σ1>σ2\sigma_1>\sigma_2σ1>σ2时，显然n1n_1n1取111，n2n_2n2 取000时，X12+X22X_1^2+X_2^2X12+X22取最大值。

由此计算出的X12+X22X_1^2+X_2^2X12+X22的最大值为σ1\sigma_1σ1即为协方差矩阵的最大特征值，其所对应的特征向量的方向就是原有数据的第一主成分方向。

A⊤e1=(10)⇒e1=A(10)\mathbf{A^\top e_1}=\binom{1}{0}\Rightarrow \mathbf{e_1}=\mathbf{A}\binom{1}{0}A⊤e1=(01)⇒e1=A(01)

上式表明的含义就是正交矩阵A\mathbf{A}A中σ1\sigma_1σ1所对应的特征向量方向就是数据的第一主成分方向。

同理则有：

A⊤e2=(01)⇒e2=A(01)\mathbf{A^\top e_2}=\binom{0}{1}\Rightarrow \mathbf{e_2}=\mathbf{A}\binom{0}{1}A⊤e2=(10)⇒e2=A(10)

向量e2\mathbf{e_2}e2是正交矩阵A\mathbf{A}A中特征值σ2\sigma_2σ2所对应的特征向量方向，该方向就是数据的第二主成分方向。

以上便是主成分分析原理的全部推导过程。

现在回到上面的例子：

555款汽车的两个特征“去中心化”后的数据显示为：

读取到的两个向量是：

X=(5.4−1.6−3.62.4−2.6),Y=(6.4−0.6−5.62.4−2.6)\mathbf{X}=\begin{pmatrix} 5.4\\ -1.6\\ -3.6\\ 2.4\\ -2.6 \end{pmatrix}, \mathbf{Y}=\begin{pmatrix} 6.4\\ -0.6\\ -5.6\\ 2.4\\ -2.6 \end{pmatrix} X=⎝⎜⎜⎜⎜⎛5.4−1.6−3.62.4−2.6⎠⎟⎟⎟⎟⎞,Y=⎝⎜⎜⎜⎜⎛6.4−0.6−5.62.4−2.6⎠⎟⎟⎟⎟⎞

协方差矩阵为：

Cov(X,Y)=14(X⋅XX⋅YY⋅XY⋅Y)=(14.317.0517.0521.3)\mathbf{\text{Cov}(X,Y)}=\frac{1}{4} \begin{pmatrix} \mathbf{X\cdot{X}} & \mathbf{X\cdot{Y}}\\ \mathbf{Y\cdot{X}} &\mathbf{Y\cdot{Y}} \end{pmatrix}= \begin{pmatrix} 14.3 & 17.05\\ 17.05 & 21.3 \end{pmatrix} Cov(X,Y)=41(X⋅XY⋅XX⋅YY⋅Y)=(14.317.0517.0521.3)

对协方差矩阵进行特征值分解：

Cov(X,Y)=(−0.63−0.77−0.770.63)(35.21000.39)(−0.63−0.77−0.770.63)\mathbf{\text{Cov}(X,Y)}= \begin{pmatrix} -0.63 & -0.77\\ -0.77 & 0.63 \end{pmatrix} \begin{pmatrix} 35.21 & 0\\ 0 & 0.39 \end{pmatrix} \begin{pmatrix} -0.63 & -0.77\\ -0.77 & 0.63 \end{pmatrix} Cov(X,Y)=(−0.63−0.77−0.770.63)(35.21000.39)(−0.63−0.77−0.770.63)

两个向量方向所能涵盖的信息量正如协方差矩阵的特征值35.2135.2135.21和0.390.390.39的比重：

如上图所示，第一主成分所能涵盖原始数据99%99\%99%的信息量，第二主成分所涵盖数据的1%1\%1%的信息量。

上述协方差特征值分解的公式表明，第一主成分对应的特征向量e1\mathbf{e_1}e1和第二主成分对应的特征向量e2\mathbf{e_2}e2分别是：

e1=(−0.63−0.77),e2=(−0.770.63)\mathbf{e_1}=\binom{-0.63}{-0.77},\mathbf{e_2}=\binom{-0.77}{0.63}e1=(−0.77−0.63),e2=(0.63−0.77)

以这两个正交基的方向为坐标系画出来的图如下所示：

经过运算，将这555个点投影到单位向量e1\mathbf{e_1}e1所在的方向后的坐标表示为：

如上图所示，将二维空间的点投影到一维空间上去并且这个方向就是原有样本的方差最大化方向，也是投影误差最小化方向，即第一主成分方向。需要注意的是，上图投影到e1\mathbf{e_1}e1上的点依旧是二维坐标的表示形式，点的坐标的二维表示如下：

X=(5.33−0.93−4.212.16−2.33),Y=(6.46−1.13−5.112.61−2.82)X= \begin{pmatrix} 5.33 \\ -0.93 \\ -4.21 \\ 2.16 \\ -2.33 \end{pmatrix}, Y= \begin{pmatrix} 6.46 \\ -1.13 \\ -5.11 \\ 2.61 \\ -2.82 \end{pmatrix} X=⎝⎜⎜⎜⎜⎛5.33−0.93−4.212.16−2.33⎠⎟⎟⎟⎟⎞,Y=⎝⎜⎜⎜⎜⎛6.46−1.13−5.112.61−2.82⎠⎟⎟⎟⎟⎞

仍需注意的是，上述公式中的X,YX,YX,Y是“去中心化”后的数据，要得到原始数据的第一主成分则需要加上这两个维度的均值，即X+4.6X+4.6X+4.6，Y+9.6Y+9.6Y+9.6。所以原始数据的第一主成分数值的二维空间表示为：

X=(9.933.960.396.762.27),Y=(16.068.474.4912.216.78)X= \begin{pmatrix} 9.93 \\ 3.96 \\ 0.39 \\ 6.76 \\ 2.27 \end{pmatrix}, Y= \begin{pmatrix} 16.06 \\ 8.47 \\ 4.49 \\ 12.21 \\ 6.78 \end{pmatrix} X=⎝⎜⎜⎜⎜⎛9.933.960.396.762.27⎠⎟⎟⎟⎟⎞,Y=⎝⎜⎜⎜⎜⎛16.068.474.4912.216.78⎠⎟⎟⎟⎟⎞

下面是将二维数据投影到一维数据的坐标表示，以e1\mathbf{e_1}e1为投影方向将原始222维数据降维到111维后的点的坐标是：

(XandY)⋅e1⇒(5.46.4−1.6−0.6−3.6−5.62.42.4−2.6−2.6)⋅(−0.63−0.77)=(−8.371.486.61−3.383.66)\mathbf{(X\text{and}Y)\cdot e_1}\Rightarrow \begin{pmatrix} 5.4 &6.4 \\ -1.6 &-0.6 \\ -3.6 & -5.6\\ 2.4 & 2.4\\ -2.6 & -2.6 \end{pmatrix}\cdot \begin{pmatrix} -0.63 \\ -0.77 \end{pmatrix}=\begin{pmatrix} -8.37 \\ 1.48 \\ 6.61 \\ -3.38 \\ 3.66 \end{pmatrix} (XandY)⋅e1⇒⎝⎜⎜⎜⎜⎛5.4−1.6−3.62.4−2.66.4−0.6−5.62.4−2.6⎠⎟⎟⎟⎟⎞⋅(−0.63−0.77)=⎝⎜⎜⎜⎜⎛−8.371.486.61−3.383.66⎠⎟⎟⎟⎟⎞

如下图所示：

同理，第二主成分的投影坐标系表示为：

降维后的第二主成分在二维空间的点的坐标表示为（“去中心化”数据）：

X=(0.1−0.660.590.26−0.28),Y=(−0.880.54−0.48−0.220.23)X= \begin{pmatrix} 0.1 \\ -0.66 \\ 0.59 \\ 0.26 \\ -0.28 \end{pmatrix}, Y= \begin{pmatrix} -0.88 \\ 0.54 \\ -0.48 \\ -0.22 \\ 0.23 \end{pmatrix} X=⎝⎜⎜⎜⎜⎛0.1−0.660.590.26−0.28⎠⎟⎟⎟⎟⎞,Y=⎝⎜⎜⎜⎜⎛−0.880.54−0.48−0.220.23⎠⎟⎟⎟⎟⎞

加上均值后的原始数据在第二主成分的投影在二维坐标系的表示为：

X=(4.713.935.184.874.31),Y=(9.5110.149.129.389.83)X= \begin{pmatrix} 4.71 \\ 3.93 \\ 5.18 \\ 4.87 \\ 4.31 \end{pmatrix}, Y= \begin{pmatrix} 9.51 \\ 10.14 \\ 9.12 \\ 9.38 \\ 9.83 \end{pmatrix} X=⎝⎜⎜⎜⎜⎛4.713.935.184.874.31⎠⎟⎟⎟⎟⎞,Y=⎝⎜⎜⎜⎜⎛9.5110.149.129.389.83⎠⎟⎟⎟⎟⎞

同理，投影到e2\mathbf{e_2}e2方向后降维到111维的坐标表示为：

(XandY)⋅e2=(5.46.4−1.6−0.6−3.6−5.62.42.4−2.6−2.6)⋅(−0.770.63)=(−0.130.85−0.76−0.340.36)\mathbf{(X\text{and}Y)\cdot{e_2}}= \begin{pmatrix} 5.4 &6.4 \\ -1.6 &-0.6 \\ -3.6 & -5.6\\ 2.4 & 2.4\\ -2.6 & -2.6 \end{pmatrix}\cdot \begin{pmatrix} -0.77 \\ 0.63 \end{pmatrix}=\begin{pmatrix} -0.13 \\ 0.85 \\ -0.76 \\ -0.34 \\ 0.36 \end{pmatrix} (XandY)⋅e2=⎝⎜⎜⎜⎜⎛5.4−1.6−3.62.4−2.66.4−0.6−5.62.4−2.6⎠⎟⎟⎟⎟⎞⋅(−0.770.63)=⎝⎜⎜⎜⎜⎛−0.130.85−0.76−0.340.36⎠⎟⎟⎟⎟⎞

由上式和上图可见，降维后的第二主成分的元素值明显小于第一主成分的元素值，自然涵盖的信息量就少，直接去掉并不会对数据的分析结果造成太大的影响。

七、上述例子的Python代码实现

主成分分析的代码实现可以分为两种：自定义pca算法和调用sklearn工具包。

上述例子中的运算结果是利用博主本人的手算和借助geogebra计算器完成，下面再看看利用代码实现的结果如何。

import numpy as npimport matplotlib.pyplot as plt# 实验数据Data = np.array([[10, 16], [3, 9], [1, 4], [7, 12], [2, 7]])print(Data)'''[[10 16][ 3 9][ 1 4][ 7 12][ 2 7]]'''# 计算均值meanVals = np.mean(Data, axis=0)print(meanVals) # [4.6 9.6]# 数据“去中心化”DataMeanRemoved = Data - meanValsprint(DataMeanRemoved)'''[[ 5.4 6.4][-1.6 -0.6][-3.6 -5.6][ 2.4 2.4][-2.6 -2.6]]'''# 计算数据的协方差矩阵CovData = np.cov(DataMeanRemoved, rowvar=False)'''rowvar：布尔值，可选如果rowvar为True（默认值），则每行代表一个变量X，另一个行为变量Y。否则，转换关系：每列代表一个变量X，另一个列为变量Y。'''CovData = np.mat(CovData) # 必须转化为矩阵形式才能做后续的运算print("协方差矩阵是：\n", CovData)'''[[14.3 17.05][17.05 21.3 ]]'''# 计算特征值和特征向量eigVals, eigVects = np.linalg.eig(CovData)# 必须是矩阵形式print("特征值是：\n", eigVals)print("特征向量是：\n", eigVects)'''特征值是：[ 0.39446927 35.20553073]特征向量是：[[-0.77494694 -0.6320263 ][ 0.6320263 -0.77494694]]'''maxEigVect = eigVects[:, 1] # 打印最大的特征值所对应的特征向量print("第一主成分所对应的特征向量为：\n", maxEigVect)'''第一主成分所对应的特征向量为：[[-0.6320263 ][-0.77494694]]'''print("去中心化后的数据是：\n", DataMeanRemoved)'''去中心化后的数据是：[[ 5.4 6.4][-1.6 -0.6][-3.6 -5.6][ 2.4 2.4][-2.6 -2.6]]'''# 第一主成分向量与去中心化后的数据进行点积得到降维后的数据LowDData = DataMeanRemoved * maxEigVectprint("降维后的数值为：\n", LowDData)'''降维后的数值为：[[-8.37260242][ 1.47621024][ 6.61499753][-3.37673577][ 3.65813042]]'''# 原始数据reconMat = (LowDData * maxEigVect.T) + meanValsprint("原始数据降维后的二维坐标表示为：\n", reconMat)'''原始数据降维后的二维坐标表示为：[[ 9.89170494 16.0883226 ][ 3.6669963 8.45601539][ 0.41914758 4.47372793][ 6.73418582 12.21679104][ 2.28796536 6.76514304]]'''

下面的结果就是使用python中的numpy模块实现。对照第六部分中的笔算结果，由此可见，两种结果相差无几。

[[-8.37260242][ 1.47621024][ 6.61499753][-3.37673577][ 3.65813042]]

下述代码是直接使用sklearn中的PCA模块实现主成分分析数据降维，运行结果无异。

from sklearn.decomposition import PCAimport numpy as npX = np.array([[10, 16], [3, 9], [1, 4], [7, 12], [2, 7]])pca = PCA(n_components = 1)pca.fit(X)print(pca.transform(X))'''[[ 8.37260242][-1.47621024][-6.61499753][ 3.37673577][-3.65813042]]'''

八、利用主成分分析对`MNIST`数据集进行降维

MNIST数据集是0-9这10个数字的手写形式的数据集，涵盖6万张数据集和1万张测试集，每张图片是28 * 28维度的图片。具体可以查看博主本人的另外一篇博客：MNIST数据集简介与使用。

使用PCA对MNIST数据集进行降维，首先将数据集中的标签“0”提取100张，并使用图片拼接技术将这100转给你手写数字0拼成一张10×1010\times 1010×10的图片，然后是降维前后的图片对比。

代码如下所述：

import tensorflow.examples.tutorials.mnist.input_data as input_dataimport numpy as npimport matplotlib.pyplot as pltfrom PIL import Imagefrom sklearn.decomposition import PCAfrom sklearn.svm import LinearSVCfrom sklearn.metrics import classification_reportmnist = input_data.read_data_sets("MNIST_data/", one_hot=False) #下载MNIST数据集imgs = mnist.train.imageslabels = mnist.train.labelsorign_imgs = []for i in range(1500):if labels[i] == 0:orign_imgs .append(imgs[i])if len(orign_imgs) == 100:breakdef array_to_img(array):array = array * 255new_img = Image.fromarray(array.astype(np.uint8))return new_imgdef comb_imgs(origin_imgs, col, row, each_width, each_height, new_type):new_img = Image.new(new_type, (col * each_width, row * each_height))for i in range(len(origin_imgs)):each_img = array_to_img(np.array(origin_imgs[i]).reshape(each_width, each_width))new_img.paste(each_img, ((i % col) * each_width, (i // col) * each_width))return new_imgdef pca(data_mat, top_n_feat=1):mean_vals = data_mat.mean(axis=0)mean_removed = data_mat - mean_vals cov_mat = np.cov(mean_removed, rowvar=False)eig_vals, eig_vects = np.linalg.eig(np.mat(cov_mat)) eig_val_index = np.argsort(eig_vals) eig_val_index = eig_val_index[:-(top_n_feat + 1): -1] reg_eig_vects = eig_vects[:, eig_val_index]low_d_data_mat = mean_removed * reg_eig_vects recon_mat = (low_d_data_mat * reg_eig_vects.T) + mean_vals return low_d_data_mat, recon_matif __name__ == "__main__":ten_origin_imgs = comb_imgs(orign_imgs, 10, 10, 28, 28, 'L')ten_origin_imgs.show()low_d_feat_for_imgs, recon_mat_for_imgs = pca(np.array(orign_imgs), 1)low_d_img = comb_imgs(recon_mat_for_imgs, 10, 10, 28, 28, 'L')low_d_img.show()

降维之前：

降维之后：

上面两幅图是降维之前和降维之后的手写数字000。利用主成分分析只取了第111个特征。显然降维之后的图片去除了原图中大量的噪声，效果表现更为规范清晰。