PCA

来源：互联网发布：编写软件程序步骤编辑：程序博客网时间：2024/06/06 02:41

PCA 经常用来进行数据的降维，压缩，在一定程度上也减少了随机噪声和一些微小误差产生的影响，一来可以方便进行数据可视化，二来可以减少原始数据的冗余。

假设原始数据的维度是 n, 我们要把它投影到 k 维空间，通过 PCA 方法选择的投影方向的直观解释有很多，比如:

使得投影后的样本点在 k 维空间中有最大分散度，也就是方差最大
使得投影误差最小
投影后模的平方的和最大
投影后的协方差矩阵是对角阵，这样就大大减少了数据的冗余

也可以说 PCA 的关键就是找一个投影方向 (也可以叫做主方向，变化幅度最大的方向)，使得在这个投影方向下，满足以上条件.

PCA 的步骤：
假设原始数据是 Xm×n ，要降到　k 维，Bm×k, 要找的投影方向矩阵为 Qn×k

对原始数据 X 做中心化处理，减掉其列均值，这样使得所有数数据点的中心就是原点，也方便了以后进行求协方差矩阵.
$X - X 的列均值＝ A$
对 A 的协方差矩阵 ATA 作奇异值分解(SVD):
$A T A = U Σ U T$
其中　对角矩阵 Σ 对角线的值由大到小排列
选取特征向量矩阵 U 的前 k 列,组成投影方向矩阵 Q
数据的投影: AQ=B
数据的恢复: BQT+X的列均值

这里有个衡量降维后数据的信息损失度的问题，k 越大，损失度自然就越低，可以用　S=∑i=1kSii∑j=1nSjj 来衡量，其中　Sii 指矩阵Σ 的第 i 个对角线元素，一般当　S>0.99 时，我们可以认为数据信息几乎没有损失。

下面通过实验进一步了解PCA

首先定义 PCA 函数

import numpy as npdef my_pca(x,k):    colmean=x.mean(axis=0)    A=x-colmean                         # 中心化                U,sigma,V=np.linalg.svd(A.T*A)      # 奇异值分解    S=np.sum(sigma[:k])/np.sum(sigma)   # 计算 S    Q=U[:,0:k]                          # 选择前 k 个特征向量    B=A*Q                               # 投影    rec=B*Q.T + colmean                 # 数据恢复    return B,S,Q,rec

例 1 :
找出主方向
这里写图片描述

可以看出，如果投影到一维空间，那么投影到红色直线所在的方向是最好的保持了数据的离散度。

例 2 :
k 的大小，也就是 S 的大小，对数据恢复的影响，原始图像的大小是 220 x 220

这里写图片描述

可以看出，随着 k 的减小，数据的恢复效果在逐渐下降。

例 3 :
在数据恢复时， A=BQT，然后在加上以前减掉的列均值，也就是说原始数据可以由 B 的列线性表示，相当于选取了另一组基，在人脸识别中，B 的列又称为特征脸.

以下是来自人脸数据库 AT&T http://www.cl.cam.ac.uk/research/dtg/attarchive/facedatabase.html 的 40 个人脸图像，每个图像的大小是 112 x 92
这里写图片描述

通过PCA可提取特征脸，特征脸如下
这里写图片描述

也就是认为，上面的每一个人脸图像，都是由这些特征脸线性组合成的

例1代码

import numpy as npimport matplotlib.pyplot as pltx1=np.random.normal(1,3,size=(100,2))x2=np.random.normal(10,3,size=(100,2))x=np.mat(np.concatenate((x1,x2),axis=0))A=x-x.mean(axis=0)b,s,q,*_=my_pca(x,2)fig,ax=plt.subplots(figsize=(5,5))ax.scatter(A[:,0],A[:,1])ax.plot([0,q[0,0]*10],[0,q[1,0]*10],color='r')ax.plot([0,-q[0,0]*10],[0,-q[1,0]*10],color='r')ax.plot([0,q[0,1]*5],[0,q[1,1]*5],color='g')ax.plot([0,-q[0,1]*5],[0,-q[1,1]*5],color='g')ax.axis='equal'

例2代码

import numpy as npimport matplotlib.pyplot as pltfrom scipy import misca=misc.imread('lena.png')img=np.mat(a[:,:,1])fig,axes=plt.subplots(1,4)ax1,ax2,ax3,ax4 = axes.flatax1.imshow(img,cmap=plt.cm.gray)ax1.set_title('original')ax1.set_xticks([])ax1.set_yticks([])B,S,Q,rec = my_pca(img,60)ax2.imshow(rec,cmap=plt.cm.gray)ax2.set_title('k={} S={:0.2f}'.format(60,S))ax2.set_xticks([])ax2.set_yticks([])B,S,Q,rec = my_pca(img,20)ax3.imshow(rec,cmap=plt.cm.gray)ax3.set_title('k={} S={:0.2f}'.format(20,S))ax3.set_xticks([])ax3.set_yticks([])B,S,Q,rec = my_pca(img,5)ax4.imshow(rec,cmap=plt.cm.gray)ax4.set_title('k={} S={:0.2f}'.format(5,S))ax4.set_xticks([])ax4.set_yticks([])

例3代码

import numpy as npimport matplotlib.pyplot as pltfrom scipy import miscfig,axes=plt.subplots(5,8,figsize=(6,4))for i in range(40):    exec('a=misc.imread(\'{}.pgm\')'.format(i+1))    if i == 0:        xx=a.reshape(-1,1)    else:        xx=np.concatenate((xx, a.reshape(-1,1)),axis=1)    img=np.mat(a)    axes.flat[i].imshow(img,cmap=plt.cm.gray)    axes.flat[i].set_xticks([])    axes.flat[i].set_yticks([])x=np.mat(xx)B,S,Q,rec=my_pca(x,40)fig,axes=plt.subplots(5,8,figsize=(6,4))for i in range(40):    a=B[:,i].reshape(112,92)    img=np.mat(a)    axes.flat[i].imshow(img,cmap=plt.cm.gray)    axes.flat[i].set_xticks([])    axes.flat[i].set_yticks([])

0 0