lle_pd.m

best routing protocol

% LLE ALGORITHM with pairwise distances (using K nearest neighbors)
%
% [Y] = lle_pd(X,K,dmax)
%
% X = pairwise distances as N x N matrix (N = #points)
% K = number of neighbors
% dmax = max embedding dimensionality
% Y = embedding as dmax x N matrix

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%function [Y] = lle(X,K,d)
function [Y] = lle_pd(X,K,d)

[N,tmp] = size(X);
fprintf(1,'LLE running on %d points\n',N);


% STEP1: FIND NEIGHBORS 
fprintf(1,'-->Finding %d nearest neighbours.\n',K);

%distance = X;
[sorted,index] = sort(X);
neighborhood = index(2:(1+K),:);



% STEP2: SOLVE FOR RECONSTRUCTION WEIGHTS
fprintf(1,'-->Solving for reconstruction weights.\n');

tol = 1e-3;
W = zeros(K,N);
for ii=1:N
  for j=1:K
    for k=1:K
      D(j,k) = X(neighborhood(j,ii), neighborhood(k,ii));
      D(j,k) = D(j,k)*D(j,k);
    end
  end
  D0 = sum(sum(D));
  for j=1:K
    for k=1:K
      %C(j,k) = (sum(D(j,:))/K+sum(D(k,:))/K-D(j,k)-D0/(K*K))/2;
      C(j,k) = (sum(D(j,:))/N+sum(D(k,:))/N-D(j,k)-D0/(N*N))/2;
      %C(j,k) = (sum(D(j,:))+sum(D(k,:))-D(j,k)-D0)/2;
    end
  end
  C = C + eye(K,K)*tol*trace(C);
  W(:,ii) = C\ones(K,1);                           % solve Cw=1
  W(:,ii) = W(:,ii)/sum(W(:,ii));                  % enforce sum(w)=1
end;


% STEP 3: COMPUTE EMBEDDING FROM EIGENVECTS OF COST MATRIX M=(I-W)'(I-W)
fprintf(1,'-->Computing embedding.\n');

% M=eye(N,N); % use a sparse matrix with storage for 4KN nonzero elements
M = sparse(1:N,1:N,ones(1,N),N,N,4*K*N); 
for ii=1:N
  w = W(:,ii);
  jj = neighborhood(:,ii);
  M(ii,jj) = M(ii,jj) - w';
  M(jj,ii) = M(jj,ii) - w;
  M(jj,jj) = M(jj,jj) + w*w';
end;

% CALCULATION OF EMBEDDING
options.disp = 0; options.isreal = 1; options.issym = 1; 
%[Y,eigenvals] = eigs(M,d+1,0,options);
[Y,eigenvals] = eig(full(M));
Y = Y(:,2:d+1)*sqrt(N); % bottom evect is [1,1,1,1...] with eval 0

%DY = allpairdistances(N, Y);
%DY = squareform(pdist(Y));
DY = euclidean_distance(Y,Y);

Y(:,:) = Y(:,:) * (sum(sum(X))/sum(sum(DY)));

fprintf(1,'Done.\n');