主成分分析理论与MATLAB实现示例

主成分分析理论与MATLAB实现示例

2017.9.3

关于主成分分析，其实还是比较简单的，但是感觉网上的资料对于理论上的分析有些欠缺，这里自己整理一下，依然是追求严谨完整，另外再加上一个自己练习的示例。

1.基本思想

现实中的数据常常会以高维数据的形式出现。举个简单的例子：对于平面上一条过原点斜线上的一些点，一般地我们需要它的(x,y)坐标来描述，但是实际上对于这个数据结构，我们只要知道点在斜线上的位置，比如到原点的距离，就可以间接地确定点的坐标。事实上，上面的过程就把二维的数据(x,y)转变为了一维（点到原点的距离）。对于有几百上千个分量的数据，我们同样可以通过选取某些“主成分”来代表这些高维数据，自然地我们希望这些主成分能够尽可能多的表示出原本数据的信息。

2.方法

接着上面的思想，我们将面对三个问题：

1.如何描述变换后前数据承载的信息？

2.将高维坐标变换成主成分坐标应该是怎样一种变换？

3.如何使得主成分能尽可能多的代表原来数据的信息？

1：当数据很集中时，比如说都在一个点的小邻域内，数据之间差别很小，此时的数据除了体现出一个总体趋势（均值）外，没有产生其他信息，而如果数据是分散开的，那么就认为它带来了较多的信息。因此在主成分分析中，我们选取数据的方差来表示所携带的信息量。此外方差还有一个性质：数据的方差等于数据每个坐标上的方差之和。

2：无疑，对一组坐标进行变换有很多种方法，主成分分析提供的便是最常见的线性变换，各个主成分也就是变换后的坐标。进一步地，我们希望不同的“主成分”代表不同的信息，或者说，他们之间是“无关的”，即各个坐标的变换之间是相互正交的。到这里，我们很容易想到正交变换，正交变换除了能满足上面的要求时，还保持了变换前后方差不变（这是因为方差依赖于距离，而正交变换是保距离的），因此正交变换是一个十分合理的选择。

3：在第一个问题中我们就指出了，我们选取方差作为信息的度量，因此，要尽可能多的获得信息，换成数学语言，就是如下的优化问题：

设我们有数据Xi=(Xi(0),⋯,Xi(m))∈Rm, i=1,⋯,n，X=(X1;⋯;Xn)∈Rn∗m，最大化方差：

max|P|=1 ∑i=1n(P⋅XTi−P⋅X¯¯¯¯¯¯¯¯)2/(n−1)

其中|P|=1是正交变换的要求，P就是从原坐标变换到新的第一个坐标的变换系数。如果我们提前为了便捷起见，将数据减去均值Xi:=Xi−X¯¯¯处理，这样就成了max|P|=1 ∑i=1n(P⋅XTi)2/(n−1)=∑i=1nP⋅(XTi⋅Xi)⋅PT/(n−1)=P⋅(∑i=1nXTi⋅Xi)⋅PT/(n−1)=P⋅A⋅PT 其中A=(∑ni=1XTi⋅Xi)/(n−1)，其实它就是(X(:,1),⋯,X(:,m))的协方差矩阵（X中一行代表一个样本此处是取各列），这个结果用定义就可以验证。于是我们的优化最终化为：

max|P|=1 PAPTA矩阵是上述的协方差矩阵。

这是一个很常规的条件极值问题，我们可以用Lagrange乘子法进行求解：

首先，由于|P|=1约束，我们有maxp∈D PAPT其中D={P|PPT=1}是一个有界闭集，因此根据连续性我们知道优化的最优解存在且最优值可达，记解为P1，并且P1APT1=infp∈D PAPT。

我们写出其拉格朗日函数L(P,λ)=PAPT+λ(1−PPT)对P1，应该有：∇L(P1,λ)=0得到2APT1−2λPT1=0APT1=λPT1

由于P1PT1=1，所以P1APT1=λ并且λ是矩阵A的特征值，P1是其对应的特征向量。而由于矩阵A是一个协方差矩阵，所以它是对称的，所以它的特征值是实数，不妨设全体特征值为λ1<⋯<λk，分别对应重述l1,⋯,lk，于是我们此时应该取最大的特征值λk。也就是P1APT1=infp∈DPAPT=λk。P1是λk对应的特征向量。

我们已经完成了前lk个主成分的选择（λk对应lk个相互正交的特征值），接下来便是在已经选择的主成分构成子空间外垂直的空间中进行优化（这是正交的要求）maxP∈E⊥λk |P|=1 PAPT同样开可以得到此时最优值为λk−1，以此类推。最终得到P1,⋯,Pm我们最终的变换就是Fj=∑i=1mPj(i)x(i) j=1,⋯,m也就是从坐标(x(1),⋯,x(m))变换为(F(1),⋯,F(m))

现在为止，我们已经回答了之前提出的三个问题，此时实际上我们已经得到了主成分分析的算法流程。

但是有一点需要注意，主成分分析对数据的尺度是敏感的，换句话说，我们如果对数据进行放缩会显著影响到我们最后的结果，因此，一般说法是，为了消除单位带来的差异，我们可以选择在对数据做预处理是不仅减去均值，还除以每个坐标上的标准差已消除单位带来的偏差：Xi(j):=Xi(j)−X¯¯¯(j)std(X(:,j))

MATLAB实例与测试

a=imread('1.jpg');
b=imread('2.jpg');
c=imread('3.jpg');
d=imread('4.jpg');
a=imresize(a,0.01);
b=imresize(b,0.01);
c=imresize(c,0.01);
d=imresize(d,0.01);
e=rgb2gray(a);
f=rgb2gray(b);
g=rgb2gray(c);
h=rgb2gray(d);
g1=e(:);
g2=f(:);
g3=g(:);
g4=h(:);
g1=double(g1);
g2=double(g2);
g3=double(g3);
g4=double(g4);
av=(g1+g2+g3+g4)/4;
G=[g1,g2,g3,g4]-repmat(av,1,4);
mat=cov(G');
[V,D]=eig(mat);
subplot(2,2,2);imshow(uint8(reshape(av+V(:,end)*    (g1-av)'*V(:,end),42,24)));
title('第一主成分')
subplot(2,2,3);imshow(uint8(reshape(av+V(:,end-1)*(g1-av)'*V(:,end-1),42,24)));
title('第二主成分')
subplot(2,2,4);imshow(uint8(reshape(av+V(:,end-2)*(g1-av)'*V(:,end-2),42,24)));
title('第三主成分')
subplot(2,2,1);imshow(e)
title('原图')

（你没看错这就是我的手2333）

上面第一个的是输入的图片，下面的是前三个主成分，可以看到这从生成的三幅图里仍然能看出原图的大体样子，他们的方差贡献分别为0.5836,0.2858和 0.1306。

事实上，因为这个实例中样本数据很少（4个），因此用PCA意义不大，仅仅作为一个例子展示。可以证明，对于落在m维空间中的n(n≤m)个点（减去均值），那么仅有n−1个主成分对应的λ（或方差）不为零，即用n−1个坐标就可以完全描述原来的数据，此时PCA可以说是没有什么意义了。，例如三维空间中的两个点，我们在做减去均值处理后，就可以选择一条空间直线穿过这两个点。我们面对的应当是在空间中大量的存在一定结构的数据。

主成分全部计入后就可以得到原来数据

subplot(2,1,2);imshow(uint8(reshape(av+V(:,end-1)*(g1-av)'*V(:,end-1)+V(:,end-2)*(g1-av)'*V(:,end-2)+V(:,end)*(g1-av)'*V(:,end),42,24)));
title('三个主成分的叠加')
subplot(2,1,1);imshow(e)
title('原图')

可以看到和原图是完全相同了。

我用的图片是压缩了100倍后得到的，所以图片不是很清晰，否则太高的清晰度我是等不及电脑去计算太大的矩阵的特征值和特征向量。

个人博客，请勿转载