mds_cs.m

这个为模式识别工具箱

% MDS_CS Classical scaling
%
% 	W = MDS_CS(D,N)
% 
% INPUT
% 	D 	Square dissimilarity matrix of the size M x M
% 	N 	Desired output dimensionality (optional; default: 2)
%
% OUTPUT
% 	W		Classical scaling mapping 
%
% DESCRIPTION  
% A linear mapping of objects given by a symmetric distance matrix D with
% a zero diagonal onto an N-dimensional Euclidean space such that the square 
% distances are preserved as much as possible. 
%
% D is assumed to approximate the Euclidean distances, i.e. 
% D_{ij} = sqrt(sum_k (x_i-x_j)^2). 
%
%
% REFERENCES
% 1. I. Borg and P. Groenen, Modern Multidimensional Scaling, Springer Verlag, Berlin, 
% 	 1997.
% 2. T.F. Cox and M.A.A. Cox, Multidimensional Scaling, Chapman and Hall, London, 1994.
%
% SEE ALSO
% MAPPINGS, MDS_INIT, MDS_CS


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

function W = mds_cs (D,n)

	if nargin < 2, n = 2; end
	if nargin < 1 | isempty(D)
		W = mapping(mfilename,n);
		W = setname(W,'Classical MDS');
		return
	end

	[m,mm] = size(D);


	if (isa(D,'dataset') | isa(D,'double')) 
		if isa(n,'mapping')
			% The projection is already defined; apply it to new data
			w = getdata(n);
			[m,n] = size(D);
	      if isa(D,'dataset'),
				lab = getlab(D); 
				D = +D;
			else
				lab = [];
			end
			D = remove_nan(D);  
			D     = D.^2;
			Q     = w{1}; 
			me    = w{2};
			L     = w{3};
	
			Z = -repmat(1/n,n,n);                 
			Z(1:n+1:end) = Z(1:n+1:end) + 1;    	 % Z = eye(n) - ones(n,n)/n
			data = -0.5 * (D - me(ones(m,1),:)) * Z * Q * diag(sqrt(abs(L))./L);          
			if ~isempty(lab), 
				W = dataset(data,lab); 
			else
				W = data;
			end
			return; 
		end
	end


	% Definition of the projection; we are here only if there is no projection yet
	if ~issym(D,1e-12),
		error('Matrix should be symmetric.')
	end


  D = remove_nan(D);  
	D = +D.^2;
	B = -repmat(1/m,m,m);					
	B(1:m+1:end) = B(1:m+1:end) + 1;		% B = eye(m) - ones(m,m)/m
	B = -0.5 * B * D * B; 							% B is now the matrix of inner products 

	ok = 0;
	p  = 1;
	k  = min (p*n,m);
	options.disp   = 0; 
	options.isreal = 1; 
	options.issym  = 1; 

	while ~ok & k <= m, 								% Find k largest eigenvalues
		k = min (p*n,m);

		[V,S,U] = eigs(B,k,'lm',options);
		s = diag(real(S));
		[ss,I] = sort(-s);
		ok = sum(s>0) >= n;
		p  = p+1;
	end
	if ~ok & k == m,
		error ('The wanted configuration cannot be found.');
	end
	J = I(1:n);
	W = mapping(mfilename,'trained',{V(:,J), mean(D,2)', s(J)},[],m,length(J));
	W = setname(W,'Classical MDS');
return









function D = remove_nan(D)
% Remove all appearing NaN's by replacing them by the nearest neighbors.
%
[m,mm] = size(D);

nanindex = find(isnan(D(:)));  
if ~isempty(nanindex),
  for i=1:m
    K = find(isnan(D(i,:)));
    I = 1:mm; 
    [md,kk] = min(D(i,:)); 
    if md < eps, 
      I(kk) = [];
    end
    D(i,K) = min(D(i,I));
  end
  % check whether the diagonal is of all zeros
  if length(intersect(find(D(:) < eps), 1:m+1:(m*m))) >= m,
    D = (D+D')/2;
  end
end