目录
1 问题描述
利用PCA或者RPCA技术实现对数据的提取或处理(数据类型不限制,例如可以是图像),给出详细过程和可视化结果。作业内容可以包括但不限于原理介绍、数据处理过程、代码(关键部分即可)、结果展示等。
2 PCA算法(参考PCA论文和PCA数学原理博客介绍)
一般我们提到降维最容易想到的算法就是主成分分析,下面我们就对其原理做一个总结。
主成分分析(Principal components analysis,以下简称PCA)主要有纯特征值分解(EVD)和奇异值分解(SVD)两种方法,它提供了一个数据驱动的分层坐标系统来表示高维相关数据。PCA是最重要的降维方法之一,在数据压缩消除冗余和数据噪音消除等领域都有广泛的应用,现在也广泛用于深度学习的数据白化操作中。
我们不仅要对EVD方法有了解,并且要会运用SVD方法进行运算减负。实际上SVD方法在实际生活中更常用,我也会重点介绍SVD方法。
一般来说,我们感兴趣的是分析一个大数据集 :
列可能是一次模拟或实验的所有测量值,即一个样本。
指标k是表示第k个不同的测量集的标签。X将由数据的时间序列组成,并且 。通常状态维度n非常大,在数百万或数十亿个自由度的数量级上。列通常被称为采样点数(snapshots),m是X中的采样样本数量。在许多系统中 ,得到一个又高又细的矩阵。通常,在单个实验中收集许多测量值,测量值相当于变量,并将这些测量值排列成行向量,即特征。测量值可以是可观察到的特征,例如特定人类个体的人口统计学特征。
PCA降维问题的优化目标:将一组n维向量降为k维(k大于0,小于n),其目标是选择k个单位(模为1)正交基,使得原始数据变换到这组基上后,各字段两两间协方差为0,而字段的方差则尽可能大(在正交的约束下,取最大的k个方差)。
而协方差矩阵统一了变量(特征)内方差及变量间协方差关系:
因此我们的很多运算是针对协方差进行的。
2.1 纯特征值分解(EVD)原理
如果一个向量v是矩阵A的特征向量,将一定可以表示成下面的形式:
其中,λ是特征向量v对应的特征值,一个矩阵的一组特征向量是一组正交向量。
将所有特征值λ排列成对角矩阵D;对应特征向量v排列成矩阵V,正交单位化后,也可写成:
对于矩阵A,有一组特征向量v,将这组向量进行正交化单位化,就能得到一组正交单位向量。特征值分解,就是将矩阵A分解为如下式:
其中,V是矩阵A的特征向量组成的矩阵,则是一个对角阵,对角线上的元素就是特征值。
在线性代数中,对称矩阵一定可以对角化。
2.2 奇异值分解(SVD)原理
SVD是特征分解的推广,是存在于每一个复值矩阵的唯一矩阵分解:
其中,和是具有正交列的酉矩阵,是一个在对角线上有实数、非负元素且对角线外有零元素的矩阵。*表示复共轭转置。
当 ,矩阵在对角线上最多有m个非零元素,因此,可以使用经济型奇异值分解(economy SVD)精确地表示X:
的列向量张成一个向量空间,它与张成的向量空间互补且正交。U的列称为X的左奇异向量,V的列称为右奇异向量。的对角线元素称为奇异值,它们按从大到小的顺序排列。X的秩等于非零奇异值的个数。
SVD与涉及相关矩阵和的特征值问题密切相关,我们发现:
由于U和V是酉矩阵,U, ,V是以下特征值问题的解:
换句话说,X的每个非零奇异值是和的特征值的正平方根,它们具有相同的非零特征值。
这提供了SVD的直观解释,其中U的列是相关矩阵的特征向量,V的列是的特征向量。我们选择按大小降序排列奇异值,因此U的列是根据它们在X的列中捕获的相关性进行分层排序的;V同样捕获了X行中的相关性。
2.3 PCA步骤
假设有m条n维数据。
输入:维样本集 ,m个样本要降维从特征维数n降为k。
输出:降维后的样本集 。
为了保持不同算法所得结果一致,首先对数据进行中心化:
现在我们计算平均列(每一行平均,即一个变量零均值化):
则平均矩阵为:
处理后的矩阵为B( ):
2.3.1 EVD实现PCA
①计算协方差矩阵
B的行的协方差矩阵C(度量特征间的相关性):
(特征维度作内积,衡量相关性)
②计算特征值和特征向量
通过计算C的特征分解可以得到主成分:
它是肯定存在的,因为C是Hermite矩阵,是一个对称矩阵。
③将特征向量按对应特征值大小从上到下按行排列成矩阵,取出最大的k个特征值对应的特征向量,按列构成矩阵 ,标准化后,变成特征向量矩阵()。
④ 即为降到k维特征后的数据(此时为)。
2.3.2 SVD实现PCA
2.3.2.1 SVD表示结果(方法一)
①计算矩阵:
②特征分解:
,V是特征矩阵
当V正交标准化:
V是SVD分解中的V
奇异值:
③将奇异按对应特征值从大到小排列,取出最大的k个奇异值,构成矩阵()。
④即为降到k维特征后的数据(此时为)
2.3.2.2 SVD加速求解特征向量(方法二)
①计算矩阵:
②SVD分解加速特征分解:
,V是特征矩阵
当V正交标准化:
V是SVD分解中的V
奇异值:
则:
③将奇异值按对应特征值大小从上到下按行排列成矩阵,取出最大的k个特征值对应的特征向量,按列构成矩阵,标准化后,变成特征向量矩阵 ()。
④ 即为降到k维特征后的数据(此时为)
2.3.3 降维特征数的选取
有时候,我们不指定降维后的k的值,而是换种方式,指定一个降维到的主成分比重阈值t。这个阈值t在(0,1]之间。假如我们的K个特征值为,则k可以通过下式得到:
容易得到f(k)为一单调递增的函数,且递增速度随着k增加而逐渐降低。对于SVD分解法,将特征值换成奇异值即可。
2.3.4 总结
对于EVD和SVD两种算法结果,其实有:
它本质上就是奇异值分解:
保留k个就有:
可以看到两种算法效果是一样的,但是SVD的两种实现都比EVD方法更高效:
当输入矩阵为 ,EVD方法计算 (),SVD方法计算(),在 (特征维度大于样本数)时,大大降低了计算复杂度。
2.4.2 实现步骤
为了得到训练集的特征脸,需要用SVD加速得到特征向量的方法来进行PCA人脸识别。(具体步骤见2.3.2.2)
15个人每个人选取4 张图片共60张作为训练样本。
(1)将每一张图像()写成列向量形式 ,并排列成数据矩阵(ReadFace.m):
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)将中心化后的训练人脸和测试人脸都投影到特征脸空间,采用欧式距离计算识别图片与每个人脸的距离,找到距离最小的图片(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:
L是低秩的,S是稀疏的扰动。
PCA和RPCA两者的区别在于对于误差的假设不同,PCA假设数据误差是服从高斯分布的,即数据噪声较小,因此PCA的S是小的扰动矩阵;RPCA假设数据噪声是稀疏的,并且可能是强的噪声,它的S可以具有任意大的量级。
L的主成分对S中的异常值和损坏数据具有鲁棒性。
这种分解对许多感兴趣的现代问题具有深远的影响,包括视频监控(背景对象出现在L中,前景对象出现在S中),人脸识别(特征脸在L中,阴影,遮挡等在S中),自然语言处理和潜在语义索引等等。
3.1 原理和实现
数学上,目标是找到满足以下条件的L和S:
subject to
然而,rank(L)和项都不是凸项,这不是一个可处理的优化问题。因此我们用L矩阵的核范数来近似L矩阵的秩,核范数指的是L的奇异值值和,所以当核范数越小时,可近似认为L的秩越低,而用S的1范数来近似矩阵S的0范数,是因为1范数指的是当矩阵中某一列的的元素的绝对值的和取最大的时候,这个某列的元素的绝对值的和则为1范数,那么如果这个和已经很小了,趋向于0,也就说明其他列的元素的绝对值的和也很小,也就可以近似认为非0元素的个数很少,所以问题就转换为这个式子了与压缩感知问题类似,可以使用的凸松弛来求解高概率的最优L和S:
subject to
值得注意的是,若X的维数为 当时,式子的解更大可能收敛。
上面式子的凸问题被称为主成分追踪(PCP),可以使用增广拉格朗日乘子(ALM)算法来解决。具体地说,增广拉格朗日可以构造为:
Y是指拉格朗日乘子。
将有约束问题转化为无约束,这需要付出代价,所以我们加入了惩罚项:
矩阵A的Frobenius范数定义为矩阵A各项元素的绝对值平方的总和
一个通解将求使最小的和,更新拉格朗日乘子,并迭代直到解收敛。
然而,对于这个特定的系统,交替方向法(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雷达扫描的图像
输出:
L:
S:
本例中,在低秩分量L中,噪点被去除并用特征面中最一致的低秩特征填充。
我们可以进一步用得到的L进行Gazebo建模:
工程分享:
通过网盘分享的文件:PCA_RPCA.zip
链接: https://pan.baidu.com/s/1lofH1Z5VKNiwvecNvwVgnA?pwd=sbk4 提取码: sbk4