主成分分析法（PCA）解析与MATLAB实践

PCA(Principal Components Analysis)即主成分分析，也称主分量分析或主成分回归分析法，是一种无监督的数据降维方法，在机器学习中常用于特征降维提取主要特征以减少计算量。PCA主要原理是将高维原数据通过一个转换矩阵，映射到另一组低维坐标系下，从而实现数据降维。举个简单的例子，设X1，X2为两组数据，将他们以坐标的形式画在坐标轴中，如下图所示，
图中点的横纵坐标分别为X1，X2的值，如果我们把数据做一个变换，使其顺时针旋转一个角度得到下图所示结果：
再向横轴方向与纵轴方向做一个投影得到以下结果：


我们可以发现数据在横轴方向的投影差异性较大，也可以说在横轴的投影包含的信息或特征较多，如此我们便可以丢弃对纵轴的投影而保留对横轴的投影，这样便把数据从二维降到了一维。以上例子用的数据只有两个维度，那么对于多维的数据我们也可以用PCA进行处理，这里需要用到一个基本概念那就是协方差（covariance）以及协方差矩阵，上例中有两组数据，那么对应的协方差矩阵就是2*2，两组数据的协方差通常表示为cov（x1，x2）。上述例子中有两组数据，则协方差矩阵可表示为：

                      cov（x1，x1）    cov（x1，x2）
                      cov（x2，x1）    cov（x2，x2）

MATLAB中求协方差的函数为cov()，调用此函数可以得到以下结果：

由于协方差矩阵通常为实对称矩阵，因此可以将其对角化得到对角矩阵：

对角元素其实就是我们所需要投影的维度中信息占比，第一个维度中的主要成分占1.284，第二个占0.0491，显然往第一个维度做投影所得到的信息较多，即第一个维度得到的是信号的主要成分。对于多维信号，其主要成分可能有多个，假设有五个信号，如下图所示：

我们同样可以先求五个信号的协方差矩阵：

  cov（x1，x1） cov（x1，x2）  cov（x1，x3）  cov（x1，x4）  cov（x1，x5）
  cov（x2，x1） cov（x2，x2）  cov（x2，x3）  cov（x2，x5）  cov（x2，x5）
  cov（x3，x1） cov（x3，x2）  cov（x3，x3）  cov（x3，x4）  cov（x3，x5）
  cov（x4，x1） cov（x4，x2）  cov（x4，x3）  cov（x4，x4）  cov（x4，x5）
  cov（x5，x1） cov（x5，x2）  cov（x5，x3）  cov（x5，x4）  cov（x5，x5）

接着再将其对角化，得到的五个成分占比如图所示（五个成分占比分别为第一个除以五个之和，第一个加第二个除以五个之和，第一个加第二个加第三个除以五个之和，以此类推）

可以看到，前两个维度成分加起来大约占到了总成分的95%，说明将信号投影到前两个维度就可以得到信号的主要信息，在MATLAB中可以用pca()函数来实现主成分提取，处理后的结果如下图所示，前两个波形即为五个信号中的主要成分。

上述第二个例子的代码如下：

clear, close all
clc

%% initialize parameters
samplerate=500; % in Hz
N=1000; % data length

freq1=5; % in Hz
freq2=7; % in Hz
taxis=[1:N]/samplerate;

PCnum=2;  % the number of PC used to reconstruct signals

%% generate test signals 
C1 = 0.75*sin(2*pi*freq1*taxis);        % 1st component: a sine wave
C2 = sawtooth(2*pi*freq2*taxis,0.5);    % 2nd component: a triangular wave

% Combine data in different proportions
X(1,:) = 0.5*C1 + 0.5*C2 + 0.1*rand(1,N);
X(2,:) = 0.7*C1 + 0.2*C2 + 0.1*rand(1,N);
X(3,:) = 0.2*C1 + 0.7*C2 + 0.1*rand(1,N);
X(4,:) = -0.3*C1 - 0.6*C2 + 0.3*rand(1,N);
X(5,:) = 0.6*rand(1,N);    % Noise only

% Center data by subtracting mean
X = X - mean(X,2)*ones(1,N);

figure, 
for i=1:size(X,1)
    subplot(size(X,1),1,i)
    plot(taxis,X(i,:)),xlim([taxis(1) taxis(end)])
end

%% Principal Components Analysis using princomp funcion (using eig and svd)
[U,PC,eigenVal]=pca(X')

for i=1:size(X,1)
    eigen_perc(i)=sum(eigenVal(1:i))/sum(eigenVal)*100;  % calculate accumulated percentage of eigenvalues
end

figure, 
for i=1:size(PC,2)
    subplot(size(PC,2),1,i)
    plot(taxis,PC(:,i)),xlim([taxis(1) taxis(end)])
end

figure,plot(eigen_perc,'-o')
xlabel('dimension'),ylabel('percentage of accumulate eigenvalues')

%% Check the covariance of principal components (PC)
cov(PC)  % make sure if the PCs are uncorrelated !
% the off-diagonal terms should be 0

%% Reconstruct the Signal only use the first PCnum PCs
newX = U(:,1:PCnum)*PC(:,1:PCnum)';

figure, 
for i=1:size(newX,1)
    subplot(size(newX,1),1,i)
    plot(taxis,X(i,:)),hold on
    plot(taxis,newX(i,:),'r'),xlim([taxis(1) taxis(end)])
end

以上内容主要参考了卢家峰老师课程的讲义，课程中的程序用作上述示例代码时有改动，上述讲解若有错误或不足之处还望各位读者批评指正！