-
Notifications
You must be signed in to change notification settings - Fork 0
/
pcaInit.m
46 lines (36 loc) · 1.22 KB
/
pcaInit.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
function [z,model] = pcaInit(y,r)
% mypca z (scaled to normal)
[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