本文主要基于同名的两篇外文参考文献 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 T B/(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 T AE/(n-1)=E T EDE T E/(n-1),因为E T =E -1 ，所以C Y =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)理论分析及应用
