PCA算法 原理与实现

  本文主要基于同名的两篇外文参考文献A Tutorial on Principal Component Analysis

  PCA,亦即主成分分析,主要用于对特征进行降维。如果数据的特征数非常多,我们可以认为其中只有一部分特征是真正我们感兴趣和有意义的,而其他特征或者是噪音,或者和别的特征有冗余。从所有的特征中找出有意义的特征的过程就是降维,而PCA是降维的两个主要方法之一(另一个是LDA).

  Jonathon Shlens的论文中举了一个物理学中测试理想情况下弹簧振动的例子,非常生动,详见[1](中文翻译见[5])。

  我们首先看一下给定一个代表数据记录的矩阵A,如果计算其主成分P,并如何利用P得到降维后的数据矩阵B,然后介绍一下这个计算过程背后的原理,最后会有在Python中实现PCA和在Weka中调用PCA算法的实例。

  1. 计算过程:

  假设我们有n条数据记录,每条记录都是m维,我们可以把这些数据记录表示成一个n*m矩阵A。

  对矩阵A的每一列,求出其平均值,对于A中的每一个元素,减去该元素所在列的平均值,得到一个新的矩阵B。

  计算矩阵Z=BTB/(n-1)。其实m*m维矩阵Z就是A矩阵的协方差矩阵。

  计算矩阵Z的特征值D和特征向量V,其中D是1*m矩阵,V是一个m*m矩阵,D中的每个元素都是Z的特征值,V中的第i列是第i个特征值对应的特征向量。

  下面,就可以进行降维了,假设我们需要把数据由m维降到k维,则我们只需要从D中挑选出k个最大的特征向量,然后从V中挑选出k个相应的特征向量,组成一个新的m*k矩阵N。

  N中的每一列就是A的主成分(Principal Component). 计算A*N得到n*k维矩阵C,就是对源数据进行降维后的结果,没条数据记录的维数从m降到了k。

  2. 原理

  要对数据进行降维的主要原因是数据有噪音,数据的轴(基)需要旋转,数据有冗余。

  (1) 噪音

  上图是一个记录弹簧振动的二维图。我们发现沿正45度方向数据的方差比较大,而沿负45度方向数据的方差比较小。通常情况下,我们都认为方差最大的方向记录着我们感兴趣的信息,所以我们可以保留正45度方向的数据信息,而负45度方向的数据信息可以认为是噪音。

  (2) 旋转

  在线性代数中我们知道,同一组数据,在不同的基下其坐标是不一样的,而我们一般认为基的方向应该与数据方差最大的方向一致,亦即上图中的基不应该是X,Y轴,而该逆时针旋转45度。

  (3) 冗余

  上图中的a,c分别代表没有冗余和高度冗余的数据。在a中,某个数据点的X,Y轴坐标值基本上是完全独立的,我们不可能通过其中一个来推测另外一个,而在c中,数据基本上都分布在一条直线上,我们很容易从一个坐标值来推测出另外一个坐标值,所以我们完全没有必要记录数据在X,Y两个坐标轴上的值,而只需要记录一个即可。数据冗余的情况跟噪音比较相似,我们只需要记录方差比较大的方向上的坐标值,方差比较小的方向上的坐标值可以看做是冗余(或噪音).

  上面三种情况,归结到最后都是要求出方差比较大的方向(基),然后在求数据在这个基下的坐标,这个过程可以表示为:

  PX=Y。

  其中k*m矩阵P是一个正交矩阵,P的每一行都对应一个方差比较大的基。m*n矩阵X和k*n矩阵Y的每一列都是一个数据(这一点和1中不同,因为这是两篇不同的论文,表示方式不一样,本质上是一样的)。

  X是原数据,P是一个新的基,Y是X在P这个新基下的坐标,注意Y中数据记录的维数从m降到了k,亦即完成了降维。

  但是我们希望得到的是一个什么样的Y矩阵呢?我们希望Y中每个基下的坐标的方差尽量大,而不同基下坐标的方差尽量小,亦即我们希望CY=YYT/(n-1)是一个对角线矩阵。

  CY  =YYT/(n-1)=P(XXT)PT/(n-1)

  令A=XXT,我们对A进行分解:A=EDET

  我们取P=ET,则CY=ETAE/(n-1)=ETEDETE/(n-1),因为ET=E-1,所以CY=D/(n-1)是一个对角矩阵。

  所以我们应该取P的每一行为A的特征向量,得到的Y才会有以上性质。

  3. 实现

  1. PCA的Python实现

  需要使用Python的科学计算模块numpy

 1   import numpy as np
 2 
 3   mat=[(2.5,0.5,2.2,1.9,3.1,2.3,2,1,1.5,1.1),    (2.4,0.7,2.9,2.2,3.0,2.7,1.6,1.1,1.6,0.9)]
 4   #转置,每一行是一条数据
 5   data=np.matrix(np.transpose(mat))
 6   data_adjust=data-mean
 7   #求协方差矩阵
 8   covariance=np.transpose(data_adjust)*data_adjust/9
 9   #求得协方差矩阵的特征值和特征向量
10   eigenvalues,eigenvectors=np.linalg.eig(covariance)         
11   feature_vectors=np.transpose(eigenvectors)
12   #转换后的数据
13   final_data=feature_vectors*np.transpose(data_adjust)

  2. 在Weka中调用PCA:

import java.io.FileReader as FileReader
import java.io.File as File
import weka.filters.unsupervised.attribute.PrincipalComponents as PCA
import weka.core.Instances as Instances
import weka.core.converters.CSVLoader as CSVLoader
import weka.filters.Filter as Filter


def main():
    #使用Weka自带的数据集cpu.arff
    reader=FileReader('DATA/cpu.arff')
    data=Instances(reader)
    pca=PCA()
    pca.setInputFormat(data)
    pca.setMaximumAttributes(5)
    newData=Filter.useFilter(data,pca)
    for n in range(newData.numInstances()):
        print newData.instance(n)
if __name__=='__main__':
    main()

  参考文献:

  [1]. Jonathon Shlens. A Tutorial on Principal Component Analysis.

  [2]. Lindsay I Smith. A Tutorial on Principal Component Analysis.

  [3]. 关于PCA算法的一点学习总结

  [4]. 主成分分析PCA算法  原理解析

  [5]. 主元分析(PCA)理论分析及应用

posted on 2013-01-13 14:26  潘的博客  阅读(5638)  评论(1编辑  收藏

导航

统计