-
Notifications
You must be signed in to change notification settings - Fork 8
/
Offline_ReProCS_PCA.m
135 lines (116 loc) · 5.04 KB
/
Offline_ReProCS_PCA.m
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
function [L_hat_offline, P_hat, S_hat_offline, T_hat, t_hat, ...
P_track_full, P_track_new, T_calc]= Offline_ReProCS_PCA(M, ...
P_init, ev_thresh, alpha, K, omega, outc, K_CS)
%%%This is the automatic version of ReProCS algorithm under the latest
%%%subspace change and support model.
%This folder contains the code accompanying pre-print.
%
%[1] "Provable Dynamic Robust PCA or Robust Subspace Tracking", Praneeth Narayanamurthy and Namrata Vaswani, IEEE Trans. Info. Theory., 2019.
%
%If you use this code please also cite the following papers
%[2] "An online algorithm for separating sparse and low-dimensional signal sequences from their sum", Han Guo, Chenlu Qiu, and Namrata Vaswani, IEEE Trans. Sig. Proc., 2014.
%[3] "Recursive Robust PCA or Recursive Sparse Recovery in Large but Structure Noise", Chenlu Qiu, Namrata Vaswani, Brain Lois, and Leslie Hogben, IEEE Trans. Info. Theory., 2014.
%[4] "Real-time Robust Principal Components' Pursuit", Chenlu Qiu, and Namrata Vaswani, Allerton, 2010.
%%% Inputs %%%
%%% M - measurement matrix %%%
%%% ev_thres - threshold for subspace change detection %%%
%%% P_init - an initial estimate of the subspace %%%
%%% t_train - the dimension of the training data %%%
%%% Algorithm parameters %%%
%%% alpha - frame length %%%
%%% K - number of projection PCA steps %%%
%%% omega - threshold for non-zero value in S %%%
%%% K_CS - number of CoSaMP iterations %%%
%%% outc - an upper bound on estimate of fraction of outliers per column
%%% Outputs %%%
%%% L_hat - estimate of the low rank matrix %%%
%%% P_hat - estimate of the subspace in which L lies %%%
%%% S_hat - estimate of the sparse signal %%%
%%% t_hat - estimate of subspace change times %%%
%% Initializations
%thresh = ev_thresh / 2;
[~, r_init] = size(P_init);
P_hat_old = P_init;
P_hat_new = [];
P_hat = [P_hat_old, P_hat_new];
[n, t_max] = size(M);
T_hat = zeros(n, t_max);
S_hat = zeros(n, t_max);
L_hat = zeros(n, t_max);
S_hat_offline = zeros(n, t_max);
L_hat_offline = zeros(n, t_max);
t_hat = [1];
% L_hat(:, 1 : t_train) = M(:, 1 : t_train);
k = 0;
cnt = 1;
phi_t = (eye(n) - P_hat * P_hat');
% In = eye(n);
ph = 1; %ph - 0 => detect, 1 => ppca
%% Main Algorithm Loop
for ii = 1 : t_max
%% Estimate support
y_t = M(:, ii) - (P_hat * (P_hat' * M(:, ii)));
DecayRate = 0.9; %other values work, may make it slower
x_t_hat_cs = cosamp_cgls(phi_t, ...
y_t, outc, DecayRate, K_CS, 1e-6);
t_hat_temp = find(abs(x_t_hat_cs) > omega);
% T_hat(t_hat_temp, ii) = 1;
%% Estimate signal components
% % [S_hat(t_hat_temp, ii), ~] = ...
% % lsqr(phi_t(:, t_hat_temp), y_t, 1e-6, 50);
% S_hat(t_hat_temp, ii) = phi_t(:, t_hat_temp) \ y_t;
S_hat(t_hat_temp, ii) = cgls(phi_t(:, t_hat_temp), y_t, ...
0, 1e-10, 10);
L_hat(:, ii) = M(:, ii) - S_hat(:, ii);
%% Subspace update
if(~mod(ii + 1 , alpha))
u = (ii + 1) / alpha;
idx = (u-1) * alpha + 1 : u * alpha ;
L_temp = L_hat(:, idx);
MM = L_temp - (P_hat_old *(P_hat_old' * L_temp));
if(~ph) %%detect phase
% phi_t = eye(n) - P_hat * P_hat';
if(svds(MM, 1) >= sqrt(alpha * ev_thresh))
ph = 1;
t_hat = [t_hat, ii];
k = 0;
end
else %%ppca phase
P_hat = simpleEVD((L_hat(:, max(1, ii - alpha + 1) : ii)), r_init);
phi_t = speye(n) - P_hat * P_hat';
k = k + 1;
if(k==K + 1)
P_hat_old = P_hat;
if(length(t_hat) == 1)
idx_offline = t_hat(end) : ii;
else
idx_offline = t_hat(end - 1) + (K) * alpha : ii;
end
%idx_offline = t_hat(end) : ii;
for kk = idx_offline
yy = M(:, kk) - (P_hat * (P_hat' * M(:, kk)));
DecayRate = 0.9;
x_t_hat_cs = cosamp_cgls(phi_t, ...
yy, outc, DecayRate, K_CS, 1e-6);
t_hat_temp = find(abs(x_t_hat_cs) > omega);
S_hat_offline(t_hat_temp, kk) = ...
cgls(phi_t(:, t_hat_temp), yy, ...
0, 1e-10, 10);
L_hat_offline(:, kk) = M(:, kk) - ...
S_hat_offline(:, kk);
end
k = 1;
ph = 0;
phi_t = speye(n) - P_hat * P_hat';
end
end
end
%% Return subspace
if((ii == 1) || ~(mod(ii + 1, alpha)))
P_track_full{cnt} = P_hat;
P_track_new{cnt} = P_hat_new;
T_calc(cnt) = ii;
cnt = cnt + 1;
end
end
end