FastMMD: Ensemble of Circular Discrepancy for Efficient Two-Sample Test

Introduction

The two-sample test is one of the most fundamental tests in statistics and has a wide range of applications. It uses samples drawn from two distributions to test whether to accept or reject the null hypothesis that they are the same or different. This task, however, is very difficult and challenging in practice since the underneath distribution information are generally unknown apriori (Bickel, 1969; Friedman & Rafsky, 1979; Hall & Tajvidi, 2002; Biau & Gyofi, 2005). The maximum mean discrepancy (MMD) is the latest test statistic designed for this task by measuring the discrepancy of two distributions by embedding them in a reproducing kernel Hilbert space (Gretton et al., 2012a). The MMD has been attracting much attention in recent two-sample test research due to its solid theoretical fundament (Smola et al., 2008; Sriperumbudur et al., 2010, 2011; Sejdinovic et al., 2013) and successful applications including biological data test, data integration and attribute matching (Gretton et al., 2009), outlier detection, data classifiability (Sriperumbudur et al., 2009), domain adaption, etc. By generalizing the MMD to kernel families as the supremum of MMDs on a class of kernels, it has also been effectively used for some basic machine learning problems such as kernel selection (Sriperumbudur et al., 2009).

Albeit its various applications, the exact MMD needs $O(N^{2}d)$ computational cost, where $N$ and $d$ denote the size and dimension of samples, respectively, to calculate the kernel values between all pairs from the assessed two-sample sets. This quadratic computational complexity greatly hampers its further application to large-scale practical problems. How to speedup the computation of MMD has thus become a hot issue in statistics and machine learning in recent years.

There are mainly two approaches proposed for this problem by approximating MMD on a subsampling set of all sample pairs. The first is MMD-linear, which is the extremely simplified MMD calculation by only using possibly fewest interactions of sample pairs (Gretton et al., 2012a). While this strategy significantly accelerates the MMD calculation to $O(Nd)$, it also brings very high variance due to its evident loss of sample pair information. To better leverage the computation cost and calculation accuracy of MMD, B-test is recently proposed (Zaremba et al., 2013). The main idea is to split two-sample sets into corresponding subsets, construct block correspondence between them, and then compute the exact MMD inner each block while omit the inter-block pair information. By changing the block size, it can vary smoothly from MMD-linear with linear complexity to exact MMD with quadratic complexity. In practice, the block size is generally set as a modest value $\sqrt{N}$ by experience. Thus the time complexity of B-test is $O(N^{3/2}d)$, correspondingly.

Actually, as the coming of the big data era, it has become a hot trend to enhancing the efficiency of kernel-based learning methods, such as support vector machines and Gaussian process, throughout machine learning, computer vision and data mining. Many efforts have been made to speedup the establishment of the kernel information and accelerate the implementation of kernel techniques (Smola & Schölkopf, 2000; Williams & Seeger, 2000; Fine & K.Scheinberg, 2001). Two of the representative developments are Random Kitchen Sinks (Rahimi & Recht, 2007, 2008) and Fastfood (Le et al., 2013), which can significantly speed up the computation for a large range of kernel functions by mapping data into a relatively low-dimensional randomized feature space. These developments inspire us for this MMD-acceleration research topic, which constitutes an important branch along this line of research.

The main difficulty and challenge of MMD calculation lie in the fact that it needs to compute the kernel values between all sample pairs of two sets. MMD-linear and B-test attain this task by only utilizing a subsampling pair subset from all. Such simplification, however, also decreases the accuracy of MMD calculation due to their neglectness of the entire sample pair information. To this aim, this paper proposes a new efficient MMD calculation strategy, which can implement this task in a more efficient and accurate way. In summary, this paper mainly contains the following four-fold contributions:

1. Through employing Bochner’s theorem and Fastfood technique (Le et al., 2013) for kernels that are spherically invariant, we reduce the MMD computation cost to $O(LN\log d)$, which has lower time complexity than current MMD approximation methods and facilitates MMD’s application in large-scale data. Moreover, our method is easy to be sequentially computed and parallelized.

2. The proposed method utilizes the interacted kernel values between all pairs from two sample sets to calculate MMD, which naturally leads to very accurate MMD result. Our experimental results substantiate that our method is with similar accuracy as exact MMD, and with significantly smaller variance than MMD-linear and B-test.

3. We have theoretically proved the uniform convergence of our method in both unbiased and biased cases. Comparatively, both MMD-linear and B-test are only feasible in unbiased cases.

4. We provide a geometrical explanation of our method in calculating MMD with shift-invariant kernels. Under this viewpoint, it is potentially useful for arousing more extensive metrics for two-sample test.

The code of our FastMMD method is available at http://gr.xjtu.edu.cn/web/dymeng/2.

Efficient MMD for Shift-Invariant Kernels

Firstly we give a brief review of MMD (Gretton et al., 2012a) and introduce some important properties of shift-invariant kernel (Rahimi & Recht, 2007). Then we propose an efficient MMD approximation method.

Substituting the empirical estimates $\mu(X_{1}):=\frac{1}{|I_{1}|}\sum_{i\in I_{1}}\phi(\mathbf{x}_{i})$ and $\mu(X_{2}):=\frac{1}{|I_{2}|}\sum_{i\in I_{2}}\phi(\mathbf{x}_{i})$ of the feature space means based on respective samples, an empirical biased estimate of MMD can then be obtained as:

where $a_{i}=\frac{1}{|I_{1}|}$ if $i\in I_{1}$, and $a_{i}=\frac{-1}{|I_{2}|}$ if $i\in I_{2}$. We can see that the time complexity of such MMD estimate is $O(N^{2}d)$. We will investigate how to accelerate its computation to $O(LN\log d)$, especially for shift-invariant kernels.

2 Efficient Approximation of MMD

The following classical theorem from harmonic analysis provides the main fundament underlying our approximation method (Genton, 2001).

Theorem 1 (Bochner) Every bounded continuous positive definite function is Fourier transform of a non-negative finite Borel measure. This means that for any bounded shift-invariant kernel $K(\mathbf{x},\mathbf{y})$, there exists a non-negative finite Borel measure $\mu$ satisfying

We assume that the discussed positive definite kernel is real valued. According to Bochner’s theorem, we have that if a shift-invariant kernel $K(\mathbf{x},\mathbf{y})$ is positive definite, there exists a proper scaled probability measure $p(\boldsymbol{\omega})$ satisfying

Since both the probability measure $p(\boldsymbol{\omega})$ and the kernel $K(\Delta)$ are real, the integrand $e^{j\boldsymbol{\omega}^{\prime}(\mathbf{x}-\mathbf{y})}$ can be replaced by $\cos(\boldsymbol{\omega}^{\prime}\mathbf{x}-\boldsymbol{\omega}^{\prime}\mathbf{y})$ in the above equation. Taking the Gaussian kernel $K(\mathbf{x},\mathbf{y};\sigma)=e^{-\frac{\|\mathbf{x}-\mathbf{y}\|^{2}}{2{\sigma}^{2}}}$ as an example, we can rewrite it as $K(\boldsymbol{\Delta};\sigma)=e^{-\frac{\|\boldsymbol{\Delta}\|^{2}}{2{\sigma}^{2}}}$, where $\boldsymbol{\Delta}=\mathbf{x}-\mathbf{y}$. Its Fourier transform is $p(\boldsymbol{\omega};\sigma)=(2\pi)^{-\frac{d}{2}}e^{-\frac{\sigma^{2}\|\boldsymbol{\omega}\|^{2}}{2}}$. By proper scaling, $\boldsymbol{\omega}$ can be viewed as a multivariate Gaussian distribution $\boldsymbol{\omega}\sim\mathcal{N}(\mathbf{0},\frac{1}{\sigma^{2}}\mathbf{I})$, where $\mathbf{I}$ is the $d\times d$ identity matrix.

Claim 2 For a shift-invariant kernel $K(\mathbf{x},\mathbf{y})=\left\langle\phi(\mathbf{x}),\phi(\mathbf{y})\right\rangle$, suppose $p(\boldsymbol{\omega})$ is its corresponding normalized measure in Bochner’s theorem, then

Claim 2 can be easily proved by substituting Eqn. (2) into kernel $K(\mathbf{x}_{i},\mathbf{x}_{j})$. A very interesting thing is that we can fortunately calculate $\sum_{i=1}^{N}\sum_{j=1}^{N}a_{i}a_{j}\cos(\boldsymbol{\omega}^{\prime}\mathbf{x}_{i}-\boldsymbol{\omega}^{\prime}\mathbf{x}_{j})$ in linear time by applying the Harmonic Addition Theorem (Nahin, 1995), as shown in Fig. 1. First, this expression can be viewed as a squared amplitude of combined sinusoids. Suppose a linear combination of $N$ sinusoids is $\sum_{i=1}^{N}a_{i}\sin(x-\boldsymbol{\omega}^{\prime}\mathbf{x}_{i})=A\sin(x-\theta)$, then its amplitude has a closed form $A^{2}=\sum_{i=1}^{N}\sum_{j=1}^{N}a_{i}a_{j}\cos(\boldsymbol{\omega}^{\prime}\mathbf{x}_{i}-\boldsymbol{\omega}^{\prime}\mathbf{x}_{j})$, see Fig. 1(a). Second, the amplitude of sinusoids with the same frequency can be calculated in a sequential way in linear time, see Fig. 1(b). By combining the above two observations, we can calculate the expression with linear time complexity. If we set $a_{i}$ in Eqn. (3) as that in empirical MMD, it turns out to be the biased estimate of MMD. As a result, Claim 2 finely implies a novel methodology to efficiently approximate MMD.

By the same spirit of Random Fourier features (Rahimi & Recht, 2007), we draw i.i.d samples from distribution $p(\boldsymbol{\omega})$, and then take the average of $A^{2}(\boldsymbol{\omega};X_{1},X_{2})$ to approximate the empirical estimate of squared MMD. The Quasi-Monte Carlo sampling techniques can also be used to generate sample sequence which has lower discrepancy than i.i.d sampling (Yang et al., 2014). Formally, we sample $L$ samples $\{\boldsymbol{\omega}_{k}\}_{k=1}^{L}$ from the distribution $p(\boldsymbol{\omega})$, and then use them to approximate $\text{\rm MMD}_{\text{b}}$:

where $A(\boldsymbol{\omega}_{k};X_{1},X_{2})$ is the amplitude of linear combination of sinusoids $\frac{1}{|I_{1}|}\sum_{i\in I_{1}}\sin(x-\boldsymbol{\omega}^{\prime}\mathbf{x}_{i})-\frac{1}{|I_{2}|}\sum_{i\in I_{2}}\sin(x-\boldsymbol{\omega}^{\prime}\mathbf{x}_{i})$. We call $L$ the number of basis functions (Neal, 1994; Le et al., 2013). The efficient calculation of $\overline{\text{\rm MMD}}_{\text{u}}$ for unbiased MMD is similar since we can also rewrite unbiased estimate of MMD as the form in Eqn. (3). The proof is provided in Appendix A.1. The aforementioned procedure for approximating MMD is described in Algorithm 1. We also provide an equivalent implementation of Algorithm 1 in Appendix A.2.

For kernels that are spherically invariant, the Fastfood technique can be employed to further speedup $\boldsymbol{\omega}$ sampling in step 2 and $\boldsymbol{\omega}^{\prime}\mathbf{x}_{u}$ calculation in steps 4 - 5 of Algorithm 1 (Le et al., 2013). This can bring further efficiency gain for MMD calculation from $O(LNd)$ to $O(LN\log d)$. In the rest of this paper, we call our original algorithm as FastMMD-Fourier and its variant using Fastfood as FastMMD-Fastfood, respectively.

3 Computational Complexity

As aforementioned, $N$ is the number of samples, $d$ is the dimension number of samples, and $L$ is the number of basis functions for approximating $p(\boldsymbol{\omega})$. Given a sampling of $\boldsymbol{\omega}$, the time complexity for calculating $A(\boldsymbol{\omega};X_{1},X_{2})$ is $O(Nd)$. The overall computational complexity of the entire FastMMD-Fourier is thus $O(LNd)$. For FastMMD-Fastfood, the computation speed is further enhanced to $O(LN\log d)$. Usually, the basis number $L$ can preset as a fixed number and are thus independent of the sample scale. As compared with the complexities of the previous MMD methods, such as $O(N^{2}d)$ for exact MMD (Gretton et al., 2012a), and $O(N^{3/2}d)$ for B-test (Zaremba et al., 2013), the proposed FastMMD methods evidently get a speed gain. Furthermore, instead of only utilizing a subsampling pair subset from all by the current MMD approximation methods, FastMMD takes into consideration all the interacted information between sample pairs. Our methods are thus expected to be more accurate.

Another interesting thing is that the calculation of $A(\boldsymbol{\omega};X_{1},X_{2})$ can be computed in a sequential way, and thus our method can be naturally implemented in stream computations. Also our method is easy to be parallelized. This further implies the potential usefulness of the proposed FastMMD methods in real large-scaled applications.

4 Approximation Guarantees

In this section, we will prove the approximation ability of the proposed FastMMD methods.

Theorem 3 (Uniform Convergence of FastMMD-Fourier) Let $\mathcal{M}$ be a compact subset of ${\bbbr}^{d}$ with diameter $\operatorname{\text{\rm diam}}(\mathcal{M})$. Then, for the biased estimate of MMD in Algorithm 1, we have:

Theorem 4 (Uniform Convergence of FastMMD-Fastfood) If we use Fastfood method (Le et al., 2013) to calculate $\{\boldsymbol{\omega}^{\prime}_{k}\mathbf{x}_{i}\}_{k=1}^{L}$ in Algorithm 1, suppose the kernel is Gaussian kernel with bandwidth $\sigma$ and $\widehat{\text{\rm MMD}}_{\text{\rm b}}$ is the biased estimate of MMD that arises from a $d\times d$ block of Fastfood In Fastfood method, the kernel expansions are calculated by constructing several $d\times d$ blocks. For the asymptotic analysis of error bound, we can treat $L$ as $d$ by padding the data with zeros if $d<L$. This allow us to consider only one block., then we have:

This bound also holds for the approximation of unbiased MMD.

From Theorems 3 and 4, we can see that the approximation of FastMMD is unbiased. The proof is provided in Appendix A.2, A.3 and A.4.

5 Tests Based on the Asymptotic Distribution of the Unbiased Statistic

The principle for constructing our FastMMD approximation complies with the spirit of random kitchen and sinks (Rahimi & Recht, 2007) and Fastfood (Le et al., 2013), both using the kernel expansion to approximate shift-invariant kernels. Specifically, we use the linear combination of Dirac-delta functions to approximate continuous function $p(\boldsymbol{\omega})$, which is uniquely determined by kernel $K(\cdot,\cdot)$. It means that we implicitly introduce a new kernel which is an approximation for the original kernel in MMD calculation.

Based on the aforementioned analysis, the uniform convergence bounds and the asymptotic distribution for general kernels in Gretton et al. (2012a), including the Theorem and Corollary 7 - 13, still hold for FastMMD. Given the asymptotic distribution of the unbiased statistic $\text{\rm MMD}_{\text{\rm u}}^{2}$, the goal of the two-sample test is to determine whether the empirical test statistic $\text{\rm MMD}_{\text{\rm u}}^{2}$ is so large as to be outside the $1-\alpha$ quantile of the null distribution. In Gretton et al. (2012a), two approaches based on asymptotic distributions are proposed. One method is using the bootstrap (Arcones & Giné, 1992) on aggregated data, and the other is approximating the null distribution by fitting Pearson curves to its first three moments. Both the two methods can be incorporated into FastMMD.

Ensemble of Circular Discrepancy

In the previous section, we proposed an efficient approximation for MMD. In this section, we give a geometric explanation for our methods by using random projection (Blum, 2006) on a circle and circular discrepancy. This explanation is expected to help us more insightfully understand such approximation and inspire more extensive metrics for the two-sample test other than MMD.

If $\boldsymbol{\omega}$ in Eqn. (2) is fixed, we can see that the positions of points projected on a unit circle sufficiently determine the kernel. In other words, the random variables $\boldsymbol{\omega}^{\prime}\mathbf{x}$ and $\boldsymbol{\omega}^{\prime}\mathbf{y}$ can be wrapped around the circumference of a unit circle without changing the value of kernel function. We first investigate the circular distribution under fixed $\boldsymbol{\omega}$ in the following, and later will discuss the cases when $\boldsymbol{\omega}$ is sampled from a multivariate distribution.

Given a fixed $\boldsymbol{\omega}$, we wrap two classes of samples on a unit circle separately. The probability density functions (PDFs) of the wrapped two random variables $X_{1}(\boldsymbol{\omega})$ and $X_{2}(\boldsymbol{\omega})$ can be mathematically expressed as:

where $\operatorname{mod(\cdot,\cdot)}$ is the modular arithmetic, and $\delta(\cdot)$ is the Dirac delta function. Distributions $p_{1}(x;\boldsymbol{\omega})$ and $p_{2}(x;\boldsymbol{\omega})$ are zero when $x\in(-\infty,0)\cup[2\pi,+\infty)$. Such distributions are called circular distributions or polar distributions (Fisher, 1993).

2 Circular Discrepancy

We now define a metric for measuring the discrepancy between $X_{1}(\boldsymbol{\omega})$ and $X_{2}(\boldsymbol{\omega})$. Later we will show that this definition is closely related to MMD.

Definition 2 Given two independent circular distributions $X_{1}\sim P_{1}$ and $X_{2}\sim P_{2}$, we define the circular discrepancy as:

where $Y\sim Q$ is also a circular distribution.

In this definition, we choose sine function as the measure for assessing the distance between two circular distributions.

Claim 5 The circular discrepancy as defined in (7) is equal to:

If $p_{1}(x)$ and $p_{2}(x)$ are probability mass functions (linear combination of Dirac delta functions), let $p_{1}(x)-p_{2}(x)=\sum_{i=1}^{N}a_{i}\delta(x-x_{i})$, and then the circular discrepancy is equal to

From Eqn. (8) in Claim 5, it can be seen that the optimal distribution of random variable $Q$ in the definition of circular discrepancy is a Dirac delta function. For the integral in Eqn. (8), if we change $y$ to $\operatorname{mod}(y+\pi,2\pi)$, the value of this expression will change its sign. Based on this observation, it is clear that there are two distributions of $Q$ that can maximize and minimize the object function, respectively. The difference of their non-zero positions is $\pi$. These two distributions construct a “optimal decision diameter” for the projected samples on a unit circle.

We then give a geometric explanation for problem (7). Note that the sine function is a distance measure with sign for two points on a circle (suppose positive angles are counterclockwise), so the definition of circular discrepancy aims to find a diameter that possibly largely separates the projected samples of two different classes. An example is shown in Fig. 2, where the orientation of the diameter corresponds to the non-zero elements of distribution $Q$ in the definition of circular discrepancy. We can see that this diameter maximizes the mean margin of two sample classes.

To have a large ensemble of circular discrepancy, the spread parameter $\sigma$ in Gaussian kernel should be neither too small nor two large. If $\sigma\rightarrow 0$, its Fourier transform $p(\boldsymbol{\omega})$ tend to be the uniform distribution. Intuition suggests that if a continuous distribution is spread out over a large region of the line then the corresponding wrapped distribution will be almost uniform on the circle. If $\sigma\rightarrow+\infty$, its Fourier transform $p(\boldsymbol{\omega})$ tend to be Dirac delta function, and all the projected points on the circle would be near zero. In both of the two cases, the ensemble of circular discrepancy is small. This observation is consistent with the asymptotic behavior of Gaussian kernel SVM (Keerthi & Lin, 2003).

3 Ensemble of Circular Discrepancy

We have discussed the circular discrepancy for a given random projection $\boldsymbol{\omega}$ on a unit circle. For shift-invariant kernels, $\boldsymbol{\omega}$ is not a fixed value, but randomly sampled from a distribution, see Eqn. (2). We use sampling method for ensemble of circular discrepancy under different random projections. This ensemble turns out to be an efficient approximation of empirical estimate of MMD, so Algorithm 1 can be explained as ensemble of circular discrepancy.

Definition 3 Suppose $p(\boldsymbol{\omega})$ is the normalized distribution in Eqn. (2) for kernel $K$; $X_{1}$ and $X_{2}$ are two circular distributions depending on $\boldsymbol{\omega}$ according to Eqn. (5)(6). The ensemble of circular discrepancy for $X_{1}$ and $X_{2}$ with shift-invariant kernel $K$ is defined as:

Claim 6 For a shift-invariant kernel $K(\mathbf{x},\mathbf{y})=\left\langle\phi(\mathbf{x}),\phi(\mathbf{y})\right\rangle$, $K(\mathbf{0})=1$; denote $\mathcal{H}$ as the associated Hilbert space with kernel $K$; $p(\boldsymbol{\omega})$ as the normalized distribution in Eqn. (2); $X_{1}$ and $X_{2}$ as two circular distributions depending on $\boldsymbol{\omega}$ defined by Eqn. (5)(6). Then the ensemble of circular discrepancy is

Proof: By substituting Eqn. (9) into Eqn. (10) and utilizing Bochner’s theorem, we can obtain

where $a_{i}=\frac{1}{|I_{1}|}$ if $i\in I_{1}$, and $a_{i}=\frac{-1}{|I_{2}|}$ if $i\in I_{2}$. ∎

Fig. 3 demonstrates the relationship between MMD, circular discrepancy and our approximation. The blue and red contours are two distributions, and $p(\boldsymbol{\omega})$ is the distribution determined by Fourier transform of the kernel. We generate i.i.d. samples $\boldsymbol{\omega}$ from distribution $p(w)$. For each generated $p(\boldsymbol{\omega})$, we project samples on a unit circle and calculate the circular discrepancy. The ensemble of discrepancy then corresponds to the MMD. We can see that the circular discrepancy constructs classifiers implicitly.

It is interesting that if we use other similarity measurements such as Mallows Distance (earth mover’s distance) (Werman et al., 1986), other than the sine function utilized in Definition 2, more extensive metrics for two-sample test can be naturally obtained. Furthermore, it should be noted that in Definition 3, we use $L_{2}$-norm average for ensemble. If we use other norms, we can also generalize more measures other than MMD for this task. All these extensions are hopeful directions for future investigation.

Experimental Results

We begin by investigating how well our methods can approximate the exact MMD as the sampling number $L$ increases. Following previous work on kernel hypothesis testing (Gretton et al., 2012b; Zaremba et al., 2013), our synthetic distribution is designed as $5\times 5$ grids of 2D Gaussian blobs. We specify two distributions, $P$ and $Q$. For distribution $P$ each Gaussian has identity covariance matrix, while for distribution $Q$ the covariance is non-spherical with a ratio $\epsilon$ of large to small covariance eigenvalues. Samples drawn from $P$ and $Q$ are presented in Fig. 4. The ratio $\epsilon$ is set as $4$, and the sample number for each distribution is set as $1000$.

We used a Gaussian kernel with $\sigma=1$, which approximately matches the scale of the variance of each Gaussian in mixture $P$. The MMD approximation results are shown in Fig. 5. We use the relative difference between exact (biased and unbiased) MMD and the approximation to quantify the error. The absolute difference also exhibits similar behavior and is thus not shown due to space limit. The results are presented as averages from $1000$ trials. As can be seen, as $L$ increases, both FastMMD-Fourier and FastMMD-Fastfood converge quickly to exact MMD in both biased and unbiased cases. Their performance are indistinguishable. It can be observed that for both methods, the good approximation can be obtained even from a modest number of basis.

2 Accuracy Test for MMD with Kernel Family

In some situations, MMD with a kernel family is preferred (Sriperumbudur et al., 2009). Here we present an experiment that illuminates the superiority of our FastMMD methods on accuracy of MMD with kernel family. The synthesis data are generated as follows. All samples are constrained into a two-dimensional rectangle: $-5\leq x_{1},x_{2}\leq 5$. The points, which are located within a circular ring in between by $x_{1}^{2}+x_{2}^{2}=1$ and $x_{1}^{2}+x_{2}^{2}=16$ are labeled as $+1$, while other points are labeled as $-1$. We generate $200$ samples for each distribution randomly as the test set. Then we use these samples to calculate MMD with a kernel family. Here the kernel family is composed of multivariate isotropic Gaussians with bandwidth varying between $0.1$ and $100$, with a multiplicative step-size of $10^{1/5}$.

We compare our method with exact MMD, MMD-linear (Gretton et al., 2012a), and B-test (Zaremba et al., 2013). Note that the latter two methods are only valid for unbiased MMD. In our methods, the number of basis function $L$ is set as $1024$. The block size in B-test is set to the default choice, i.e, the square of sample size $\sqrt{N}$. For our method, we repeat $1000$ times and use the curves and error bars to represent means and standard deviations of MMD, respectively. Since both MMD-linear and B-test depend on the permutation of data samples, we make $1000$ permutations of the samples. From Fig. 6, it can be seen that the means of all these methods are consistent with the true values. Also it can be seen that our FastMMD-Fourier and FastMMD-Fastfood have similar accuracy, and their deviations are much smaller than those of MMD-linear and B-test.

For some applications, we need to find the kernel that has the maximal MMD. Since our methods have lower variance, they incline to find the correct $\sigma$ with higher probability than MMD-linear and B-test.

3 Efficiency Test on Synthetic Data

In order to evaluate the efficiency of our methods, we generate samples uniformly from $[0,0.95]^{d}$ and $[0.95,1]^{d}$. The efficiency of different methods is shown in Fig. 7. In Fig. 7(a), the number of samples is varied from $10^{3}$ to $10^{5}$, and the data dimension $d$ is set as $16$. In Fig. 7(b), the number of samples is set as $10^{4}$, and the dimension $d$ is varied from $8$ to $1024$.

All competing methods are implemented in Matlab except that we use Spiral WHT Packagehttp://www.spiral.net/ to perform the Fast Walsh-Hadamard transform in Fastfood. The comparison methods include exact MMD, MMD-linear and B-testhttps://github.com/wojzaremba/btest. We run all the codes on a PC with AMD Athlon X2 250 ($800$ MHz) CPU and $4$ GB RAM memory.

From Fig. 7, we can see that when $N$ or $d$ varies from small to large, our methods gradually become more efficient than exact MMD and B-test. When $d=16$ and $N=10^{5}$, FastMMD-Fourier and FastMMD-Fastfood are $2000/5$x and $4000/10$x faster than exact MMD / B-test, respectively. As for MMD-linear, since it is the extreme simplified subsampling version of MMD calculation, it always runs very fast. However, as the sample size or dimension increases, the computation times of our methods still depict a more slowly increase trend than that of MMD-linear. This empirically confirms the efficiency of the proposed FaseMMD methods.

4 Efficiency Test on Large-Scale Real Data

In order to validate the efficiency of FastMMD on real world data, we also perform experiments on PASCAL VOC 2007 http://pascallin.ecs.soton.ac.uk/challenges/VOC/voc2007/, which is popular in computer vision community. The dataset consists of $9963$ images, and there are 20 object categories in this dataset with some images containing multiple objects. In our experiments, we choose $4015$ images which contain persons as one sample set, and the remaining $5948$ images as another sample set. We use the VLFeat toolboxhttp://www.vlfeat.org/ to extract features. For each image, the feature is constructed by bag-of-words representation and spatial pyramid. The codebook size for bag-of-words is $1024$, and the spatial pyramid contains $1\times 1$, $2\times 2$ and $4\times 4$ levels. Thus the feature dimension for each image is $1024\times(1+2\times 2+4\times 4)=21504$. Finally, the feature vectors are normalized by $L_{1}$ norm.

We use MMD to measure the discrepancy between these two image set. In this MMD calculation, the sample number is $N=9963$, and data dimension is $d=21504$. The number of basis function $L$ is fixed as $1024$. The bandwidth $\sigma$ is set as $10^{-2.2}=0.0631$. We perform this experiment on a PC with Intel Core i7-3770 $@3.4$ GHz CPU and $32$ GB RAM memory.

The results of efficiency and accuracy are shown in Tab. 1. We compare our method with exact MMD, MMD-linear and B-test. Except for exact MMD, we repeat $10$ times for each method, and report the mean and standard deviation of MMD, execution time, and relative error. We can see that our FastMMD methods have smaller approximation error and smaller deviation, and meanwhile it has three orders of speedup compared to exact MMD.

5 Type II Error

Given MMD, several strategies can be employed to calculate the test threshold (Gretton et al., 2009, 2012a; Chwialkowski et al., 2014). The bootstrap strategy (Arcones & Giné, 1992) is utilized in our experiments since it can be easily integrated into our FastMMD method. Also the bootstrap is preferred for large-scale datasets since it costs $O(N^{2}d)$, faster than most other methods for this task such as Person, with cost $O(N^{3}d)$ (Gretton et al., 2012a).

The data used for this experiment is generated from Gaussian blob distributions as described in Section 4.1. The sample size is set as $1000$ for two distributions. The bandwidth is selected by maximizing MMD. The selected bandwidth is $\sigma=1$, and it approximately matches the scale of variance of each Gaussian blob in distribution $P$. We find that all of biased/unbiased FastMMD methods have similar good performance, and we thus only demonstrate the result of FastMMD-Fourier for biased MMD.

The level $\alpha$ for Type I error is set as $0.05$, and the number of bootstrap shuffles is $1000$. The Type II error is shown in Fig. 8(a). We can see that the Type II error drops quickly when increasing number of basis. It demonstrates empirically that increasing number of basis can decrease the Type II error.

We also compare our method with Bootstrap (Bootstrap approach with exact MMD), Pearson (moment matching to Pearson curves), spectrum (Gram matrix eigenspectrum) and Gamma (two-parameter Gamma approximation), where the former two approaches are presented in Gretton et al. (2012a) and the latter two are introduced in Gretton et al. (2012b). The basis number utilized in FastMMD is fixed as $256$. The type II errors of all competing methods with respect to $\epsilon$ are shown in Fig. 8(b). The execution times for different approaches are $88.7$s (Bootstrap), $37.9$s (Pearson), $33.2$s (spectrum), $0.73$s (Gamma), and $4.5$ (our FastMMD+Bootstrap). We can see that the proposed FastMMD method can achieve comparable results with other methods, while it is evidently more efficient than Bootstrap, Pearson and spectrum. This actually substantiates the efficiency of FastMMD since our method and Bootstrap are the same except with different strategies for calculating MMD. Our method is less efficient than Gamma method because bootstrap methods need to calculate MMD for many times (it is determined by the number of bootstrap shuffles which is $1000$ in this experiment). How to design effective statistic based on FastMMD and optimal threshold remains a challenging issue for future investigation.

Conclusions

In this paper, we propose a method, called FastMMD, for efficiently calculating the maximum mean discrepancy (MMD) for shift-invariant kernels. Taking advantage of Fourier transform of shift-invariant kernels, we get a linear time complexity approximation method. We prove the theoretical convergence of the proposed method in both unbiased and biased cases, and further present a geometric explanation for it. This explanation on one side delivers new insight for intrinsic MMD mechanism, and on the other side is hopeful to inspire more extensive new metrics for two-sample test, which will be investigated in our future research. Future work also includes finding the significance threshold using more efficient and effective strategies other than bootstrap.

Acknowledgement

We thank Alex Smola for valuable suggestions. We also thank anonymous reviewers for their constructive comments, which helped us to improve the manuscript. This research was supported by the National Basic Research Program of China (973 Program) under Grant No. 2013CB329404, the China NSFC project under contract 61373114, 11131006, and the Civil Aviation Administration of China jointly funded project No. U1233110.

Appendix

The unbiased estimate of MMD is (See Lemma 6 in Gretton et al. (2012a))

Denote $S_{1}=\frac{1}{m^{2}}\sum_{i\in I_{1}}\sum_{j\in I_{1}}K(\mathbf{x}_{i},\mathbf{x}_{j})$, $S_{2}=\frac{1}{n^{2}}\sum_{i\in I_{2}}\sum_{j\in I_{2}}K(\mathbf{x}_{i},\mathbf{x}_{j})$. If the kernel $K$ is shift-invariant, let $k_{0}=K(\mathbf{0},\mathbf{0})=K(\mathbf{x}_{i},\mathbf{x}_{i})$ for any $\mathbf{x}_{i}$, and then we have

Denote $\mu_{1}=\frac{1}{m}\sum_{i\in I_{1}}\phi(\mathbf{x}_{i})$ and $\mu_{2}=\frac{1}{n}\sum_{i\in I_{2}}\phi(\mathbf{x}_{i})$. The biased and unbiased estimate of MMD can be reformulated as:

A.2 Preliminary for the Proof of Theorem 3 and Theorem 4

The following two claims and an equivalent form of Algorithm 1 provide basic tools to analyze the MMD approximation.

Claim 7 (Uniform Convergence for Linear Combination of Approximation) Let $\mathcal{M}$ be a compact subset of ${\bbbr}^{d}$, $f(\cdot;\boldsymbol{\omega})$ be a function with parameter $\boldsymbol{\omega}$. If $\boldsymbol{\omega}$ is drawn from a distribution, and the difference of two functions is bounded as

where $C(\epsilon)$ is function about $\epsilon$. Then

Proof: For any fixed $\boldsymbol{\omega}$, we have

Based on Eqn. (14), we know that with probability more than $1-C(\epsilon)$, it hold that $\sup_{\mathbf{x}_{i},\mathbf{y}_{j}\in\mathcal{M}}\left|f(\mathbf{x}_{i},\mathbf{y}_{j};\boldsymbol{\omega})-g(\mathbf{x}_{i},\mathbf{y}_{j})\right|<\epsilon$. So with probability more than $1-C(\epsilon)$, it holds that $\sup_{\{\mathbf{x}_{i}\},\{\mathbf{y}_{j}\}\subset\mathcal{M}}\left|\sum_{(i,j)}a_{i}a_{j}f(\mathbf{x}_{i},\mathbf{y}_{j};\boldsymbol{\omega})-\sum_{(i,j)}a_{i}a_{j}g(\mathbf{x}_{i},\mathbf{y}_{j})\right|$ $<\sum_{(i,j)}\left|a_{i}a_{j}\right|\cdot\epsilon$. The proof is then completed. ∎

and then make a variable substitution $x\leftarrow x^{\prime}/b$.

Define a function $f(x)=2\log x+\sqrt{a^{2}+x^{-2}}+\log(\sqrt{2}a)-a$. Next we prove three properties of $f(x)$.

(i) $\lim_{x\rightarrow 0^{+}}f(x)=+\infty$, because $x^{-2}$ dominates the other when $x\rightarrow 0^{+}$. Also it is obvious that $\lim_{x\rightarrow+\infty}f(x)=+\infty$.

(ii) $f^{\prime}(x)=\frac{1}{x}\left(2-\frac{1}{x^{2}\sqrt{a^{2}+x^{-2}}}\right)$. Let $f^{\prime}(x^{*})=0$, and then we obtain the only solution $x^{*}=\sqrt{\frac{-1+\sqrt{1+a^{2}}}{2a^{2}}}$. Since $f(x)$ is continuous and differentiable, $f(0^{+})=+\infty$, and $f(+\infty)=+\infty$, we have that $f(x^{*})$ is the minimum of this function.

Denote $g(a)=\log(-1+\sqrt{1+a^{2}})-\log(\sqrt{2}a)+(1+\sqrt{1+a^{2}})-a$. We can verify that $g^{\prime}(a)=\sqrt{1+a^{-2}}-1>0$ and $g(1)>0$. So $f(x^{*})=g(a)>0$. ∎

We provide an equivalent form of FastMMD. It is easy to verify that

where $\mathbf{v}(\mathbf{x}_{i})=\left[\cos(\boldsymbol{\omega}^{\prime}_{1}\mathbf{x}_{i}),\cdots,\cos(\boldsymbol{\omega}^{\prime}_{L}\mathbf{x}_{i}),\sin(\boldsymbol{\omega}^{\prime}_{1}\mathbf{x}_{i}),\cdots,\sin(\boldsymbol{\omega}^{\prime}_{L}\mathbf{x}_{i})\right]$.

The time complexity of this method is also linear. Note that the expression on the right hand is the square length for weighted addition of random Fourier features (Rahimi & Recht, 2007), which means this method uses the random Fourier features to calculate MMD. The procedure for approximating MMD is described in Algorithm 2. Note that this algorithm is equivalent to Algorithm 1 as described in the maintext. Algorithm 1 provides us a geometric explanation of MMD, and Algorithm 2 is convenient for the following uniform convergence analysis.

A.3 Proof of Theorem 3

From Algorithm 2, it is easy to see that:

where $a_{i}=\frac{1}{|I_{1}|}$ if $i\in I_{1}$ and $a_{i}=\frac{1}{|I_{2}|}$ if $i\in I_{2}$, $\boldsymbol{\omega}=\{\boldsymbol{\omega}_{i}\}_{i=1}^{L}$. $K(\cdot)$ and $\mathbf{z}(\mathbf{x}_{i};\boldsymbol{\omega})$ are the same as that in Algorithm 2.

According to Eqn. (5) in Gretton et al. (2012a):

Based on the above Eqn. (18)(19), we have:

According to Claim 1 in Random Features (Rahimi & Recht, 2007):

Denote the right hand of the inequality (21) as $C(\epsilon)$. Based on Claim 7, the Eqn. (20) is bounded by $C\left(\frac{\epsilon}{\sum_{i=1}^{N}\sum_{j=1}^{N}|a_{i}a_{j}|}\right)=C\left(\frac{\epsilon}{4}\right)$, then we can obtain the uniform convergence for biased estimate of MMD.

Next we prove the case for unbiased estimate of MMD. According to Algorithm 2 and by virtue of certain algebraically equivalent transformation, we obtain:

According to Eqn. (3) in Gretton et al. (2012a), we know that

Based on Eqn. (22)(23) and Claim 7, $\text{\rm Pr}\left[\sup_{\mathbf{x}_{1},\cdots\mathbf{x}_{N}\in\mathcal{M}}\left|\overline{\text{\rm MMD}}_{\text{b}}^{2}-\text{\rm MMD}_{\text{b}}^{2}\right|\geq\epsilon\right]$ is bounded by $C\left(\frac{\epsilon}{\sum_{(i,j)}|a_{i}a_{j}|}\right)=C\left(\frac{\epsilon}{4}\right)$, and then we obtain the uniform convergence for unbiased estimate of MMD. ∎

A.4 Proof of Theorem 4

$\mathcal{M}_{\boldsymbol{\Delta}}$ is compact and with diameter at most twice $\operatorname{\text{\rm diam}}(\mathcal{M})$. And then we can find an $\epsilon$-net that covers $\mathcal{M}_{\boldsymbol{\Delta}}$ using at most $T=\left(4\operatorname{\text{\rm diam}}(\mathcal{M})/r\right)^{d}$ balls of radius $r$ (Cucker & Smale, 2001). Let $\left\{\boldsymbol{\Delta}_{k}\right\}_{k=1}^{T}$ denote the centers of these balls. We have $\left|f(\boldsymbol{\Delta})\right|<\epsilon$ for all $\boldsymbol{\Delta}\in\mathcal{M}_{\boldsymbol{\Delta}}$ if the following two conditions hold for all $k$: (i) $\left|f(\boldsymbol{\Delta}_{k})\right|<\epsilon/2$; (ii) $\left|f(\mathbf{x}_{i},\mathbf{x}_{j})\right|<\epsilon/2$, if $\mathbf{x}_{i}-\mathbf{x}_{j}$ belongs to the ball $k$ of the $\epsilon$-net. Next we bound the probability of these two events:

(i) The union bound followed by Hoeffding’s inequality applied to the anchors in the $\epsilon$-net gives

(ii) If we use Fastfood for FastMMD in Algorithm 2, and suppose the estimate of kernel arises from a $d\times d$ block of Fastfood, then according to Theorem 6 in Fastfood literature (Le et al., 2013) we have:

Since $d$ usually is large in Fastfood method, we suppose that $(\log d)/2\geq 1$. Considering Claim 8 and the fact $\|\mathbf{x}_{i}-\mathbf{x}_{j}\|\leq 2r$, this inequality can be further reformulated as

Combining (24) and (25) gives a bound in terms of the free variable $r$:

This has the form $1-\kappa_{1}r^{-d}-\kappa_{2}r^{2}$. Setting $r=\left(\frac{\kappa_{1}}{\kappa_{2}}\right)^{\frac{1}{d+2}}$ turns this to $1-\kappa_{2}^{\frac{d}{d+2}}\kappa_{1}^{\frac{2}{d+2}}$, and assuming that $\frac{\log d\operatorname{\text{\rm diam}}(\mathcal{M})}{d\sigma^{2}\epsilon^{2}}\geq 1$ and $\operatorname{\text{\rm diam}}(\mathcal{M})\geq 1$, we obtain that

Denote the right hand of the inequality (26) as $C(\epsilon)$. Similar to the proof in Theorem 3, we have

Then we complete the proof of the uniform convergence for biased estimate of MMD. The convergence for unbiased estimate of MMD can be proved in the similar way. ∎

A.5 Proof of Claim 5

Let $p_{1}(x)$, $p_{2}(x)$ and $q(y)$ be the PDFs for random variables $X_{1}$, $X_{2}$ and $Y$, respectively. We then have

If $p_{1}(x)-p_{2}(x)=\sum_{i=1}^{N}a_{i}\delta(x-x_{i})$, then

The second equation holds because the sum of sinusoids with the same frequency is also a sinusoid with that frequency. According to the trigonometric identity, it holds that

where $\operatorname{atan2}$ is the four-quadrant arctangent function. The supremum of this problem is $A$ when $y=\theta$. ∎

References