-
Notifications
You must be signed in to change notification settings - Fork 2
/
GA.m
126 lines (102 loc) · 3.02 KB
/
GA.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
function B = GA(Jp,cplmat,W1,S1,W2,S2,nscale,cscale1,cscale2,a1,a2,tol)
% computing optimal edge weigth matrix B using stochastic gradient ascent algorithm
% Yu Hang, Jan. 2015, NTU
% initialize B mattrix
p = size(Jp,1);
n = size(cplmat,1);
[edgerow,edgecol] = find(triu(Jp,1));
nedge = length(edgerow);
W0 = log(a2*(ones(nedge,1))/sqrt(nedge));
W = W0;
% initialize parameters for ADADELTA
eps = 1e-6;
rho = 0.95;
d2W = 0;
g2W = 0;
gW_array = zeros(nedge,n);
% eta = 1e-3;
% compute matrix E
E = sparse(nedge,p);
for i = 1:nedge
if edgecol(i)<p
E(i,edgerow(i)) = 1;
E(i,edgecol(i)) = -1;
else
E(i,edgerow(i)) = 1;
end
end
% begin stochastic gradient ascent algorithm
for nt = 1:5e5
Bv = exp(W);
% ne = ceil(n*rand(1));
% compute low rank approximation matrix
% L = LowRankApp(4,W1,S1,W2,S2,nscale,cscale1,cscale2);
% RC = QBC\L;
% gW = (n*cplmat(ne,:).'.*sum((E*RC).*(E*L),2)-4*a1*(norm(Bv)^2-a2^2)*Bv).*Bv;
parfor ne = 1:n
BC = sparse(edgerow,edgecol,Bv.*cplmat(ne,:).',p,p);
BC = BC+BC.';
QBC = spdiags(sum(BC,2),0,-BC);
QBC = [QBC(1:p-1,1:p-1),sparse(p-1,1);sparse(1,p-1),1];
gW_array(:,ne) = cplmat(ne,:).'.*sum(E*(QBC\speye(p)).*E,2);
end
gW = (mean(gW_array,2)-4*a1*(norm(Bv)^2-a2^2)*Bv).*Bv;
g2W = rho*g2W+(1-rho)*gW.^2;
dW = sqrt(d2W+eps)./sqrt(g2W+eps).*gW;
d2W = rho*d2W+(1-rho)*dW.^2;
W = W+dW;%eta*gW; %
if rem(nt,10) == 0
Wh = W;
save resulttmptg2;
mean(abs(Wh-W0))
if mean(abs(Wh-W0)) < 1e-4
break;
else
W0 = Wh;
% eta = eta*0.995;
end
end
end
W0 = log(a2*exp(W)/norm(exp(W)));
W = W0;
d2W = 0;
g2W = 0;
% eta = 1e-3;
for nt = 1:5e5
Bv = exp(W);
B = sparse(edgerow,edgecol,Bv,p,p);
B = B+B.';
QB = spdiags(sum(B,2),0,-B);
QB = [QB(1:p-1,1:p-1),sparse(p-1,1);sparse(1,p-1),1];
% compute low rank approximation matrix
% L = LowRankApp(4,W1,S1,W2,S2,nscale,cscale1,cscale2);
% R = QB\L;
% ne = ceil(n*rand(1));
% RC = QBC\L;
% gW = (n*sum((spdiags(cplmat(ne,:).',0,nedge,nedge)*(E*RC)-E*R).*(E*L),2)-4*a1*(norm(Bv)^2-a2^2)*Bv).*Bv;
parfor ne = 1:n
BC = sparse(edgerow,edgecol,Bv.*cplmat(ne,:).',p,p);
BC = BC+BC.';
QBC = spdiags(sum(BC,2),0,-BC);
QBC = [QBC(1:p-1,1:p-1),sparse(p-1,1);sparse(1,p-1),1];
gW_array(:,ne) = cplmat(ne,:).'.*sum(E*(QBC\speye(p)).*E,2);
end
gW = (mean(gW_array,2)-sum(E*(QB\speye(p)).*E,2)-4*a1*(norm(Bv)^2-a2^2)*Bv).*Bv;
g2W = rho*g2W+(1-rho)*gW.^2;
dW = sqrt(d2W+eps)./sqrt(g2W+eps).*gW;
d2W = rho*d2W+(1-rho)*dW.^2;
W = W+dW; %eta*gW; %
if rem(nt,10) == 0
Wh = W;
save resulttmptg2;
mean(abs(Wh-W0))
if mean(abs(Wh-W0)) < 1e-6
break;
else
W0 = Wh;
% eta = eta*0.995;
end
end
end
B = sparse(edgerow,edgecol,a2/norm(exp(Wh))*exp(Wh),p,p);
B = B+B.';