% LLE ALGORITHM (using K nearest neighbors) with Low-Dimensional
% Neighborhood Representation
%
% [Y] = lle_ldnr(X,K,dmax)
%
% X = data as D x N matrix (D = dimensionality, N = #points)
% K = number of neighbors
% dmax = max embedding dimensionality
% Y = embedding as dmax x N matrix
%
% This is a modification of the original LLE code as appears in LLE's
% homepage http://www.cs.toronto.edu/~roweis/lle/code/lle.m
%
% The computation of the weight vectors in this version is based on first
% projection the local neighborhood and then find the best weight vector
% with respecet to the projected neighborhood.
% (see LLE with low-dimensional neighborhood representation,
% http://arxiv.org/abs/0808.0780 )
%
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function [Y] = lle_ldnr(X,K,d)
[D,N] = size(X);
fprintf(1,'Local PCA LLE running on %d points in %d dimensions\n',N,D);
% STEP1: COMPUTE PAIRWISE DISTANCES & FIND NEIGHBORS
fprintf(1,'-->Finding %d nearest neighbours.\n',K);
X2 = sum(X.^2,1);
distance = repmat(X2,N,1)+repmat(X2',1,N)-2*X'*X;
[sorted,index] = sort(distance);
neighborhood = index(2:(1+K),:);
% STEP2: SOLVE FOR RECONSTRUCTION WEIGHTS
fprintf(1,'-->Solving for reconstruction weights.\n');
W = zeros(K,N);
for ii=1:N
z = X(:,neighborhood(:,ii))'-repmat(X(:,ii),1,K)'; % shift ith pt to origin
[U,L,V]=svd(z);
U=U(:,d+1:K);
a=U'*ones(K,1);
a=a/(a'*a);
W(:,ii) =U*a;
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 = Y(:,1:d)'*sqrt(N); % bottom evect is [1,1,1,1...] with eval 0
fprintf(1,'Done.\n');