zoukankan      html  css  js  c++  java
  • 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)理论分析及应用

  • 相关阅读:
    gThumb 3.1.2 发布,支持 WebP 图像
    航空例行天气预报解析 metaf2xml
    Baruwa 1.1.2 发布,邮件监控系统
    Bisect 1.3 发布,Caml 代码覆盖测试
    MoonScript 0.2.2 发布,基于 Lua 的脚本语言
    Varnish 入门
    快速增量备份程序 DeltaCopy
    恢复模糊的图像 SmartDeblur
    Cairo 1.12.8 发布,向量图形会图库
    iText 5.3.4 发布,Java 的 PDF 开发包
  • 原文地址:https://www.cnblogs.com/kemaswill/p/2858387.html
Copyright © 2011-2022 走看看