mani.m

包含了大多数流形学习方法的代码,有PCA,ISOMAP,LLE,HLLE

tic;
options.dims = handles.d;
[Y, R] = Isomap(D, handles.K, options);
handles.Y = Y.coords{1}';
handles.runTime = toc;
guidata(hObject, handles);
PlotEmbedding(hObject,eventdata,handles,3);
end;
% Run LLE.
try
tic;
handles.Y = lle(handles.X',handles.K,handles.d)';
handles.runTime = toc;
guidata(hObject, handles);
PlotEmbedding(hObject,eventdata,handles,4);
end;
% Run HLLE.
try
tic;
[Y, mse] = HLLE(handles.X',handles.K,handles.d);
handles.Y = Y';
handles.runTime = toc;
guidata(hObject, handles);
PlotEmbedding(hObject,eventdata,handles,5);
end;
% Run Laplacian.
try
tic;
[E,V] = leigs(handles.X, 'nn', handles.K, handles.d+1);
handles.Y = E(:,1:handles.d);
handles.runTime = toc;
guidata(hObject, handles);
PlotEmbedding(hObject,eventdata,handles,6);
end;
% Run Diffusion Map
try
tic;
handles.Y = diffusionKernel(handles.X,handles.sigma,handles.alpha,handles.d);
handles.runTime = toc;
guidata(hObject, handles);
PlotEmbedding(hObject,eventdata,handles,7);
end;
% Run LTSA, if possible.
try
tic;
[T,NI] = LTSA (handles.X', handles.d, handles.K);
handles.Y = T';
handles.runTime = toc;
guidata(hObject, handles);
PlotEmbedding(hObject,eventdata,handles,8);
end;
% Finished running all algorithms.
updateString{1} = 'Ran all 8 algorithms.';
updateString{2} = 'Was it worth the wait?';
set(handles.UpdatesText,'String',updateString);


% --- Executes on button press in ExampleButton.
function ExampleButton_Callback(hObject, eventdata, handles)
exampleValue = get(handles.ExampleMenu,'Value');
N = round(str2double(get(handles.PointsEdit,'String')));   % Number of points to use.
N = max(1,N);
ExParam = str2double(get(handles.ParamEdit,'String'));  % Second parameter for example
ExParam = max(0.0,ExParam);
switch exampleValue
    case 1  % Swiss Roll
        tt = (3*pi/2)*(1+2*rand(1,N));  
        height = 21*rand(1,N);
        handles.X = [tt.*cos(tt); height; ExParam*tt.*sin(tt)]';
        handles.ColorVector = tt';
        updateString{1} = 'Swiss Roll example loaded.';
    case 2  % Swiss Hole
        % Swiss Roll w/ hole example taken from Donoho & Grimes
        tt = (3*pi/2)*(1+2*rand(1,2*N));  
        height = 21*rand(1,2*N);
        kl = repmat(0,1,2*N);
        for ii = 1:2*N
            if ( (tt(ii) > 9)&(tt(ii) < 12))
                if ((height(ii) > 9) & (height(ii) <14))
                    kl(ii) = 1;
                end;
            end;
        end;
        kkz = find(kl==0);
        tt = tt(kkz(1:N));
        height = height(kkz(1:N));
        handles.X = [tt.*cos(tt); height; ExParam*tt.*sin(tt)]';     
        handles.ColorVector = tt';
        updateString{1} = 'Swiss Hole example loaded.';
    case 3  % Corner Planes
        k = 1;
        xMax = floor(sqrt(N));
        yMax = ceil(N/xMax);
        cornerPoint = floor(yMax/2);
        for x = 0:xMax
            for y = 0:yMax
                if y <= cornerPoint
                    X(k,:) = [x,y,0];
                    ColorVector(k) = y;
                else
                    X(k,:) = [x,cornerPoint+(y-cornerPoint)*cos(pi*ExParam/180),(y-cornerPoint)*sin(pi*ExParam/180)];
                    ColorVector(k) = y;
                end;
                k = k+1;
            end;
        end;
        handles.X = X;
        handles.ColorVector = ColorVector';
        updateString{1} = 'Corner Planes example loaded.';
    case 4  % Punctured Sphere by Saul & Roweis
        inc = 9/sqrt(N);   %inc = 1/4;
        [xx,yy] = meshgrid(-5:inc:5);
        rr2 = xx(:).^2 + yy(:).^2;
        [tmp ii] = sort(rr2);
        Y = [xx(ii(1:N))'; yy(ii(1:N))'];
        a = 4./(4+sum(Y.^2));
        handles.X = [a.*Y(1,:); a.*Y(2,:); ExParam*2*(1-a)]';
        handles.ColorVector = handles.X(:,3);
        updateString{1} = 'Punctured Sphere example loaded.';
    case 5  % Twin Peaks by Saul & Roweis
        inc = 1.5 / sqrt(N);  % inc = 0.1;
        [xx2,yy2] = meshgrid(-1:inc:1);
        zz2 = sin(pi*xx2).*tanh(3*yy2);
        xy = 1-2*rand(2,N);
        handles.X = [xy; sin(pi*xy(1,:)).*tanh(3*xy(2,:))]';
        handles.X(:,3) = ExParam * handles.X(:,3);
        handles.ColorVector = handles.X(:,3);
        updateString{1} = 'Twin Peaks example loaded.';
    case 6  % 3D Clusters
        numClusters = ExParam;
        numClusters = max(1,numClusters);
        Centers = 10*rand(numClusters,3);
        D = L2_distance(Centers',Centers',1);
        minDistance = min(D(find(D>0)));
        k = 1;
        N2 = N - (numClusters-1)*9;
        for i = 1:numClusters
            for j = 1:ceil(N2/numClusters)
               X(k,1:3) = Centers(i,1:3)+(rand(1,3)-0.5)*minDistance/sqrt(12);
               ColorVector(k) = i;
               k = k + 1;
           end;
           % Connect clusters with straight line.
           if i < numClusters
               for t = 0.1:0.1:0.9
                    X(k,1:3) = Centers(i,1:3) + (Centers(i+1,1:3)-Centers(i,1:3))*t;
                    ColorVector(k) = 0;
                    k = k+1;
                end;
           end;
        end;
        handles.X = X;
        handles.ColorVector = ColorVector;
        updateString{1} = '3D Clusters example loaded.';
    case 7  % Toroidal Helix by Coifman & Lafon
        noiseSigma=0.05;   %noise parameter
        t = (1:N)'/N;
        t = t.^(ExParam)*2*pi;
        handles.X = [(2+cos(8*t)).*cos(t) (2+cos(8*t)).*sin(t) sin(8*t)]+noiseSigma*randn(N,3);
        handles.ColorVector = t;
        updateString{1} = 'Toroidal Helix example loaded.';
    case 8  % Gaussian randomly sampled
        X = ExParam * randn(N,3);
        X(:,3) = 1 / (ExParam^2 * 2 * pi) * exp ( (-X(:,1).^2 - X(:,2).^2) / (2*ExParam^2) );
        handles.X = X;
        handles.ColorVector = X(:,3);
        updateString{1} = 'Gaussian example loaded.';
    case 9  % Occluded disks
        m = 20;   % Image size m x m.
        Rd = 3;   % Disk radius.
        Center = (m+1)/2;
        u0 = zeros(m,m);
        for i = 1:m
            for j = 1:m
                if (Center - i)^2 + (Center - j)^2 < ExParam^2
                    u0(i,j) = 1;
                end;
            end;
        end;
        for diskNum = 1:N
            DiskX(diskNum) = (m-1)*rand+1;
            DiskY(diskNum) = (m-1)*rand+1;
            u = u0;
            for i = 1:m
                for j = 1:m
                    if (DiskY(diskNum) - i)^2 + (DiskX(diskNum) - j)^2 < Rd^2
                        u(i,j) = 1;
                    end;
                end;
            end;
            X(diskNum,:) = reshape(u,1,m*m);
        end;
        % Since this is a special manifold, plot separately.
        axes(handles.maniAXES);
        t = 0:0.1:2*pi+0.1;
        plot(ExParam*cos(t)+Center,ExParam*sin(t)+Center);
        axis([0.5 m+0.5 0.5 m+0.5]);
        hold on;
        handles.ColorVector = (sqrt((DiskX-Center).^2+(DiskY-Center).^2))';
        scatter(DiskX,DiskY,12,handles.ColorVector');
        hold off;
        handles.X = X;
        updateString{1} = 'Occluded Disks example loaded.';
        updateString{3} = 'Warning: Embedding may be slow.';
end;
assignin ('base','maniX',handles.X);
updateString{2} = 'Manifold data written to matrix "maniX"';
set(handles.UpdatesText,'String',updateString);
handles.isExample = 1;
guidata(hObject, handles);
if exampleValue ~= 9    % Already plotted Occluded Disks example.
    PlotManifold(hObject,eventdata,handles,0);
end;


% --- Plot the input manifold.
function PlotManifold (hObject, eventdata, handles, plotNumber)
[m,n] = size(handles.X);
if get(handles.ColorCheck,'Value') == 1 & handles.isExample == 0
    handles.ColorVector = evalin ('base',get(handles.ColorEdit,'String'));
    Clength = length(handles.ColorVector);
    if Clength < m
        handles.ColorVector(Clength+1:m) = handles.ColorVector(Clength);
    elseif Clength > m
        handles.ColorVector = handles.ColorVector(1:m);
    end;
end;
guidata(hObject, handles);
if plotNumber == 0
    axes(handles.maniAXES);
else
    figure;
    subplot(331);
end;
if n == 2
    if get(handles.ColorCheck,'Value') == 1 | handles.isExample == 1
        scatter(handles.X(:,1),handles.X(:,2),12,handles.ColorVector,'filled');
    else
        scatter(handles.X(:,1),handles.X(:,2),12,'filled');
    end;
    axis tight;
elseif n == 3
    if get(handles.ColorCheck,'Value') == 1 | handles.isExample == 1
        scatter3(handles.X(:,1),handles.X(:,2),handles.X(:,3),12,handles.ColorVector,'filled');
    else
        scatter3(handles.X(:,1),handles.X(:,2),handles.X(:,3),12,'filled');
    end;
    axis tight;
elseif n == 1
    if get(handles.ColorCheck,'Value') == 1 | handles.isExample == 1
        scatter(handles.X(:,1),ones(m,1),12,handles.ColorVector,'filled');
    else
        scatter(handles.X(:,1),ones(m,1),12,'filled');
    end;
    axis tight;
else
    cla;
    axis([-1 1 -1 1]);
    text(-0.7,0,'Only plots 2D or 3D data');
    axis off;
end;


% --- Plot the output embedding.
function PlotEmbedding (hObject, eventdata, handles, plotNumber)
[m,n] = size(handles.Y);
if get(handles.ColorCheck,'Value') == 1 & handles.isExample == 0
    handles.ColorVector = evalin ('base',get(handles.ColorEdit,'String'));
    Clength = length(handles.ColorVector);
    if Clength < m
        handles.ColorVector(Clength+1:m) = handles.ColorVector(Clength);
    elseif Clength > m
        handles.ColorVector = handles.ColorVector(1:m);
    end;
end;
guidata(hObject, handles);
if plotNumber == 0
    axes(handles.embedAXES);
else
    subplot(3,3,plotNumber+1);
end;
if n == 2
    if get(handles.ColorCheck,'Value') == 1 | handles.isExample == 1
        scatter(handles.Y(:,1),handles.Y(:,2),12,handles.ColorVector,'filled');
    else
        scatter(handles.Y(:,1),handles.Y(:,2),12,'filled');
    end;
    axis tight;
elseif n == 3
    if get(handles.ColorCheck,'Value') == 1 | handles.isExample == 1
        scatter3(handles.Y(:,1),handles.Y(:,2),handles.Y(:,3),12,handles.ColorVector,'filled');
    else
        scatter3(handles.Y(:,1),handles.Y(:,2),handles.Y(:,3),12,'filled');
    end;
    axis tight;
elseif n == 1
    if get(handles.ColorCheck,'Value') == 1 | handles.isExample == 1
        scatter(handles.Y(:,1),ones(m,1),12,handles.ColorVector,'filled');
    else
        scatter(handles.Y(:,1),ones(m,1),12,'filled');
    end;
    axis tight;
else
    cla;
    axis([-1 1 -1 1]);
    text(-0.7,0,'Only plots 2D or 3D data');
    axis off;
end;
timeunits = 's';
if handles.runTime >= 60
    handles.runTime = handles.runTime / 60;
    timeunits = 'm';
end;
if handles.runTime >= 60
    handles.runTime = handles.runTime / 60;
    timeunits = 'h';
end;
switch plotNumber
    case 1
        title(['MDS: ',num2str(handles.runTime),timeunits]);
        pause(0.01);
    case 2
        title(['PCA: ',num2str(handles.runTime),timeunits]);
        pause(0.01);
    case 3
        title(['ISOMAP: ',num2str(handles.runTime),timeunits]);
        pause(0.01);
    case 4
        title(['LLE: ',num2str(handles.runTime),timeunits]);
        pause(0.01);
    case 5
        title(['Hessian LLE: ',num2str(handles.runTime),timeunits]);
        pause(0.01);
    case 6
        title(['Laplacian: ',num2str(handles.runTime),timeunits]);
        xlabel(['KNN = ',num2str(handles.K)]);
        pause(0.01);
    case 7
        title(['Diffusion Map: ',num2str(handles.runTime),timeunits]);
        xlabel(['Alpha = ',num2str(handles.alpha)]);
        pause(0.01);
    case 8
        title(['LTSA: ',num2str(handles.runTime),timeunits]);
        xlabel(['Sigma = ',num2str(handles.sigma)]);
end;

% --- LLE ALGORITHM (using K nearest neighbors)
% Written by Sam Roweis & Lawrence Saul
function [Y] = lle(X,K,d)
warning off;
[D,N] = size(X);
% STEP1: COMPUTE PAIRWISE DISTANCES & FIND NEIGHBORS 
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
if(K>D) 
  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)
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


% --- ISOMAP Algorithm
% Written by Tenenbaum, de Silva, and Langford (2000). 
function [Y, R] = Isomap(D, K, options); 
%%%%% Step 0: Initialization and Parameters %%%%%
N = size(D,1); 
INF =  1000*max(max(D))*N;  %% effectively infinite distance
dims = options.dims; 
comp = 1;  % Assume one component.
Y.coords = cell(length(dims),1); 
R = zeros(1,length(dims)); 
%%%%% Step 1: Construct neighborhood graph %%%%%
[tmp, ind] = sort(D); 
for i=1:N
    D(i,ind((2+K):end,i)) = INF; 
end;
D = min(D,D');    %% Make sure distance matrix is symmetric
%%%%% Step 2: Compute shortest paths %%%%%
for k=1:N
     D = min(D,repmat(D(:,k),[1 N])+repmat(D(k,:),[N 1])); 
end
%%%%% Remove outliers from graph %%%%%
n_connect = sum(~(D==INF));        %% number of points each point connects to
[tmp, firsts] = min(D==INF);       %% first point each point connects to
[comps, I, J] = unique(firsts);    %% represent each connected component once
size_comps = n_connect(comps);     %% size of each connected component
[tmp, comp_order] = sort(size_comps);  %% sort connected components by size
comps = comps(comp_order(end:-1:1));    
size_comps = size_comps(comp_order(end:-1:1)); 
n_comps = length(comps);               %% number of connected components
if (comp>n_comps)                
     comp=1;                              %% default: use largest component
end
Y.index = find(firsts==comps(comp)); 
D = D(Y.index, Y.index); 
N = length(Y.index); 
%%%%% Step 3: Construct low-dimensional embeddings (Classical MDS) %%%%%
opt.disp = 0;