PCA和RPCA算法流程以及在人脸识别和噪声分离中的应用

新出炉的烤鸡

已于 2024-09-04 16:58:48 修改

阅读量1.6k

点赞数 5

文章标签： matlab 机器学习

于 2023-08-29 16:23:42 首次发布

本文链接：https://blog.csdn.net/nih2221/article/details/132548583

版权

1 问题描述

2 PCA算法(参考PCA论文和PCA数学原理博客介绍）

2.3.2.2 SVD加速求解特征向量（方法二）

1 问题描述

利用PCA或者RPCA技术实现对数据的提取或处理（数据类型不限制，例如可以是图像），给出详细过程和可视化结果。作业内容可以包括但不限于原理介绍、数据处理过程、代码（关键部分即可）、结果展示等。

2 PCA算法(参考PCA论文和PCA数学原理博客介绍）

一般我们提到降维最容易想到的算法就是主成分分析，下面我们就对其原理做一个总结。

主成分分析（Principal components analysis，以下简称PCA）主要有纯特征值分解（EVD）和奇异值分解（SVD）两种方法，它提供了一个数据驱动的分层坐标系统来表示高维相关数据。PCA是最重要的降维方法之一，在数据压缩消除冗余和数据噪音消除等领域都有广泛的应用，现在也广泛用于深度学习的数据白化操作中。

我们不仅要对EVD方法有了解，并且要会运用SVD方法进行运算减负。实际上SVD方法在实际生活中更常用，我也会重点介绍SVD方法。

一般来说，我们感兴趣的是分析一个大数据集 $X\in C^{n\times m}$ :

$X=\begin{bmatrix} |&|& &| \\ x_1&x_2&\cdots&x_m\\ |&|& &| \end{bmatrix}$

列 $x_k \in C^{n}$ 可能是一次模拟或实验的所有测量值，即一个样本。

指标k是表示第k个不同的测量集的标签。X将由数据的时间序列组成，并且 $x_k = x(k\Delta t )$ 。通常状态维度n非常大，在数百万或数十亿个自由度的数量级上。列通常被称为采样点数（snapshots），m是X中的采样样本数量。在许多系统中 $n\gg m$ ，得到一个又高又细的矩阵。通常，在单个实验中收集许多测量值，测量值相当于变量，并将这些测量值排列成行向量，即特征。测量值可以是可观察到的特征，例如特定人类个体的人口统计学特征。

PCA降维问题的优化目标：将一组n维向量降为k维（k大于0，小于n），其目标是选择k个单位（模为1）正交基，使得原始数据变换到这组基上后，各字段两两间协方差为0，而字段的方差则尽可能大（在正交的约束下，取最大的k个方差）。

而协方差矩阵统一了变量（特征）内方差及变量间协方差关系：

$\tfrac{1}{m-1}XX^{*}=\tfrac{1}{m-1}\begin{bmatrix} \sum_{j=1}^{m}x^{2}_{1j} & \cdots &\sum_{j=1}^{m}x_{j1}x_{nj} \\ \vdots & \ddots &\vdots \\ \sum_{j=1}^{m}x_{nj}x_{j1} & \cdots & \sum_{j=1}^{m}x^{2}_{nj} \end{bmatrix}$

因此我们的很多运算是针对协方差进行的。

2.1 纯特征值分解（EVD）原理

如果一个向量v是矩阵A的特征向量，将一定可以表示成下面的形式：

$Av=\lambda v$

其中，λ是特征向量v对应的特征值，一个矩阵的一组特征向量是一组正交向量。

将所有特征值λ排列成对角矩阵D；对应特征向量v排列成矩阵V，正交单位化后，也可写成：

$AV=\Sigma V$

对于矩阵A，有一组特征向量v，将这组向量进行正交化单位化，就能得到一组正交单位向量。特征值分解，就是将矩阵A分解为如下式：

$A=V^{-1}\Sigma V$

其中，V是矩阵A的特征向量组成的矩阵， $\Sigma$ 则是一个对角阵，对角线上的元素就是特征值。

在线性代数中，对称矩阵一定可以对角化。

2.2 奇异值分解（SVD）原理

SVD是特征分解的推广，是存在于每一个复值矩阵 $X\in C^{n\times m}$ 的唯一矩阵分解:

$X=U\Sigma ^{'}V^{*}$

其中， $U\in C^{n\times n}$ 和 $V\in C^{m\times m}$ 是具有正交列的酉矩阵， $\Sigma^{'}\in R^{n\times m}$ 是一个在对角线上有实数、非负元素且对角线外有零元素的矩阵。*表示复共轭转置。

当 $n\geqslant m$ ，矩阵 $\Sigma ^{'}$ 在对角线上最多有m个非零元素，因此，可以使用经济型奇异值分解（economy SVD）精确地表示X:

$X=U\Sigma ^{'}V^{*}=\begin{bmatrix} \hat{U}& \hat{U}^{\perp } \end{bmatrix}\begin{bmatrix} \hat{\Sigma } \\ 0 \end{bmatrix}V^{*}$

$\hat{U}^{\perp }$ 的列向量张成一个向量空间，它与 $U^{\perp }$ 张成的向量空间互补且正交。U的列称为X的左奇异向量，V的列称为右奇异向量。 $\hat{\Sigma}\in C^{m\times m}$ 的对角线元素称为奇异值，它们按从大到小的顺序排列。X的秩等于非零奇异值的个数。

SVD与涉及相关矩阵 $XX^{*}$ 和 $X^{*}X$ 的特征值问题密切相关，我们发现：

$XX^{*}=U\begin{bmatrix} \hat{\Sigma }\\ 0 \end{bmatrix}V^{*}V\begin{bmatrix} \hat{\Sigma }& 0 \end{bmatrix}U^{*}=U\begin{bmatrix} \hat{\Sigma }^{2}&0 \\ 0 & 0 \end{bmatrix}U^{*}$

$X^{*}X=V\begin{bmatrix} \hat{\Sigma }& 0 \end{bmatrix}U^{*}U\begin{bmatrix} \hat{\Sigma}\\ 0 \end{bmatrix}V^{*}=V\hat{\Sigma}^{2}V^{*}$

由于U和V是酉矩阵，U， $\hat{\Sigma}$ ，V是以下特征值问题的解:

$XX^{*}U=U\begin{bmatrix} \hat{\Sigma }^{2}&0 \\ 0 & 0 \end{bmatrix}$

$X^{*}XV=V\hat{\Sigma}^{2}$

换句话说，X的每个非零奇异值是 $X^{*}X$ 和 $XX^{*}$ 的特征值的正平方根，它们具有相同的非零特征值。

这提供了SVD的直观解释，其中U的列是相关矩阵 $XX^{*}$ 的特征向量，V的列是 $X^{*}X$ 的特征向量。我们选择按大小降序排列奇异值，因此U的列是根据它们在X的列中捕获的相关性进行分层排序的;V同样捕获了X行中的相关性。

2.3 PCA步骤

假设有m条n维数据。

输入： $n\times m$ 维样本集，m个样本要降维从特征维数n降为k。

输出：降维后的样本集 $X^{'}$ 。

为了保持不同算法所得结果一致，首先对数据进行中心化：

现在我们计算平均列(每一行平均，即一个变量零均值化)：

$\bar{x}_i=\tfrac{1}{m}\sum_{j=1}^{n}X_{ij}$

则平均矩阵为：

$\bar{X}=\begin{bmatrix} \bar{x}_1\\ \vdots \\ \bar{x}_2 \end{bmatrix}\begin{bmatrix} 1 & \cdots & 1 \end{bmatrix}$

处理后的矩阵为B（ $n\times m$ ）:

$B=X-\bar{X}$

2.3.1 EVD实现PCA

①计算协方差矩阵

B的行的协方差矩阵C（度量特征间的相关性）：

$C=\tfrac{1}{m-1}BB^{*}$ （特征维度作内积，衡量相关性）

②计算特征值和特征向量

通过计算C的特征分解可以得到主成分:

$C=U\Sigma U^{-1}$

它是肯定存在的，因为C是Hermite矩阵，是一个对称矩阵。

③将特征向量按对应特征值大小从上到下按行排列成矩阵，取出最大的k个特征值对应的特征向量，按列构成矩阵 $\begin{bmatrix} u_1 & \cdots & u_k \end{bmatrix}$ ,标准化后，变成特征向量矩阵 $U_k$ （ $n\times k$ ）。

④ $X^{'}=U_k^{*}X$ 即为降到k维特征后的数据（此时为 $k\times m$ ）。

2.3.2 SVD实现PCA

2.3.2.1 SVD表示结果（方法一）

①计算矩阵：

$D=B^{*}B$

②特征分解：

$D=V\Sigma V^{-1}$ ，V是特征矩阵

当V正交标准化：

V是SVD分解 $B=U\Sigma ^{'}V^{*}$ 中的V

奇异值：

$\hat{\Sigma}=\sqrt{\Sigma}$

③将奇异按对应特征值从大到小排列，取出最大的k个奇异值，构成矩阵 $\Sigma^{'}_k$ （ $k\times m$ ）。

④ $X^{'}=\Sigma_k^{'}V^{*}$ 即为降到k维特征后的数据（此时为 $k\times m$ ）

2.3.2.2 SVD加速求解特征向量（方法二）

①计算矩阵：

$D=B^{*}B$

②SVD分解加速特征分解：

$D=V\Sigma V^{-1}$ ，V是特征矩阵

当V正交标准化：

V是SVD分解 $B=U\Sigma ^{'}V^{*}$ 中的V

奇异值：

$\hat{\Sigma}=\sqrt{\Sigma}$

则：

$U=B(\hat{\Sigma}V^{*})^{-1}=BV(\hat{\Sigma})^{-1}$

③将奇异值 $\hat{\Sigma}$ 按对应特征值大小从上到下按行排列成矩阵，取出最大的k个特征值对应的特征向量，按列构成矩阵 $\begin{bmatrix} u_1 & \cdots & u_k \end{bmatrix}$ ,标准化后，变成特征向量矩阵 $U_k$ （ $n\times k$ )。

④ $X^{'}=U_k^{*}X$ 即为降到k维特征后的数据（此时为 $k\times m$ ）

2.3.3 降维特征数的选取

有时候，我们不指定降维后的k的值，而是换种方式，指定一个降维到的主成分比重阈值t。这个阈值t在(0,1]之间。假如我们的K个特征值为 $\lambda_1\geqslant \lambda_2\geqslant \cdots \geqslant \lambda_k$ ,则k可以通过下式得到：

$f(k)=\tfrac{\sum_{i=1}^{k}\lambda_i}{\sum_{i=1}^{K}\lambda_i}\geqslant t$

容易得到f(k)为一单调递增的函数，且递增速度随着k增加而逐渐降低。对于SVD分解法，将特征值换成奇异值即可。

2.3.4 总结

对于EVD和SVD两种算法结果，其实有：

$U^{*}X=U^{*}U\Sigma^{'}V^{*}=\Sigma^{'}V^{*}$

它本质上就是奇异值分解：

$X=U\Sigma^{'}V^{*}$

保留k个就有：

$X^{'}=U_k^{*}X=\Sigma_k^{'}V^{*}$

可以看到两种算法效果是一样的，但是SVD的两种实现都比EVD方法更高效：

当输入矩阵为，EVD方法计算 $C=\tfrac{1}{m-1}BB^{*}$ （ $n\times n$ ），SVD方法计算 $D=B^{*}B$ （ $m \times m$ ），在 $n\gg m$ (特征维度大于样本数)时，大大降低了计算复杂度。

2.4.2 实现步骤

为了得到训练集的特征脸，需要用SVD加速得到特征向量的方法来进行PCA人脸识别。（具体步骤见2.3.2.2）

15个人每个人选取4 张图片共60张作为训练样本。

(1)将每一张图像（ $128\times128$ ）写成列向量形式 $x_i(16384\times 1)$ ，并排列成数据矩阵（ReadFace.m)：

$X=\begin{bmatrix} |&|& &| \\ x_1&x_2&\cdots&x_m\\ |&|& &| \end{bmatrix}(m=60)$

function Training_Data = ReadFace(Training_Path)
% ---------- Construct 2D matrix from all of the 1D image vectors in the training data file ------------
flist = dir(strcat(Training_Path,'\*.jpg'));
Training_Data = [];
for imidx = 1:length(flist)  %flist为60，即六十张图片
    fprintf('Constructing Training Image Data Space [%d] \n', imidx); 
    path = strcat(Training_Path,strcat('\',int2str(imidx),'.jpg'));
    img = imread(path);
    [irow icol] = size(img);
    temp = reshape(img',irow*icol,1);   % 将二维图片转换成一维向量
    Training_Data = [Training_Data temp]; %遍历结束为16384*60的矩阵
end
fprintf('\n');

(2)求均值向量，对矩阵进行中心化处理,利用SVD分解加速求A的特征向量U，按奇异值贡献率大于80%选取降维维数，得到最终的特征脸(EigenfaceCore.m)：

function [m, A, Eigenfaces] = EigenfaceCore(Training_Data)

%----------------------------Calculate the mean image ------------------------
% ---------------------compute the covariance matrix --------------------------
m = mean(Training_Data,2);  %列平均

Train_Number = size(Training_Data,2);
temp_m = [];  
%平均矩阵
for i = 1 : Train_Number
    temp_m = [temp_m m];
end
A = double(Training_Data) - temp_m;

%SVD加速得到特征向量U
%-------------------------------------------------
L = A'*A;                %降低计算复杂度
[V,S] = eig(L);          %V是A的右奇异向量
V = orth(V);  %正交化
[S_sort,index] = sort(diag(S),'descend');  %降序排列
S = diag(sqrt(S_sort));  %求奇异值
V = V(:,index);          %对应特征向量
U = A*V*diag(1./diag(S));%利用SVD分解求AA’的特征向量U
s = sum(diag(S))
latent = 100*diag(S)/s;
pareto(latent);          %画方差贡献率图
%奇异值排序并选择贡献率阈值为90%进行降维
U_k = [];
percents = 0;
for k = 1 : size(S,2) 
    if( percents > 90 ) % 设置阈值
        U_k = U(:,1:k);
        k
        break
    end
    percents = percents + 100*S(k,k)/s;
end
Eigenfaces = U_k;  %特征空间

(3)将中心化后的训练人脸和测试人脸都投影到特征脸空间,采用欧式距离计算识别图片与每个人脸的距离 $D_i$ ，找到距离最小的图片(Recognition.m)：

function OutputName = Recognition(TestImage, m, A, Eigenfaces)
%-------------Project the selected test image and all of the training
%图片投影到特征空间
%投影处理后的训练集图片
ProjectedImages = Eigenfaces'*A;
Train_Number = size(A,2);
%投影测试集图片
InputImage = imread(TestImage);
temp = InputImage(:,:,1);

[irow icol] = size(temp)
InImage = reshape(temp',irow*icol,1);
Difference = double(InImage)-m; 
Projected_TestImage = Eigenfaces'*Difference;

%  找到训练图片中与识别图片欧氏距离最小的
Euc_dist = [];
for i = 1 : Train_Number
    q = ProjectedImages(:,i);
    temp = ( norm( Projected_TestImage - q ) )^2;
    Euc_dist = [Euc_dist temp];
end

[Euc_dist_min , Recognized_index] = min(Euc_dist);
OutputName = strcat(int2str(Recognized_index),'.jpg');   %拼接图片名称

(4)其它代码：

Visualize_faces.m:

function Visualize_faces(face,imgrow, imgcol)
%-------------------Show the maxmum nine pictures of Eigenfaces---------------
    Num_Eigenvalue = size(face,2);
    img = zeros(imgrow, imgcol);
    for i=1:min(Num_Eigenvalue,9)
        img(:) = face(:,i);
        subplot(3,3,i);
        imshow(img',[]);
    end
end

main.m:

clear,clc;

Training_Path = 'D:\matlab\bin\Face-Recognition-Using-PCA-master\Face-Recognition-Using-PCA-master\TrainDatabase';  %Set your directory for training data file
Testing_Path = 'D:\matlab\bin\Face-Recognition-Using-PCA-master\Face-Recognition-Using-PCA-master\TestDatabase';    %Set your directory for testing data file

disp('Pick a Testing Photo From TestDatabase please')
[filename, pathname] = uigetfile({'*.jpg'},'Pick a Testing Photo From TestDatabase please');   %选择测试的识别图片
disp('Hold a second for computing')
TestImage = [pathname, filename];
im = imread(TestImage);
[imgrow, imgcol] = size(im);

Training_Data = ReadFace(Training_Path);
[m, A, Eigenfaces] = EigenfaceCore(Training_Data);  %得到平均脸，处理后的数据，特征脸
OutputName = Recognition(TestImage, m, A, Eigenfaces);
%可视化
figure('Name','Eigenface')
Visualize_faces(Eigenfaces, imgrow, imgcol)
figure('Name','Meanface')
Visualize_faces(m, imgrow, imgcol)
SelectedImage = strcat(Training_Path,'\',OutputName);
SelectedImage = imread(SelectedImage);

figure('name','Recognition Result')
subplot(1,2,1);
imshow(im)
title('Test Image');

subplot(1,2,2);
imshow(SelectedImage);
title('Recognition Result');
disp('Done')

2.4.3 结果

方差贡献率的帕累托图（只显示前10）：

平均脸：

特征脸（只画前9个）：

识别结果：

直接用特征值分解得到特征脸Matlab计算了十分钟远没算完；用SVD加速得到特征脸只用了0.01s，效率大大提升。

3 RPCA算法

PCA可以说是目前应用最广泛的数据分析和降维的统计学工具。然而实际应用中，当数据错误或者缺失时，PCA不能很好的抓住数据的真实子空间结构，因此效果比较差，特别是错误或缺失较大时，效果更差。不幸的是，在图像处理、web数据分析和生物信息学等现代应用程序中，严重的数据缺失无处不在，其中一些测量数据可能会被任意破坏(由于闭塞、恶意篡改或传感器故障)，或者与我们试图识别的低维结构完全无关。

主成分分析(PCA) 非常容易受到异常值和损坏数据的影响，在异常值方面很脆弱。为了改善这种敏感性，cand等人开发了一种鲁棒主成分分析(RPCA)，试图将数据矩阵X分解为结构化低秩矩阵L和包含异常值和损坏数据的稀疏矩阵S:

$X=L+S$

L是低秩的，S是稀疏的扰动。

PCA和RPCA两者的区别在于对于误差的假设不同，PCA假设数据误差是服从高斯分布的，即数据噪声较小，因此PCA的S是小的扰动矩阵；RPCA假设数据噪声是稀疏的，并且可能是强的噪声，它的S可以具有任意大的量级。

L的主成分对S中的异常值和损坏数据具有鲁棒性。

这种分解对许多感兴趣的现代问题具有深远的影响，包括视频监控(背景对象出现在L中，前景对象出现在S中)，人脸识别(特征脸在L中，阴影，遮挡等在S中)，自然语言处理和潜在语义索引等等。

3.1 原理和实现

数学上，目标是找到满足以下条件的L和S:

$min_{L,S}rank(L)+||S||_0$ subject to $L+S=X$

然而，rank(L)和 $||S||_0$ 项都不是凸项，这不是一个可处理的优化问题。因此我们用L矩阵的核范数来近似L矩阵的秩，核范数指的是L的奇异值值和，所以当核范数越小时，可近似认为L的秩越低，而用S的1范数来近似矩阵S的0范数，是因为1范数指的是当矩阵中某一列的的元素的绝对值的和取最大的时候，这个某列的元素的绝对值的和则为1范数，那么如果这个和已经很小了，趋向于0，也就说明其他列的元素的绝对值的和也很小，也就可以近似认为非0元素的个数很少，所以问题就转换为这个式子了与压缩感知问题类似，可以使用的凸松弛来求解高概率的最优L和S:

$min_{L,S}||L||_*+\lambda||S||_1$ subject to $L+S=X$

值得注意的是，若X的维数为 $n \times m$ 当 $\lambda=\tfrac{1}{\sqrt{max(n,m)}}$ 时，式子的解更大可能收敛。

上面式子的凸问题被称为主成分追踪(PCP)，可以使用增广拉格朗日乘子(ALM)算法来解决。具体地说，增广拉格朗日可以构造为:

$\Gamma (L,S,Y)=||L||_*+\lambda||S||_1+<Y,X-L-S>$

Y是指拉格朗日乘子。

将有约束问题转化为无约束，这需要付出代价，所以我们加入了惩罚项：

$\Gamma (L,S,Y)=||L||_*+\lambda||S||_1+<Y,X-L-S>+\tfrac{u}{2}||X-L-S||^{2}_F$

矩阵A的Frobenius范数定义为矩阵A各项元素的绝对值平方的总和

一个通解将求使 $\Gamma$ 最小的 $L_k$ 和 $S_k$ ，更新拉格朗日乘子 $Y_{k+1}=Y_k+u(X-L_k-S_k)$ ，并迭代直到解收敛。

然而，对于这个特定的系统，交替方向法(ADM)提供了一个简单的过程来找到L和S：

①收缩操运算为(shrink.m)：

1.	function out = shrink(X,tau)
2.	    out = sign(X).*max(abs(X)-tau,0);
3.	end

②构造奇异值阈值运算(SVT.m) ：

1.	function out = SVT(X,tau)
2.	    [U,S,V] = svd(X,'econ');
3.	    out = U*shrink(S,tau)*V';
4.	end

③可以使用和SVT运算迭代求解L和S（使用交替方向法ADM）(RPCA.m):

1.	function [L,S] = RPCA(X)
2.	[n1,n2] = size(X); mu = n1*n2/(4*sum(abs(X(:))));%此时解更有可能收敛
3.	lambda = 1/sqrt(max(n1,n2));
4.	thresh = 1e-7*norm(X,'fro');%收敛阈值
5.	L = zeros(size(X)); S = zeros(size(X)); Y = zeros(size(X));
6.	count = 0;
7.	%ADM收敛
8.	while((norm(X-L-S,'fro')>thresh)&&(count<1000)) 
9.	    L = SVT(X-S+(1/mu)*Y,1/mu);
10.	    S = shrink(X-L+(1/mu)*Y,lambda/mu);
11.	    Y = Y + mu*(X-L-S);
12.	    count = count + 1
13.	end

3.2 结果

我使用一张我用激光雷达扫描的二维图片来进行处理。

主函数(main.m)：

1.	c;
2.	X = double(imread('B315.png'))
3.	[L,S] = RPCA(X)
4.	figure(1);
5.	imshow(L/256);
6.	figure(2);
7.	imshow(S/256);

输入X：rplindar雷达扫描的图像

输出：

S：

本例中，在低秩分量L中，噪点被去除并用特征面中最一致的低秩特征填充。

我们可以进一步用得到的L进行Gazebo建模：

工程分享：
通过网盘分享的文件：PCA_RPCA.zip
链接: https://pan.baidu.com/s/1lofH1Z5VKNiwvecNvwVgnA?pwd=sbk4 提取码: sbk4

新出炉的烤鸡

关注

5
点赞
踩
21

收藏

觉得还不错? 一键收藏
3
评论
复制链接

分享到 QQ

分享到新浪微博

扫一扫