Dissimilarity embedding for indefinite inner product space

Given a (symmetric) dissimilarity matrix, the function computes the embedding of the matrix into the potentially indefinite inner product space defined through the matrix, which in general is given by metrics between vectors.

The embedding is generated in a manner such that it keeps the distances encoded within the entries of the matrix, which can be generated with different metrics such as Ln, maximum or even fractional ones.

References to literature are given.

Follow-ups (to come): Projection of the resulting Xnk - configuration Matrix (see source code documentation) onto a Pontryagin space, allowing to represent new vectors in the indefinite inner product space given by the matrix.

The language is Octave and can easily be changed to Matlab: Just substitute all endXXXX by end

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