annlle2.m

一个lle的matlab代码

% LLE ALGORITHM (using K nearest neighbors)
% [Y] = lle(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
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function [Y] = annlle(X1,X2,X3,K,d)
X=[ X1 X2 X3];
[D,N] = size(X);
u1=mean(X1,2);
u2=mean(X2,2);
u3=mean(X3,2);
u=(1/3)*(u1+u2+u3);
%%%%%%%%%%%%%%%%%%%%%%%%
%注意：当维数大于样本点个数的时候，协方差可能出现奇异
%为此：我们使用样本点的散布矩阵来代替协方差阵
A=zeros(D,D);

for i=1:N
    sw=(X(:,i)-u)*(X(:,i)-u)';
    A=A+sw;
end        
%%%%%%%%%%%%%%%%%%%%%%%散布矩阵计算结束
%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%计算类内、类间矩阵
sb=(1/3)*((u1-u)*(u1-u)'+(u2-u)*(u2-u)'+(u3-u)*(u3-u)');
n=size(sw);
I=eye(n);
s=inv(sw^(1/2));
swb=s*(s*sb*s+I)*s;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%
fprintf(1,'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),:);
%结束
%基于自适应最近邻的算法
% STEP1: COMPUTE PAIRWISE DISTANCES & FIND NEIGHBORS 
fprintf(1,'-->Finding %d nearest neighbours.\n',K);
distance=[];
for i=1:N     
    for j=1:N; 
    distance(i,j)=(X(:,i)-X(:,j))'*swb*(X(:,i)-X(:,j));
    end
end
[sorted,index] = sort(distance);
neighborhood = index(2:(1+K),:);

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

if(K>D) 
  fprintf(1,'   [note: K>D; regularization will be used]\n'); 
  tol=1e-3; % regularlizer in case constrained fits are ill conditioned
else
  tol=0;
end

W = zeros(K,N);
for ii=1:N
   z = X(:,neighborhood(:,ii))-repmat(X(:,ii),1,K); % shift ith pt to origin
   C = z'*z;                                        % local covariance
   C = C + eye(K,K)*tol*trace(C);                   % regularlization (K>D)
   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 = Y(:,2:d+1)'*sqrt(N); % bottom evect is [1,1,1,1...] with eval 0

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

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


% other possible regularizers for K>D
%   C = C + tol*diag(diag(C));                       % regularlization
%   C = C + eye(K,K)*tol*trace(C)*K;                 % regularlization