-
Notifications
You must be signed in to change notification settings - Fork 2
/
NIPALS.m
78 lines (69 loc) · 1.67 KB
/
NIPALS.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
function [T, P, U, Q, W, B, B2] = NIPALS (X,Y,a)
%
% Inputs:
% X X matrix normaliz
% Y Y matrix normaliz
% a low reduction value
% Outputs:
% T score for X
% P loading for X
% U score for Y
% Q loading for Y
% B regression coefficients
% B2 regression Matrix
%Initialization
Ej=X;
Fj=Y;
[nX,mX] = size(Ej);
[nY,pY] = size(Fj);
% a value controll, until j_max=min(m,n)-> rank(X), Pag.81
j_max = min([mX,nX]);
if a>j_max
disp('STOPPED: a value not valid')
return
end
% preallocating variables for speed
T= zeros(nX,a);
P= zeros(mX,a);
U= zeros(nX,a);
Q= zeros(pY,a);
W= zeros(mX,a);
B= zeros(a,a);
B2= zeros(mX,pY);
%max number of itertion to convergence
n_max_it = 1e3;
for j = 1 : a
tj_prec = Ej(:,1);
uj = Fj(:,1);
% compute Lv and SV until converfence
for i = 1 : n_max_it
wj = Ej' * uj/ norm(Ej' * uj);
tj = Ej * wj;
qj = Fj' * tj/ norm (Fj' * tj);
uj = Fj * qj;
if norm(tj_prec - tj) < 10e-10
%for convergence 10e-10, example pag, 79
break
end
% updatig tj_prec for next iteration
tj_prec = tj;
end
pj = Ej' * tj/(tj' * tj);
% scaling not necessary (6.15) to (6.17)
bj = uj'*tj/(tj'*tj);
% updatig Ej e Fj for next iteration
Ej = Ej - tj * pj';
Fj = Fj - bj * tj * qj';
% save iteration result:
T(:, j)= tj;
P(:, j)= pj;
U(:, j)= uj;
Q(:, j)= qj;
W(:, j)= wj;
B(j,j)= bj;
B2= W*pinv(P'*W)*pinv(T'*T)*T'*Y;
% check for residual to stop:
if (norm(Ej) ==0 && norm(Fj) ==0)
break
end
end