基于特征向量的主成分分析(PCA)原理解释

引子

首先看一下如何对一维向量的进行分解，我们知道，一个 $n$ 维向量 $a$ 可以由 $n$ 个正交向量线性 $v_i,i=1,2,...,n$ 组合而来，即
$$a=\sum _{k=1}^n{k}_i{v}_i\text{\hspace{1em}(1)}$$
在图像与信号处理领域常用的傅里叶变换（Fourier Transform）、小波变换（Wavelet Transform）等，都是基于此理论。

我们知道，随着傅里叶级数的阶数的增加，其结果将越来越接近原始信号，因此实际上我们只需要记录能量较高的几个分量及其幅值，就可以基本复原原始信号。那么对于图像这样的二维数据，要怎样分解呢？

实际上，答案是对每一行向量分别进行分解，它们使用同一组正交基。

如何获取正交基

我们知道，对于矩阵 $A$ ，如果有一组标准正交基 $v_i$ ，那么矩阵 $A$ 的每一行都可以由向量 $v_i$ 线性组合得到。

然而正交基有很多，如何选取正交基，可以从中挑选尽可能少的向量，而能表示尽可能多信息？（如何提高压缩效率）

假设我们已经有了一组正交向量 $u_i$ ，那么对于向量 $u$ ，如何衡量它所表示的信息的多少呢？我们知道，向量 $a$ 与向量 $u$ 做内积，即
$$a\cdot u=\sum _{j=1}^n{a}_i{u}_i\text{\hspace{1em}(2)}$$
表征了向量 $a$ 中向量 $u$ 占据的分量大小，对于数据 $A_{m\times n}=[a_1^T,a_2^T,...,a_m^T{]}^T$ ，其每一行 $a_i$ 中向量 $u$ 的分量大小 $e_i$ 为
$$Au=\left[\begin{array}{c}{a}_1\\ a_2\\ ⋮\\ a_m\end{array}\right]u=\left[\begin{array}{c}{e}_1\\ e_2\\ ⋮\\ e_m\end{array}\right]\text{\hspace{1em}(3)}$$
为了使向量 $u$ 能够表示尽可能多信息，即 $e_i$尽可能大，我们计算其平方和，并希望其尽可能大
$$\begin{array}{lc}E& ={\sum }_{j=1}^m{e}_j^2\\ & =(Au{)}^T(Au)\\ & =u^T{A}^{T}Au\end{array}\text{\hspace{1em}(4)}$$
这里我们发现，有一个熟悉的结构出现了！回忆矩阵的特征向量，对于矩阵 $R$ ，其特征向量组成的矩阵 $Q=[v_1,v_2,\cdots ,v_n]$ ，其特征值组成的对角矩阵 Σ = d i a g { λ 1 , λ 2 , ⋯   , λ n } \Sigma=diag\{\lambda_1,\lambda_2,\cdots,\lambda_n \} Σ=diag{λ1​,λ2​,⋯,λn​} ，有
$$R=Q\Sigma Q^{-1}\text{\hspace{1em}(5)}$$
变换一下形式
$$\Sigma =Q^{-1}RQ\text{\hspace{1em}(6)}$$
如果我们把 $A^{T}A$ 视作 $R_1$ ，把 $[u_1,u_2,\cdots ,u_n]$ 视作 $Q_1$ ，我们有
$${\Sigma }_1=Q_1^T{R}_1{Q}_1\text{\hspace{1em}(7)}$$
然而我们发现式(7)与式(6)还存在微小差别，那么让我们继续回忆线性代数的知识，什么情况下矩阵 $Q$ 的逆等于其自身的转置呢？答案是如果方阵 $Q$ 的每一列(行)向量都是单位向量，且两两之间相互正交，则 $Q^{T}Q=I$ ， $Q{Q}^T=I$ ，即 $Q^{-1}=Q^T$ 。

也就是说，如果我们将能够使 $R_1$ 的特征向量单位正交，那么问题将转化为：选择一特征向量，使其对应的特征值尽可能大。

什么样的矩阵的特征向量单位正交？

• 满足什么特征的矩阵的特征向量相互正交？

我们知道，实对称矩阵的特征向量是自然正交的。

由于矩阵 $(A^{T}A{)}^T=A^{T}A$ ，因此矩阵 $A^{T}A$ 是实对称阵，即其特征向量是一组正交向量。

• 特征向量单位化

至于单位化，由于 $Ax=\lambda x$，所以
$$A\frac{x}{||x||}=\lambda \frac{x}{||x||}\text{\hspace{1em}(8)}$$
所以 $x/||x||$ 也是 $A$ 的特征向量。

至此，我们得到了 $R_1$ 的一组单位正交的特征向量组。

求矩阵 $R_1$ 的特征值与特征向量，然后对特征值进行排序，从最大特征值对应的特征向量开始挑选，直到达到所要求的信息保留率为止，则挑选出的向量组 $[v_{k_1},v_{k_2},...,v_{k_p}]$ 即所谓的“主成分”，记为 $V$。

如何计算主成分得分

对于矩阵 $A=[b_1^T,b_2^T,...,b_m^T{]}^T$ 的每一行向量 $b_i$，与向量 $v_j$ 做内积，即向量 $b_i^T$ 中向量 $v_j$ 的分量大小，因此主成分得分为
$$F_{m\times p}=AV=\left[\begin{array}{ccccc}{b}_1^T\cdot v_1& b_1^T\cdot v_2& \cdots & b_1^T\cdot v_p& \\ b_2^T\cdot v_1& b_2^T\cdot v_2& \cdots & b_2^T\cdot v_p& \\ ⋮& ⋮& \ddots & ⋮& \\ b_m^T\cdot v_1& b_m^T\cdot v_1& \cdots & b_m^T\cdot v_p& \end{array}\right]\text{\hspace{1em}(9)}$$
其中 $F_{ij}$ 表示 $A$ 的第 $i$ 行数据中，$v_j$ 分量的大小。

如何恢复数据

我们根据
$$A_2=F{V}^T\text{\hspace{1em}(10)}$$
进行数据恢复（ 因为 $v_i$ 单位正交， $V{V}^T=I$ ，由(9)式可推知(10)式 ） ，$A_2$ 即 $A$ 的近似。

标准PCA流程

好了，现在我们已经了解了如何用特征向量来进行数据压缩，让我们来看一下标准的PCA是如何做的：

1. 将数据 $A_{m\times n}$ 的每一列进行标准化

$$A_2=\left[\frac{a_1-mean(a_1)}{std(a_1)},\frac{a_2-mean(a_2)}{std(a_2)},...,\frac{a_n-mean(a_n)}{std(a_n)}\right]\text{\hspace{1em}(11)}$$

其中
$$\overline{X}=mean(X)=\frac{1}{m}\sum _{j=1}^m{x}_j\text{\hspace{1em}(12)}$$

$$std(X)=\frac{1}{m-1}\sqrt{\sum _{j=1}^m(x_j-\overline{X}{)}^2}\text{\hspace{1em}(13)}$$

2. 计算 $A_2$ 的协方差矩阵，矩阵 $A_2=[X_1,X_2,\cdots ,X_n]$ 的协方差矩阵为

$$R=\sigma (X_1,X_2,...,X_n)=\frac{1}{n-1}\left[\begin{array}{ccccc}\sigma (X_1,X_1)& \sigma (X_1,X_2)& \cdots & \sigma (X_1,X_n)& \\ \sigma (X_2,X_1)& \sigma (X_2,X_2)& \cdots & \sigma (X_2,X_n)& \\ ⋮& ⋮& \ddots & ⋮& \\ \sigma (X_n,X_1)& \sigma (X_n,X_2)& \cdots & \sigma (X_n,X_n)& \end{array}\right]\text{\hspace{1em}(14)}$$

其中方差计算公式为
$$\sigma (X,Y)=\frac{1}{m-1}\sum _{j=1}^m(x_j-\overline{X})(y_j-\overline{Y})\text{\hspace{1em}(15)}$$

由于经过变换后，$A_2$ 的每一列向量均值为0，因此
$$R=\frac{1}{n-1}{A}_2^T{A}_2\text{\hspace{1em}(16)}$$
根据之前的讨论，这里的协方差矩阵 $R$ 的特征向量，将是一组单位正交向量。

3. 计算 $R$ 的特征值与对应的特征向量，并按照特征值大小降序排列。
4. 依据主成分信息保留率 $T$ ，计算出保留的特征向量个数 $k$ ，使

$$\sum _{i=1}^k{\lambda }_k>T\text{\hspace{1em}(17)}$$

5. 获得主成分 $V$

$$V=[v_1,v_2,...,v_k]\text{\hspace{1em}(18)}$$

6. 计算主成分得分 $F$

$$F=A_{2}V\text{\hspace{1em}(19)}$$

7. 算法结束，输出主成分 $V_{n\times k}$ 、主成分得分 $F_{m\times k}$ 、$A$ 的列向量的 均值 $M_{1\times n}$ 与 标准差 $D_{1\times n}$

PCA数据恢复

1. 计算 $B=[b_1,b_2,...,b_n]=F{V}^T$ ；
2. 对每一列进行还原，得到复原数据 $C=[c_1,c_2,...,c_n]$ ，其中

$$c_j=st{d}_j{b}_j+mea{n}_j\text{\hspace{1em}(20)}$$

matlab实验

第一种，使用标准的PCA方法：

%----------主成分分析(PCA)---------
% author: 今朝无言
function [V,F,M,D]=PCA(m,T)
    % m为二维数据，T为信息保留率
    M=mean(m,1);
    D=std(m,1);
    s=size(m);
    
    m2=(m-ones(s(1),1)*M)./(ones(s(1),1)*D); %标准化
    R=cov(m2); %协方差矩阵
    [x,y]=eig(R); %求特征向量x、特征值y
    
    y=diag(y);
    [y,index]=sort(y);
    y=y(end:-1:1); index=index(end:-1:1); %降序排列
    x=x(:,index); %将x按y的次序排列
    
    DS(:,1)=y; %特征值
    DS(:,2)=DS(:,1)/sum(DS(:,1)); %贡献率
    DS(:,3)=cumsum(DS(:,2)); %累积贡献率
    
    k=find(DS(:,3)>=T); k=k(1); %选取的特征向量个数
    V=x(:,1:k); %k个特征向量
    
    F=m2*V; %主成分得分
end

第二种，不归一化:

%----------主成分分析(PCA)---------
% author: 今朝无言
function [V,F]=PCA2(m,T)
    % m为二维数据，T为信息保留率
    R=m'*m;
    [x,y]=eig(R); %求特征向量x、特征值y
    
    y=diag(y);
    [y,index]=sort(y);
    y=y(end:-1:1); index=index(end:-1:1); %降序排列
    x=x(:,index); %将x按y的次序排列
    
    DS(:,1)=y; %特征值
    DS(:,2)=DS(:,1)/sum(DS(:,1)); %贡献率
    DS(:,3)=cumsum(DS(:,2)); %累积贡献率
    
    k=find(DS(:,3)>=T); k=k(1); %选取的特征向量个数
    V=x(:,1:k); %k个特征向量
    
    F=m*V; %主成分得分   Fij表示m第i行数据对特征向量j的得分（数据i对向量j的投影长度）
end

实验验证：

% test PCA
clc,clear,close all

%% 
m=imread('1.jpg');
m=double(rgb2gray(m));
s=size(m);

%% 标准PCA
T=0.95;
[V,F,M,D]=PCA(m,T);

m2=F*V';
m2=m2.*(ones(s(1),1)*D)+ones(s(1),1)*M;

figure
set(gcf,'unit','normalized','position',[0.2,0.2,0.64,0.32]);
subplot(1,2,1); imshow(uint8(m))
subplot(1,2,2); imshow(uint8(m2))
title('PCA')

%% 不归一化的PCA
T2=0.999;
%由于第一个向量比重太大，因此T要设置的大一些，然而在更小的压缩率下，恢复效果比上述PCA强很多！
[V2,F2]=PCA2(m,T2);

m3=F2*V2';

figure
set(gcf,'unit','normalized','position',[0.2,0.2,0.64,0.32]);
subplot(1,2,1); imshow(uint8(m))
subplot(1,2,2); imshow(uint8(m3))
title('PCA2')

%% 
fprintf('PCA 压缩率=%f\t\t',...
    (size(V(:),1)+size(F(:),1)+size(M(:),1)+size(D(:),1))/size(m(:),1))
fprintf('PCA MSE=%f\n',sum(sqrt(m2(:)-m(:)))/(s(1)*s(2)))

fprintf('PCA2 压缩率=%f\t\t',...
    (size(V2(:),1)+size(F2(:),1))/size(m(:),1))
fprintf('PCA2 MSE=%f\n',sum(sqrt(m3(:)-m(:)))/(s(1)*s(2)))

结果：

可以看到，第二种方法在获得更小的压缩率的同时，同时数据恢复效果也更好。单从视觉上看，尽管均方误差差不多，但前者的噪声肉眼明显可见，而后者几乎完全看不到噪声。