Robust principal component analysis?(RPCA简单理解)

参考文献：Candès, E.J., Li, X., Ma, Y., and Wright, J.: ‘Robust principal component analysis?’, J. ACM, 2011, 58, (3), pp. 11

作者主页有很多关于low-rank的代码：http://perception.csl.illinois.edu/matrix-rank/sample_code.html

主要算法公式如下：

关于这个低秩分解的代码，其实相当简单：

例如：公式（5-2）对应的代码为：

S = wthresh(M-L+Y/miu,'s',lambda/miu);

公式（5-3）对应的代码：

% ---- update L --- % 
[U,D,V] = svd(M-S+Y/miu,'econ'); 
D1 = wthresh(D,'s',1/miu);    
L = U*D1*V';

我自己根据公式编写的代码如下：

% 求解 argmin rank(L) + ||S||_1  s.t. M = L+S
% Reference: Wright, J., Ganesh, A., Rao, S., Peng, Y. and Ma, Y. (2009) 
% Robust principal component analysis: Exact recovery of corrupted low-rank matrices via convex optimization. 
% In: Proceedings of Advances in neural information processing systems. 2080-2088.
L = zeros(size(M));
S = L;
Y = L;
norm_two = lansvd(M, 1, 'L'); % computes the K largest singular values.
miu = 1.25/norm_two;
% miu = 0.1;
max_miu = 1e7;
lambda = 0.005;
rho = 1.5;
max_iter = 200;
for iter = 1:max_iter
    % --- update S ----%
    S = wthresh(M-L+Y/miu,'s',lambda/miu);
    % ---- update L --- %
    [U,D,V] = svd(M-S+Y/miu,'econ');
    D1 = wthresh(D,'s',1/miu);    
    L = U*D1*V';
    Y = Y+miu*(M-L-S);
    miu = min(max_miu,rho*miu);
    obj(iter) = norm(M-L-S,'fro')^2;
    if iter > 2 && abs(obj(iter) - obj(iter-1)) < 1e-7
        break;
    end
end
figure;
imshow(reshape(L(:,80),50,40),[]);  title('低秩部分');

figure;
imshow(reshape(S(:,80),50,40),[]); title('稀疏部分');
figure;
imshow(reshape(M(:,80),50,40),[]); title('原始图像');

对于AR数据库，我们对遮挡脸进行了试验，结果如下：

而同时，我们调用作者主页编写的权威代码得到的结果为：（这里原始文章的代码下载地址为：链接: http://pan.baidu.com/s/1i5m4QrV 密码: fduh）

 

这里我们小小总结下调这篇文章算法参数的心得：

1、lambda值越大时，对于lambda约束的矩阵，其值就越小，几乎为0，矩阵越稀疏，甚至与稀疏到0~~；相反的，另外的那个矩阵则占据主导位置。

     所以，在作者原始代码中，lambda=1/sqrt(m), m为图像特征维数。经验：lambda常可取值[0.001,0.01]；

2、关于miu这个值，作者用的是miu=1.25/S_L, S_L表示原始矩阵M的最大特征值。miu一般可取0.15等

3、总感觉尽管低秩矩阵可以去掉遮挡，但是是以丢失细节特征为代价。

4、一般而言，遮挡矩阵M中，若全为遮挡脸，那么恢复效果一般比较差，而当存在一些未遮挡脸时，恢复效果会比较好！！