demo_scurveandwave.m

一种半监督流形学习算法

function demo_scurveandwave
%
% function demo_scurveandwave.m
% demonstrates the use of graph laplacian for aligning two manifold data
% sets.
% Main functions are regongraph_inverse and regongraph_inverse.
% Please refer to  < www.seas.upenn.edu/~jhham/papers/AISTATS05.pdf >
% for details on the algorithm.
%
% Written by Jihun Ham, 2005, GRASP lab, University of Pennsylvania

load testdata.mat
%clear; close all;

% na=200; nb=400; nc=600; nx=na+nb; ny=na+nc;
% a=1:na; b=na+1:na+nb; c=na+nb+1:na+nb+nc;
% Data=genrectangle(na,nb,nc,2E-4);
% 
% %Data=genuniform(n); n=size(Data,2); 
% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% % unsupervised feature extraction
% % 
% %Data=[Data(1,:)-0.5;2/3*Data(2,:)-1/3];
% 
% clf; subplot(221);
% X=genscurve(Data(:,[a,b])+0.5,1,0); axis square;
% subplot(221);
% Y=genpeaks(Data(:,[a,c])+0.5,1,0); axis square;
% Data=(Data-repmat(min(Data,[],2),[1 na+nb+nc]))./repmat(max(Data,[],2)-min(Data,[],2),[1 na+nb+nc]);

opt.disp = 0; opt.isreal = 1; opt.issym = 1; 

[neighborx,mdistx]=findneighbor(X,6);
[Lx,Wx]=graphlaplacian(X,neighborx,'affine');
[Ex,Dx]=eigs(Lx'*Lx,3,'SM',opt);[dummy,index]=sort(diag(Dx));
Ex=Ex(:,index(2:3)); ddx=1./dummy(2:end);
Ex=(Ex-repmat(min(Ex,[],1),[nx 1]))./repmat(max(Ex,[],1)-min(Ex,[],1),[nx 1]);

[neighbory,mdisty]=findneighbor(Y,10);
[Ly,Wy]=graphlaplacian(Y,neighbory,'affine'); 
[Ey,Dy]=eigs(Ly'*Ly,3,'SM',opt);[dummy,index]=sort(diag(Dy));
Ey=Ey(:,index(2:3)); ddy=1./dummy(2:end);
Ey=(Ey-repmat(min(Ey,[],1),[ny 1]))./repmat(max(Ey,[],1)-min(Ey,[],1),[ny 1]);

clf; colormap default;
subplot(221),cla,scatter3(X(1,:),X(2,:),X(3,:),8,Data(1,[a b]),'filled');title('X : s-curve')
subplot(222),cla,scatter3(Y(1,:),Y(2,:),Y(3,:),8,Data(1,[a c]),'filled');title('Y : wave')
subplot(223),cla,scatter(Ex(:,1),Ex(:,2),8,Data(1,[a,b]),'filled'); axis off square; title('embedding of X');
subplot(224),cla,scatter(Ey(:,1),Ey(:,2),8,Data(1,[a c]),'filled'); axis off square; title('embedding of Y');


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Aligning by parameter maps 

na1=50;

lx=randperm(nx); lx=lx(1:na1); ux=setdiff(1:nx,lx); Fl=Data(:,lx)'; 
Fu=regongraph_inverse(Wx,Fl,lx);
F=zeros(nx,2); F(lx,:)=Fl; F(ux,:)=Fu;

ly=randperm(ny); ly=ly(1:na1); uy=setdiff(1:ny,ly); 
tData=Data(:,[a c]); Gl=tData(:,ly)'; 
Gu=regongraph_inverse(Wy,Gl,ly);
G=zeros(ny,2); G(ly,:)=Gl; G(uy,:)=Gu;

figure;
subplot(121) ;cla, axis off; grid off; hold on; title('unaligned')
view([15 15]); axis([-0.15 1.25 -0.15 1.25 -0.05 1.25]); axis vis3d; 
line([-0.1 1.2],[-0.1 -0.1],[0 0],'Color',[0 0 0],'LineWidth',1);
line([-0.1 -0.1],[-0.1 1.2],[0 0],'Color',[0 0 0],'LineWidth',1);
line([-0.1 1.2],[-0.1 -0.1],[1 1],'Color',[0 0 0],'LineWidth',1);
line([-0.1 -0.1],[-0.1 1.2],[1 1],'Color',[0 0 0],'LineWidth',1);
line([-0.1 1.2],[-0.1 -0.1],[0.5 0.5],'Color',[0 0 0],'LineWidth',1);
line([-0.1 -0.1],[-0.1 1.2],[0.5 0.5],'Color',[0 0 0],'LineWidth',1);
line([-0.1 -0.1],[-0.1 -0.1],[0 1.2],'Color',[0 0 0],'LineWidth',1);
scatter3(Ey(:,1),Ey(:,2),zeros(ny,1),6,Data(1,[a,c]),'filled');
for n=1:50
    line([Gl(n,1) Ey(ly(n),1)],[Gl(n,2) Ey(ly(n),2)],[0.5 0],'Color',[.3 .3 .3],'LineWidth',0.5,'LineStyle','-');
end
scatter3(Gl(:,1),Gl(:,2),0.5*ones(na1,1),6,tData(1,ly),'filled');
scatter3(Fl(:,1),Fl(:,2),0.5*ones(na1,1),6,Data(1,lx),'filled');
for n=1:50
    line([Ex(lx(n),1) Fl(n,1)],[Ex(lx(n),2) Fl(n,2)],[1 0.5],'Color',[.3 .3 .3],'LineWidth',0.5,'LineStyle','-');
end
scatter3(Ex(:,1),Ex(:,2), ones(nx,1),6,Data(1,[a,b]),'filled');

subplot(122);cla,axis off; grid off; hold on; title('aligned by semi-supervision');
view([15 15]); axis([-0.15 1.25 -0.15 1.25 -0.05 1.25]); axis vis3d; 
line([-0.1 1.2],[-0.1 -0.1],[0 0],'Color',[0 0 0],'LineWidth',1);
line([-0.1 -0.1],[-0.1 1.2],[0 0],'Color',[0 0 0],'LineWidth',1);
line([-0.1 1.2],[-0.1 -0.1],[1 1],'Color',[0 0 0],'LineWidth',1);
line([-0.1 -0.1],[-0.1 1.2],[1 1],'Color',[0 0 0],'LineWidth',1);
line([-0.1 1.2],[-0.1 -0.1],[0.5 0.5],'Color',[0 0 0],'LineWidth',1);
line([-0.1 -0.1],[-0.1 1.2],[0.5 0.5],'Color',[0 0 0],'LineWidth',1);
line([-0.1 -0.1],[-0.1 -0.1],[0 1.2],'Color',[0 0 0],'LineWidth',1);
scatter3(G(:,1),G(:,2),zeros(ny,1),6,Data(1,[a,c]),'filled');
for n=1:na1
    line([Gl(n,1) G(ly(n),1)],[Gl(n,2) G(ly(n),2)],[0.5 0],'Color',[.3 .3 .3],'LineWidth',0.5,'LineStyle','-');
end
scatter3(Gl(:,1),Gl(:,2),0.5*ones(na1,1),6,tData(1,ly),'filled');
scatter3(Fl(:,1),Fl(:,2),0.5*ones(na1,1),6,Data(1,lx),'filled');
for n=1:na1
    line([F(lx(n),1) Fl(n,1)],[F(lx(n),2) Fl(n,2)],[1 0.5],'Color',[.3 .3 .3],'LineWidth',0.5,'LineStyle','-');
end
scatter3(F(:,1),F(:,2), ones(nx,1),6,Data(1,[a,b]),'filled');

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Aligning by pairwise correspondence

na2=200; 

[Fc,Gc]=regongraph_coupled(Lx,Ly,na2,2);

Fc=(Fc-repmat(min(Fc,[],1),[nx 1]))./repmat(max(Fc,[],1)-min(Fc,[],1),[nx 1]);
Gc=(Gc-repmat(min(Gc,[],1),[ny 1]))./repmat(max(Gc,[],1)-min(Gc,[],1),[ny 1]);

figure;
subplot(121);cla,axis off; grid off; hold on; title('unaligned')
view([15 15]); axis([-0.15 1.25 -0.15 1.25 -0.05 0.75]); axis vis3d; 
line([-0.1 1.2],[-0.1 -0.1],[0 0],'Color',[0 0 0],'LineWidth',1);
line([-0.1 -0.1],[-0.1 1.2],[0 0],'Color',[0 0 0],'LineWidth',1);
line([-0.1 1.2],[-0.1 -0.1],[0.5 0.5],'Color',[0 0 0],'LineWidth',1);
line([-0.1 -0.1],[-0.1 1.2],[0.5 0.5],'Color',[0 0 0],'LineWidth',1);
line([-0.1 -0.1],[-0.1 -0.1],[0 0.7],'Color',[0 0 0],'LineWidth',1);
scatter3(Ey(:,1),Ey(:,2),zeros(ny,1),6,Data(1,[a,c]),'filled');
for n=1:50%na2
    line([Ex(n,1) Ey(n,1)],[Ex(n,2) Ey(n,2)],[0.5 0],'Color',[.3 .3 .3],'LineWidth',0.5,'LineStyle','-');
end
scatter3(Ex(:,1),Ex(:,2),0.5*ones(nx,1),6,Data(1,[a,b]),'filled');

subplot(122); cla,axis off; grid off; hold on; title('aligned by correspondence');
view([15 15]); axis([-0.15 1.25 -0.15 1.25 -0.05 0.75]); axis vis3d; 
line([-0.1 1.2],[-0.1 -0.1],[0 0],'Color',[0 0 0],'LineWidth',1);
line([-0.1 -0.1],[-0.1 1.2],[0 0],'Color',[0 0 0],'LineWidth',1);
line([-0.1 1.2],[-0.1 -0.1],[0.5 0.5],'Color',[0 0 0],'LineWidth',1);
line([-0.1 -0.1],[-0.1 1.2],[0.5 0.5],'Color',[0 0 0],'LineWidth',1);
line([-0.1 -0.1],[-0.1 -0.1],[0 0.7],'Color',[0 0 0],'LineWidth',1);
scatter3(Gc(:,1),Gc(:,2),zeros(ny,1),6, Data(1,[a,c]),'filled');
for n=1:50%na2
    line([Fc(n,1) Gc(n,1)],[Fc(n,2) Gc(n,2)],[0.5 0],'Color',[.3 .3 .3],'LineWidth',0.5,'LineStyle','-');
end
scatter3(Fc(:,1),Fc(:,2),0.5*ones(nx,1),6,Data(1,[a,b]),'filled');


function [L,W,E] = graphlaplacian(X,neighbor,method,param,normalized,mdist)
%function [L,W,E] = graphlaplacian(X,neighbor,method,param,normalized,mdist)

% W'1=1, L=(D-W)(D-W');

normalized=0;
%if nargin<5 normalized=0; end
[d,N]=size(X);
W=zeros(N);
E=zeros(N,1);
switch method
case 'adjacency'
    for ii=1:N
        W(ii,neighbor(:,ii))=1;
    end
case 'affine'
    K=size(neighbor,1); tol=1E-2;%tol=param;
    for ii=1:N
  		z = X(:,neighbor(:,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)
		t = C\ones(K,1);                           % solve Cw=1
   		t = t/sum(t);                  % enforce sum(w)=1
       	W(ii,neighbor(:,ii))=t';
        E(ii)=norm(X(:,ii)-X(:,neighbor(:,ii))*t);
    end;
    index=find(W<1E-8 & W>-1E-8); W(index)=0;
    
case 'nonneg'
    K=size(neighbor,1);
    if (K>d) disp('Possibly ill-conditioned..'); end
    for ii=1:N
   		c=lsqnonneg(X(:,neighbor(:,ii)),X(:,ii));
       	W(ii,neighbor(:,ii))=c';
        E(ii)=norm(X(:,ii)-X(:,neighbor(:,ii))*c);
    end;
    index=find(W<1E-8); W(index)=0;
    
case 'convex'
    options=optimset('Display','off','LargeScale','off');
    %Display - Level of display [ off | iter | notify | final ] 
    K=size(neighbor,1);
    if (K>d) disp('Possibly ill-conditioned..'); end
	for ii=1:N
        Xn=X(:,neighbor(:,ii));
        c=quadprog(Xn'*Xn,-Xn'*X(:,ii),[],[],ones(1,K),1,zeros(K,1),ones(K,1),[],options); 
       	W(ii,neighbor(:,ii))=c';
        E(ii)=norm(X(:,ii)-X(:,neighbor(:,ii))*c);
    end;

case 'gaussian'
	Xsq = sum(X.^2,1);
	distance = repmat(Xsq,N,1)+repmat(Xsq',1,N)-2*X'*X;
	for ii=1:N
		jj=neighbor(:,ii);
		W(ii,jj)=exp(-distance(ii,jj)/2*(mdist(ii)*param)^2);
	end;
    index=find(W<1E-8); W(index)=0;
otherwise
    disp('method not supported.');
end

%W=(W+W')/2;
D=diag(sum(W,2));
L=D-W;

if (normalized==1)
    d=1./diag(D);
    L=diag(d)*L;
    W=diag(diag(L))-L;
end
L=sparse(L); W=sparse(W);

function [neighbor,nndist] = findneighbor(X,K)
%function [neighbor,nndist] = findneighbor(X,K);

N=size(X,2);
dist_e=euclidist(X);
[sorted,index]=sort(dist_e);
neighbor=index(2:(K+1),:);

nndist=sparse(N,N);
for i=1:N
    nndist(index(2:K+1),i)=sorted(2:K+1,i);
end
nndist=(nndist+nndist')/2;

function dist=euclidist(X1,X2,M);
% function dist=euclidist(X1,X2,M);

[D,N]=size(X1);
if (nargin==1 | isempty(X2)), X2=X1; end

if nargin==3
    X1=X1.*repmat(M(:),[1 N]);
    X2=X2.*repmat(M(:),[1 N]);
end

[D1,N1] = size(X1); [D2,N2]=size(X2);
X1sq = sum(X1.^2,1); X2sq=sum(X2.^2,1);
dist_squared = repmat(X1sq',1,N2)+repmat(X2sq,N1,1)-2*X1'*X2;
ind=find(dist_squared<0);
dist_squared(ind)=0;
dist=sqrt(dist_squared);


function Fu=regongraph_inverse(W,Fl,l)
% function Fu=regongraph_inverse(W,Fl,l)
% find harmonic functions (Fu) of unlabled nodes (u) given the values (Fl) at
% the labled nodes (l)

n=length(W); 
D=diag(sum(W,2));
u=setdiff(1:n,l);
Fu=inv(D(u,u)-W(u,u))*W(u,l)*Fl;


function [F,G,dd]=regongraph_coupled(Lx,Ly,na,ndim)
% function [F,G,dd]=regongraph_coupled(Lx,Ly,na,ndim)
% Regression on graphs: solves eigenvalue problem

opt.disp=0;
opt.issym=1;
opt.isreal=1;

nx=length(Lx); ny=length(Ly);
nb=nx-na; nc=ny-na;

L=[Lx,zeros(na+nb,nc);zeros(nc,na+nb+nc)] ...
 +[Ly(1:na,1:na),zeros(na,nb),Ly(1:na,na+1:na+nc); ...
   zeros(nb,na+nb+nc);Ly(na+1:na+nc,1:na),zeros(nc,nb),Ly(na+1:na+nc,na+1:na+nc)];
[E,D]=eigs(L'*L,ndim+1,'SM',opt);[dummy,index]=sort(diag(D));
E=E(:,index(2:end)); dd=dummy(2:end);
F=E(1:nx,:); G=E([1:na, nx+1:nx+nc],:);