您的位置：首页 > 编程语言 > MATLAB

主成分分析（PCA）及其MATLAB的实现方法

2020-03-06 14:58 2511 查看

文章目录

  • 原理与步骤简述
  • MATLAB的实现方法
  • 应用

    • 说明：下文中，粗斜体字母均表示矩阵（如A\boldsymbol AA）；为不引起歧义，列向量也均加箭头表示（如a⃗\vec aa）

    概述

    PCA的目的

    假设现在有这样一个情景：现在要统计并可视化分析男大学生体测成绩，如果只参考立定跳远和1000m成绩两项指标，我们可以以立定跳远成绩作为xxx轴，1000m成绩作为yyy轴做出散点图，每个点代表一个学生；若统计三项指标，我们也可以在三维空间中做出散点图；但如果要统计四项乃至更多的指标，我们就无法再以此方法进行数据的可视化。

    而主成分分析（Principal Component Analysis，PCA）的方法，可以将具有多个观测变量的高维数据集降维，使人们可以从事物之间错综复杂的关系中找出一些主要的方面，从而能更加有效地利用大量统计数据进行定量分析，并可以更好地进行可视化、回归等后续处理。

    PCA的几何意义

    先将问题简化为二维情形。有NNN个样品，具有两个观测变量X1, X2X_1,\ X_2X1​, X2​，做出散点图（如下图中的蓝色点），这样，在由X1, X2X_1,\ X_2X1​, X2​组成的坐标空间中，NNN个样品的分布情况如带状。现在问：如果现在要将两个观测变量缩减为一个，应该如何选取？

    可以直观地看出，这NNN个样品无论沿X1X_1X1​轴还是沿X2X_2X2​轴方向，均有较大的离散性（其离散程度可以分别用观测变量X1X_1X1​的方差和X2X_2X2​的方差定量表示），也就是说，只考虑X1, X2X_1,\ X_2X1​, X2​的其中一个，原始数据均会有较大损失。

    现在考虑以下线性组合，变换坐标空间：
    {T1=X1cos⁡θ+X2sin⁡θT2=−X1sin⁡θ+X2cos⁡θ \begin{cases} T_1 = X_1 \cos \theta + X_2 \sin \theta \\ T_2 = -X_1 \sin \theta + X_2 \cos \theta \end{cases} {T1​=X1​cosθ+X2​sinθT2​=−X1​sinθ+X2​cosθ​

    即

    [T1T2]=[X1X2][cos⁡θ−sin⁡θsin⁡θcos⁡θ]=[X1X2]W(1) \begin{bmatrix} T_1 & T_2 \end{bmatrix}=\begin{bmatrix} X_1 & X_2 \end{bmatrix} \begin{bmatrix} \cos \theta & -\sin \theta \\ \sin \theta & \cos \theta \end{bmatrix} = \begin{bmatrix} X_1 & X_2 \end{bmatrix} \boldsymbol W \tag{1} [T1​​T2​​]=[X1​​X2​​][cosθsinθ​−sinθcosθ​]=[X1​​X2​​]W(1)

    式中W\boldsymbol WW为旋转变换矩阵，有如下性质：

    • W\boldsymbol WW的第iii列组成的列向量（也就是后面将要提到的特征向量），代表的就是新坐标空间的基底TiT_iTi​在原坐标空间中的坐标，且其为单位向量；
    • W\boldsymbol WW为正交阵，即满足：WTW=I\boldsymbol W^{\rm T} \boldsymbol W=\boldsymbol IWTW=I。

    经过旋转，NNN个数据点在T1T_1T1​轴上的离散程度最大，因而变量T1T_1T1​代表了原始数据的绝大部分信息，这样，即使不考虑变量T2T_2T2​也不会损失太多数据信息。这个T1T_1T1​即为第一主成分（Principal Component 1，PC1），如图中箭头所示。若将所有数据点投影到T1T_1T1​轴上（图中橙色点），就得到了降维后的数据。若有多个主成分，则：

    • 这些主成分之间相互独立，即没有重叠的信息，亦即这些特征向量之间正交，cov(Ti,Tj)=0,i≠j{\rm cov} \left( T_i, T_j \right) = 0,\quad i \ne jcov(Ti​,Tj​)=0,i​=j；
    • 主成分的方差依次递减，var(T1)≥var(T2)≥⋯{\rm var}\left( T_1 \right) \ge {\rm var}\left( T_2 \right) \ge \cdotsvar(T1​)≥var(T2​)≥⋯

    也就是说，PCA并不会对原有数据做任何的改变，而只是将“观看”原有数据的视角转换了，即，在原有数据空间中的数据的相对位置，与在主成分空间（Principal Component Space）中的相对位置是完全相同的，相当于只是更换了原有数据的基底。

    原理与步骤简述

    算法一：特征分解（Eigen Decomposition）

    假设有一n×mn\times mn×m维的数据矩阵A=[a⃗1Ta⃗2T⋮a⃗nT]\boldsymbol A = \begin{bmatrix} \vec a_1^{\rm T} \\ \vec a_2^{\rm T} \\ \vdots \\ \vec a_n^{\rm T} \end{bmatrix}A=⎣⎢⎢⎢⎡​a1T​a2T​⋮anT​​⎦⎥⎥⎥⎤​，其中nnn为样本量，mmm为观测变量的数量。PCA的步骤如下：

    1. 先对A\boldsymbol AA进行中心化（整体平移数据，使数据中心在(0,0)(0,0)(0,0)）：

        对A\boldsymbol AA求列上的平均值：a⃗T‾=1n∑i=1na⃗iT\overline {\vec a^{\rm T}} = \dfrac 1 n \sum _{i=1}^n \vec a_i ^{\rm T}aT=n1​∑i=1n​aiT​（结果为一行向量）；
      • 记Aˉ=[11⋮1]a⃗T‾=[a⃗T‾a⃗T‾⋮a⃗T‾]\boldsymbol {\bar A} = \begin{bmatrix} 1 \\ 1 \\ \vdots \\ 1 \end{bmatrix}\overline {\vec a^{\rm T}} = \begin{bmatrix} \overline {\vec a^{\rm T}} \\ \overline {\vec a^{\rm T}} \\ \vdots \\ \overline {\vec a^{\rm T}} \end{bmatrix}Aˉ=⎣⎢⎢⎢⎡​11⋮1​⎦⎥⎥⎥⎤​aT=⎣⎢⎢⎢⎡​aTaT⋮aT​⎦⎥⎥⎥⎤​；
      • 中心化后的数据矩阵X=A−Aˉ\boldsymbol X = \boldsymbol A - \boldsymbol {\bar A}X=A−Aˉ（n×mn\times mn×m维）。

    2. 计算X\boldsymbol XX的协方差矩阵C\boldsymbol CC：
      C=XTX(2) \boldsymbol C = \boldsymbol X^{\rm T}\boldsymbol X \tag{2} C=XTX(2)

    3. 对X\boldsymbol XX做特征分解，即求解特征方程
      ∣C−λI∣=0(3) \left| \boldsymbol C - \lambda \boldsymbol I\right|=0 \tag{3} ∣C−λI∣=0(3)
      可得到mmm个特征值（Eigenvalues）λi\lambda_iλi​。再解方程
      (C−λiI)w⃗i=0(4) \left( \boldsymbol C - \lambda _i \boldsymbol I\right)\vec w_i =0 \tag{4} (C−λi​I)wi​=0(4)
      其中w⃗i=[w1iw2i⋮wmi],i=1,2,⋯ ,m\vec w_i= \begin{bmatrix} w_{1i} \\ w_{2i} \\ \vdots \\ w_{mi}\end{bmatrix},\quad i=1,2,\cdots ,mwi​=⎣⎢⎢⎢⎡​w1i​w2i​⋮wmi​​⎦⎥⎥⎥⎤​,i=1,2,⋯,m，得到mmm个特征向量（Eigenvectors）w⃗i\vec w_iwi​，将它们组成矩阵 W\boldsymbol WW。可以验证，∑j=1mwji=1\sum_{j=1}^m w_{ji}=1∑j=1m​wji​=1。

    4. 将特征值降序排列，其对应的特征向量也排列到对应位置（调换W\boldsymbol WW的列）。我们这样做的原因是，var(Ti)=λi{\rm var} \left( T_i\right)=\lambda_ivar(Ti​)=λi​。

      进行特征还原：
      T=XW(5) \boldsymbol T = \boldsymbol X \boldsymbol W \tag{5} T=XW(5)
      其中：

        T\boldsymbol TT，n×mn\times mn×m维，称为主成分得分（principal component scores），即为新坐标空间中的数据点
      • W\boldsymbol WW，m×mm\times mm×m维，为特征向量组成的矩阵（称为loadings）

    5. 我们可以只取W\boldsymbol WW的前rrr列，即将m×mm\times mm×m维矩阵缩减为m×rm\times rm×r维矩阵，记作Wr\boldsymbol W_rWr​，则有：
      Tr=XWr(6) \boldsymbol T_r = \boldsymbol X \boldsymbol W_r \tag{6} Tr​=XWr​(6)
      Tr\boldsymbol T_rTr​同样为T\boldsymbol TT从n×mn\times mn×m维缩减为n×rn\times rn×r维的结果，相当于将原有的mmm个观测变量缩减为最主要的rrr个，即达到了我们的目的——降维。

    算法二：奇异值分解（Singular Value Decomposition，SVD）

    对矩阵X\boldsymbol XX进行奇异值分解：
    X=UΣV∗=UΣVT \boldsymbol X = \boldsymbol U \boldsymbol \Sigma \boldsymbol V^* = \boldsymbol U \boldsymbol \Sigma \boldsymbol V^{\rm T} X=UΣV∗=UΣVT
    其中：

    • U\boldsymbol UU（n×nn\times nn×n 维），V\boldsymbol VV（m×mm\times mm×m 维）分别称为左奇异向量和右奇异向量；
    • V∗\boldsymbol V^*V∗表示矩阵V\boldsymbol VV的共轭转置矩阵，因为X\boldsymbol XX为实数矩阵，所以可写为V∗=VT\boldsymbol V^*=\boldsymbol V^{\rm T}V∗=VT；
    • Σ\boldsymbol \SigmaΣ为一矩形对角矩阵（n×mn\times mn×m维），其对角线元素称为奇异值（singular value）。

    在PCA的问题中恒有：W=V\boldsymbol W = \boldsymbol VW=V。与W\boldsymbol WW相似，V\boldsymbol VV具有以下性质：

    • Σ\boldsymbol \SigmaΣ的对角线元素也满足σ1>σ2>⋯\sigma_1> \sigma_2 > \cdotsσ1​>σ2​>⋯，即，也是降序排列的；
    • U\boldsymbol UU和V\boldsymbol VV的第iii列，对应于Σ\boldsymbol \SigmaΣ的第iii个元素（第iii大元素）σi\sigma_iσi​；
    • U\boldsymbol UU和V\boldsymbol VV满足：U∗U=I\boldsymbol U^* \boldsymbol U = \boldsymbol IU∗U=I，V∗V=I\boldsymbol V^* \boldsymbol V = \boldsymbol IV∗V=I。

    故，公式(5)(5)(5)可写为：

    T=XW=XV=UΣV∗V=UΣ \begin{aligned} \boldsymbol T & = \boldsymbol X \boldsymbol W\\ & = \boldsymbol X \boldsymbol V\\ & = \boldsymbol U \boldsymbol \Sigma \boldsymbol V^*\boldsymbol V \\ & = \boldsymbol U \boldsymbol \Sigma \\ \end{aligned} T​=XW=XV=UΣV∗V=UΣ​

    即：

    T=UΣ(7) \boldsymbol T = \boldsymbol U \boldsymbol \Sigma \tag{7} T=UΣ(7)

    我们同样可以对Σ\boldsymbol \SigmaΣ只取第一个r×rr \times rr×r的块，记作Σr\boldsymbol \Sigma _ rΣr​，相应地U\boldsymbol UU也只取前rrr列，记作Ur\boldsymbol U_rUr​，则有

    Tr=UrΣr(8) \boldsymbol T_r = \boldsymbol U_r \boldsymbol \Sigma _r \tag{8} Tr​=Ur​Σr​(8)

    rrr的选取标准

    计算方差的累积贡献率：
    f(k)=∑i=1iλk∑i=1mλi,k=1,2,⋯(9) f(k)=\dfrac{\sum _{i=1}^i \lambda_k}{\sum_{i=1}^m \lambda_i},\quad k = 1,2,\cdots \tag{9} f(k)=∑i=1m​λi​∑i=1i​λk​​,k=1,2,⋯(9)
    作出其图像。因为λ1>λ2>⋯\lambda_1> \lambda_2 > \cdotsλ1​>λ2​>⋯，故f(k)f(k)f(k)为一单调递增的函数，且其递增速度随着kkk增加逐渐降低。

    一般地，我们可以取使得f(r)≥f(r) \gef(r)≥某一阈值（如95%95\%95%）的最小的rrr，这样最多只会损失掉5%的方差。

    对于SVD法，将公式(9)(9)(9)中的λ\lambdaλ换为σ\sigmaσ即可。

    两种算法的比较

    在采用特征分解法时，我们无法避免计算XTX\boldsymbol X^{\rm T}\boldsymbol XXTX，而在观测变量数mmm非常大时，这一算法的劣势将被无限放大（协方差矩阵为m×mm\times mm×m维）。

    而采用SVD算法，则只需要计算Tr=UrΣr\boldsymbol T_r = \boldsymbol U_r \boldsymbol \Sigma _rTr​=Ur​Σr​，而Σr\boldsymbol \Sigma _rΣr​为对角阵，显然这一算法的计算量要小很多（这一点类似于DFT与FFT之间的比较）。默认情况下，MATLAB中的

    pca
    函数也会使用SVD算法。

    MATLAB的实现方法

    我们先载入MATLAB自带的数据集

    fisheriris
    （该数据集统计了三种鸢尾花的花萼长、花萼宽、花瓣长、花瓣宽），然后进行中心化处理，并计算协方差矩阵：

    load fisheriris;
X = meas;     % n = 150, m = 4

% 中心化
meanX = ones(size(X,1), 1) * mean(X);
centredX = X - meanX;

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

    特征分解法：利用
    eig
    函数

    [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)    % 可以验证每个特征向量各元素的平方和均为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 = 2即可

figure;
scatter(Tr(:,1), Tr(:,2), 130, categorical(species), '.');
colormap(winter);
xlabel('Principal Component 1');
ylabel('Principal Component 2');


    SVD法：利用
    svd
    函数

    [U, Sigma, V] = svd(X);    % 可以检验，W和V完全相同（向量的正负号不影响）
Vr = V(:, 1:2);    % 提取前两个主成分的特征向量
Tr = X * Vr;  %  新坐标空间的数据点
% 画图部分同上，略

    利用
    pca
    函数

    pca
    的常用调用格式如下：

    [loadings, scores] = pca(X, 'NumComponents', r);

    其中：

    • loadings
      为Wr\boldsymbol W_rWr​矩阵（m×rm\times rm×r维），即主成分系数；
    • scores
      为Tr\boldsymbol T_rTr​矩阵（n×rn\times rn×r维）；
    • eigenvalues
      为所有特征值组成的列向量。

    且默认情况下，

    pca
    会自动将数据
    X
    中心化。

    在本例中，我们可以略去中心化的步骤，直接调用该函数：

    [Wr, Tr, ev] = pca(X, 'NumComponents',2);
% 画图部分略

    应用

    聚类分析

    如上面鸢尾花的例子中，降维后的数据仍可以清晰地分为三类。这样，当我们拿到一种鸢尾花，计算相应的T1T_1T1​和T2T_2T2​，将结果画在散点图中，我们就可以判断出其属于哪一种鸢尾花。

    例如，我们在电商平台浏览并购买商品时，平台就会收集你的年龄、性别、购买商品平均价格、购买频率、最初浏览商品直到最终购买之间的时间间隔等等大量、多维度的信息，然后进行降维，将你归于“大学生”“白领”“一家三口”等类别，然后定向为你推送促销商品的通知。

    图像压缩

    假设有一张n×mn \times mn×m的图片X\boldsymbol XX，根据(6)(6)(6)式，我们可以利用矩阵Wr\boldsymbol W_rWr​将其降至n×rn\times rn×r维。如果我们只传输Tr\boldsymbol T_rTr​和Wr\boldsymbol W_rWr​，就可以反推还原出X\boldsymbol XX，而压缩比可达nmr(n+m)≤nm2r\dfrac {nm}{r(n+m)} \le \dfrac {\sqrt{nm}} {2r}r(n+m)nm​≤2rnm​​

    人脸检测与匹配

    假设有nnn个人脸训练样本，每个样本共mmm个像素，每个样本由其像素灰度值展开组成一个行向量，按列组成矩阵X\boldsymbol XX。

    同样先求其平均向量（称为“平均脸”），中心化后求协方差矩阵，并进行特征分解。任何一幅人脸图像都可以变换到主成分空间，得到“特征脸”（Eigenfaces）。

    这样，若将待识别的人脸做同样的变换，遍历已有的特征脸中，寻找最为接近的特征脸，即完成了匹配。

    • 点赞
    • 收藏
    • 分享
    • 文章举报
    jz8_AWarmohb 发布了5 篇原创文章 · 获赞 0 · 访问量 507 私信 关注
    内容来自用户分享和网络整理，不保证内容的准确性，如有侵权内容，可联系管理员处理 点击这里给我发消息
    标签： 
    相关文章推荐
    新的分享
    章节导航