mds.m

The pattern recognition matlab toolbox

        % Plot progress if requested.

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

        if (success)
          sigma = sigma0/sqrt(pnorm2);      % SIGMA: a small step from y to yy
          yy    = y + sigma .*p;
          [e,a] = mds_samstress(opt.q,yy,yrep,Ds,[],opt.nanid,opt.isratio); 
          g2    = mds_gradients(opt.q,yy,yrep,a*Ds,[],opt.nanid); % Gradient for yy
          s     = (g2-g1)/sigma;        % Approximation of Hessian*P directly, instead of computing
          delta = p(:)' * s(:);         % the Hessian, since only DELTA = P'*Hessian*P is in fact needed
        end

        % Regularize the Hessian indirectly to make it positive definite.
        % DELTA is now computed as P'* (Hessian + regularization * Identity) * P
        delta = delta + (lambda1 - lambda) * pnorm2;  

        % Indicate if the Hessian is negative definite; the regularization above was too small.
        if (delta < 0)
          lambda1 = 2 * (lambda - delta/pnorm2);  % Now the Hessian will be positive definite
          delta   = -delta + lambda * pnorm2;     % This is obtained after plugging lambda1 
          lambda  = lambda1;                      % into the formulation a few lines above.
        end

        mi = - p(:)' * g1(:);
        yy = y + (mi/delta) .*p;                  % mi/delta is a step size 
        [ee,a] = mds_samstress(opt.q,yy,yrep,Ds,[],opt.nanid,opt.isratio); 

        % This minimization procedure is based on the second order approximation of the  
        % stress function by using the gradient and the Hessian approximation. The Hessian 
        % is regularized, but maybe not sufficiently. The ratio Dc (<=1) below indicates 
        % a proper approximation, if Dc is close to 1.

        Dc = 2 * delta/mi^2 * (e - ee); 
        e  = ee;

        % If Dc > 0, then the stress can be successfully decreased.
        success = (Dc >= 0);
        if (success)
          y       = yy;
          [ee,a]  = mds_samstress(opt.q,yy,yrep,Ds,[],opt.nanid,opt.isratio); 
          g1      = mds_gradients(opt.q,yy,yrep,a*Ds,[],opt.nanid); 
          gnorm2  = g1(:)' * g1(:);
          lambda1 = 0;

          beta = max(0,(g1(:)'*(g1(:)-g0(:)))/mi);  
          p   = -g1 + beta .* p;                 % P is a new conjugate direction  

          if (g1(:)'*p(:) >= 0 | mod(it-1,n*m) == 0),
            p = -g1;                            % No much improvement, restart
          end   

          if (Dc >= 0.75)
            lambda = 0.25 * lambda;             % Make the regularization smaller
          end

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

        else        % Dc < 0  
          % Note that for Dc < 0, the iteration number IT is not increased, 
          % so the Hessian is further regularized until SUCCESS equals 1.
          lambda1 = lambda;
        end

      % The approximation of the Hessian was poor or the stress was not 
      % decreased (Dc < 0), hence the regularization lambda is enlarged.
      if (Dc < 0.25)
        lambda = lambda + delta * (1 - Dc)/pnorm2;
      end
    end
  end

return;



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

%MDS_STRESS - Calculate Sammon stress during optimization
%
%   E = MDS_SAMSTRESS(Q,Y,YREP,Ds,D,NANINDEX)
%
% INPUT
%   Q          Indicator of the Sammon stress; Q = -2,-1,0,1,2
%   Y         Current lower-dimensional configuration
%   YREP      Configuration of the representation objects; it should be
%             empty when no representation set is considered
%   Ds        Original distance matrix
%   D         Approximate distance matrix (optional; otherwise computed from Y)
%   NANINDEX  Index of the missing values; marked in Ds by NaN (optional; to
%               be found in Ds)
%
% OUTPUT
%   E         Sammon stress
%
% DESCRIPTION
% Computes the Sammon stress between the original distance matrix Ds and the
% approximated distance matrix D between the mapped configuration Y and the
% configuration of the representation set YREP, expressed as follows:
%
% E = 1/(sum_{i<j} Ds_{ij}^(q+2)) sum_{i<j} (Ds_{ij} - D_{ij})^2 * Ds_{ij}^q
%
% It is directly used in the MDS_SAMMAP routine.

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

function [e,alpha] = mds_samstress (q,y,yrep,Ds,D,nanindex,isratio)

  % If D is not given, calculate it from Y, the current mapped points.

  if (nargin < 5) | (isempty(D))
    if (~isempty(yrep))
      D = sqrt(distm(y,yrep));     
    else
      D = sqrt(distm(y));     
    end
  end

  if (nargin < 6)
    nanindex = [];    % Not given, so calculate below.
  end

  if (nargin < 7)
    isratio = 0;      % Assume this is meant.
  end

  todefine = isempty(yrep);

  [m,n]  = size(y); [mm,k] = size(Ds);
  if (m ~= mm)
    error('The sizes of Y and Ds do not match.');
  end
  if (any(size(D) ~= size(Ds)))
    error ('The sizes of D and Ds do not match.');
  end
  m2 = m*k;

  % Convert to double.
  D  = +D; Ds = +Ds;

  % I is the index of non-NaN, non-zero (> eps) values to be included 
  % for the computation of the stress.

  I = 1:m2; 
  if (~isempty(nanindex))
    I(nanindex) = [];
  end

  O = [];
  if (todefine), O = 1:m+1:m2; end                               
  I = setdiff(I,O);

  % If OPTIONS.isratio is set, calculate optimal ALPHA to scale with.

  if (isratio)
    alpha = sum((Ds(I).^q).*D(I).^2)/sum((Ds(I).^(q+1)).*D(I));
    Ds    = alpha*Ds;
  else
    alpha = 1; 
  end
    
  % C is the normalization factor.
  c = sum(Ds(I).^(q+2)); 

  % If Q == 0, prevent unnecessary calculation (X^0 == 1).
  if (q ~= 0)
    e = sum(Ds(I).^q .*((Ds(I)-D(I)).^2))/c;
  else
    e = sum(((Ds(I)-D(I)).^2))/c;
  end

return; 



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

%MDS_GRADIENTS - Gradients for variants of the Sammon stress
%
%   [G1,G2,CC] = MDS_GRADIENTS(Q,Y,YREP,Ds,D,NANINDEX)
%
% INPUT
%   Q          Indicator of the Sammon stress; Q = -2,-1,0,1,2
%   Y         Current lower-dimensional configuration
%   YREP      Configuration of the representation objects; it should be
%             empty when no representation set is considered
%   Ds        Original distance matrix
%   D         Approximate distance matrix (optional; otherwise computed from Y)
%   nanindex  Index of missing values; marked in Ds by NaN (optional;
%             otherwise found from Ds)
%
% OUTPUT
%   G1        Gradient direction
%   G2        Approximation of the Hessian by its diagonal  
%
% DESCRIPTION  
% This is a routine used directly in the MDS_SAMMAP routine.
%
% SEE ALSO
% MDS, MDS_INIT, MDS_SAMMAP, MDS_SAMSTRESS

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

function [g1,g2,c] = mds_gradients(q,y,yrep,Ds,D,nanindex)

  % If D is not given, calculate it from Y, the current mapped points.

  if (nargin < 5) | (isempty(D))
    if (~isempty(yrep))
      D = sqrt(distm(y,yrep));     
    else
      D = sqrt(distm(y));     
    end
  end

  % If NANINDEX is not given, find it from Ds.

  if (nargin < 6)
    nanindex = find(isnan(Ds(:))==1);  
  end

  % YREP is empty if no representation set is defined yet.
  % This happens when the mapping should be defined.

  todefine = (isempty(yrep));
  if (todefine) 
    yrep  = y;
  end

  [m,n]  = size(y); [mm,k] = size(Ds);
  if (m ~= mm)
    error('The sizes of Y and Ds do not match.');
  end
  if (any(size(D) ~= size(Ds)))
    error ('The sizes of D and Ds do not match.');
  end
  m2 = m*k;

  % Convert to doubles.

  D  = +D; Ds = +Ds;

  % I is the index of non-NaN, non-zero (> eps) values to be included 
  % for the computation of the gradient and the Hessians diagonal.

  I = (1:m2)';
  if (~isempty(nanindex))
    I(nanindex) = [];
  end
  K    = find(Ds(I) <= eps | D(I) <= eps);
  I(K) = [];

  % C is a normalization factor.

  c = -2/sum(Ds(I).^(q+2)); 

  % Compute G1, the gradient.

  h1 = zeros(m2,1);
  if (q == 0)                    % Prevent unnecessary computation when Q = 0.
    h1(I) = (Ds(I)-D(I)) ./ D(I);  
  else  
    h1(I) = (Ds(I)-D(I)) ./ (D(I).*(Ds(I).^(-q)));
  end
  h1 = reshape (h1',m,k);
  g2 = h1 * ones(k,n);          % Here G2 is assigned only temporarily,
  g1 = c * (g2.*y - h1*yrep);    % for the computation of G1.     

  % Compute G2, the diagonal of the Hessian, if requested.

  if (nargout > 1)
    h2 = zeros(m2,1);
    switch (q)
      case -2, 
        h2(I) = -1./(Ds(I).*D(I).^3);
      case -1, 
        h2(I) = -1./(D(I).^3);
      case 0, 
        h2(I) = - Ds(I)./(D(I).^3);
      case 1, 
        h2(I) = - Ds(I).^2./(D(I).^3);
      case 2, 
        h2(I) = -(Ds(I)./D(I)).^3;
      end 
    h2 = reshape (h2',m,k);
    g2 = c * (g2 + (h2*ones(k,n)).*y.^2 + h2*yrep.^2 - 2*(h2*yrep).*y);  
  end

return;



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

%MDS_PLOT Plots the results of the MDS mapping in a 2D or 3D figure
% 
%   MDS_PLOT(Y,YREP,LAB,REPLAB,E)
% 
% INPUT
%   Y       Configuration in 2D or 3D space
%   YREP    Configuration of the representation points in 2D or 3D space
%   LAB     Labels of Y
%   REPLAB  Labels of YREP
%   E       Stress value
%
% DESCRIPTION  
% Used directly in MDS_SAMMAP routine.


function mds_plot (y,yrep,lab,replab,e)

  % This is done in a rather complicated way in order to speed up drawing.

  K = min(size(y,2),3);
  if (K > 1)
    y   = +y;
    col = 'brmk';
    sym = ['+*xosdv^<>p']';
    ii  = [1:44];
    s   = [col(ii-floor((ii-1)/4)*4)' sym(ii-floor((ii-1)/11)*11)];
    [nlab,lablist] = renumlab([replab;lab]); c = max(nlab);
    [ms,ns] = size(s); if (ms == 1), s = setstr(ones(m,1)*s); end

    yy = [yrep; y];  

    cla; 
    if (K == 2)
      hold on;
      for i = 1:c
        J = find(nlab==i); plot(yy(J,1),yy(J,2),s(i,:));
      end
    else 
      for i = 1:c
        J = find(nlab==i); plot3(yy(J,1),yy(J,2),yy(J,3),s(i,:));
        if (i == 1), hold on; grid on; end
      end
    end

    title(sprintf('STRESS: %f', e));
    axis equal; 
     drawnow;
  end

return;