mds.m

The pattern recognition matlab toolbox

function printinfo(opt)

  fprintf(opt.st,'Sammon mapping, error function with the parameter q=%d\n',opt.q);

  switch (opt.optim)
    case 'pn',  
      fprintf(opt.st,'Minimization by Pseudo-Newton algorithm\n');
    case 'scg',  
      fprintf(opt.st,'Minimization by Scaled Conjugate Gradients algorithm\n');
    otherwise 
      error(strcat('Possible initialization methods: pn (Pseudo-Newton), ',...
                   'or scg (Scaled Conjugate Gradients).'));
  end

return

% **********************************************************************************

% [I,J,P] = MDS_REPPOINTS(D)
%
% Finds the indices of repeated/left points. J contains the indices of
% repeated points in D. This means that for each j in J, there is a point
% k ~= j such that D(J(j),k) = 0. I contains the indices of the remaining 
% points, and P those of the points in I that were identical to those in J.
% Directly used in MDS routine.

function [I,J,P] = mds_reppoints(D)

  epsilon = 1e-20;      % Differences smaller than this are assumed to be zero.

  [m,mm] = size(D);

  I = 1:m; J = []; P = [];

  if (m == mm) & (all(abs(D(1:m+1:end)) <= epsilon))
    K = intersect (find (triu(ones(m),1)), find(D < epsilon));
    if (~isempty(K))
      P = mod(K,m);
      J = fix(K./m) + 1;           
      I(J) = [];                    
    end
  else
    [J,P] = find(D<=epsilon); 
    I(J)  = [];
  end

return;

% **********************************************************************************


% MDS_SETOPT  Sets parameters for the MDS mapping
%
%   OPT = MDS_SETOPT(OPT_GIVEN)
% 
% INPUT
%   OPT_GIVEN Parameters for the MDS mapping, described below (default: [])
%
% OUTPUT
%   OPT       Structure of chosen options; if OPT_GIVEN is empty, then 
%               OPT.q       = 0
%               OPT.optim   = 'pn'
%               OPT.init    = 'cs'
%               OPT.etol    = 1e-6 (the precise value depends on q)
%               OPT.maxiter = inf 
%               OPT.itmap   = 'yes' 
%               OPT.isratio = 0 
%               OPT.st      = 1 
%               OPT.inspect = 2 
%
% DESCRIPTION  
% Parameters for the MDS mapping can be set or changed. OPT_GIVEN consists
% of the following fields: 'q', 'init', 'optim', 'etol', 'maxiter','itmap',
% 'isratio', 'st' and 'inspect'. OPTIONS can include all the fields or some
% of them only. The fields of OPT have some default values, which can be
% changed by the OPT_GIVEN field values. If OPT_GIVEN is empty, then OPT
% contains all default values. For a description of the fields, see MDS.

%
% Copyright: Elzbieta Pekalska, [e-mail], 2000-2003
% Faculty of Applied Sciences, Delft University of Technology
%

function opt = mds_setopt (opt_given)

  opt.q       = 0;
  opt.init    = 'cs'; 
  opt.optim   = 'pn';
  opt.st      = 1; 
  opt.itmap   = 'yes'; 
  opt.maxiter = inf; 
  opt.inspect = 2; 
  opt.etol    = inf; 
  opt.isratio = 0; 
  opt.nanid   = [];    % Here are some extra values; set up in the MDS routine.
  opt.plotlab = [];    % Not to be changed by the user.

  if (~isempty(opt_given))
    if (~isstruct(opt_given))
      error('OPTIONS should be a structure with at least one of the following fields: q, init, etol, optim, maxiter, itmap, isratio, st or inspect.');
    end
    fn = fieldnames(opt_given);
    if (~all(ismember(fn,fieldnames(opt))))
      error('Wrong field names; valid field names are: q, init, optim, etol, maxiter, itmap, isratio, st or inspect.')
    end 
    for i = 1:length(fn)
      opt = setfield(opt,fn{i},getfield(opt_given,fn{i}));
    end
  end

  if (isempty(intersect(opt.q,-2:1:2)))
    error ('OPTIONS.q should be -2, -1, 0, 1 or 2.');
  end

  if (opt.maxiter < 2)
    error ('OPTIONS.iter should be at least 1.');
  end

  if (isinf(opt.etol))
    switch (opt.q)                 % Different defaults for different stresses.
      case -2, 
        opt.etol = 1e-6;
      case {-1,0}
        opt.etol = 10*sqrt(eps); 
      case {1,2}
        opt.etol = sqrt(eps); 
      end 
  elseif (opt.etol <= 0) | (opt.etol >= 0.1)
    error ('OPTIONS.etol should be positive and smaller than 0.1.');
  end

  if (~ismember(opt.optim, {'pn','scg'}))
    error('OPTIONS.optim should be pn or scg.');
  end

  if (~ismember(opt.itmap, {'yes','no'}))
    error('OPTIONS.itmap should be yes or no.');
  end

return



% **********************************************************************************

% MDS_SAMMAP Sammon iterative nonlinear mapping for MDS
%
%   [YMAP,ERR] = MDS_SAMMAP(D,Y,YREP,OPTIONS)
% 
% INPUT
%   D       Square (M x M) dissimilarity matrix
%   Y       M x N matrix containing starting configuration, or
%   YREP    Configuration of the representation points
%   OPTIONS Various parameters of the  minimization procedure put into 
%           a structure consisting of the following fields: 'q', 'optim',
%           'init','etol','maxiter', 'isratio', 'itmap', 'st' and 'inspect'
%           (default:
%             OPTIONS.q       = 0          
%             OPTIONS.optim   = 'pn'      
%             OPTIONS.init    = 'cs'
%             OPTIONS.etol    = 1e-6 (the precise value depends on q)
%             OPTIONS.maxiter = inf 
%             OPTIONS.isratio = 0 
%             OPTIONS.itmap   = 'yes' 
%             OPTIONS.st      = 1 
%             OPTIONS.inspect = 2).
% 
% OUTPUT
%   YMAP     Mapped configuration
%   ERR     Sammon stress 
%
% DESCRIPTION  
% Maps the objects given by a symmetric distance matrix D (with a zero
% diagonal) onto, say, an N-dimensional configuration YMAP by an iterative
% minimization of a variant of the Sammon stress. The minimization starts
% from the initial configuration Y; see MDS_INIT.
%
% YREP is the Sammon configuration of the representation set. It is used
% when new points have to be projected. In other words, if D is an M x M
% symmetric distance matrix, then YREP is empty; if D is an M x N matrix,
% then YMAP is sought such that D can approximate the distances between YMAP
% and YREP.
%
% Missing values can be handled by marking them by NaN in D.
%
% SEE ALSO
% MAPPINGS, MDS, MDS_CS, MDS_INIT, MDS_SETOPT

%
% Copyright: Elzbieta Pekalska, [e-mail], 2000-2003
% Faculty of Applied Sciences, Delft University of Technology
%

function [y,err] = mds_sammap(Ds,y,yrep,opt)

  if (nargin < 4)
    opt = [];                 % Will be filled by MDS_SETOPT, below.
  end
  if isempty(opt), opt = mds_setopt(opt); end
  if (nargin < 3)
    yrep = []; 
  end

  % Extract labels and calculate distance matrix.

  [m,n] = size(y);
  if (~isempty(yrep))
    replab = getlab(yrep);
    yrep = +yrep;
    D = sqrt(distm(y,yrep)); 
  else
    D = sqrt(distm(y));     
    yrep = [];
    replab = [];
  end

  if (isempty(opt.plotlab))
    opt.plotlab = ones(m,1);
  end

  it   = 0;           % Iteration number.
  eold = inf;         % Previous stress (error).

  % Calculate initial stress.

  [e,a]= mds_samstress(opt.q,y,yrep,Ds,D,opt.nanid,opt.isratio); err = e;

  fprintf(opt.st,'iteration: %4i   stress: %3.8d\n',it,e);  

  switch (opt.optim)  
    case 'pn',        % Pseudo-Newton minimization.

      epr   = e;      % Previous values of E, EOLD for line search algorithm.
      eepr  = e;
      add   = 1e-4;   % Avoid G2 equal 0; see below.
      
      BETA   = 1e-3;   % Parameter for line search algorithm.
      lam1   = 0;     
      lam2   = 1;      % LAM (the step factor) lies in (LAM1, LAM2).
      LMIN   = 1e-10;  % Minimum acceptable value for LAM in line search.
      lambest = 1e-4; % LAM = LAMBEST, if LAM was not found.

      % Loop until the error change falls below the tolerance or until
      % the maximum number of iterations is reached.

      while (abs(e-eold) >= opt.etol*(1 + abs(e))) & (it < opt.maxiter)

        % Plot progress if requested.

        if (opt.st > 0) & (opt.inspect > 0) & (mod(it,opt.inspect) == 0)
%          if (it == 0), figure(1); clf; end
          mds_plot(y,yrep,opt.plotlab,replab,e); 
        end

        yold = y; eold = e;

        % Calculate the stress.

        [e,a]      = mds_samstress (opt.q,yold,yrep,Ds,[],opt.nanid,opt.isratio);

        % Calculate gradients and pseudo-Newton update P.

        [g1,g2,cc] = mds_gradients (opt.q,yold,yrep,a*Ds,[],opt.nanid);
        p        = (g1/cc)./(abs(g2/cc) + add);
        slope    = g1(:)' * p(:);

        lam = 1;                 

        % Search for a suitable step LAM using a line search.

        while (1)

          % Take a step and calculate the delta error DE.

          y      = yold + lam .* p;  
          [e,a]  = mds_samstress(opt.q,y,yrep,Ds,[],opt.nanid,opt.isratio);
          de     = e - eold;

          % Stop if the condition for a suitable lam is fulfilled.

          if (de < BETA * lam * slope)
            break;
          end
   
          % Try to find a suitable step LAM.

          if (lam ~= 1)
            r1 = de - lam*slope; 
            r2 = epr - eepr - lam2*slope; 
            aa = (2*de + r1/lam^2 - r2/lam2^2)/(lam2-lam1)^3;
            bb = ((-lam2*r1)/lam^2 +(lam*r2)/lam2^2)/(lam2-lam1)^2;

            if (abs(aa) <= eps)
              lamtmp = -0.5 * slope/bb;
            else
              lamtmp = (-bb + sqrt(max(0,(bb^2 - 3*aa*slope))))/(3*aa);  
            end

            % Prevent LAM from becoming too large.

            if (lamtmp > 0.5 * lam), lamtmp = 0.5 * lam; end

          else
            lamtmp = -0.5 * slope/(e - eold - slope);
          end

          % Prevent LAM from becoming too small.

          lam2 = lam;    
          lam  = max(lamtmp, 0.1*lam);  
          epr  = e;
          eepr = eold;   

          if (lam < LMIN)
            y = yold + lambest .* p;
            [e,a] = mds_samstress(opt.q,y,yrep,Ds,[],opt.nanid,opt.isratio);
            break;
          end
        end

        it  = it + 1; 
        err = [err; e];
        fprintf(opt.st,'iteration: %4i   stress: %3.8d\n',it,e); 

      end


    case 'scg',              % Scaled Conjugate Gradient minimization.

      sigma0  = 1e-8;                        
      lambda  = 1e-8;       % Regularization parameter.
      lambda1 = 0;                                                             

      % Calculate initial stress and direction of decrease P.

      [e,a]   = mds_samstress(opt.q,y,yrep,Ds,D,opt.nanid,opt.isratio);   
      g1      = mds_gradients(opt.q,y,yrep,a*Ds,D,opt.nanid);   % Gradient
      p       = -g1;                                            % Direction of decrease
      gnorm2  = g1(:)' * g1(:);
      success = 1;                                              

      % Sanity check.

      if (gnorm2 < 1e-15)
        prwarning(2,['Gradient is nearly zero: ' gnorm2 ' (unlikely); initial configuration is the right one.']);  
        return
      end

      % Loop until the error change falls below the tolerance or until
      % the maximum number of iterations is reached.

      while (abs(eold - e) > opt.etol * (1 + e)) & (it < opt.maxiter)

        g0     = g1;            % Previous gradient
        pnorm2 = p(:)' * p(:);
        eold   = e;