Reproducing results with Shogun

[jeslago]

Hi,

I was one of the attendant to the ds3 summer school. I was trying to go over the things that we learnt and repeat them using shogun (my goal is to use hypothesis testing via shogun for my own research). However, I cannot get the same results in shogun and in the code developed at DS3.

Something as simple as computing the MMD metric outputs different results using shogun w.r.t. using the MMD implementation of the summer school. I can show that with the following example:

import numpy as np
from tqdm import tqdm_notebook as tqdm
from scipy.spatial.distance import squareform, pdist, cdist
import shogun as sg

data = np.load("blobs.npz")
X_learn = data["X"]
Y_learn = data["Y"]

def two_sample_permutation_test(test_statistic, X, Y, num_permutations, prog_bar=True):
    assert X.ndim == Y.ndim
    
    statistics = np.zeros(num_permutations)
    
    range_ = range(num_permutations)
    if prog_bar:
        range_ = tqdm(range_)
    for i in range_:
        # concatenate samples
        if X.ndim == 1:
            Z = np.hstack((X,Y))
        elif X.ndim == 2:
            Z = np.vstack((X,Y))
            
        # IMPLEMENT: permute samples and compute test statistic
        perm_inds = np.random.permutation(len(Z))
        Z = Z[perm_inds]
        X_ = Z[:len(X)]
        Y_ = Z[len(X):]
        my_test_statistic = test_statistic(X_, Y_)
        statistics[i] = my_test_statistic
    return statistics

def quadratic_time_mmd(X,Y,kernel):
    assert X.ndim == Y.ndim == 2
    K_XX = kernel(X,X)
    K_XY = kernel(X,Y)
    K_YY = kernel(Y,Y)
       
    n = len(K_XX)
    m = len(K_YY)
    
    # IMPLEMENT: unbiased MMD statistic (could also use biased, doesn't matter if we use permutation tests)
    np.fill_diagonal(K_XX, 0)
    np.fill_diagonal(K_YY, 0)
    mmd = np.sum(K_XX) / (n*(n-1))  + np.sum(K_YY) / (m*(m-1))  - 2*np.sum(K_XY)/(n*m)
    return mmd


def gauss_kernel(X, Y=None, sigma=1.0):
    """
    Computes the standard Gaussian kernel k(x,y)=exp(- ||x-y||**2 / (2 * sigma**2))

    X - 2d array, samples on left hand side
    Y - 2d array, samples on right hand side, can be None in which case they are replaced by X
    
    returns: kernel matrix
    """

    # IMPLEMENT: compute squared distances and kernel matrix
    sq_dists = sq_distances(X,Y)
    K = np.exp(-sq_dists / (2 * sigma**2))
    return K

def sq_distances(X,Y=None):
    assert(X.ndim==2)

    # IMPLEMENT: compute pairwise distance matrix. Don't use explicit loops, but the above scipy functions
    # if X=Y, use more efficient pdist call which exploits symmetry
    if Y is None:
        sq_dists = squareform(pdist(X, 'sqeuclidean'))
    else:
        assert(Y.ndim==2)
        assert(X.shape[1]==Y.shape[1])
        sq_dists = cdist(X, Y, 'sqeuclidean')

    return sq_dists


log_sigma=-2
num_permutations=200

# Shogun implementation
feat_p=sg.RealFeatures(X_learn.T.astype(np.float64))
feat_q=sg.RealFeatures(Y_learn.T.astype(np.float64))
kernel=sg.GaussianKernel(2 * (10**log_sigma)**2)

mmd=sg.QuadraticTimeMMD(feat_p,feat_q)
mmd.set_kernel(kernel)

mmd.set_statistic_type(sg.ST_UNBIASED_FULL)
statistic=mmd.compute_statistic()

mmd.set_null_approximation_method(sg.NAM_PERMUTATION)
mmd.set_num_null_samples(num_permutations)


# DS3 summer school implementation
my_kernel = lambda X,Y : gauss_kernel(X,Y,sigma=10**log_sigma)
my_mmd = lambda X,Y : quadratic_time_mmd(X,Y, my_kernel)
my_statistic = my_mmd(X_learn, Y_learn)
statistics = two_sample_permutation_test(my_mmd, X_learn, Y_learn, num_permutations, prog_bar=False)
p_value = np.mean(my_statistic <= np.sort(statistics))

print(statistic)
print(my_statistic)
print(mmd.compute_p_value(statistic))
print(p_value)

Any guess on what might be happening? The MMD implementation should be the same in the toolbox and in the code. Might it be the kernel?

Edit: I have tried to use linear kernels and I still get different MMD values. I used the linear_kernel method from the summer school).

[jeslago]

I have also compared the Shogun implementation with the implementation that you provided in the summer school to be used with Tensorflow. I still get inconsistent results:

import tensorflow as tf
import mmd
data = np.load("blobs.npz")
X = data["X"]
Y = data["Y"]

sigma_median=5.079463545341823

with tf.Session() as sess:
    init = tf.global_variables_initializer()
    sess.run(init)
    print("MMD and ratio:", sess.run(mmd.rbf_mmd2_and_ratio(X,Y, sigma=sigma_median)))

feat_p=sg.RealFeatures(X.T.astype(np.float64))
feat_q=sg.RealFeatures(Y.T.astype(np.float64))
kernel=sg.GaussianKernel(2 * (sigma_median)**2)

mmd=sg.QuadraticTimeMMD(feat_p,feat_q)
mmd.set_kernel(kernel)

mmd.set_statistic_type(sg.ST_BIASED_FULL)
print(mmd.compute_statistic())

[karlnapf]

I think Shogun returns a scaled version of the MMD statistic (Note that this does not matter if you do the test as either the permutation test is used, or the spectral test is scaled appropriately).
A better way to compare would be to compute the p-value of the statistic and compare that ... and it should be similar (more similar for more permutations)
@lambday what does the compute_statistic call compute exactly? The docs do not reflect it http://shogun-toolbox.org/api/latest/classshogun_1_1CQuadraticTimeMMD.html

[jeslago]

Thanks for the answer. Indeed, the p-values are actually more similar, so I guess the MMD statistic from shogun is somehow scaled. It would be good though to know the scaling of it!

[lambday]

hey @jeslago, the MMD statistic in Shogun (the one that the method compute_statistic returns) is \frac{n_x\timex n_y}{n_x + n_y}\times MMD^2 estimate, where n_x is the number of samples from P and n_y is the number of samples from Q. You can check the API doc for CMMD class (link below) to see for more details.

HTH :)

[CMMD class api doc] http://www.shogun-toolbox.org/api/latest/classshogun_1_1CMMD.html