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=B^TB/(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中每个基下的坐标的方差尽量大，而不同基下坐标的方差尽量小，亦即我们希望C_Y=YY^T/(n-1)是一个对角线矩阵。

  C_Y  =YY^T/(n-1)=P(XX^T)P^T/(n-1)

  令A=XX^T，我们对A进行分解：A=EDE^T

  我们取P=E^T，则C_Y=E^TAE/(n-1)=E^TEDE^TE/(n-1),因为E^T=E^-1，所以C_Y=D/(n-1)是一个对角矩阵。

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

  3. 实现

  1. PCA的Python实现

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

  import numpy as np

  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)]
  #转置,每一行是一条数据
  data=np.matrix(np.transpose(mat))
  data_adjust=data-mean
  #求协方差矩阵
  covariance=np.transpose(data_adjust)*data_adjust/9
  #求得协方差矩阵的特征值和特征向量
  eigenvalues,eigenvectors=np.linalg.eig(covariance)         
  feature_vectors=np.transpose(eigenvectors)
  #转换后的数据
  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)理论分析及应用