-
Notifications
You must be signed in to change notification settings - Fork 7
/
LRPR_track_large.m
133 lines (126 loc) · 5.62 KB
/
LRPR_track_large.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
function [B_hat, Uo, Uo_track, t_calc] = LRPR_track_large(Params, Paramsrwf, Y, Ysqrt, A, X)
[~,k_max] = size(Y);
det_mode = 0; %0 means not-detect and 1 means detect
khat = [1];
l=0;
ctr = 1;
t_calc = [];
% %need to define alpha, L, thresh
for kk = 2 : k_max
if((~mod(kk-1, Params.alpha)) &&(kk == khat(end) + Params.alpha) &&(det_mode == 0)) %% this means alpha frames after previous dectection
[~,Y_init,Ai] = ...
Generate_Mes_temp(X(:, kk - Params.alpha:kk), Params,Params.m_init, Params.alpha);
Yu = zeros(Params.n, Params.n);
for nh = 1 : Params.alpha
normest = sqrt((9/Params.m) * Y_init(:,nh)' * Y_init(:, nh));
Ytr = Y_init(:,nh) .* (abs(Y_init(:, nh)) > normest);
Yu = Yu + Ai(:,:,nh) * diag(Ytr) * Ai(:,:,nh)';
end
Yu = Yu / Params.q / Params.m;
[U_hat,~,~] = svds(Yu, Params.r);
Uo = U_hat;
Uo_track{ctr} = Uo;
t_calc = [t_calc, kk];
ctr = ctr + 1;
%%estimate the coefficients bhat's
[Ysqrt1,~,Ab] = ...
Generate_Mes_temp(X(:, kk - Params.alpha:kk), Params,Params.m_u, Params.alpha);
[~,~,Au] = ...
Generate_Mes_temp(X(:, kk - Params.alpha:kk), Params,Params.m_b, Params.alpha);
B_hat = zeros(Params.r, Params.alpha);
Chat = zeros(Params.m, Params.alpha);
for ni = 1 : Params.alpha
Amatrix = Ab(:,:,ni)' * Uo;% Design matrices for coefficients
A1 = @(I) Amatrix * I;
At = @(Y) Amatrix' * Y;
Paramsrwf.Tb_LRPRnew = Params.Tb_LRPRnew(1); %%hacky soln -- T_a iterations of RWF!
[bhat] = RWFsimple(Ysqrt1(:,ni), Paramsrwf, A1, At);
B_hat(:,ni) = bhat;
x_k = Uo * B_hat(:,ni);
Chat(:, ni) = (Au(:,:,ni)'* x_k >= 0) - (Au(:,:,ni)'* x_k < 0);
end
%fprintf('subspace initialized at %d and ctr is %d\n', kk, ctr);
elseif((~mod(kk-1, Params.alpha)) && (kk ~= khat(end) + Params.alpha) &&(det_mode == 0)) %%not initializing
l = l + 1;
[Ysqrt1,~,Ab] = ...
Generate_Mes_temp(X(:, kk - Params.alpha:kk), Params,Params.m_b, Params.alpha);
[Ysqrt_u,~,Au] = ...
Generate_Mes_temp(X(:, kk - Params.alpha:kk), Params,Params.m_u, Params.alpha);
B_hat = zeros(Params.r, Params.alpha);
Chat = zeros(Params.m, Params.alpha);
for ni = 1 : Params.alpha
Amatrix = Ab(:,:,ni)' * Uo;% Design matrices for coefficients
A1 = @(I) Amatrix * I;
At = @(Y) Amatrix' * Y;
Paramsrwf.Tb_LRPRnew = Params.Tb_LRPRnew(1); %%hacky soln -- T_a iterations of RWF!
[bhat] = RWFsimple(Ysqrt1(:,ni), Paramsrwf, A1, At);
B_hat(:,ni) = bhat;
x_k = Uo * B_hat(:,ni);
Chat(:, ni) = (Au(:,:,ni)'* x_k >= 0) - (Au(:,:,ni)'* x_k < 0);
end
[Qb,~] = qr(B_hat');
Bo = Qb(:,1:Params.r)';
%[Ysqrt_u,~,Au] = Generate_Mes(X,Params,Params.m_u);
Zvec = zeros(Params.m*Params.alpha, 1);
for nt = 1 : Params.alpha
strt_idx = Params.m*(nt-1) + 1;
end_idx = strt_idx + Params.m - 1;
TempVec = Chat(:, nt) .* Ysqrt_u(:,nt);
Zvec(strt_idx:end_idx, 1) = TempVec;
end
Uvec = cgls_new(@mult_H, @mult_Ht , Zvec, 0, 1e-16,30);
U_hat = reshape(Uvec, Params.n, Params.r);
[Qu,~] = qr(U_hat);
Uo = Qu(:, 1:Params.r);
Uo_track{ctr} = Uo;
t_calc = [t_calc, kk];
ctr = ctr+1;
%fprintf('subspace updated at %d and ctr is %d\n', kk, ctr);
%X_hat = Uo * B_hat;
end
if((~mod(kk-1, Params.alpha)) &&(l >= Params.L)) %%means epsilon accurate subspace is obtained and need to check for ss change
%fprintf('checking for change at %d\n', kk)
det_mode = 1;
Uo_track{ctr} = Uo;
t_calc = [t_calc, kk];
ctr = ctr + 1;
end
if((~mod(kk-1, Params.alpha)) && (det_mode == 1))
[~,Y_init,Ai] = ...
Generate_Mes_temp(X(:, kk - Params.alpha:kk), Params,Params.m_init, Params.alpha);
Yu = zeros(Params.n, Params.n);
for nh = 1 : Params.alpha
normest = sqrt((9/Params.m_init) * Y_init(:,nh)' * Y_init(:, nh));
Ytr = Y_init(:,nh) .* (abs(Y_init(:, nh)) > normest);
Yu = Yu + Ai(:,:,nh) * diag(Ytr) * Ai(:,:,nh)';
end
Yu = Yu / Params.alpha / Params.m_init;
Y_u_det = Yu - 2 * (Uo * (Uo' * Yu)) + Uo * (Uo' * ((Yu * Uo) * Uo')) ;
[~,sig_det] = svds(Y_u_det);
sig_det = diag(sig_det);
fprintf('the detection criterion value at %d is %f, \t %f\n',kk, ...
max(sig_det) - min(sig_det), Params.thresh );
if(max(sig_det) - min(sig_det) >= Params.thresh )
det_mode = 0;
l = 0;
khat = [khat, kk];
fprintf('change detected at %d\n', kk);
%fprintf('Max: %f \t Min: %f\n', max(sig_det(:)), min(sig_det(:)));
end
end
end
function x_out = mult_H(x_in)
X_mat = reshape(x_in, Params.n, Params.r);
x_out = zeros(Params.alpha*Params.m, 1);
for na = 1: Params.alpha
x_out((na-1)*Params.m + 1 : na*Params.m) = Au(:,:,na)' * X_mat * Bo(:,na);
end
end
function w_out = mult_Ht(w_in)
w_out = zeros(Params.n*Params.r, 1);
for na = 1: Params.alpha
tmp_vec = Au(:,:,na) * w_in((na-1)*Params.m+1:na*Params.m);
w_out = w_out + kron(Bo(:,na), tmp_vec);
end
end
end