PCA matlab实现
2014-03-10 18:28
190 查看
PCA 流程如下:
1、去均值 2、计算协方差矩阵 3、计算协方差特征值和特征向量 4、降序排列特征值选取较大的特征值,选择相应的特征值和特征向量
以下按照步骤编写matlab代码。
Matlab函数mean可得:如下
协方差定义:
具体求解:
*注意分母为(n-1)而不是n,因为这样定义的协方差方差是总体方差的无偏估计(具体可见:
http://en.wikipedia.org/wiki/Unbiased_estimator#Sample_variancehttp://www.zhihu.com/question/20099757)
协方差矩阵如下:
其元素aij表示变量i,j之间的协方差cov(i,j);
(关于协方差矩阵的含义,可见blog:/article/7692386.html)
其中Xj为去中心化之后的特征向量。
去中心化后,
可表示为如下:
求和之后,可以得到协方差矩阵。
B.
同理,上面的向量相加我们可以直接用矩阵相乘的形式得到。
设X为去中心后的特征矩阵,那么
C.
直接用matlab自带的协方差的函数cov()计算,不过注意cov按行计算,实际运用时要转置。
我们使用矩阵形式,得到如下代码:
根据PCA的原理,我们需要寻找使协方差矩阵对角化的变换矩阵
(可见blog:/article/7692386.html)。
一个方阵可以写成如下形式:
其中Q为其特征向量组成的矩阵,
为其特征值组成的对角矩阵,转化一下式子,得到:
于是,我们现在只要得到协方差矩阵的归一化特征向量,组成转化矩阵Q即可。
使用matlab自带的函数eig(),
代码:
SVD(奇异值分解):
A.奇异值
设A为m*n阶实矩阵,则存在m阶正交阵U和n阶正交阵V,使得
A = U*S*V’
其中S=diag(σi,σ2,……,σr),σi>0 (i=1,…,r),r=rank(A)。
其中:
对任意矩阵A,它的奇异值就是AA'或A'A的非零特征值的开方(它们有相同的非零特征值),这些特征值都是正数。
U 为AA'单位特征向量矩阵。
V为A’A单位特征向量矩阵。
B.奇异值与特征值的联系
奇异值有类似于特征值的性质,当矩阵为共轭对称矩阵时,特征值=奇异值。不过一般情况是不相同的。
如果把矩阵看做一个线性变换,那么特征值表征了其特征向量方向的能量大小。根据定义我们可以看出,奇异值也有类似的性质
C.奇异值分解与PCA的关系
通过变换,我们可以得到:
也就是,已知A,V,我们可以求得U。
如上,如果AA’维数过大,计算机不好求解其特征向量U,那么我们可以转而求A’*A的特征向量V。
求解PCA的过程中,对于小样本问题,样本维数M>>样本个数N,那么X*X’得到的协方差矩阵为M*M,不好特征分解。如果我们根据SVD的原理,解X’*X(N*N)的特征向量,最后再变化,也能达到同样的目的。
(SVD具体可见:
/article/5178349.html
http://szshdy.blog.163.com/blog/static/1322012512010511156587/)
通过奇异值分解,得到协方差矩阵特征向量,代码如下:
最终,整合上述的代码,得到如下完整的PCA代码:
1、去均值 2、计算协方差矩阵 3、计算协方差特征值和特征向量 4、降序排列特征值选取较大的特征值,选择相应的特征值和特征向量
以下按照步骤编写matlab代码。
1.去均值
Matlab函数mean可得:如下Mean_Image=mean(Train_SET,2); Train_SET=Train_SET-Mean_Image*ones(1,Train_NUM);
2.计算协方差矩阵
协方差定义:
具体求解:
*注意分母为(n-1)而不是n,因为这样定义的协方差方差是总体方差的无偏估计(具体可见:http://en.wikipedia.org/wiki/Unbiased_estimator#Sample_variancehttp://www.zhihu.com/question/20099757)
协方差矩阵如下:
其元素aij表示变量i,j之间的协方差cov(i,j);
(关于协方差矩阵的含义,可见blog:/article/7692386.html)
计算方法:
A.其中Xj为去中心化之后的特征向量。
去中心化后,
可表示为如下:
求和之后,可以得到协方差矩阵。
B.
同理,上面的向量相加我们可以直接用矩阵相乘的形式得到。
设X为去中心后的特征矩阵,那么
C.
直接用matlab自带的协方差的函数cov()计算,不过注意cov按行计算,实际运用时要转置。
我们使用矩阵形式,得到如下代码:
R=Train_SET*Train_SET'/(Train_NUM-1);
3.计算特征值与特征向量
根据PCA的原理,我们需要寻找使协方差矩阵对角化的变换矩阵(可见blog:/article/7692386.html)。
一个方阵可以写成如下形式:
其中Q为其特征向量组成的矩阵,
为其特征值组成的对角矩阵,转化一下式子,得到:
于是,我们现在只要得到协方差矩阵的归一化特征向量,组成转化矩阵Q即可。
使用matlab自带的函数eig(),
代码:
[V,S]=eig(R);
**小样本问题:
上面的代码存在一个问题,就是常见的小样本问题。样本维数>>样本个数,这样得到的协方差矩阵很大,直接求解时间复杂度过高。于是我们通过另一种方式来求解。SVD(奇异值分解):
A.奇异值
设A为m*n阶实矩阵,则存在m阶正交阵U和n阶正交阵V,使得
A = U*S*V’
其中S=diag(σi,σ2,……,σr),σi>0 (i=1,…,r),r=rank(A)。
其中:
对任意矩阵A,它的奇异值就是AA'或A'A的非零特征值的开方(它们有相同的非零特征值),这些特征值都是正数。
U 为AA'单位特征向量矩阵。
V为A’A单位特征向量矩阵。
B.奇异值与特征值的联系
奇异值有类似于特征值的性质,当矩阵为共轭对称矩阵时,特征值=奇异值。不过一般情况是不相同的。
如果把矩阵看做一个线性变换,那么特征值表征了其特征向量方向的能量大小。根据定义我们可以看出,奇异值也有类似的性质
C.奇异值分解与PCA的关系
通过变换,我们可以得到:
也就是,已知A,V,我们可以求得U。
如上,如果AA’维数过大,计算机不好求解其特征向量U,那么我们可以转而求A’*A的特征向量V。
求解PCA的过程中,对于小样本问题,样本维数M>>样本个数N,那么X*X’得到的协方差矩阵为M*M,不好特征分解。如果我们根据SVD的原理,解X’*X(N*N)的特征向量,最后再变化,也能达到同样的目的。
(SVD具体可见:
/article/5178349.html
http://szshdy.blog.163.com/blog/static/1322012512010511156587/)
通过奇异值分解,得到协方差矩阵特征向量,代码如下:
R=Train_SET'*Train_SET/(Train_NUM-1); [V,S]=Find_K_Max_Eigen(R,Eigen_NUM); disc_value=S; disc_set=zeros(NN,Eigen_NUM); Train_SET=Train_SET/sqrt(Train_NUM-1); for k=1:Eigen_NUM disc_set(:,k)=(1/sqrt(disc_value(k)))*Train_SET*V(:,k); end
4.完整代码
最终,整合上述的代码,得到如下完整的PCA代码:function [disc_set,disc_value,Mean_Image]=Eigenface_f(Train_SET,Eigen_NUM)
[NN,Train_NUM]=size(Train_SET);
if NN<=Train_NUM
Mean_Image=mean(Train_SET,2); Train_SET=Train_SET-Mean_Image*ones(1,Train_NUM);
R=Train_SET*Train_SET'/(Train_NUM-1);
[V,S]=Find_K_Max_Eigen(R,Eigen_NUM);
disc_value=S;
disc_set=V;
else % 小样本问题,svd
Mean_Image=mean(Train_SET,2); Train_SET=Train_SET-Mean_Image*ones(1,Train_NUM);
R=Train_SET'*Train_SET/(Train_NUM-1); [V,S]=Find_K_Max_Eigen(R,Eigen_NUM); disc_value=S; disc_set=zeros(NN,Eigen_NUM); Train_SET=Train_SET/sqrt(Train_NUM-1); for k=1:Eigen_NUM disc_set(:,k)=(1/sqrt(disc_value(k)))*Train_SET*V(:,k); end
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function [Eigen_Vector,Eigen_Value]=Find_K_Max_Eigen(Matrix,Eigen_NUM)
[NN,NN]=size(Matrix);
[V,S]=eig(Matrix); %Note this is equivalent to; [V,S]=eig(St,SL); also equivalent to [V,S]=eig(Sn,St); %
S=diag(S);
[S,index]=sort(S);
Eigen_Vector=zeros(NN,Eigen_NUM);
Eigen_Value=zeros(1,Eigen_NUM);
p=NN;
for t=1:Eigen_NUM
Eigen_Vector(:,t)=V(:,index(p));
Eigen_Value(t)=S(p);
p=p-1;
end
相关文章推荐
- PCA分析以及MATLAB实现
- 主成分分析(PCA)原理与故障诊断(SPE、T^2以及结合二者的综合指标)-MATLAB实现
- PCA(主成分)分析及MATLAB实现
- 主成分分析PCA的matlab实现
- PCA的MATLAB实现
- PCA-MATLAB 实现
- PCA降维算法总结以及matlab实现PCA(个人的一点理解)
- PCA算法学习_2(PCA理论的matlab实现)
- 基于PCA的人脸识别_Matlab实现(个人研读之后的一些总结)
- PCA降维算法总结以及matlab实现PCA(个人的一点理解)
- 运用PCA(主成分分析法)进行人脸识别的MATLAB 代码实现
- [置顶] PCA原理分析和Matlab实现方法(三)
- PCA算法学习_2(PCA理论的matlab实现)
- PCA (主成分分析)matlab实现
- PCA降维算法总结以及matlab实现PCA(个人的一点理解)
- PCA降维算法总结以及matlab实现PCA(个人的一点理解)
- PCA降维算法总结以及matlab实现PCA(个人的一点理解)
- PCA降维算法总结以及matlab实现PCA
- matlab实现基于PCA的人脸识别算法
- 【模式识别】PCA检测人脸的简单示例MATLAB实现