引子

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

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

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

如何获取正交基

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

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

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

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

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

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

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

满足什么特征的矩阵的特征向量为单位向量？

如果一个矩阵的列向量的范数相等，则其特征向量是单位化的。（至于为什么，留个坑，以后来填 oVo）

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

因此要使满足特征向量为单位向量这一条，我们要对矩阵 $R_1=A^T A$ 做处理，我们对矩阵 $A$ 的列向量进行处理，使其拥有相同的范数，（在PCA中一般是进行标准化， $a_i=(a_i-mean(a_i))/std(a_i)$ ，其范数为 $\sqrt{n-1}$ ，因此效果一样）：
$A_2=AP= \left[ a_1,a_2,...,a_n \right] \left[ \begin{array}{c} \frac{1}{||a_1||}& \quad 0& \quad \cdots& \quad 0 \\ 0& \quad \frac{1}{||a_2||}& \quad \cdots& \quad 0 \\ \vdots& \quad \vdots& \quad \ddots& \quad \vdots& \\ 0& \quad 0& \quad \cdots& \quad \frac{1}{||a_2||} \\ \end{array} \right] \tag {8}$

则矩阵 $R_2=A_2^T A_2$ 的列向量范数也相同（这个坑之后一块儿填 oVo）。至此，我们已经找到了特征向量自然单位正交的矩阵 $R_2$ 。

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

如何计算主成分得分

对于矩阵 $A_2=[b_1^T,b_2^T,...,b_m^T]^T$ 的每一行向量 $b_i$ ，与向量 $v_j$ 做内积，即向量 $b_i^T$ 中向量 $v_j$ 的分量大小，因此主成分得分为
$F_{m\times p}=A_2 V= \left[ \begin{array}{c} b_1^T\cdot v_1& \quad b_1^T\cdot v_2& \quad \cdots& \quad b_1^T\cdot v_p \\ b_2^T\cdot v_1& \quad b_2^T\cdot v_2& \quad \cdots& \quad b_2^T\cdot v_p \\ \vdots& \quad \vdots& \quad \ddots& \quad \vdots& \\ b_m^T\cdot v_1& \quad b_m^T\cdot v_1& \quad \cdots& \quad b_m^T\cdot v_p \\ \end{array} \right] \tag {9}$
其中 $F_{ij}$ 表示 $A_2$ 的第 $i$ 行数据中， $v_j$ 分量的大小。

如何恢复数据

我们根据
$A_3=FV^T \tag {10}$
进行数据恢复（因为 $v_i$ 单位正交， $VV^T=I$ ，由(9)式可推知(10)式） $A_2$ 的近似，其列向量是经过标准化处理的，因此要得到数据 $A$ 的近似，还要进行反变换。如果是使用式(8)进行的标准化，那么反变换为
$A_4=A_3P^{-1}= A_3 \left[ \begin{array}{c} ||a_1||& \quad 0& \quad \cdots& \quad 0 \\ 0& \quad ||a_2||& \quad \cdots& \quad 0 \\ \vdots& \quad \vdots& \quad \ddots& \quad \vdots& \\ 0& \quad 0& \quad \cdots& \quad ||a_2|| \\ \end{array} \right] \tag {11}$
即将每一列进行复原。

标准PCA流程

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

将数据 $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] \tag {12}$

其中
$\overline X=mean(X)=\frac{1}{m}\sum_{j=1}^{m}x_j \tag {13}$

$std(X)=\frac{1}{m-1}\sqrt{\sum_{j=1}^{m}(x_j-\overline X)^2} \tag {14}$

经过此步， $A_2$ 的每一列的范数均为 $\sqrt{n-1}$ ；

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

$\sigma(X_1,X_2,...,X_n)= \frac{1}{n-1} \left[ \begin{array}{c} \sigma(X_1,X_1)& \quad \sigma(X_1,X_2)& \quad \cdots& \quad \sigma(X_1,X_n) \\ \sigma(X_2,X_1)& \quad \sigma(X_2,X_2)& \quad \cdots& \quad \sigma(X_2,X_n) \\ \vdots& \quad \vdots& \quad \ddots& \quad \vdots& \\ \sigma(X_n,X_1)& \quad \sigma(X_n,X_2)& \quad \cdots& \quad \sigma(X_n,X_n) \\ \end{array} \right] \tag {15}$
其中方差计算公式为
$\sigma (X,Y)=\frac{1}{m-1}\sum_{j=1}^{m}(x_j-\overline X)(y_j-\overline Y) \tag {16}$

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

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

$\sum_{i=1}^{k}\lambda_k>T \tag {18}$

获得主成分 $V$

$V=[v_1,v_2,...,v_k] \tag {19}$

计算主成分得分 $F$

$F=A_2V \tag {20}$

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

PCA数据恢复

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

$c_j=std_jb_j+mean_j \tag {21}$

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

第二种，采用式（8）所示的归一化方法:

%----------主成分分析(PCA)---------
% author: 今朝无言
function [V,F,A]=PCA2(m,T)
    % m为二维数据，T为信息保留率
    A=sqrt(diag(m'*m))';
    s=size(m);
    
    m2=m./(ones(s(1),1)*A); %归一化
    R=m2'*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; %主成分得分   Fij表示m第i行数据对特征向量j的得分（数据i对向量j的投影长度）
end

实验验证：

% test PCA
clc,clear,close all

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

%% 
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')

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

m3=F2*V2';
m3=m3.*(ones(s(1),1)*A);

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(A(:),1))/size(m(:),1))
fprintf('PCA2 MSE=%f\n',sum(sqrt(m3(:)-m(:)))/(s(1)*s(2)))

结果：

在这里插入图片描述