RPCA

RPCA

低秩矩阵分解（PRCA）

数学模型

首先给出低秩矩阵分解数学模型
min ⁡ A , E ∥ A ∥ ∗ + ∥ E ∥ 0 s . t . D = A + E \min_{A,E} \left\|A\right\|_* +\left\|E\right\|_0\quad s.t.D=A+E A,Emin​∥A∥∗​+∥E∥0​s.t.D=A+E
采用增广拉格朗日乘子法求解
L ( A , E , Y , μ ) = ∥ A ∥ ∗ + λ ∥ E ∥ 1 + ⟨ Y , D − A − E ⟩ + μ 2 ∥ D − A − E ∥ F 2 L(A,E,Y,\mu)=\left\|A\right\|_*+\lambda \left\|E\right\|_1\ +\left \langle Y,D-A-E\right \rangle+\frac{\mu}{2}\left\|D-A-E\right\|_F^2 L(A,E,Y,μ)=∥A∥∗​+λ∥E∥1​ +⟨Y,D−A−E⟩+2μ​∥D−A−E∥F2​
因此，原问题就变为求解
min ⁡ A , E L ( A , E , Y , μ ) \min_{A,E}L(A,E,Y,\mu) A,Emin​L(A,E,Y,μ)
很容易解得
E ^ = arg min ⁡ E { λ ∥ E ∥ 1 + μ 2 ∥ E − ( D − A + Y μ ) ∥ F 2 } \hat{E}=\argmin_E\left \{ \lambda \left\|E\right\|_1+\frac{\mu}{2}\left\|E-(D-A+\frac{Y}{\mu})\right\|_F^2\right \} E^=Eargmin​{λ∥E∥1​+2μ​∥∥∥∥​E−(D−A+μY​)∥∥∥∥​F2​}
A A A的解就不赘述

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

采用一共120帧的matlab自带视频，其中每一帧的图片大小为 120 × 160 120\times160 120×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_FrameNmov(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 < 2lambda = 1 / sqrt(m);
endif nargin < 3tol = 1e-7;
elseif tol == -1tol = 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 ~convergediter = iter + 1temp_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)*YA_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 < tolconverged = true;endif ~converged && iter >= maxIterdisp('Maximum iterations reached') ;converged = 1 ;end
endend