Code

Dissimilarity embedding for indefinite inner product space

Andreas Beschorner November 29, 2010 Coded in Octave
function  [Xnk, lambda, calX, sig, T] = IIPSembedd(d)
% syntax: [Xnk, lambda, calX, sig, T] = IIPSembedd(d)
% creates indefinite inneer product space embedding for square dissimilarity matrix d
% Xnk will be the embedding/ mapping matrix, sig the signature of the 
% resulting (in)definite inner product space.
% calX is the matrix of ordered eigenvectors, lambda the associated vector
% of (sortted) eigenvalues.
% T is the gramian in the indefinite inner product space.
% Reference:
%   Pekalska, Elzbieta et al: A Generalized Kernel Approach to Dissimilarity-based Classification
%		Journal of Machine Learning Research 2 (2001), pp. 175-211
%	Pekalska, Elzbieta et al: The Dissimilarity Representation for Pattern Recognition
%		World Scientific, Singapore, December 2005
%	Goldfarb, Lev: A New Approach to Pattern Recogntion, 1985
%		
	n = size(d, 1);
	D2 = d.^2;                          % elementwise squaring
	J = zeros(n) + 1/n;                 % 1/n * [1] - Matrix
	J = eye(n) - J;                     % get centering matrix ...
	T = -0.5 *J * D2 * J;               % ... and Gramian
	T = 0.5 * (T + T');                 % symmetrize T (e.g. numerical errors)

% get Eigenvalue diagonal matrix and corresponding Eigenvector-matrix
	[calX, lambda] = eig(T);            % todo: lambda=real(lambda)?
	lambda=real(lambda);                % just in case. should always be real for a symmetric dissimilarity matrix 'd'.
	evs = diag(lambda);                 % get vector of eigenvalues from diagonal matrix
	[sev, idx] = sort(-evs);            % sort!
	sgn = sign(sev);                    % get dimension information for sorting
	ll = find(sgn>=0, 1);
	if (isempty(ll))
		ll = length(sgn);
	endif
	sig = ([length(sgn)-ll, ll]);

% Sort eigenvalues and associated eigenvectors:
% first positive EVs (descending), follwed by zero and negatives
% (descending w.r.t. absolute value).
	evs = [-sev(1:ll-1); -sev(end:-1:ll)];
	lambda=evs;
	idx = [idx(1:ll-1); idx(end:-1:ll)];
	calX = calX(:, idx);

% Selection of important EVs, for instance by threshold 0.00001, can be inserted here

% Compute final configuration
	Xnk=calX * diag(sqrt(abs(evs)));
endfunction