主成分分析（PCA）算法模型实现及分析（源码在文章后附录）

• 1 引言
• 2 关于PCA原理和算法实现
• 2.1 PCA基本原理
• • 2.2 协方差计算
  • 2.3 PCA实现步骤
  • •   （1）PCA算法实现步骤
    •   （2）基于特征值分解协方差矩阵实现PCA算法
    • （3）基于SVD分解协方差矩阵实现PCA
  • 2.4 简单的总结一下
• MATLAB代码附页：

Author： Nirvana Of Phoenixl
Proverbs for you：There is no doubt that good things will always come, and when it comes late, it can be a surprise.

1 引言

  主成分分析（PCA）是一种能够极大提升无监督特征学习速度的数据降维算法。在许多领域研究与应用当中，通常需要对含有多个变量的数据进行观测并分析规律，而许多变量之间可能存在相关性，从而增加了问题分析的复杂性。如果分别对每个指标进行分析，分析往往是孤立的，不能完全利用数据中的信息，因此盲目减少指标会损失很多有用的信息，从而产生错误的结论。
  因而需要减少分析指标的同时，确保指标包含信息的相对完整，以达到对所收集数据进行全面分析的目的。由于各变量之间存在一定的相关关系，因此可以考虑将关系紧密的变量变成尽可能少的新变量，使这些新变量是两两不相关的，那么就可以用较少的综合指标分别代表存在于各个变量中的各类信息。
  主成分分析（Principal Component Analysis，PCA）的方法，可以将具有多个观测变量的高维数据集降维，使人们可以从事物之间错综复杂的关系中找出一些主要的方面，从而能更加有效地利用大量统计数据进行定量分析，并可以更好地进行可视化、回归等后续处理。

2 关于PCA原理和算法实现

2.1 PCA基本原理

  PCA主要思想：是将n维特征映射到k维上，即在原有n维特征的基础上重新构造出来的k维特征，此过程产生的新的正交特征成为主成分。
  PCA的工作：就是从原始的空间中顺序地找一组相互正交的坐标轴，新的坐标轴的选择与数据本身是密切相关的。其中，（1）第一个新坐标轴选择是原始数据中方差最大的方向，（2）第二个新坐标轴选取是与第一个坐标轴正交的平面中使得方差最大的，（3）第三个轴是与第1，2个轴正交的平面中方差最大的。依次类推，可以得到n个这样的坐标轴ID。通过这种方式获得的新的坐标轴，我们发现，大部分方差都包含在前面k个坐标轴中，后面的坐标轴所含的方差几乎为0。于是，我们可以忽略余下的坐标轴，只保留前面k个含有绝大部分方差的坐标轴。事实上，这相当于只保留包含绝大部分方差的维度特征，而忽略包含方差几乎为0的特征维度，实现对数据特征的降维处理。

2.2 协方差计算

  通过计算数据矩阵的协方差矩阵，然后得到协方差矩阵的特征值特征向量，选择特征值最大(即方差最大)的k个特征所对应的特征向量组成的矩阵，即将数据矩阵转换到新的空间当中，实现数据特征的降维。如式（2-1）到（2-3）所示实现样本协方差的计算。
  

（1）计算样本均值：

其中x表示样本，n表示样本总数。
  
（2）样本方差：

（3）样本X和样本Y的协方差：

  其中cov（X,Y）表示两样本的协方差，用记为E。
  方差的计算公式是针对一维特征，即就是针对同一特征不同样本的取值计算得到的；而协方差至少要满足二维特征，实际上方差是协方差的特殊情况。

  由于得到协方差矩阵的特征值特征向量有两种方法：特征值分解协方差矩阵、奇异值分解协方差矩阵，所以PCA算法有两种实现方法：基于特征值分解协方差矩阵实现PCA算法、基于SVD分解协方差矩阵实现PCA算法。本文选择MATLAB实现PCA算法。

2.3 PCA实现步骤

  （1）PCA算法实现步骤

  选定MATLAB作为实验平台实现主成分分析达到降维目的，采用SVD算法是该平台默认的PCA实现算法。在完成明确2.2节的基本知识后，下面就可以确定具体实现降维的操作。
  第一步：原始样本数据获取。
  第二步：分别取求每个特征的平均值，针对所有样本，减去相应的均值。
  第三步：求解协方差矩阵，如2-1所示步骤。
  第四步: 奇异值分解，求取协方差矩阵的特征值和特征向量。
  第五步: 倒序排列特征值（从大到小），选择最大特征值作为主成分，成为新的样本。
  第六步：将特征值最大的d个向量作为投影向量，构成D*d维的投影矩阵W，对于任意维样本，将其投影选取的特征向量（主成分方向）上。

  （2）基于特征值分解协方差矩阵实现PCA算法

在MATLAB中需要载入其自带的数据集fisheriris，该数据集总共统计了三种鸢尾花的花萼长、花萼宽、花瓣长和花瓣宽，然后进行中心化处理，并计算协方差矩阵，如图2.1所示。

            图2.1导入数据集并中心化处理

  利用特征值分解法：使用eig函数，如图2.2所示，实现主成分分析，主要包括特征值矩阵的提取、按升序排列特征值等。

              图2.2特征分解及协方差计算

  通过计算方差的累积贡献率，如式2-4所示，结合数据模型可以实现数据的降维，基于方差贡献率可以确定最终降维的维数，一般来说是通过对数据的观察来确定主元个数，而利用此方法可以简单的确定PCA主成分分析中的主元个数。

其中， 是递增的，因此f (k)为单调递增的函数，且递增速度随着k增加逐渐降低。

  一般来说，取f（r）≥ 某一阈值（如95%）的最小的r，这样最多损失5%的方差，不同算法取的 不同或者替换为 即可，通过画作图实现，如图2.3所示。

                图2.3方差贡献率确定维数

  实现效果如2.4所示，可以判断降维数，用以分析。

                图2.4方差贡献率

  结合MATLAB确定的降维维数，依据数据集完成降维，如图2.5所示。

  通过上述中心化、计算协方差和特征分解并进行作图实现降维处理，如图2.6实现。

                 图2.6 PCA降维

（3）基于SVD分解协方差矩阵实现PCA

  实际上MATLAB中本身就集成了SVD算法，在一般情况下，MATLAB中的PCA算法也会使用SVD，所以在也可以通过奇异值分解（SVD）实现。同样地，也可以利用pca函数实现，两者的调用格式如图2.7所示，其作图部分与（2）中所示作图一致。

                图2.7 不同PCA实现

2.4 简单的总结一下

  PCA主成分分析算法，通过对高维数据的降维处理，将原本数据集使用合理的方法采用主成分代替，也就是利用主成分来代表高维数据尽可能的不失真的表示原本数据集。

  通过上面鸢尾花的例子，降维后的数据仍然可以清晰地分为三类。当需要确定一种鸢尾花是，计算相应的T1和T2主成分得分（Principal component scores），即为新空间中的数据点，将其结果画在散点图中，就可以判断出其属于哪一种鸢尾花，同样的道理应用到更多的场合也可以实现。

  PCA可以应用到很多场合比如聚类分析，然后应用到的电商场合的推送；图像的压缩和人脸检测与匹配等。

MATLAB代码附页：

load fisheriris;              %导入数据集
X = meas;                 % n = 150, m = 4
meanX = ones(size(X,1), 1) * mean(X);  % 中心化处理

centredX = X - meanX;

C = cov(centredX);	        % 直接调用cov直接计算协方差矩阵即可

[W, Lambda] = eig(C);       % W是特征向量组成的矩阵（4×4），Lambda是特征值组成的对角矩阵
ev = (diag(Lambda))';		% 提取特征值
ev = ev(:, end:-1:1);		    % eig计算出的特征值是升序的，这里手动倒序（W同理）
W = W(:, end:-1:1);
sum(W.*W, 1)             % 可以验证每个特征向量各元素的平方和均为

Wr = W(:, 1:2);             % 提取前两个主成分的特征向量
Tr = centredX * Wr;         %  新坐标空间的数据点
% 作图
figure;
    stairs(cumsum(ev)/sum(ev), 'LineWidth',1.5);
    axis([1 4 0 1]);
    xlabel('$ k $', 'Interpreter', 'latex');
    ylabel('$ f(k)=\frac{\sum _{i=1}^i \lambda_k}{\sum_{i=1}^m \lambda_i} $',...
        'Interpreter', 'latex');
    hold on;
    plot([1 4], [0.95 0.95], '--');  % 从图中可以看出，r为方差贡献率，取r = 2
figure;
    scatter(Tr(:,1), Tr(:,2), 130, categorical(species), '.');
    colormap(winter);
    xlabel('Principal Component 1');
    ylabel('Principal Component 2');
[U, Sigma, V] = svd(X);              % 可以检验，W和V完全相同（向量的正负号不影响）
Vr = V(:, 1:2);                      % 提取前两个主成分的特征向量
Tr = X * Vr;                        % 新坐标空间的数据点
% 画图部分同上
[loadings, scores] = pca(X, 'NumComponents', r);
[Wr, Tr, ev] = pca(X, 'NumComponents',2);   % 画图部分