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

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

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

一、引言

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

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

主成分分析（Principal Component Analysis，简称PCA）是一种通过降维技术把多变量化为少数几个主成分的统计分析方法。

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

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

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

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

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

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

  如下面两幅图所示：

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

四、算法流程

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

五、所需的数学公式

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

样本均值：

样本方差：

样本和样本的协方差：

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

代表的是无偏估计。

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

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

能够使得这些特征向量只发生伸缩变换，而变换的效果等价于特征向量与某个常量相乘。

矩阵的特征值分解：

是矩阵所对应的特征向量，是矩阵的特征值组成的对角矩阵。

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

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

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

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

款车的均值为：

为方便后续方差的计算，首先将数据进行“去中心化”，即，“去中心化”后的数据显示为：

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

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

款车的方差为：

款车的使用年限属性，如下表所示：

  相应的坐标系表示为：

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

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

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

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

的线性组合：

是平面直角坐标系的两个标准正交基：

为中心进行旋转则和的值也会发生变化。但是，无论坐标系旋转到何种角度，是不变的。当轴旋转到与向量共线时，点则无需向量的表示，降成一维的点，如下图所示：

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

直接去掉向量对点并不影响：

上图点可以表示成：

，：

则点和点的向量表示为：

和点在轴上投影的长度要多一写。为了降维，应该多分配给，少分配给。

，，在轴的投影如下图所示：

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

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

，。设向量为，向量为。根据向量点积的公式有：

  则有：

  上式转变为二次型形式：

是一个二次型可以进行特征值（奇异值）分解：

其中，为正交矩阵，即，或者并且。

是对角矩阵，

其中，和是矩阵的特征值（奇异值），。

回代：

，由于是正交矩阵，所得到的也是单位向量，即：

则有

的最大值问题就转化成了如下的最优化问题。

上述问题可以使用拉格朗日乘数法来解答。其实当时，显然取， 取时，取最大值。

的最大值为即为协方差矩阵的最大特征值，其所对应的特征向量的方向就是原有数据的第一主成分方向。

中所对应的特征向量方向就是数据的第一主成分方向。

  同理则有：

是正交矩阵中特征值所对应的特征向量方向，该方向就是数据的第二主成分方向。

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

现在回到上面的例子：

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

  读取到的两个向量是：

  协方差矩阵为：

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

和的比重：

的信息量，第二主成分所涵盖数据的的信息量。

和第二主成分对应的特征向量分别是：

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

个点投影到单位向量所在的方向后的坐标表示为：

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

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

将二维数据投影到一维数据的坐标表示，以为投影方向将原始维数据降维到维后的点的坐标是：

  如下图所示：

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

“去中心化”数据）：

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

投影到方向后降维到维的坐标表示为：

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

七、上述例子的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拼成一张的图片，然后是降维前后的图片对比。

  代码如下所述：

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()

  降维之前：

  降维之后：

。利用主成分分析只取了第个特征。显然降维之后的图片去除了原图中大量的噪声，效果表现更为规范清晰。

九、总结

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