laplacian_eigen.m

The goal of SPID is to provide the user with tools capable to simulate, preprocess, process and clas

function mappedX = laplacian_eigen(X, no_dims, k, sigma, eig_impl)
%LAPLACIAN_EIGEN Performs non-linear dimensionality reduction using Laplacian Eigenmaps
%
%   mappedX = laplacian_eigen(X, no_dims, k, sigma, eig_impl)
%
% Performs non-linear dimensionality reduction using Laplacian Eigenmaps.
% The data is in matrix X, in which the rows are the observations and the
% columns the dimensions. The variable dim indicates the preferred amount
% of dimensions to retain (default = 2). The variable k is the number of 
% neighbours in the graph (default = 12).
% The reduced data is returned in the matrix mappedX.
%
%

% This file is part of the Matlab Toolbox for Dimensionality Reduction v0.1b.
% The toolbox can be obtained from http://www.cs.unimaas.nl/l.vandermaaten
% You are free to use, change, or redistribute this code in any way you
% want. However, it is appreciated if you maintain the name of the original
% author.
%
% (C) Laurens van der Maaten
% Maastricht University, 2007

    if ~exist('no_dims', 'var')
        no_dims = 2;
    end
    if ~exist('k', 'var')
        k = 12;
    end
	if ~exist('sigma', 'var')
		sigma = 1;
    end
    if ~exist('eig_impl', 'var')
        eig_impl = 'Matlab';
    end
    
    % Construct neighborhood graph
    disp('Constructing neighborhood graph...');
    if size(X, 1) < 4000
        G = squareform(pdist(X, 'euclidean'));
        % Compute neighbourhood graph
        [tmp, ind] = sort(G); 
        for i=1:size(G, 1)
            G(i, ind((2 + k):end, i)) = 0; 
        end
        G = sparse(G);
        G = max(G, G');             % Make sure distance matrix is symmetric
    else
        G = find_nn(X, k);
    end
    G = G .^ 2;
	G = G ./ max(max(G));
    
    % Compute weights (W = G)
    disp('Computing weight matrices...');
    
    % Compute Gaussian kernel (heat kernel-based weights)
    G(G ~= 0) = exp(-G(G ~= 0) / (2 * sigma ^ 2));
        
    % Construct diagonal weight matrix
    D = diag(sum(G, 2));
    
    % Compute Laplacian
    L = D - G;
    L(isnan(L)) = 0; D(isnan(D)) = 0;
	L(isinf(L)) = 0; D(isinf(D)) = 0;
    
    % Construct eigenmaps (solve Ly = labda*Dy)
    disp('Constructing Eigenmaps...');
	tol = 0;
    if strcmp(eig_impl, 'JDQR')
        options.Disp = 0;
        options.LSolver = 'bicgstab';
        [mappedX, lambda] = jdqz(L, D, no_dims + 1, tol, options);			% only need bottom (no_dims + 1) eigenvectors
    else
        options.disp = 0;
        options.isreal = 1;
        options.issym = 1;
        [mappedX, lambda] = eigs(L, D, no_dims + 1, tol, options);			% only need bottom (no_dims + 1) eigenvectors
    end
    
    % Sort eigenvectors in ascending order
    lambda = diag(lambda);
    [tmp, ind] = sort(lambda, 'ascend');
    
    % Final embedding
	mappedX = mappedX(:,ind(2:no_dims + 1));