转载于http://blog.csdn.net/guyuealian/article/details/68487833
网上关于PCA(主成分分析)原理和分析的博客很多,本博客并不打算长篇大论推论PCA理论,而是用最精简的语言说明鄙人对PCA的理解,并在最后给出用Matlab计算PCA过程的三种方法,方便大家对PCA的理解。

    源代码和附件下载地址: http://download.csdn.net/detail/guyuealian/9799160

一、PCA原理简要说明

PCA算法主要用于降维,就是将样本数据从高维空间投影到低维空间中,并尽可能的在低维空间中表示原始数据。PCA的几何意义可简单解释为:

0维-PCA:将所有样本信息都投影到一个点,因此无法反应样本之间的差异;要想用一个点来尽可能的表示所有样本数据,则这个点必定是样本的均值。
      1维-PCA:相当于将所有样本信息向样本均值的直线投影;
      2维-PCA:将样本的平面分布看作椭圆形分布,求出椭圆形的长短轴方向,然后将样本信息投影到这两条长短轴方向上,就是二维PCA。(投影方向就是平面上椭圆的长短轴方向);
      3维-PCA:样本的平面分布看作椭圆形分布,投影方法分别是椭圆球的赤道半径a和b,以及是极半径c(沿着z轴);

PCA简而言之就是根据输入数据的分布给输入数据重新找到更能描述这组数据的正交的坐标轴,比如下面一幅图,对于那个椭圆状的分布,最方便表示这个分布的坐标轴肯定是椭圆的长轴短轴而不是原来的x ,y轴。

那么如何求出这个长轴和短轴呢?于是线性代数就来了:我们需要先求出这堆样本数据的协方差矩阵,然后再求出这个协方差矩阵的特征值和特征向量,对应最大特征值的那个特征向量的方向就是长轴(也就是主元)的方向,次大特征值的就是第二主元的方向,以此类推。

二、PCA实现方法

      实现PCA的方法, 可【1】直接调用Matlab工具箱princomp( )函数实现,也可【2】 自己实现PCA的过程,当然也可以【3】使用快速PCA算法的方法。

(1)方法一:[COEFF SCORE latent]=princomp(X)

参数说明:
    1)COEFF 是主成分分量,即样本协方差矩阵的特征向量;
    2)SCORE主成分,是样本X在低维空间的表示形式,即样本X在主成份分量COEFF上的投影 ,若需要降k维,则只需要取前k列主成分分量即可
    3)latent:一个包含样本协方差矩阵特征值的向量;

      实例:假设有8个样本,每个样本有4个特征(属性),使用PCA方法实现降维(k维,k小于特征个数4),并提取前2个主成份的特征,即将原始数据从4维空间降维到2维空间。
  1. %% 样本矩阵X,有8个样本,每个样本有4个特征,使用PCA降维提取k个主要特征(k<4)
  2. k=2;                            %将样本降到k维参数设置
  3. X=[1 2 1 1;                     %样本矩阵
  4. 3 3 1 2;
  5. 3 5 4 3;
  6. 5 4 5 4;
  7. 5 6 1 5;
  8. 6 5 2 6;
  9. 8 7 1 2;
  10. 9 8 3 7]
  11. %% 使用Matlab工具箱princomp函数实现PCA
  12. [COEFF SCORE latent]=princomp(X)
  13. pcaData1=SCORE(:,1:k)            %取前k个主成分
%% 样本矩阵X,有8个样本,每个样本有4个特征,使用PCA降维提取k个主要特征(k<4)
k=2; %将样本降到k维参数设置
X=[1 2 1 1; %样本矩阵
3 3 1 2;
3 5 4 3;
5 4 5 4;
5 6 1 5;
6 5 2 6;
8 7 1 2;
9 8 3 7]
%% 使用Matlab工具箱princomp函数实现PCA
[COEFF SCORE latent]=princomp(X)
pcaData1=SCORE(:,1:k) %取前k个主成分

运行结果:

  1. X =
  2. 1     2     1     1
  3. 3     3     1     2
  4. 3     5     4     3
  5. 5     4     5     4
  6. 5     6     1     5
  7. 6     5     2     6
  8. 8     7     1     2
  9. 9     8     3     7
  10. COEFF =
  11. 0.7084   -0.2826   -0.2766   -0.5846
  12. 0.5157   -0.2114   -0.1776    0.8111
  13. 0.0894    0.7882   -0.6086    0.0153
  14. 0.4735    0.5041    0.7222   -0.0116
  15. SCORE =
  16. -5.7947   -0.6071    0.4140   -0.0823
  17. -3.3886   -0.8795    0.4054   -0.4519
  18. -1.6155    1.5665   -1.0535    1.2047
  19. -0.1513    2.5051   -1.3157   -0.7718
  20. 0.9958   -0.5665    1.4859    0.7775
  21. 1.7515    0.6546    1.5004   -0.6144
  22. 2.2162   -3.1381   -1.6879   -0.1305
  23. 5.9867    0.4650    0.2514    0.0689
  24. latent =
  25. 13.2151
  26. 2.9550
  27. 1.5069
  28. 0.4660
  29. pcaData1 =
  30. -5.7947   -0.6071
  31. -3.3886   -0.8795
  32. -1.6155    1.5665
  33. -0.1513    2.5051
  34. 0.9958   -0.5665
  35. 1.7515    0.6546
  36. 2.2162   -3.1381
  37. 5.9867    0.4650
X =
1 2 1 1
3 3 1 2
3 5 4 3
5 4 5 4
5 6 1 5
6 5 2 6
8 7 1 2
9 8 3 7
COEFF =
0.7084 -0.2826 -0.2766 -0.5846
0.5157 -0.2114 -0.1776 0.8111
0.0894 0.7882 -0.6086 0.0153
0.4735 0.5041 0.7222 -0.0116
SCORE =
-5.7947 -0.6071 0.4140 -0.0823
-3.3886 -0.8795 0.4054 -0.4519
-1.6155 1.5665 -1.0535 1.2047
-0.1513 2.5051 -1.3157 -0.7718
0.9958 -0.5665 1.4859 0.7775
1.7515 0.6546 1.5004 -0.6144
2.2162 -3.1381 -1.6879 -0.1305
5.9867 0.4650 0.2514 0.0689
latent =
13.2151
2.9550
1.5069
0.4660
pcaData1 =
-5.7947 -0.6071
-3.3886 -0.8795
-1.6155 1.5665
-0.1513 2.5051
0.9958 -0.5665
1.7515 0.6546
2.2162 -3.1381
5.9867 0.4650

(2)方法二:自己编程实现

PCA的算法过程,用一句话来说,就是“将所有样本X减去样本均值m,再乘以样本的协方差矩阵C的特征向量V,即为PCA主成分分析”,其计算过程如下:
    [1].将原始数据按行组成m行n列样本矩阵X(每行一个样本,每列为一维特征)
    [2].求出样本X的协方差矩阵C和样本均值m;(Matlab可使用cov()函数求样本的协方差矩阵C,均值用mean函数)
    [3].求出协方差矩阵的特征值D及对应的特征向量V;(Matlab可使用eigs()函数求矩阵的特征值D和特征向量V)
    [4].将特征向量按对应特征值大小从上到下按行排列成矩阵,取前k行组成矩阵P;(eigs()返回特征值构成的向量本身就是从大到小排序的)
    [5].Y=(X-m)×P即为降维到k维后的数据;

    PS:关于协方差矩阵,很多人有点郁闷,有些教程用协方差矩阵,而有些资料是用散步矩阵计算,其实协方差矩阵和散步矩阵就是一个倍数关系:协方差矩阵C×(n-1)=散步矩阵S。我们可以使用Matlab工具箱的cov函数求协方差矩阵C。
  1. %% 自己实现PCA的方法
  2. [Row Col]=size(X);
  3. covX=cov(X);                                    %求样本的协方差矩阵(散步矩阵除以(n-1)即为协方差矩阵)
  4. [V D]=eigs(covX);                               %求协方差矩阵的特征值D和特征向量V
  5. meanX=mean(X);                                  %样本均值m
  6. %所有样本X减去样本均值m,再乘以协方差矩阵(散步矩阵)的特征向量V,即为样本的主成份SCORE
  7. tempX= repmat(meanX,Row,1);
  8. SCORE2=(X-tempX)*V                              %主成份:SCORE
  9. pcaData2=SCORE2(:,1:k)
 %% 自己实现PCA的方法
[Row Col]=size(X);
covX=cov(X); %求样本的协方差矩阵(散步矩阵除以(n-1)即为协方差矩阵)
[V D]=eigs(covX); %求协方差矩阵的特征值D和特征向量V
meanX=mean(X); %样本均值m
%所有样本X减去样本均值m,再乘以协方差矩阵(散步矩阵)的特征向量V,即为样本的主成份SCORE
tempX= repmat(meanX,Row,1);
SCORE2=(X-tempX)*V %主成份:SCORE
pcaData2=SCORE2(:,1:k)
运行结果:
  1. SCORE2 =
  2. -5.7947    0.6071   -0.4140    0.0823
  3. -3.3886    0.8795   -0.4054    0.4519
  4. -1.6155   -1.5665    1.0535   -1.2047
  5. -0.1513   -2.5051    1.3157    0.7718
  6. 0.9958    0.5665   -1.4859   -0.7775
  7. 1.7515   -0.6546   -1.5004    0.6144
  8. 2.2162    3.1381    1.6879    0.1305
  9. 5.9867   -0.4650   -0.2514   -0.0689
  10. pcaData2 =
  11. -5.7947    0.6071
  12. -3.3886    0.8795
  13. -1.6155   -1.5665
  14. -0.1513   -2.5051
  15. 0.9958    0.5665
  16. 1.7515   -0.6546
  17. 2.2162    3.1381
  18. 5.9867   -0.4650
SCORE2 =
-5.7947 0.6071 -0.4140 0.0823
-3.3886 0.8795 -0.4054 0.4519
-1.6155 -1.5665 1.0535 -1.2047
-0.1513 -2.5051 1.3157 0.7718
0.9958 0.5665 -1.4859 -0.7775
1.7515 -0.6546 -1.5004 0.6144
2.2162 3.1381 1.6879 0.1305
5.9867 -0.4650 -0.2514 -0.0689 pcaData2 =
-5.7947 0.6071
-3.3886 0.8795
-1.6155 -1.5665
-0.1513 -2.5051
0.9958 0.5665
1.7515 -0.6546
2.2162 3.1381
5.9867 -0.4650

对比方法一和方法可知,主成份分量SCORE和SCORE2的绝对值是一样的(符号只是相反方向投影而已,不影响分析结果),其中pcaData是从SCORE提取前2列的数据,这pcaData就是PCA从4维降到2维空间的数据表示形式,pcaData可以理解为:通过PCA降维,每个样本可以用2个特征来表示原来4个特征了。

(3)方法三:使用快速PCA方法

     PCA的计算中最主要的工作量是计算样本协方差矩阵的本征值和本征向量。假设样本矩阵X的大小为n ×d (nd 维样本特征向量),则样本散布矩阵(协方差矩阵) S 将是一个d×d的方阵,故当维数d较大时计算复杂度会非常高。例如当维数d=10000,S是一个10 000 ×10 000的矩阵,此时如果采用上面的princomp函数计算主成份,Matlab通常会出现内存耗尽的问题(out of memory), 即使有足够多的内存,要得到S的全部本征值可能也要花费数小时的时间。
    快速PCA的方法相关理论,可以参考张铮的《精通Matlab数字图像处理与识别》第12章(P307),PDF可在附件下载http://download.csdn.net/detail/guyuealian/9799160

fastPCA函数用来对样本矩阵A进行快速主成分分析和降维(降至k维),其输出pcaA为维后的k维样本特征向量组成的矩阵,每行一个样本,列数k为降维后的样本特征维数,相当于princomp函数中的输出SCORE, 而输出V为主成分分量,相当于princomp函数中的输出COEFF。

  1. %% 使用快速PCA算法实现的方法
  2. [pcaData3 COEFF3] = fastPCA(X, k )
%% 使用快速PCA算法实现的方法
[pcaData3 COEFF3] = fastPCA(X, k )

其中fastPCA函数的代码实现如下:

  1. function [pcaA V] = fastPCA( A, k )
  2. % 快速PCA
  3. % 输入:A --- 样本矩阵,每行为一个样本
  4. %      k --- 降维至 k 维
  5. % 输出:pcaA --- 降维后的 k 维样本特征向量组成的矩阵,每行一个样本,列数 k 为降维后的样本特征维数
  6. %      V --- 主成分向量
  7. [r c] = size(A);
  8. % 样本均值
  9. meanVec = mean(A);
  10. % 计算协方差矩阵的转置 covMatT
  11. Z = (A-repmat(meanVec, r, 1));
  12. covMatT = Z * Z';
  13. % 计算 covMatT 的前 k 个本征值和本征向量
  14. [V D] = eigs(covMatT, k);
  15. % 得到协方差矩阵 (covMatT)' 的本征向量
  16. V = Z' * V;
  17. % 本征向量归一化为单位本征向量
  18. for i=1:k
  19. V(:,i)=V(:,i)/norm(V(:,i));
  20. end
  21. % 线性变换(投影)降维至 k 维
  22. pcaA = Z * V;
  23. % 保存变换矩阵 V 和变换原点 meanVec
function [pcaA V] = fastPCA( A, k )
% 快速PCA
% 输入:A --- 样本矩阵,每行为一个样本
% k --- 降维至 k 维
% 输出:pcaA --- 降维后的 k 维样本特征向量组成的矩阵,每行一个样本,列数 k 为降维后的样本特征维数
% V --- 主成分向量
[r c] = size(A);
% 样本均值
meanVec = mean(A);
% 计算协方差矩阵的转置 covMatT
Z = (A-repmat(meanVec, r, 1));
covMatT = Z * Z';
% 计算 covMatT 的前 k 个本征值和本征向量
[V D] = eigs(covMatT, k);
% 得到协方差矩阵 (covMatT)' 的本征向量
V = Z' * V;
% 本征向量归一化为单位本征向量
for i=1:k
V(:,i)=V(:,i)/norm(V(:,i));
end
% 线性变换(投影)降维至 k 维
pcaA = Z * V;
% 保存变换矩阵 V 和变换原点 meanVec

运行结果为:

  1. pcaData3 =
  2. -5.7947   -0.6071
  3. -3.3886   -0.8795
  4. -1.6155    1.5665
  5. -0.1513    2.5051
  6. 0.9958   -0.5665
  7. 1.7515    0.6546
  8. 2.2162   -3.1381
  9. 5.9867    0.4650
  10. COEFF3 =
  11. 0.7084   -0.2826
  12. 0.5157   -0.2114
  13. 0.0894    0.7882
  14. 0.4735    0.5041

最新文章

  1. Qt5.7中使用MySQL Driver
  2. 解决Myeclipse下Debug出现Source not found以及sql server中导入数据报错
  3. Java基础-JVM
  4. 验证码点击刷新 this.src=this.src+&#39;?&#39;+Math.random()
  5. 转载:windows的mysql提权方式
  6. SqlServer扩展存储过程
  7. layer.js定制弹窗
  8. Android 代码监控apk安装,卸载,替换
  9. (转)CentOS5.5 下搭建 PHP 环境(最佳的LAMP环境)
  10. php根据经纬度计算距离和方向--摘录自http://haotushu.sinaapp.com/post-520.html
  11. C#总结项目《影院售票系统》编写总结二
  12. HDU 3729 I&amp;#39;m Telling the Truth(二部图最大匹配+结果输出)
  13. ConcurrentHashMap源码分析(一)
  14. hive上mysql元数据库配置
  15. Java输出流FileOutputStream使用详解
  16. VS2017 ASP.NET MVC 5.0 开部署问题汇总
  17. vue 不常见操作
  18. virtualenv 的使用
  19. 浅谈移动端 View 的显示过程
  20. 《JavaScript》字符转义

热门文章

  1. @Transactional注解真的有必要声明rollbackFor属性吗?
  2. 使用request对象进行简单的注册以及信息显示
  3. 文本转语音TTS(文本阅读和视频配音制作)MP3
  4. 知识图谱顶会论文(SIGIR-2022) MorsE:归纳知识图嵌入的元知识迁移
  5. 题解 P3395 路障
  6. 说说switch关键字
  7. 【题解】CF45I TCMCF+++
  8. 倍福Ads协议通信测试
  9. hyperworks2021位安装教程
  10. kotlin的suspend对比csharp的async&amp;await