RPCA_1

低秩矩阵分解（PRCA）

数学模型

首先给出低秩矩阵分解数学模型
$$\underset{A,E}{\min }{\Vert A\Vert }_{\ast }+{\Vert E\Vert }_{0}s.t.D=A+E$$
采用增广拉格朗日乘子法求解
$$L(A,E,Y,\mu )={\Vert A\Vert }_{\ast }+\lambda {\Vert E\Vert }_1+⟨Y,D-A-E⟩+\frac{\mu }{2}{\Vert D-A-E\Vert }_F^2$$
因此，原问题就变为求解
$$\underset{A,E}{\min }L(A,E,Y,\mu )$$
很容易解得
$$\hat{E}=\underset{E}{\mathrm{argmin}}\left\{\lambda {\Vert E\Vert }_1+\frac{\mu }{2}{\Vert E-(D-A+\frac{Y}{\mu })\Vert }_F^2\right\}$$
$A$的解就不赘述

Matlab采用RPCA实现视频背景建模

采用一共120帧的matlab自带视频，其中每一帧的图片大小为$120\times 160$，选取其中的第110帧，结果如下：



Matlab代码如下：

%% this to read avi by using mmread to get every frame
video = VideoReader('traffic.avi');
nFrames = video.NumberOfFrames;   %得到帧数
H = video.Height;     %得到高度
W = video.Width;      %得到宽度
Rate = video.FrameRate;
Cal_FrameN=nFrames;
% Preallocate movie structure.

 mov(1:Cal_FrameN) = struct('cdata',zeros(H,W,3,'uint8'),'colormap',[]);
 
%read one frame every time
for i = 1:Cal_FrameN
    mov(i).cdata = read(video,i);
    P = mov(i).cdata;
    double_P=im2double(P);
    %imshow( double_P),title('原始图片3');
    Frame=rgb2gray( double_P);%当前灰度帧
    %下采样
    %Down_Frame=Frame(2:4:H,2:4:W);
    M(:, i) = Frame(:);
   %  disp('当前播帧数：'),disp(i);
    % imshow(P),title('原始图片');
     
 
%     P2=rgb2gray(P);
 
 
end

 
 %[A_hat, E_hat]=inexact_alm_rpca(M);
 [A_hat, E_hat]=rpca_myself(M);
% A=reshape(A_hat(:,120), size(Down_Frame));
% E=reshape(E_hat(:,120), size(Down_Frame));
% D_real=reshape(M(:,120), size(Down_Frame));
D=reshape(M(:,110),120,160);
A=reshape(A_hat(:,110),120,160);
E=reshape(E_hat(:,110),120,160);
subplot(1,3,1)
imshow(D),title('原始图片');
subplot(1,3,2)
imshow(A),title('背景图片');
subplot(1,3,3)
imshow(E),title('运动项');

主要调用的函数为 rpcamyself.m

function [A_hat,E_hat]=rpca_myself(D,lambda,tol)
[m, n] = size(D);
if nargin < 2
    lambda = 1 / sqrt(m);
end

if nargin < 3
    tol = 1e-7;
elseif tol == -1
    tol = 1e-7;
end
maxIter = 1000;


% initialize
Y = D;
%norm_two = lansvd(Y, 1, 'L');
norm_two = svd(Y);
norm_two=norm_two(1);
norm_inf = norm( Y(:), inf) / lambda;
dual_norm = max(norm_two, norm_inf);
Y = Y / dual_norm;

A_hat = zeros( m, n);
E_hat = zeros( m, n);
mu = 1.25/norm_two % this one can be tuned
mu_bar = mu * 1e7
rho = 1.5          % this one can be tuned
d_norm = norm(D, 'fro');

iter = 0;
% total_svd = 0;
converged = false;
stopCriterion = 1;
% sv = 10;
while ~converged
    iter = iter + 1
    
    temp_T = D - A_hat + (1/mu)*Y;
    E_hat = max(temp_T - lambda/mu, 0);
    E_hat = E_hat+min(temp_T + lambda/mu, 0);
    [U, S, V] = svd(D - E_hat + (1/mu)*Y, 'econ');
    diagS = diag(S);
    svp = length(find(diagS > 1/mu));%rank of matrix D - E_hat + (1/mu)*Y
    A_hat = U(:, 1:svp) * diag(diagS(1:svp) - 1/mu) * V(:, 1:svp)';

    
    Z=D-A_hat-E_hat;
    
    Y = Y + mu*Z;
    mu = min(mu*rho, mu_bar);
    stopCriterion = norm(Z, 'fro') / d_norm;
    if stopCriterion < tol
        converged = true;
    end
    

    if ~converged && iter >= maxIter
        disp('Maximum iterations reached') ;
        converged = 1 ;
    end
end


end