[37] | 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 | |
---|