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)理論分析及應用