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

一、引言

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

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

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

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

---

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

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

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

---

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

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

  如下面两幅图所示：

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

---

四、算法流程

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

---

五、所需的数学公式

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

样本均值：
x ˉ = 1 n ∑ i = 1 n x i \bar{x}=\frac{1}{n}\sum_{i=1}^nx_i xˉ=n1​i=1∑n​xi​

样本方差：

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

样本 X X X和样本 Y Y Y的协方差：
Cov ( X , Y ) = E [ ( X − E [ X ] ) ( Y − E [ Y ] ) ] = 1 n − 1 ∑ i = 1 n ( x i − x ˉ ) ( y i − 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−11​i=1∑n​(xi​−xˉ)(yi​−yˉ​)​

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

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

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

A ⋅ v = λ ⋅ v \mathbf{A}\cdot{\mathbf{v}}=\lambda\cdot{\mathbf{v}} A⋅v=λ⋅v
  其中， v v v是矩阵 A A A的特征向量， λ \lambda λ是矩阵 A A A的特征值，每个特征值所对应的特征向量相互正交。

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

矩阵的特征值分解：

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

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

---

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

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

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

车型	价格（十万）
A	10
B	3
C	1
D	7
E	2

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

   5 5 5款车的均值为：

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

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

车型	价格（十万）
A	5.4
B	-1.6
C	-3.6
D	2.4
E	-2.6

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

显然上述图片可以看做是以 0 0 0为中心的数据分布，这样做则更具有可观测性。通过上图可以清楚地看到 A A A和 D D D的售价比较高， C C C， E E E和 B B B的售价稍微低一些。

   5 5 5款车的方差为：

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

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

车型	价格（十万）	使用年限（年）
A	10	16
B	3	9
C	1	4
D	7	12
E	2	7

  相应的坐标系表示为：

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

车型	价格（十万）	使用年限（年）
A	5.4	6.4
B	-1.6	-0.6
C	-3.6	-5.6
D	2.4	2.4
E	-2.6	-2.6

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

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

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

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

A = ( x y ) = x e 1 + y e 2 A=\binom{x}{y}=x\mathbf{e_1}+y\mathbf{e_2} A=(yx​)=xe1​+ye2​

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

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

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

A = 3.82 e 1 + 0 e 2 A=3.82\mathbf{e_1}+0\mathbf{e_2} A=3.82e1​+0e2​

直接去掉向量 e 1 \mathbf{e_1} e1​对点 A A A并不影响：

上图 A A A点可以表示成：

A = 3.82 e 1 A=3.82\mathbf{e_1} A=3.82e1​

  如果是两个点，情况就稍微复杂了点。有如下图所示的两个点 A = ( x 1 y 1 ) A=\binom{x_1}{y_1} A=(y1​x1​​)， C = ( x 2 y 2 ) C=\binom{x_2}{y_2} C=(y2​x2​​)：

则点 A A A和点 C C C的向量表示为：

w = x 1 e 1 + y 1 e 2 , v = x 2 e 1 + y 2 e 2 \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=x1​e1​+y1​e2​,v=x2​e1​+y2​e2​

  由上图可知，点 A A A和点 C C C在 x x x轴上投影的长度要多一写。为了降维，应该多分配给 x 1 , x 2 x_1,x_2 x1​,x2​，少分配给 y 1 , y 2 y_1,y_2 y1​,y2​。

  对于上面的两个点 A = ( X 1 Y 1 ) A=\binom{X_1}{Y_1} A=(Y1​X1​​)， B = ( X 2 Y 2 ) B=\binom{X_2}{Y_2} B=(Y2​X2​​)，在 x x x轴的投影如下图所示：

  随着以 O O O为原点的坐标轴的不断旋转， X 1 X_1 X1​和 X 2 X_2 X2​的值将不断发生变化，当旋转到如下图所示的情况， X 1 X_1 X1​和 X 2 X_2 X2​的值最大。

  尽量多的分配给 X 1 X_1 X1​和 X 2 X_2 X2​，则要借鉴最小二乘法的思想：

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

  假设坐标系的正交基向量为 e 1 = ( e 11 e 12 ) \mathbf{e_1}=\binom{e_{11}}{e_{12}} e1​=(e12​e11​​)， e 2 = ( e 21 e 22 ) \mathbf{e_2}=\binom{e_{21}}{e_{22}} e2​=(e22​e21​​)。设向量 O A \mathbf{OA} OA为 a \mathbf{a} a，向量 O B \mathbf{OB} OB为 b \mathbf{b} b。根据向量点积的公式有：

X 1 = a ⋅ e 1 = ( x 1 y 1 ) ⋅ ( e 11 e 12 ) = x 1 e 11 + y 1 e 12 X_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​=(y1​x1​​)⋅(e12​e11​​)=x1​e11​+y1​e12​

X 2 = b ⋅ e 1 = ( x 2 y 2 ) ⋅ ( e 11 e 12 ) = x 2 e 11 + y 2 e 12 X_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​=(y2​x2​​)⋅(e12​e11​​)=x2​e11​+y2​e12​

  则有：

X 1 2 + X 2 2 = ( x 1 e 11 + y 1 e 12 ) 2 + ( x 2 e 11 + y 2 e 12 ) 2 = x 1 2 e 11 2 + 2 x 1 y 1 e 11 e 12 + y 1 2 e 12 2 + x 2 2 e 11 2 + 2 x 2 y 2 e 11 e 12 + y 2 2 e 12 2 = ( x 1 2 + x 2 2 ) 2 e 11 2 + 2 ( x 1 y 1 + x 2 y 2 ) e 11 e 12 + ( y 1 2 + y 2 2 ) 2 e 12 2 \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​​=(x1​e11​+y1​e12​)2+(x2​e11​+y2​e12​)2=x12​e112​+2x1​y1​e11​e12​+y12​e122​+x22​e112​+2x2​y2​e11​e12​+y22​e122​=(x12​+x22​)2e112​+2(x1​y1​+x2​y2​)e11​e12​+(y12​+y22)2e122​​

  上式转变为二次型形式：

X 1 2 + X 2 2 = e 1 ⊤ ( ( x 1 2 + x 2 2 ) 2 x 1 y 1 + x 2 y 2 x 1 y 1 + x 2 y 2 ( y 1 2 + y 2 2 ) 2 ) ⏟ P e 1 = e 1 ⊤ P e 1 X_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​)2x1​y1​+x2​y2​​x1​y1​+x2​y2​(y12​+y22)2​)​P​e1​=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并且 A A ⊤ = A A − 1 = I \mathbf{AA^{\top}}=\mathbf{AA^{-1}}=\mathbf{I} AA⊤=AA−1=I。

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

Σ = ( σ 1 0 0 σ 2 ) \mathbf{\Sigma}=\begin{pmatrix} \sigma_1 & 0 \\ 0 & \sigma_2 \end{pmatrix} Σ=(σ1​0​0σ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回代：

X 1 2 + X 2 2 = e 1 ⊤ P e 2 = e 1 ⊤ A Σ A ⊤ e 1 = ( A ⊤ e 1 ) ⊤ Σ ( A ⊤ e 1 ) \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 ⊤ e 1 \mathbf{n}=\mathbf{A^\top{e_1}} n=A⊤e1​，由于 A \mathbf{A} A是正交矩阵，所得到的 n \mathbf{n} n也是单位向量，即：

n = ( n 1 n 2 ) , n 1 2 + n 2 2 = 1 \mathbf{n}=\binom{n_1}{n_2},n_1^2+n_2^2=1 n=(n2​n1​​),n12​+n22​=1

则有

X 1 2 + X 2 2 = ( A ⊤ e 1 ) ⊤ Σ ( A ⊤ e 1 ) = n ⊤ Σ n = ( n 1 n 2 ) ( σ 1 0 0 σ 2 ) ( n 1 n 2 ) = ( n 1 σ 1 n 2 σ 2 ) ( n 1 n 2 ) = σ 1 n 1 2 + σ 2 n 2 2 \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=(n1​​n2​​)(σ1​0​0σ2​​)(n1​n2​​)=(n1​σ1​​n2​σ2​​)(n1​n2​​)=σ1​n12​+σ2​n22​​

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

max ⁡ ( X 1 2 + X 2 2 ) ⟺ { max ⁡ ( σ 1 n 1 2 + σ 2 n 2 2 ) n 1 2 + n 2 2 = 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(σ1​n12​+σ2​n22​)n12​+n22​=1σ1​>σ2​​

上述问题可以使用拉格朗日乘数法来解答。其实当 σ 1 > σ 2 \sigma_1>\sigma_2 σ1​>σ2​时，显然 n 1 n_1 n1​取 1 1 1， n 2 n_2 n2​ 取 0 0 0时， X 1 2 + X 2 2 X_1^2+X_2^2 X12​+X22​取最大值。

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

A ⊤ e 1 = ( 1 0 ) ⇒ e 1 = A ( 1 0 ) \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 ⊤ e 2 = ( 0 1 ) ⇒ e 2 = A ( 0 1 ) \mathbf{A^\top e_2}=\binom{0}{1}\Rightarrow \mathbf{e_2}=\mathbf{A}\binom{0}{1} A⊤e2​=(10​)⇒e2​=A(10​)

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

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

  现在回到上面的例子：

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

车型	价格（十万）	使用年限（年）
A	5.4	6.4
B	-1.6	-0.6
C	-3.6	-5.6
D	2.4	2.4
E	-2.6	-2.6

  读取到的两个向量是：

X = ( 5.4 − 1.6 − 3.6 2.4 − 2.6 ) , Y = ( 6.4 − 0.6 − 5.6 2.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 ) = 1 4 ( X ⋅ X X ⋅ Y Y ⋅ X Y ⋅ Y ) = ( 14.3 17.05 17.05 21.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⋅X​X⋅YY⋅Y​)=(14.317.05​17.0521.3​)

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

Cov ( X , Y ) = ( − 0.63 − 0.77 − 0.77 0.63 ) ( 35.21 0 0 0.39 ) ( − 0.63 − 0.77 − 0.77 0.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.210​00.39​)(−0.63−0.77​−0.770.63​)

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

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

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

e 1 = ( − 0.63 − 0.77 ) , e 2 = ( − 0.77 0.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​)

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

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

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

X = ( 5.33 − 0.93 − 4.21 2.16 − 2.33 ) , Y = ( 6.46 − 1.13 − 5.11 2.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 , Y X,Y X,Y是“去中心化”后的数据，要得到原始数据的第一主成分则需要加上这两个维度的均值，即 X + 4.6 X+4.6 X+4.6， Y + 9.6 Y+9.6 Y+9.6。所以原始数据的第一主成分数值的二维空间表示为：

X = ( 9.93 3.96 0.39 6.76 2.27 ) , Y = ( 16.06 8.47 4.49 12.21 6.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​⎠⎟⎟⎟⎟⎞​

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

( X and Y ) ⋅ e 1 ⇒ ( 5.4 6.4 − 1.6 − 0.6 − 3.6 − 5.6 2.4 2.4 − 2.6 − 2.6 ) ⋅ ( − 0.63 − 0.77 ) = ( − 8.37 1.48 6.61 − 3.38 3.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.6​6.4−0.6−5.62.4−2.6​⎠⎟⎟⎟⎟⎞​⋅(−0.63−0.77​)=⎝⎜⎜⎜⎜⎛​−8.371.486.61−3.383.66​⎠⎟⎟⎟⎟⎞​

  如下图所示：

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

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

X = ( 0.1 − 0.66 0.59 0.26 − 0.28 ) , Y = ( − 0.88 0.54 − 0.48 − 0.22 0.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.71 3.93 5.18 4.87 4.31 ) , Y = ( 9.51 10.14 9.12 9.38 9.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​⎠⎟⎟⎟⎟⎞​

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

( X and Y ) ⋅ e 2 = ( 5.4 6.4 − 1.6 − 0.6 − 3.6 − 5.6 2.4 2.4 − 2.6 − 2.6 ) ⋅ ( − 0.77 0.63 ) = ( − 0.13 0.85 − 0.76 − 0.34 0.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.6​6.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 np
import 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 - meanVals
print(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 * maxEigVect
print("降维后的数值为：\n", LowDData)
'''
降维后的数值为：
 [[-8.37260242]
 [ 1.47621024]
 [ 6.61499753]
 [-3.37673577]
 [ 3.65813042]]
'''
# 原始数据
reconMat = (LowDData * maxEigVect.T) + meanVals
print("原始数据降维后的二维坐标表示为：\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 PCA
import numpy as np
X = 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 × 10 10\times 10 10×10的图片，然后是降维前后的图片对比。

  代码如下所述：

import tensorflow.examples.tutorials.mnist.input_data as input_data
import numpy as np
import matplotlib.pyplot as plt
from PIL import Image
from sklearn.decomposition import PCA
from sklearn.svm import LinearSVC
from sklearn.metrics import classification_report

mnist = input_data.read_data_sets("MNIST_data/", one_hot=False) #下载MNIST数据集
imgs = mnist.train.images
labels = mnist.train.labels

orign_imgs = []
for i in range(1500):
    if labels[i] == 0:
        orign_imgs .append(imgs[i])
    if len(orign_imgs) == 100:
        break

def array_to_img(array):
    array = array * 255
    new_img = Image.fromarray(array.astype(np.uint8))
    return new_img
   

def 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_img
  

def 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_mat
 
if __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()  

  降维之前：

  降维之后：

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

---

八、总结

  综上所述，主成分分析的确能够摒弃数据样本中的噪声并以一个较大的概率保留数据的有利信息，但缺点有可能是在降维的过程中损失掉有用的信息。