Indefiniteness of p-Wasserstein and Bottleneck distance based kernels

test_negative_type_simple.m

function [pos_ev_idx,neg_ev_idx,seed,data] = test_negative_type_simple(options)
% TEST_NEGATIVE_TYPE_SIMPLE evaluate the eigenvalues of Gram matrices 
%   
%   This function evaluates the number of positive/negative eigenvalues 
%   of the Gram matrix M, either constructed from the (1) q-Wasserstein 
%   distance or (2) the bottleneck distance.
%
%   Both distances lead to nonnegative (symmetric) Gram matrices and
%   the question is if the Gram matrices are (conditionally) negative
%   definite (see suppl. material). 
%
%   Input: 
%   ------
%
%   Mx1 cell array 'options', e.g.,
%
%       options{1}.q = 1;                % controls the order
%       options{1}.num_sets = 10;        % NxN Gram matrix, here: 10x10
%       options{1}.type = 'wasserstein'; % use Wasserstein, else: 'bottleneck'
%       options{1}.seed = 212575;        % RNG seed, -1 for random
%       options{1}.squared = 0;          % square the distance
%       options{1}.disp = 1;             % show results
%   
%   Calling 
%
%       >> [~,~,seed,data] = test_negative_type_simple(options);
%
%   will evaluate the eigenvalues of the 10x10 Gram matrix with elements 
%  
%       g_ij = d_{W,1}(.,.)
%
%   i.e., the 1-Wasserstein distance between two point sets.
%
%   Returns:
%   --------
%
%       pos_ev_idx  -  Positions of positive eigenvalues
%       neg_ev_idx  -  Positions of negative eigenvalues 
%       seed        -  Seed of RNG (useful, if -1 initially) 
%       data        -  (2 x 2 x num_sets) matrix of input points
% 
% Author(s): Anonymous

M = length(options);
    for m=1:M
        
        current_options = options{m};
        % if seed == -1, then generate seed and try ...
        if current_options.seed == -1
            seed = round(sum(100*clock));
            current_options.seed = seed;
        end
        seed = current_options.seed;
        rng(current_options.seed);

        % create data
        data = round(rand(2, 2, current_options.num_sets) * 10) - 1;

        % compute 
        distance_matrix = zeros(current_options.num_sets);
        edge_weights = zeros(2);
        for setA = 1 : current_options.num_sets
            for setB = 1 : current_options.num_sets
                for pointA = 1:2
                    for pointB = 1:2
                        startPoint = squeeze(data(:,pointA,setA));
                        endPoint = squeeze(data(:,pointB,setB));
                        edge_weights(pointA,pointB) = norm( startPoint -  endPoint, inf);
                    end
                end

                if strcmp(current_options.type,'wasserstein')
                    distance_matrix(setA,setB) = min( ...
                        edge_weights(1,2)^current_options.q + ...
                        edge_weights(2,1)^current_options.q, ...
                        edge_weights(1,1)^current_options.q + ...
                        edge_weights(2,2)^current_options.q)^(1/current_options.q);
                elseif strcmp(current_options.type,'bottleneck')
                    distance_matrix(setA,setB) = min( ...
                        max( edge_weights(1,2), edge_weights(2,1) ), ...
                        max( edge_weights(1,1), edge_weights(2,2)));
                else
                    error('unknown distance');
                end
             end
        end

        % square distance or not
        if current_options.squared
            gram_matrix = distance_matrix.^2;
        else
            gram_matrix = distance_matrix;
        end

        assert(length(find((gram_matrix>=0)==1)) == numel(gram_matrix));

        spectrum_of_gram_matrix = eig( gram_matrix );
        
        % # positive and negative eigenvalues
        pos_ev_idx = find(spectrum_of_gram_matrix>0);
        neg_ev_idx = find(spectrum_of_gram_matrix<0);

        if current_options.disp
            fprintf('---------------------------\n');
            fprintf('Counterexample #%d\n',m);
            if (strcmp(current_options.type,'bottleneck'))
                fprintf('Distance: d_B\n');
            else
                fprintf('Distance: d_{W,%d}\n', ...
                current_options.q);
            end
            fprintf('Gram matrix size: %d x %d\n', ...
                current_options.num_sets, ...
                current_options.num_sets);
            fprintf('---------------------------\n');
            fprintf('#pos. EVs: %d\n', length(pos_ev_idx));
            for n=1:length(pos_ev_idx)
                fprintf('\t%.10f\n', spectrum_of_gram_matrix(pos_ev_idx(n)));
            end
            fprintf('#neg. EVs: %d\n', length(neg_ev_idx));
            for n=1:length(neg_ev_idx)
                fprintf('\t%.10f\n', spectrum_of_gram_matrix(neg_ev_idx(n)));
            end
        end
    end
end