PCA数学原理及编程实现

一. 一个场景

已知一家超市，销售A，B，C，D四种产品，现对每种产品的一周之内每天的销售情况记录如下：

         A    B   C    D

周一 2     0    8    9

周二 4     0   11  13

周三 3     1   10  12

周四 2     3   11  10

周五 1     0   12    9

周六10   11    1    2

周日11   12    2    1

我们将A，B，C，D在一周内的销售情况分别记为事件A，事件B，事件C，事件D，则A，B，C，D均为一个7维的向量，即我们可以得到4条7维向量组成的数据。

不难发现，A，B事件相似度高(在周一到周五销售较少，在周六和周日销售较多)； C，D事件相似度高。于是，我们可以假设存在一个向量t，将A，B，C，D分别投射到t上，得到一个新的数据：

    A    B   C    D

t   A'   B'   C'   D'

如果在这个向量t上，A和B的距离很近，C和D的距离很近，同时AB，CD的距离很远，我们就成功的将一个7维的数据降维到了一维。  

二. 降维目标抽象化

在上述场景中，我们并不知道向量t是什么，也未明确如何表征在向量t上，A，B，C，D之间的距离，同时，如果要将数据降为二维改如何表述呢？

实际上，对于上述t向量，我们需要使所有事件在向量t上的投影所组成的新数组(A‘，B'，C'，D')的方差最大，也就是在在t向量方向上，所有事件的差异可以最大化，故t向量即为这组4条7维数组的第一主成分。而在与第一主成分垂直的六维空间里，再寻找第二向量t2，使得所有事件在向量t2上的投影所组成的新数组(A''，B''，C''，D'')的方差最大，向量t2即为这组4条7维数组的第二主成分，以此类推，可获得第一到第七主成分。

降维问题的优化目标：将一组N维向量(每一列为一个记录)降为K维（K大于0，小于N），其目标是选择K个单位（模为1）正交基，使得原始数据变换到这组基上后，在各维度上两两间协方差为0，而每个维度内的方差则尽可能大（在正交的约束下，取最大的K个方差）。

三.算法推演

1.假设原始矩阵为X(每列为一个事件，每行为一个特征维度)，X的协方差矩阵为C1；变换基矩阵为P，变换后的矩阵为Y，Y的协方差矩阵为C2；PCA的实质就是找到矩阵P，使得C2为对角阵。

2.又, Y=PX, C1=Xt(X), C2=Yt(Y)； 故C2=(PX)(t(PX))=P(X(t(X))t(P)=PC1t(P)

3.故，优化目标变成了寻找一个变换基矩阵P，满足PC1t(P)是一个对角矩阵，并且对角元素按从大到小依次排列，那么P的前K行就是要寻找的基，用P的前K行组成的矩阵乘以X就使得X从N维降到了K维并满足上述优化条件。

4.根据协方差矩阵的特征，以其特征值为行的矩阵即为变化矩阵P。

四.算法

设有m条n维数据；

1）将原始数据按列组成n行m列矩阵X（每一列代表一个事件，每一行代表一个特征值维度）

2）将X的每一行（代表一个属性字段）进行零均值化，即减去这一行的均值

3）求出协方差矩阵C=(1/m)Xt(X)

4）求出协方差矩阵的特征值及对应的特征向量(模为1)

5）将特征向量按对应特征值大小从上到下按行排列成矩阵，取前k行组成矩阵P

6）Y=PX即为降维到k维后的数据

五.算法实现：R语言

> A=c(2,4,3,2,1,10,11)

> B=c(0,0,1,3,0,11,12)
> C=c(8,11,10,11,12,1,2)

> D=c(9,13,12,10,9,2,1)

> sale=data.frame(A,B,C,D)

> rownames(sale)=c('MON','TUE','WED','THU',"FRI",'SAT','SUN')

> data.pca<-prcomp(sale)

> data.pca
Standard deviations:
[1] 9.0818425 1.6889849 1.3284286 0.6531112

Rotation:
         PC1         PC2         PC3        PC4
A  0.4302964  0.62644196 -0.01270685 -0.6498108
B  0.5787863  0.11382130 -0.62985188  0.5053095
C -0.4840149 -0.09097089 -0.77652218 -0.3930232
D -0.4955613  0.76572807  0.01176625  0.4098062
> plot(x=data.pca$rotation[,1],y=data.pca$rotation[,2],pch=19,col=c(rep(2,length(data))),cex=1.25)

> biplot(data.pca)
> plot(data.pca, type="lines")