PCA（Principal Component Analysis）

PCA通过线性变换将原始数据变换为一组各维度线性无关的表示，可用于提取数据的主要特征分量，常用于高维数据的降维。

PCA算法是如何实现的？

简单来说，就是将数据从原始的空间中转换到新的特征空间中，例如原始的空间是三维的(x,y,z)，x、y、z分别是原始空间的三个基，我们可以通过某种方法，用新的坐标系(a,b,c)来表示原始的数据，那么a、b、c就是新的基，它们组成新的特征空间。在新的特征空间中，可能所有的数据在c上的投影都接近于0，即可以忽略，那么我们就可以直接用(a,b)来表示数据，这样数据就从三维的$(x,y,z)$降到了二维的(a,b)。

问题是如何求新的基(a,b,c)?

一般步骤是这样的：先对原始数据零均值化，然后求协方差矩阵，接着对协方差矩阵求特征向量和特征值，这些特征向量组成了新的特征空间。具体的细节，推荐Andrew Ng的网页教程：Ufldl 主成分分析 ，写得很详细。

PCA本质上是将方差最大的方向作为主要特征，并且在各个正交方向上将数据“离相关”，也就是让它们在不同正交方向上没有相关性。

python算法实现

零均值化（特征中心化）

假如原始数据集为矩阵data_mat，data_mat中每一行代表一个样本，每一列代表同一个特征。零均值化就是求每一列的平均值，然后该列上的所有数都减去这个均值。也就是说，这里零均值化是对每一个特征而言的，零均值化都，每个特征的均值变成0。实现代码如下：

def zero_mean(data_mat):
    mean_val = np.mean(data_mat,axis=0)     #按列求均值，即求各个特征的均值
    data = data_mat-mean_val
    return data, mean_val

函数中用numpy中的mean方法来求均值，axis=0表示按列求均值。
该函数返回两个变量，data是零均值化后的数据，mean_val是每个特征的均值，是给后面重构数据用的。

求协方差矩阵

data, mean_val = zero_mean(data_mat)
cov_mat = np.cov(data,rowvar=False)

numpy中的cov函数用于求协方差矩阵，cov_mat即所求的协方差矩阵。

求特征值、特征向量

调用numpy中的线性代数模块linalg中的eig函数，可以直接由cov_mat求得特征值和特征向量：

eig_vals, eig_vects=np.linalg.eig(np.mat(cov_mat))

eig_vals存放特征值，行向量。
eig_vects存放特征向量，每一列带别一个特征向量。
特征值和特征向量是一一对应的

保留主要的成分[即保留值比较大的前n个特征]

第三步得到了特征值向量eig_vals，假设里面有m个特征值，我们可以对其排序，排在前面的n个特征值所对应的特征向量就是我们要保留的，它们组成了新的特征空间的一组基n_eig_vect。将零均值化后的数据乘以n_eig_vect就可以得到降维后的数据。代码如下：

eig_val_indice = np.argsort(eig_vals)            # 对特征值从小到大排序

n_eig_val_indice = eig_val_indice[-1:-(n_components+1):-1]        #最大的n_components个特征值的下标

n_eig_vect = eig_vects[:,n_eig_val_indice]        #最大的n个特征值对应的特征向量

lowD_data_mat = data * n_eig_vect                 # 低维特征空间的数据

recon_mat =(lowD_data_mat * n_eig_vect.T) + mean_val  #重构数据

选择主成分个数

文章写到这里还没有完，应用PCA的时候，对于一个1000维的数据，我们怎么知道要降到几维的数据才是合理的？即n要取多少，才能保留最多信息同时去除最多的噪声？一般，我们是通过方差百分比来确定n的，这一点在Ufldl教程中说得很清楚，并且有一条简单的公式，下面是该公式的截图：

一般而言，设$x_1,x_2,...,x_n$表示$\sum$的特征值（由大到小排序），使得对应的$x_j$为对应于特征向量$u_j$的特征值。那么如果我们保留前$k$个成分，则保留的方差百分比可计算为

$$\frac{\sum_{j=1}^kx_j}{\sum_{j=1}^nx_j}$$

根据这条公式，可以写个函数，函数传入的参数是百分比percentage和特征值向量，然后根据percentage确定n，代码如下

def percentage2n(eigVals,percentage):  
    sortArray=np.sort(eigVals)   #升序  
    sortArray=sortArray[-1::-1]  #逆转，即降序  
    arraySum=sum(sortArray)  
    tmpSum=0  
    num=0  
    for i in sortArray:  
        tmpSum+=i  
        num+=1  
        if tmpSum>=arraySum*percentage:  
            return num  

预测：

def transform(dat, n_eig_vect , mean_val ):
     data = dat - mean_val
     return data * n_eig_vect

总代码：

def zero_mean(data_mat):
    mean_val = np.mean(data_mat,axis=0)     #按列求均值，即求各个特征的均值
    data = data_mat - mean_val
    return data, mean_val


def pca(data_mat,n_components = 2):

    data, mean_val = zero_mean(data_mat)
    cov_mat = np.cov(data,rowvar=False)

    eig_vals, eig_vects=np.linalg.eig(np.mat(cov_mat))

    eig_val_indice = np.argsort(eig_vals)            # 对特征值从小到大排序

    n_eig_val_indice = eig_val_indice[-1:-(n_components+1):-1]        #最大的n_components个特征值的下标

    n_eig_vect = eig_vects[:,n_eig_val_indice]        #最大的n个特征值对应的特征向量

    lowD_data_mat = data * n_eig_vect                 # 低维特征空间的数据

    recon_mat =(lowD_data_mat * n_eig_vect.T) + mean_val  #重构数据

    return lowD_data_mat, recon_mat


def percentage2n(eig_vals, percentage):
    sortArray=np.sort(eig_vals)   #升序
    sortArray=sortArray[-1::-1]  #逆转，即降序
    arraySum=sum(sortArray)
    tmpSum=0
    num=0
    for i in sortArray:
        tmpSum+=i
        num+=1
        if tmpSum>=arraySum*percentage:
            return num

python第三方实现

# coding=utf8

import numpy as np

from sklearn.decomposition import PCA
import matplotlib.pyplot as plt

# t-distributed Stochastic Neighbor Embedding.
tsne = TSNE(perplexity=30, n_components=2, init='pca', n_iter=5000)

a = np.array([[1,2,3,4,5,6],
     [6,5,4,3,2,1],
     [2,2,2,2,2,2]])

b = [[1,2,3,4,5,6],
     [1,2,1,2,1,2],
     [6,5,4,3,2,1],
     [2,2,2,2,2,2]]

pca = PCA(2,whiten=True)

# c = pca.fit(a)    # train

c = pca.fit_transform(a)

d = pca.transform(b)

print c
print d

plt.figure(figsize=(10, 5))
plt.subplot(121)
plt.scatter(c[:, 0], c[:, 1])
plt.subplot(122)
plt.scatter(d[:, 0], d[:, 1])
plt.show()

较大规模下的PCA降维方法

http://xueshu.baidu.com/s?wd=paperuri%3A%28c866d9a88ecc4b604a9e1069b646bd7c%29&filter=sc_long_sign&tn=SE_xueshusource_2kduw22v&sc_vurl=http%3A%2F%2Fwenku.baidu.com%2Fview%2F3212768be87101f69e3195e1.html&ie=utf-8&sc_us=4083164024310701723

附

回到目录

引

(1) PCA的数学原理

(2) 原文

(3) 主成分分析(PCA)——基于python+numpy

(4) scikit-learn中PCA的使用方法