function [y, E] = sammon(x, n, opts)
%
% SAMMON - apply Sammon's nonlinear mapping
%
% Y = SAMMON(X) applies Sammon's nonlinear mapping procedure on
% multivariate data X, where each row represents a pattern and each column
% represents a feature. On completion, Y contains the corresponding
% co-ordinates of each point on the map. By default, a two-dimensional
% map is created. Note if X contains any duplicated rows, SAMMON will
% fail (ungracefully).
%
% [Y,E] = SAMMON(X) also returns the value of the cost function in E (i.e.
% the stress of the mapping).
%
% An N-dimensional output map is generated by Y = SAMMON(X,N) .
%
% A set of optimisation options can also be specified using a third
% argument, Y = SAMMON(X,N,OPTS) , where OPTS is a structure with fields:
%
% MaxIter - maximum number of iterations
% TolFun - relative tolerance on objective function
% MaxHalves - maximum number of step halvings
% Input - {'raw','distance'} if set to 'distance', X is
% interpreted as a matrix of pairwise distances.
% Display - {'off', 'on', 'iter'}
% Initialisation - {'pca', 'random'}
%
% The default options structure can be retrieved by calling SAMMON with
% no parameters.
%
% References :
%
% [1] Sammon, John W. Jr., "A Nonlinear Mapping for Data Structure
% Analysis", IEEE Transactions on Computers, vol. C-18, no. 5,
% pp 401-409, May 1969.
%
% See also : SAMMON_TEST
%
% File : sammon.m
%
% Date : Monday 12th November 2007.
%
% Author : Gavin C. Cawley and Nicola L. C. Talbot
%
% Description : Simple vectorised MATLAB implementation of Sammon's non-linear
% mapping algorithm [1].
%
% References : [1] Sammon, John W. Jr., "A Nonlinear Mapping for Data
% Structure Analysis", IEEE Transactions on Computers,
% vol. C-18, no. 5, pp 401-409, May 1969.
%
% History : 10/08/2004 - v1.00
% 11/08/2004 - v1.10 Hessian made positive semidefinite
% 13/08/2004 - v1.11 minor optimisation
% 12/11/2007 - v1.20 initialisation using the first n principal
% components.
%
% Thanks : Dr Nick Hamilton (nick@maths.uq.edu.au) for supplying the
% code for implementing initialisation using the first n
% principal components (introduced in v1.20).
%
% To do : The current version does not take advantage of the symmetry
% of the distance matrix in order to allow for easy
% vectorisation. This may not be a good choice for very large
% datasets, so perhaps one day I'll get around to doing a MEX
% version using the BLAS library etc. for very large datasets.
%
% Copyright : (c) Dr Gavin C. Cawley, November 2007.
%
% This program is free software; you can redistribute it and/or modify
% it under the terms of the GNU General Public License as published by
% the Free Software Foundation; either version 2 of the License, or
% (at your option) any later version.
%
% This program is distributed in the hope that it will be useful,
% but WITHOUT ANY WARRANTY; without even the implied warranty of
% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
% GNU General Public License for more details.
%
% You should have received a copy of the GNU General Public License
% along with this program; if not, write to the Free Software
% Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
%
% use the default options structure
if nargin < 3
opts.Display = 'iter';
opts.Input = 'raw';
opts.MaxHalves = 20;
opts.MaxIter = 500;
opts.TolFun = 1e-9;
opts.Initialisation = 'random';
end
% the user has requested the default options structure
if nargin == 0
y = opts;
return;
end
% create a two-dimensional map unless dimension is specified
if nargin < 2
n = 2;
end
% set level of verbosity
if strcmp(opts.Display, 'iter')
display = 2;
elseif strcmp(opts.Display, 'on')
display = 1;
else
display = 0;
end
% create distance matrix unless given by parameters
if strcmp(opts.Input, 'distance')
D = x;
else
D = euclid(x, x);
end
% remaining initialisation
N = size(x, 1);
scale = 0.5/sum(D(:));
D = D + eye(N);
Dinv = 1./D;
if strcmp(opts.Initialisation, 'pca')
[UU,DD] = svd(x);
y = UU(:,1:n)*DD(1:n,1:n);
else
y = randn(N, n);
end
one = ones(N,n);
d = euclid(y,y) + eye(N);
dinv = 1./d;
delta = D - d;
E = sum(sum((delta.^2).*Dinv));
% get on with it
for i=1:opts.MaxIter
% compute gradient, Hessian and search direction (note it is actually
% 1/4 of the gradient and Hessian, but the step size is just the ratio
% of the gradient and the diagonal of the Hessian so it doesn't matter).
delta = dinv - Dinv;
deltaone = delta*one;
g = delta*y - y.*deltaone;
dinv3 = dinv.^3;
y2 = y.^2;
H = dinv3*y2 - deltaone -2*y.*(dinv3*y) + y2.*(dinv3*one);
s = -g(:)./abs(H(:));
y_old = y;
% use step-halving procedure to ensure progress is made
for j=1:opts.MaxHalves
y(:) = y_old(:) + s;
d = euclid(y,y) + eye(N);
dinv = 1./d;
delta = D - d;
E_new = sum(sum((delta.^2).*Dinv));
if E_new < E
break;
else
s = 0.5*s;
end
end
% bomb out if too many halving steps are required
if j == opts.MaxHalves
if display
fprintf(1, 'Warning : MaxHalves exceeded.\n');
end
break;
end
% evaluate termination criterion
if abs((E - E_new)/E) < opts.TolFun
if display
fprintf(1, 'Optimisation terminated - TolFun exceeded.\n');
end
break;
end
% report progress
E = E_new;
if display > 1
fprintf(1, 'epoch = %d : E = %12.10f\n', i, E*scale);
end
end
% fiddle stress to match the original Sammon paper
E = E*scale;
% all done
function d = euclid(x,y)
d = sqrt(sum(x.^2,2)*ones(1,size(y,1))+ones(size(x,1),1)*sum(y.^2,2)'-2*(x*y'));
