1 | function [Y,details]=sdeIncXL(X,k,FirstSet,d,pars); |
---|
2 | % |
---|
3 | % [Y,details]=sdeInc(X,k,FirstSet,d,pars); |
---|
4 | % |
---|
5 | % |
---|
6 | % Required parameters: |
---|
7 | % |
---|
8 | % X matrix containing the input vectors as columns |
---|
9 | % k size of local neighborhood |
---|
10 | % FirstSet Size of Subsample |
---|
11 | % (only for speedup purpose |
---|
12 | % d dimension of output |
---|
13 | % |
---|
14 | % |
---|
15 | % Optional parameters: |
---|
16 | % |
---|
17 | % pars: |
---|
18 | % pars.slack=1 Use slack for first set (*DEFAULT*) |
---|
19 | % pars.slack=0 Don't use slack for first set |
---|
20 | % |
---|
21 | % pars.solver=0 Use CSDP for non-first sets (*DEFAULT*) |
---|
22 | % pars.solver=1 Use SeDuMi for non-first sets |
---|
23 | % pars.solver=2 SDPT3 (very fast - not fully tested) |
---|
24 | % |
---|
25 | % pars.verify Automatically increases K if graph is not connected |
---|
26 | % (Note: If no K is specified pars.verify=1) |
---|
27 | % pars.verify=1 (default) on |
---|
28 | % pars.verify=0 off |
---|
29 | % |
---|
30 | % |
---|
31 | % pars.factor Weight of Slack variables (only if pars.slack=1) |
---|
32 | % the very closer to 1 ther harder do the slack |
---|
33 | % variables get |
---|
34 | % pars.factor=0.9999 (default) |
---|
35 | % |
---|
36 | % pars.angles=1 preserves angles in local neighborhood (*DEFAULT*) |
---|
37 | % pars.angles=0 does not preserve angles in local neighborhood |
---|
38 | % |
---|
39 | % pars.repell Mx2 list of two points that should repell each other |
---|
40 | % pars.repell=[] (*default*) |
---|
41 | % |
---|
42 | % pars.rweight factor between 0 and 1, of how much the repelling |
---|
43 | % should weight compared to the normal objective |
---|
44 | % function |
---|
45 | % pars.rweight=0.3 (*default*) |
---|
46 | % |
---|
47 | % pars.save Filename to save temporare states |
---|
48 | % default='tempsave' |
---|
49 | % |
---|
50 | % pars.continue Filename to continue from temporare states |
---|
51 | % default='' |
---|
52 | % |
---|
53 | % pars.i initial subset (default random) |
---|
54 | % |
---|
55 | % Output: |
---|
56 | % |
---|
57 | % Y Matrix with output vectors |
---|
58 | % details |
---|
59 | % details.SY Eigenvectors of subsample |
---|
60 | % details.SD Eigenvalues of subsample |
---|
61 | % details.time time needed for computation |
---|
62 | % |
---|
63 | % (version 1.2) |
---|
64 | % (copyright 2004 by Kilian Q. Weinberger: |
---|
65 | % http://www.seas.upenn.edu/~kilianw ) |
---|
66 | [D N]=size(X); |
---|
67 | i=randperm(N); % bug in matlab |
---|
68 | |
---|
69 | if(~exist('pars') | ~isfield(pars,'factor')) pars.factor=0.9999;end; |
---|
70 | if(~isfield(pars,'slack')) pars.slack=1;end; |
---|
71 | if(~isfield(pars,'verify')) pars.verify=1;end; |
---|
72 | if(~isfield(pars,'solver')) pars.solver=0;end; |
---|
73 | if(~isfield(pars,'repell')) pars.repell=[];end; |
---|
74 | if(~isfield(pars,'rweight')) pars.rweight=0.3;end; |
---|
75 | if(~isfield(pars,'save')) pars.save='tempsave.mat';end; |
---|
76 | if(~isfield(pars,'i')) pars.i=randperm(N);end; |
---|
77 | |
---|
78 | A=sparse([],[],[],d*(d+1)/2,(d+1)^2); |
---|
79 | acounter=1;dim=d;b=[]; |
---|
80 | for i=1:dim |
---|
81 | for j=i:dim |
---|
82 | a=zeros(dim+1); |
---|
83 | a(i,j)=1; |
---|
84 | A(acounter,:)=reshape(a,1,(dim+1)^2); |
---|
85 | acounter=acounter+1; |
---|
86 | b=[b;i==j]; |
---|
87 | end; |
---|
88 | end; |
---|
89 | pars.stampA=A; |
---|
90 | pars.stampb=b; |
---|
91 | clear('A'); |
---|
92 | |
---|
93 | |
---|
94 | pars |
---|
95 | |
---|
96 | |
---|
97 | if(~isfield(pars,'continue')) |
---|
98 | |
---|
99 | fprintf('Dividing problem into subproblems ...'); |
---|
100 | % this line is executed twice due to a bug in Matlab 6 |
---|
101 | i=pars.i; |
---|
102 | if(length(pars.repell)>0) |
---|
103 | repeller=unique(pars.repell(:)); |
---|
104 | % keyboard; |
---|
105 | [temp,temp2,reppos]=intersect(repeller,i); |
---|
106 | % [temp,reppos]=ismember(repeller,i) |
---|
107 | if(length(reppos>0)) |
---|
108 | i=[repeller' setdiff(i,repeller)]; |
---|
109 | |
---|
110 | j=zeros(1,length(X)); |
---|
111 | j(repeller)=1:length(repeller); |
---|
112 | pars.repell(:)=j(pars.repell(:)); |
---|
113 | end; |
---|
114 | end; |
---|
115 | firsti=i(1:FirstSet); |
---|
116 | fprintf('done\n'); |
---|
117 | |
---|
118 | |
---|
119 | fprintf('computing distances ....'); |
---|
120 | X2 = sum(X.^2); |
---|
121 | DMfirsti = repmat(X2(firsti),FirstSet,1)+repmat((X2(firsti))',1,FirstSet)-2*(X(:,firsti)'*X(:,firsti)); |
---|
122 | %clear('X2'); |
---|
123 | fprintf('done\n'); |
---|
124 | |
---|
125 | |
---|
126 | if(pars.verify) |
---|
127 | fprintf('Checking if graph is connected ...'); |
---|
128 | k=neighborsUDD(DMfirsti,k); |
---|
129 | fprintf('done\n'); |
---|
130 | end; |
---|
131 | |
---|
132 | [V,DD]=sde(DMfirsti,k,pars); |
---|
133 | |
---|
134 | Y=V(1:d,:); |
---|
135 | |
---|
136 | index=firsti;clear('firsti'); |
---|
137 | jj=ones(1,N); |
---|
138 | jj(index)=0; |
---|
139 | jj=find(jj); |
---|
140 | |
---|
141 | |
---|
142 | if(~strcmp(pars.save,'')) |
---|
143 | fprintf('Saving temporary file...'); |
---|
144 | save(pars.save); |
---|
145 | fprintf('done\n'); |
---|
146 | end; |
---|
147 | |
---|
148 | fprintf('Sorting additional points ...\n'); |
---|
149 | distances=zeros(1,length(jj)); |
---|
150 | for tt1=1:length(distances); |
---|
151 | x=X(:,jj(tt1)); |
---|
152 | distances(tt1)=min(X2(index)'-((X(:,index))'*x).*2+x'*x); |
---|
153 | if(mod(tt1,1000)==0) fprintf('Points processed: %i/%i\n',tt1,length(distances)); end; |
---|
154 | end; |
---|
155 | [temp,tt1]=sort(distances); |
---|
156 | jj=jj(tt1); |
---|
157 | |
---|
158 | fprintf('done\n'); |
---|
159 | |
---|
160 | |
---|
161 | ScaleFactor=max(max(abs(Y))); |
---|
162 | Y=Y./ScaleFactor; |
---|
163 | %DM=DM./(ScaleFactor^2); |
---|
164 | |
---|
165 | err=[]; |
---|
166 | shift=zeros(d,1); |
---|
167 | precomputed=0; |
---|
168 | |
---|
169 | else |
---|
170 | load(pars.continue); |
---|
171 | [temp,precomputed]=size(Y); |
---|
172 | end; |
---|
173 | |
---|
174 | tic; |
---|
175 | while(length(jj)>0) |
---|
176 | |
---|
177 | [temp ntotal]=size(jj); |
---|
178 | |
---|
179 | %distances=DM(index,jj); |
---|
180 | %[temp iio]=sort(min(distances));s |
---|
181 | ii=1; |
---|
182 | |
---|
183 | x=X(:,jj(ii)); |
---|
184 | |
---|
185 | |
---|
186 | distances=X2(index)'-((X(:,index))'*x).*2+x'*x; |
---|
187 | |
---|
188 | |
---|
189 | disttone=zeros(1,k); |
---|
190 | ne=zeros(1,k); |
---|
191 | for te=1:k |
---|
192 | [tt1 tt2]=min(distances); |
---|
193 | disttone(te)=tt1; |
---|
194 | ne(te)=tt2; |
---|
195 | distances(tt2)=inf; |
---|
196 | end; |
---|
197 | % [disttone ne]=sort(distances(:,ii)); |
---|
198 | |
---|
199 | |
---|
200 | |
---|
201 | [y er]=sdeAddPoint3(Y(:,ne)-repmat(shift,1,k),disttone./ScaleFactor^2,pars); |
---|
202 | Y=[Y,y+shift]; |
---|
203 | % XX=[XX x]; |
---|
204 | err=[err er.err]; |
---|
205 | |
---|
206 | % save('temp4'); |
---|
207 | index=[index jj(ii)]; |
---|
208 | |
---|
209 | jj=jj(2:end); |
---|
210 | |
---|
211 | % Xle=X(:,jj); |
---|
212 | |
---|
213 | % centralize Y |
---|
214 | %Y=Y-repmat(y./(N-ntotal+1),1,N-ntotal+1); |
---|
215 | shift=shift+y./(N-ntotal+1); |
---|
216 | lindex=N-ntotal+1; |
---|
217 | fprintf('%i ',lindex); |
---|
218 | |
---|
219 | if(mod(lindex,100)==0) |
---|
220 | timeremaining(length(index)-precomputed,length(jj)); |
---|
221 | if(mod(lindex,1000)==0 & ~strcmp(pars.save,'')) |
---|
222 | save(pars.save); |
---|
223 | end; |
---|
224 | end; |
---|
225 | end; |
---|
226 | |
---|
227 | details.SY=Y(:,1:FirstSet).*ScaleFactor; |
---|
228 | details.index=index; |
---|
229 | details.SD=DD; |
---|
230 | details.time=toc; |
---|
231 | details.err=err; |
---|
232 | details.k=k; |
---|
233 | Y(:,index)=Y.*ScaleFactor; |
---|
234 | |
---|
235 | |
---|
236 | |
---|
237 | |
---|
238 | |
---|
239 | |
---|
240 | |
---|
241 | |
---|
242 | |
---|
243 | function [y,details]=sdeAddPoint3(Anchors,distances,pars); |
---|
244 | |
---|
245 | if(~exist('pars') | ~isfield(pars,'factor')) pars.factor=0.9999;end; |
---|
246 | if(~isfield(pars,'solver')) pars.solver=0;end; |
---|
247 | if(~isfield(pars,'rweight')) pars.rweight=0;end; |
---|
248 | |
---|
249 | |
---|
250 | |
---|
251 | [dim,k]=size(Anchors); |
---|
252 | |
---|
253 | % Construct A |
---|
254 | % First fix the identity matrix in Zxs |
---|
255 | %a=[ones(dim,1); 1]; |
---|
256 | %A=sparse([],[],[],dim*(dim+1)/2+k,(dim+1)^2); |
---|
257 | b=[]; |
---|
258 | |
---|
259 | %acounter=1; |
---|
260 | |
---|
261 | b=pars.stampb; |
---|
262 | A=pars.stampA; |
---|
263 | [temp1 temp2]=size(A); |
---|
264 | acounter=temp1+1; |
---|
265 | |
---|
266 | |
---|
267 | % Now add anchor constraints |
---|
268 | for i=1:k |
---|
269 | a=[Anchors(:,i); -1]; |
---|
270 | |
---|
271 | A(acounter,:)=[reshape(a*a',1,(length(a))^2)]; |
---|
272 | acounter=acounter+1; |
---|
273 | |
---|
274 | b=[b; distances(i)]; |
---|
275 | end; |
---|
276 | |
---|
277 | |
---|
278 | [constraints ll]=size(A); |
---|
279 | A=[[zeros(constraints-k+1,2*(k-1));[eye(k-1) -1*eye(k-1)]] A]; |
---|
280 | |
---|
281 | K.l=2*(k-1); |
---|
282 | K.s=dim+1; |
---|
283 | |
---|
284 | if(pars.rweight==1) |
---|
285 | c=[ones(1,K.l) zeros(1,(K.s)^2)]; |
---|
286 | else |
---|
287 | c=[ones(1,K.l)*pars.factor zeros(1,(K.s)^2-1) -1*(1-pars.factor) ]; |
---|
288 | end; |
---|
289 | |
---|
290 | |
---|
291 | |
---|
292 | %if(pars.solver) |
---|
293 | switch(pars.solver) |
---|
294 | case{0}, |
---|
295 | pars2.printlevel=0; |
---|
296 | [kk temp temp info]=csdp(A,b,c,K,pars2); |
---|
297 | M=reshape(kk(K.l+1:K.l+K.s^2),K.s,K.s); |
---|
298 | case{1}, |
---|
299 | pars2.fid=0; |
---|
300 | [kk temp info]=sedumi(A,b,c,K,pars2); |
---|
301 | M=reshape(kk(end-K.s^2+1:end),K.s,K.s); |
---|
302 | case{2}, |
---|
303 | % keyboard; |
---|
304 | %fprintf('converting into SQLP...'); |
---|
305 | [blk,A,c,b]=convert_sedumi(A',b,c,K); |
---|
306 | %fprintf('DONE\n'); |
---|
307 | sqlparameters; |
---|
308 | OPTIONS.verbose=0; |
---|
309 | [obj,kkt,temp,temp]=sqlp(blk,A,c,b,OPTIONS); |
---|
310 | M=kkt{length(kkt)}; |
---|
311 | %keyboard; |
---|
312 | %M=reshape(kk(K.l+1:K.l+K.s^2),K.s,K.s); |
---|
313 | info=0; |
---|
314 | end; |
---|
315 | |
---|
316 | y=M(1:dim,end); |
---|
317 | G=M(end,end); |
---|
318 | %fprintf('Trace Error: %d\n',abs(G-y'*y)); |
---|
319 | details.err=abs(G-y'*y); |
---|
320 | |
---|
321 | |
---|
322 | |
---|
323 | |
---|
324 | |
---|
325 | |
---|
326 | |
---|
327 | |
---|
328 | |
---|
329 | |
---|
330 | |
---|
331 | |
---|
332 | |
---|
333 | |
---|
334 | |
---|
335 | |
---|
336 | |
---|
337 | |
---|
338 | %%% SIMPLE BUT USEFUL UTIL FUNCTIONS |
---|
339 | |
---|
340 | |
---|
341 | function timeremaining(sofar,left); |
---|
342 | time=toc/(sofar)*left; |
---|
343 | days=floor(time/(60^2*24)); |
---|
344 | hours=floor((time-days*60^2*24)/60^2); |
---|
345 | minutes=floor((time-hours*60^2)/60); |
---|
346 | seconds=floor(time-minutes*60); |
---|
347 | |
---|
348 | fprintf('\nTime remaining: '); |
---|
349 | if(days>0) fprintf('%i Day(s) ',days);end; |
---|
350 | if(hours>0 ) fprintf('%i Hour(s) ',hours);end; |
---|
351 | if(minutes>0 & days==0) fprintf('%i Minute(s) ',minutes);end; |
---|
352 | if(seconds>0 & hours==0 & days==0) fprintf('%i Second(s) ',seconds);end; |
---|
353 | if(time==0) fprintf('You are DONE!!');end; |
---|
354 | fprintf('\n'); |
---|
355 | |
---|
356 | |
---|
357 | function neighbors=getneighborsUDD(DD,K); |
---|
358 | ne=getneighborsDD(DD,K); |
---|
359 | N=length(DD); |
---|
360 | |
---|
361 | for i=1:N |
---|
362 | neighbors{i}=[]; |
---|
363 | end; |
---|
364 | |
---|
365 | for i=1:N |
---|
366 | neighbors{i}=merge(sort(neighbors{i}),sort(ne(:,i))); |
---|
367 | for j=1:K |
---|
368 | neighbors{ne(j,i)}=merge(neighbors{ne(j,i)}, i); |
---|
369 | end; |
---|
370 | end; |
---|
371 | |
---|
372 | |
---|
373 | |
---|
374 | function result=merge(x,y); |
---|
375 | result=unique([x(:);y(:)]); |
---|
376 | |
---|
377 | |
---|
378 | function v=vec(M); |
---|
379 | v=M(:); |
---|
380 | |
---|
381 | |
---|
382 | function k=neighborsUDD(DD,K); |
---|
383 | N=length(DD); |
---|
384 | if(nargin<2) |
---|
385 | K=2; |
---|
386 | end; |
---|
387 | k=K; |
---|
388 | while(k<N & (1-connectedUDD(DD,k))) |
---|
389 | k=k+1; |
---|
390 | fprintf('Trying K=%i \n',k); |
---|
391 | end; |
---|
392 | |
---|
393 | |
---|
394 | |
---|
395 | function result=connectedUDD(DD,K); |
---|
396 | |
---|
397 | % result = connecteddfs (X,K) |
---|
398 | % |
---|
399 | % X input vector |
---|
400 | % K number of neighbors |
---|
401 | % |
---|
402 | % Returns: result = 1 if connected 0 if not connected |
---|
403 | |
---|
404 | if(nargin<2) |
---|
405 | fprintf('Number of Neighbors not specified!\nSetting K=4\n'); |
---|
406 | K=4; |
---|
407 | end; |
---|
408 | N=length(DD); |
---|
409 | ne=getneighborsUDD(DD,K); |
---|
410 | maxSize=0; |
---|
411 | for i=1:N |
---|
412 | if(length(ne{i})>maxSize) maxSize=length(ne{i});end; |
---|
413 | end; |
---|
414 | neighbors=ones(maxSize,N); |
---|
415 | for i=1:N |
---|
416 | neighbors(1:length(ne{i}),i)=ne{i}; |
---|
417 | end; |
---|
418 | oldchecked=[]; |
---|
419 | checked=[1]; |
---|
420 | |
---|
421 | while((size(checked)<N) & (length(oldchecked)~=length(checked))) |
---|
422 | next=neighbors(:,checked); |
---|
423 | next=transpose(sort(next(:))); |
---|
424 | next2=[next(2:end) 0]; |
---|
425 | k=find(next-next2~=0); |
---|
426 | next=next(k); |
---|
427 | |
---|
428 | oldchecked=checked; |
---|
429 | checked=neighbors(:,next(:)); |
---|
430 | checked=transpose(sort(checked(:))); |
---|
431 | checked2=[checked(2:end) 0]; |
---|
432 | k=find(checked-checked2~=0); |
---|
433 | checked=checked(k); |
---|
434 | % if(length(oldchecked)==length(checked)) |
---|
435 | % prod(double(checked==oldchecked)); |
---|
436 | % end; |
---|
437 | end; |
---|
438 | result=(length(checked)==N); |
---|
439 | |
---|
440 | |
---|
441 | |
---|
442 | function X = mat(x,n) |
---|
443 | if nargin < 2 |
---|
444 | n = floor(sqrt(length(x))); |
---|
445 | if (n*n) ~= length(x) |
---|
446 | error('Argument X has to be a square matrix') |
---|
447 | end |
---|
448 | end |
---|
449 | X = reshape(x,n,n); |
---|
450 | |
---|
451 | |
---|
452 | |
---|
453 | |
---|
454 | function neighbors=getneighborsDD(DD,K); |
---|
455 | % PAIRWISE DISTANCES |
---|
456 | N=length(DD); |
---|
457 | % NEIGHBORS |
---|
458 | [sorted,index] = sort(DD); |
---|
459 | neighbors = index(2:(1+K),:); |
---|
460 | |
---|
461 | |
---|
462 | |
---|
463 | |
---|
464 | |
---|
465 | |
---|
466 | |
---|
467 | |
---|
468 | |
---|