Python使用三種方法實(shí)現(xiàn)PCA算法

主成分分析，即Principal Component Analysis（PCA），是多元統(tǒng)計(jì)中的重要內(nèi)容，也廣泛應(yīng)用于機(jī)器學(xué)習(xí)和其它領(lǐng)域。它的主要作用是對(duì)高維數(shù)據(jù)進(jìn)行降維。PCA把原先的n個(gè)特征用數(shù)目更少的k個(gè)特征取代，新特征是舊特征的線性組合，這些線性組合最大化樣本方差，盡量使新的k個(gè)特征互不相關(guān)。

主成分分析（PCA） vs 多元判別式分析（MDA）

PCA和MDA都是線性變換的方法，二者關(guān)系密切。在PCA中，我們尋找數(shù)據(jù)集中最大化方差的成分，在MDA中，我們對(duì)類間最大散布的方向更感興趣。

一句話，通過PCA，我們將整個(gè)數(shù)據(jù)集（不帶類別標(biāo)簽）映射到一個(gè)子空間中，在MDA中，我們致力于找到一個(gè)能夠最好區(qū)分各類的最佳子集。粗略來講，PCA是通過尋找方差最大的軸（在一類中，因?yàn)?a href='/map/pca/' style='color:#000;font-size:inherit;'>PCA把整個(gè)數(shù)據(jù)集當(dāng)做一類），在MDA中，我們還需要最大化類間散布。

在通常的模式識(shí)別問題中，MDA往往在PCA后面。

PCA的主要算法如下：

1. 組織數(shù)據(jù)形式，以便于模型使用；
2. 計(jì)算樣本每個(gè)特征的平均值；
3. 每個(gè)樣本數(shù)據(jù)減去該特征的平均值（歸一化處理）；
4. 求協(xié)方差矩陣；
5. 找到協(xié)方差矩陣的特征值和特征向量；
6. 對(duì)特征值和特征向量重新排列（特征值從大到小排列）；
7. 對(duì)特征值求取累計(jì)貢獻(xiàn)率；
8. 對(duì)累計(jì)貢獻(xiàn)率按照某個(gè)特定比例，選取特征向量集的字跡合；
9. 對(duì)原始數(shù)據(jù)（第三步后）。

其中協(xié)方差矩陣的分解可以通過按對(duì)稱矩陣的特征向量來，也可以通過分解矩陣的SVD來實(shí)現(xiàn)，而在Scikit-learn中，也是采用SVD來實(shí)現(xiàn)PCA算法的。

本文將用三種方法來實(shí)現(xiàn)PCA算法，一種是原始算法，即上面所描述的算法過程，具體的計(jì)算方法和過程，可以參考：A tutorial on Principal Components Analysis, Lindsay I Smith. 一種是帶SVD的原始算法，在Python的Numpy模塊中已經(jīng)實(shí)現(xiàn)了SVD算法，并且將特征值從大從小排列，省去了對(duì)特征值和特征向量重新排列這一步。最后一種方法是用Python的Scikit-learn模塊實(shí)現(xiàn)的PCA類直接進(jìn)行計(jì)算，來驗(yàn)證前面兩種方法的正確性。

用以上三種方法來實(shí)現(xiàn)PCA的完整的Python如下：

import numpy as np
from sklearn.decomposition import PCA
import sys
#returns choosing how many main factors
def index_lst(lst, component=0, rate=0):
  #component: numbers of main factors
  #rate: rate of sum(main factors)/sum(all factors)
  #rate range suggest: (0.8,1)
  #if you choose rate parameter, return index = 0 or less than len(lst)
  if component and rate:
    print('Component and rate must choose only one!')
    sys.exit(0)
  if not component and not rate:
    print('Invalid parameter for numbers of components!')
    sys.exit(0)
  elif component:
    print('Choosing by component, components are %s......'%component)
    return component
  else:
    print('Choosing by rate, rate is %s ......'%rate)
    for i in range(1, len(lst)):
      if sum(lst[:i])/sum(lst) >= rate:
        return i
    return 0
 
def main():
  # test data
  mat = [[-1,-1,0,2,1],[2,0,0,-1,-1],[2,0,1,1,0]]
   
  # simple transform of test data
  Mat = np.array(mat, dtype='float64')
  print('Before PCA transforMation, data is:\n', Mat)
  print('\nMethod 1: PCA by original algorithm:')
  p,n = np.shape(Mat) # shape of Mat
  t = np.mean(Mat, 0) # mean of each column
   
  # substract the mean of each column
  for i in range(p):
    for j in range(n):
      Mat[i,j] = float(Mat[i,j]-t[j])
       
  # covariance Matrix
  cov_Mat = np.dot(Mat.T, Mat)/(p-1)
   
  # PCA by original algorithm
  # eigvalues and eigenvectors of covariance Matrix with eigvalues descending
  U,V = np.linalg.eigh(cov_Mat)
  # Rearrange the eigenvectors and eigenvalues
  U = U[::-1]
  for i in range(n):
    V[i,:] = V[i,:][::-1]
  # choose eigenvalue by component or rate, not both of them euqal to 0
  Index = index_lst(U, component=2) # choose how many main factors
  if Index:
    v = V[:,:Index] # subset of Unitary matrix
  else: # improper rate choice may return Index=0
    print('Invalid rate choice.\nPlease adjust the rate.')
    print('Rate distribute follows:')
    print([sum(U[:i])/sum(U) for i in range(1, len(U)+1)])
    sys.exit(0)
  # data transformation
  T1 = np.dot(Mat, v)
  # print the transformed data
  print('We choose %d main factors.'%Index)
  print('After PCA transformation, data becomes:\n',T1)
   
  # PCA by original algorithm using SVD
  print('\nMethod 2: PCA by original algorithm using SVD:')
  # u: Unitary matrix, eigenvectors in columns
  # d: list of the singular values, sorted in descending order
  u,d,v = np.linalg.svd(cov_Mat)
  Index = index_lst(d, rate=0.95) # choose how many main factors
  T2 = np.dot(Mat, u[:,:Index]) # transformed data
  print('We choose %d main factors.'%Index)
  print('After PCA transformation, data becomes:\n',T2)
   
  # PCA by Scikit-learn
  pca = PCA(n_components=2) # n_components can be integer or float in (0,1)
  pca.fit(mat) # fit the model
  print('\nMethod 3: PCA by Scikit-learn:')
  print('After PCA transformation, data becomes:')
  print(pca.fit_transform(mat)) # transformed data      
main()

運(yùn)行以上代碼，輸出結(jié)果為：

這說明用以上三種方法來實(shí)現(xiàn)PCA都是可行的。這樣我們就能理解PCA的具體實(shí)現(xiàn)過程啦~~有興趣的讀者可以用其它語言實(shí)現(xiàn)一下哈