主成分分析降维（MNIST数据集）
北京 | 高性能计算之GPU CUDA课程11月24-26日
刘凯欣，中国矿业大学在校学生，曾参加过ThoughtWorks举办的结对编程活动。
今天看了用主成分分析简化数据，就顺便用MNIST数据集做了下实验，想直观地看一下效果，并通过完成这个小demo深入理解下原理。
我发现“是什么、能做什么、怎么用、效果是什么、原理是什么、优缺点是什么”这样的思路能让我更好地接受一个新知识，之所以把原理放在效果后面，是因为我比较喜欢先看看它的作用，可视化意义之后能提起我对一个知识的兴趣，加深对它意义的理解，后面看数学原理会容易，所以整篇文章就以这样的思路组织整理。
主成分分析（Principal Component Analysis，PCA），一种降维方法，在PCA中，数据从原来的坐标系转换到了新的坐标系，新坐标系由数据本身决定，在新坐标系中，第一个坐标轴选择的是原始数据中方差最大的方向，第二个坐标轴选择的是和第一个坐标轴正交且具有最大方差的方向。该过程一直重复，重复次数为原始数据中特征的数目。我们会发现，大部分方差都包含在最前面的几个新坐标轴中。因此，我们可以忽略余下的坐标轴，即对数据进行了降维处理。
初看这段话感觉是抽象的。方差大意味着什么？方差是衡量源数据和期望值相差的度量值，方差越大，数据差别越大。选择方差最大的方向，就是选择数据差别最大的方向。重复特征数目次，就是说找第一个特征（第一维）方差最大的方向（即覆盖数据点最多的一条直线），做第一个轴，正交且最大方差方向做第二个轴，在此基础上再看第二个特征（第二维），找方差最大方向做第一个轴，正交且最大方差方向做第二个轴，依次类推。这样执行后会发现前几个坐标轴已经差不多囊括所有大差异了，剩下的就不要了，所以实现了降维。
上面从理论上讲了主成分分析和它是如何一步一步实现降维的，有一个感性认识。
降维，在多个指标中只取重要的几个指标，能使复杂问题简单化，就像说话说重点一样。
要做的事就是使用tensorflow里的MNIST数据集，取前100张图片中所有的手写数字7图片，对他们进行主成分分析，输出经过降维反变换回去的图片，对比差异，看看降维后的效果。
引入MNIST数据集、numpy和PIL的Image
import tensorflow.examples.tutorials.mnist.input_data as input_data import numpy as np from PIL import Image
获得MNIST数据集的所有图片和标签
mnist = input_data.read_data_sets("MNIST_data/", one_hot=False)
imgs = mnist.train.images
labels = mnist.train.labels
这里可以看看imgs和labels的type和shape，对于一个python初学者来说总是想搞清楚各个变量的类型和长相。
print(type(imgs)) # <type 'numpy.ndarray'>
print(type(labels)) # <type 'numpy.ndarray'>
print(imgs.shape) # (55000, 784)
print(labels.shape) # (55000,)
取前1000张图片里的100个数字7
origin_7_imgs = [] for i in range(1000):
if labels[i] == 7 and len(origin_7_imgs) < 100:
origin_7_imgs.append(imgs[i])
看看shape
print(np.array(origin_7_imgs).shape) # (100, 784)
把10张图片排成2x5的表格
由于一张图片是一个784维的一维数组，变成我们想看的图片就需要把它reshape成28x28的二维数组，然后再用Image里的方法，把它拼成一张2x5的大图。
由于tensorflow中MNIST都是灰度图（L），所以shape是（55000，784），每张图的dtype是float32，如果是彩色图（RGB），shape可能是（55000，784，3），图的dtype是uint8，从array转到Image需要用下面的方法：
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)) # 第二个参数为每次粘贴起始点的横纵坐标。在本例中，分别为（0，0）（28，0）（28*2，0）依次类推，第二行是（0，28）（28，28），（28*2，28）类推 new_img.paste(each_img, ((i % col) * each_width, (i / col) * each_width)) returnnew_img
效果图
ten_origin_7_imgs=comb_imgs(origin_7_imgs, 10, 10, 28, 28, 'L') ten_origin_7_imgs.show()
数字7原图
实现主成分分析算法（详细代码解析在文章后面的原理部分）
def pca(data_mat, top_n_feat=99999999):
"""
主成分分析：
输入：矩阵data_mat ，其中该矩阵中存储训练数据，每一行为一条训练数据 保留前n个特征top_n_feat，默认全保留
返回：降维后的数据集和原始数据被重构后的矩阵（即降维后反变换回矩阵）
"""
# 获取数据条数和每条的维数
num_data,dim = data_mat.shape
print(num_data) # 100
print(dim) # 784
# 数据中心化，即指变量减去它的均值
mean_vals = data_mat.mean(axis=0) #shape:(784,)
mean_removed = data_mat - mean_vals # shape:(100, 784)
# 计算协方差矩阵（Find covariance matrix）
cov_mat = np.cov(mean_removed, rowvar=0) # shape：(784, 784)
# 计算特征值(Find eigenvalues and eigenvectors)
eig_vals, eig_vects = linalg.eig(mat(cov_mat)) # 计算特征值和特征向量，shape分别为（784，）和(784, 784)
eig_val_index = argsort(eig_vals) # 对特征值进行从小到大排序，argsort返回的是索引，即下标 eig_val_index = eig_val_index[:-(top_n_feat + 1) : -1] # 最大的前top_n_feat个特征的索引
# 取前top_n_feat个特征后重构的特征向量矩阵reorganize eig vects,
# shape为(784, top_n_feat)，top_n_feat最大为特征总数
reg_eig_vects = eig_vects[:, eig_val_index]
# 将数据转到新空间
low_d_data_mat = mean_removed * reg_eig_vects # shape: (100, top_n_feat), top_n_feat最大为特征总数
recon_mat = (low_d_data_mat * reg_eig_vects.T) + mean_vals # 根据前几个特征向量重构回去的矩阵，shape:(100, 784)
return low_d_data_mat, recon_mat
调用PCA进行降维
low_d_feat_for_7_imgs, recon_mat_for_7_imgs = pca(np.array(origin_7_imgs), 1) # 只取最重要的1个特征 print(low_d_feat_for_7_imgs.shape) # (100, 1) print(recon_mat_for_7_imgs.shape) # (100, 784)
看降维后只用1个特征向量重构的效果图
low_d_img = comb_imgs(recon_mat_for_7_imgs, 10, 10, 28, 28, 'L') low_d_img.show()
数字7降维后的图
降维前后对比图
不难发现降维后数字7长得规则多了，或许降维后再用tensorflow入门教程的softmax进行分类accuracy会更高。
前面转坐标轴从理论上考虑，这里主要从数学的角度考虑。
第一个主成分是数据差异最大（方差最大）的方向，第二个主成分是数据差异次大且与第一个主成分正交的方向。通过数据集的协方差矩阵及其特征值分析，就能求得这些主成分的值。
统计学中的几个概念
平均值
这个最为熟悉最不容易忘记，描述样本集合的中间点。
标准差
描述样本集合中各个点到平均值的距离。
方差
标准差的平方。
方差
协方差
方差是用来描述一维数据的，协方差用来描述二维数据，用来描述两个随机变量之间的关系，如果是正值，则说明两变量正相关，负值，说明负相关，0，说明不相关，即相互独立。
协方差
从公式可以看出协方差的一些性质：
1、cov(X, X) = var(X)
2、cov(X,Y) = cov(Y, X)
协方差矩阵
协方差可以描述二维数据，但是对于多维数据来说，我们只能两个点两个点地计算多次协方差，一个n维数据，我们需要计算C(n, 2)=A(n,2)/2=n!/((n-2)!*2)个协方差，自然就需要用矩阵来组织这些数据。所以协方差矩阵的定义为：
协方差矩阵
比如数据集有三个维度，X，Y，Z，则协方差矩阵为
三维协方差矩阵
可见，矩阵的对角线为方差，由于cov(X,Y) = cov(Y, X)，所以是一个对称矩阵。
注意，协方差矩阵计算的是不同维度之间的协方差，不是不同样本之间的协方差。
目的就是找出差异最大的方向，也就是影响最大的几个特征，数学上通过协方差矩阵来找差异最大的特征，排序，最后找到降维后的特征矩阵。
# 数据中心化，即指变量减去它的均值 mean_vals = data_mat.mean(axis=0) #shape:(784,) mean_removed = data_mat - mean_vals # shape:(100, 784) # 计算协方差矩阵（Find covariance matrix） cov_mat = np.cov(mean_removed, rowvar=0)
协方差矩阵需要计算平均值，上面强调了计算的是不同维度的协方差，数据每行是一个样本，每列是一个维度，因此计算的是列的平均值，即axis=0，因此shape为（784，）。使用np的cov函数计算协方差矩阵，api入下：
numpy.cov(m, y=None, rowvar=True, bias=False, ddof=None, fweights=None, aweights=None)[source]
**rowvar代表是否转置。在API里，默认rowvar是True，也就是行是variable，列是observation，我们这里列是observation，行是variable。
** eig_vals, eig_vects = linalg.eig(mat(cov_mat)) # 计算特征值和特征向量
mat(cov_mat)：将输入转成矩阵。和matrix，asmatrix不同，如果输入已经是举着嗯或者ndarray，它不会制作副本。相当于matrix(data, copy=False)
详细API请点这里（https://docs.scipy.org/doc/numpy/reference/generated/numpy.mat.html）linalg.eig(a)：计算特征值和特征向量
详细API请点这里（https://docs.scipy.org/doc/numpy/reference/generated/numpy.linalg.eig.html）
矩阵乘法对应了一个变换，在这个变换的过程中，原向量主要发生旋转、伸缩的变化。如果矩阵对某一个向量或某些向量只发生伸缩变换，不对这些向量产生旋转的效果，那么这些向量就称为这个矩阵的特征向量，伸缩的比例就是特征值。
eig_val_index = argsort(eig_vals) # 对特征值进行从小到大排序，argsort返回的是索引，即下标
numpy.argsort(a, axis=-1, kind='quicksort', order=None)
详细API请点这里（https://docs.scipy.org/doc/numpy/reference/generated/numpy.argsort.html）
eig_val_index = eig_val_index[:-(top_n_feat + 1) : -1] # 最大的前top_n_feat个特征的索引 reg_eig_vects = eig_vects[:, eig_val_index]
这里有一个语法问题，[::]代表切片，[开始：结束：步长]，负号代表从后往前每隔步长个取一次，比如有一个array[1, 2, 3, 4, 5]，取[:-4:-2]，0是第一个，-1是最后一个(在这里是5的下标)，从最后一个往前数，一直数到-4（在这里是2的下标），每两个取1个数，最后得到的array是[5, 3]。
[:, eig_val_index]代表第一维不变，第二维要eig_val_index个，所以它的shape是（784，top_n_feat）
# 将数据转到新空间 low_d_data_mat = mean_removed * reg_eig_vects # shape: (100, top_n_feat), top_n_feat最大为特征总数 recon_mat = (low_d_data_mat * reg_eig_vects.T) + mean_vals # 根据前几个特征向量重构回去的矩阵，shape:(100, 784)
一个shape是（100，784）的矩阵，乘以一个shape是（784，top_n_feat）的矩阵，最后得到降维的矩阵（100， top_n_feat）
recon_mat再将矩阵变回（100，784），得到降维后再重构的矩阵。
优点：降低数据的复杂性，识别最重要的特征
缺点：不一定需要，且可能损失有用信息
适用数据类型：数值型数据
原文链接：http://www.jianshu.com/p/b9f2c92dfeaa
查阅更为简洁方便的分类文章以及最新的课程、产品信息，请移步至全新呈现的“LeadAI学院官网”：
www.leadai.org
请关注人工智能LeadAI公众号，查看更多专业文章
大家都在看
TensorFlow从1到2
| 第三章 深度学习革命的开端：卷积神经网络
点击“阅读原文”直接打开【北京站 | GPU CUDA 进阶课程】报名链接