-
Notifications
You must be signed in to change notification settings - Fork 0
/
pcaInit_v02.m
47 lines (37 loc) · 1.25 KB
/
pcaInit_v02.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
function [z,model] = pcaInit_v02(y,r)
% v_02 : z = u * s (non scaled to normal)
% mypca
[n,dim] = size(y);
y_mean = mean(y);
y = y - repmat(y_mean, n,1);
if r >= 1 %calculated first r bases.
[u,s,v] = svds(y,r);
rank = r;
eigenValue = diag(s);
CumuEnergy = cumsum(eigenValue)./sum(eigenValue);
elseif r<1
[u,s,v] = svd(y,'econ');
eigenValue = diag(s);
CumuEnergy = cumsum(eigenValue)./sum(eigenValue);
idx = find(CumuEnergy >=r);
rank = idx(1);
u = u(:,1:rank);
s = s(1:rank,1:rank);
v = v(:,1:rank);
end
% z = u;
z = u * s;
% base = s*v';
model.rank = rank;
model.energy = CumuEnergy;
model.eigenValue = eigenValue;
model.s = s;
model.v = v;
model.y_mean = y_mean;
% model.project = @(yte) (yte - repmat(y_mean, size(yte,1) ,1))*v*inv(s);
invs = diag( (diag(s)+eps).^(-1) );
% model.project = @(yte) (yte - repmat(y_mean, size(yte,1) ,1))*v*invs;
model.project = @(yte) (yte - repmat(y_mean, size(yte,1) ,1))*v;
% model.recover = @(zte) zte*s*v' + repmat(y_mean, size(zte,1) ,1);
model.recover = @(zte) zte*v' + repmat(y_mean, size(zte,1) ,1);
end