Title: POST SELECTION INFERENCE WITH INCOMPLETE MAXIMUM MEAN DISCREPANCY ESTIMATOR

Abstract: Measuring divergence between two distributions is essential in machine learning and statistics and has various applications including binary classification, change point detection, and two-sample test. Furthermore, in the era of big data, designing divergence measure that is interpretable and can handle high-dimensional and complex data becomes extremely important. In this paper, we propose a post selection inference (PSI) framework for divergence measure, which can select a set of statistically significant features that discriminate two distributions. Specifically, we employ an additive variant of maximum mean discrepancy (MMD) for features and introduce a general hypothesis test for PSI. A novel MMD estimator using the incomplete U-statistics, which has an asymptotically normal distribution (under mild assumptions) and gives high detection power in PSI, is also proposed and analyzed theoretically. Through synthetic and real-world feature selection experiments, we show that the proposed framework can successfully detect statistically significant features. Last, we propose a sample selection framework for analyzing different members in the Generative Adversarial Networks (GANs) family.

1. INTRODUCTION:

Computing the divergence between two probability distributions is fundamental to machine learning and has many important applications such as binary classification (Friedman et al., 2001) , change point detection (Yamada et al., 2013a; Liu et al., 2013) , two-sample test (Gretton et al., 2012; Yamada et al., 2013b) , and generative models such as generative adversarial networks (GANs) (Goodfellow et al., 2014; Li et al., 2015b; Nowozin et al., 2016) , to name a few. Recently, interpreting the difference between distributions has become an important task in applied machine learning (Mueller & Jaakkola, 2015; Jitkrittum et al., 2016) since it can facilitate scientific discovery. For instance, in biomedical binary classification tasks, it is common to analyze which variables or features are different between two different distributions (classes).

The simplest approach to measure divergence between two probability densities would be parametric methods. For example, t-test can be used if the distributions to be compared are known and well-defined (Anderson, 2001) . However, in many real-world problems, the property of distribution is not known a priori, and therefore model assumptions are likely to be violated. In contrast, non-parametric methods can be applied to any distributions without prior assumptions. The maximum mean discrepancy (MMD) (Gretton et al., 2012) is an example of non-parametric discrepancy measures and is defined as the difference between the mean embeddings of two distributions in a reproducing kernel Hilbert space (RKHS). For RKHSs associated with characteristic kernels (Sriperumbudur et al., 2011) , the difference in mean embeddings defines a proper metric, and thus MMD can capture potentially non-linear difference between distributions (i.e., higher order moments). However, since MMD considers the entire d dimensional vector, it is hard to interpret how individual features contribute to the discrepancy.

To deal with the interpretability issue, divergence measure with feature selection has been actively studied (Yamada et al., 2013a; Mueller & Jaakkola, 2015; Jitkrittum et al., 2016) . For instance, in Mueller & Jaakkola (2015) , the Wasserstein divergence is employed as a divergence measure and 1 regularizer is used for feature selection. However, these approaches focus on detecting a set of features that discriminate two distributions. But for scientific discovery applications such as biomarker discovery, it might be preferable to test the significance of each selected feature (e.g., one biomarker). A naive approach would be to select features from one dataset and then test the selected features using the same data. However, in such case, selection bias is included and thus the false positive rate cannot be controlled. Therefore, it is crucial in hypothesis testing that the selection event should be taken into account to correct the bias. To the best of our knowledge, there is not yet an existing framework that tests the significance of selected features that distinguish between two different distributions on the same dataset.

In this paper, we propose mmdInf, a post selection inference (PSI) algorithm for distribution comparison, which finds a set of statistically significant features that discriminate between two distributions. mmdInf enjoys several compelling properties. First, it is a non-parametric method based on kernels, and thus it can detect features that distinguish various types of distributions. Second, the proposed framework is general and can be applied to not only feature selection but also other distribution comparison problems, such as dataset comparison. In mmdInf, we employ the recently developed PSI algorithm (Lee et al., 2016) and use MMD as a divergence measure. However, the standard empirical squared MMD estimator has a degenerated null distribution, which violates the requirement of a proper Normal distribution for the current PSI algorithm. To address this issue, we apply the block estimate (Zaremba et al., 2013) and the linear estimate (Gretton et al., 2012) of MMD. Furthermore, we propose a new empirical estimate of MMD based on the incomplete U-statistics (incomplete MMD estimator) (Blom, 1976; Janson, 1984; Lee, 1990) and show that it asymptotically follows the normal distribution, and has much greater detection power (compared to the existing estimators of MMD) in PSI. Finally, we propose a framework to analyze different members in the GAN (Goodfellow et al., 2014) family based on mmdInf. We elucidate the theoretical properties of the incomplete U-statistics estimate of MMD and show that mmdInf can successfully detect significant features through feature selection experiments.

The contributions of our paper are summarized as follows:

(1) We propose a non-parametric PSI algorithm mmdInf for distribution comparison.

(2) We propose the incomplete MMD estimator and investigate its theoretical properties (see Figure 1 ).

(3) We propose a sample selection framework based on mmdInf that can be used for analyzing generative models.

2. RELATED WORK:

Divergence measures in general can be divided into two categories: f -divergence (Ali & Silvey, 1966) , such as the Kullback-Leibler divergence (Cover & Thomas, 2012) or the α-divergence (Rényi et al., 1961; Póczos & Schneider, 2011) , and integral probability metric (Müller, 1997) , including the total variation distance (Shorack & Shorack, 2000) and the Wasserstein distance in its dual (Villani, 2008) . In this work we consider the maximum mean discrepancy (MMD) (Gretton et al., 2012) , which falls into the second category. Although these divergence measures can be used for testing the discrepancy between p(x) and q(x), it is hard to test the significance of one of k selected features from the entire d features, where the setup is useful for scientific discovery tasks such as biomarker discovery. Recently, a MMD-based change point detection algorithm, which can compute the p-value of the largest score, has been proposed (Li et al., 2015a) . However, this method can only test the largest score (i.e., one-feature). Thus, it is not clear whether the approach can be extended to feature selection in general.

A novel testing framework for the post selection inference (PSI) has been recently proposed (Lee et al., 2016) , in which statistical inference after feature selection with Lasso (Tibshirani, 1996) is investigated. This work shows that statistical inference conditioned on the selection event can be done for linear regression models with Gaussian noise if the selection event can be written as a set of linear constraints. However, the PSI algorithm needs to assume a Gaussian response, which is a relatively strong assumption, and consequently it cannot be directly applied for non-Gaussian output problems such as classification. To deal with this issue, a kernel based PSI algorithm (hsicInf) for independence test has been proposed (Yamada et al., 2018) , in which the Hilbert-Schmidt Independence Criterion (HSIC) (Gretton et al., 2005) is employed to measure the independence between input and its output, and thus significance test can be performed on non-Gaussian data. In this paper, we propose an alternative kernel based inference algorithm for distribution comparison called mmdInf, which can be used for both feature selection and binary classification.

In Sections 3.4 and 5.3, we manifest how we analyze different members in the GANs (Goodfellow et al., 2014) family. Here, we discuss several evaluation metrics that have been proposed to compare generative models. For example, the inception scores (Salimans et al., 2016) and the mode scores (Che et al., 2016) measure the quality and diversity of the generated samples, but they were not able to detect overfitting and mode dropping / collapsing for generated samples. The Frechet inception distance (FID) (Heusel et al., 2017) defines a score using the first two moments of the real and generated distributions, whereas the classifier two-sample tests (Lopez-Paz & Oquab, 2016) considers the classification accuracy of a binary classifier as a statistic for two-sample test. Although the above metrics are reasonable in terms of discriminability, robustness, and efficiency, the distances between samples are required to be computed in a suitable feature space. We can also use the kernel density estimation (KDE); or more recently, Wu et al. (2016) proposed to apply the annealed importance sampling (AIS) to estimate the likelihood of the decoder-based generative models. Nevertheless, these approaches need the access to the generative model for computing the likelihood, which are less favorable comparing to the model agnostic approaches which rely only on a finite generated sample set. On the other hand, the maximum mean discrepancy (MMD) (Gretton et al., 2012) is preferred against its competitors (Sutherland et al., 2016; Huang et al., 2018) . Therefore, we propose the mmdInf based GANs analysis framework.

3. PROPOSED METHOD (M M DIN F ):

In this section, we introduce a PSI algorithm with MMD (Gretton et al., 2012) .

3.1. POST SELECTION INFERENCE (PSI):

Suppose we are given independent and identically distributed (i.i.d.) samples

X = {x i } m i=1 ∈ R d×m from a d-dimensional distribution p and i.i.d. samples Y = {y j } n j=1 ∈ R d×n from another d-dimensional distribution q.

Our goal is to find k < d features that differentiate between X and Y and test whether each selected feature is statistically significant.

Let S be a set of selected features, we consider the following hypothesis test:

• H 0 : d s=1 η s D(X (s) , Y (s) ) = 0 | S was selected, • H 1 : d s=1 η s D(X (s) , Y (s) ) = 0 | S was selected,

where D(X (s) , Y (s) ) is the estimated discrepancy measure for the selected feature s and η = [η 1 , . . . , η d ] ∈ R d is an arbitrary pre-defined parameter. To test the s-th feature, we can set η as a unit vector whose s-th position is 1 and zero otherwise.

We employ the post-selection inference (PSI) framework (Lee et al., 2016) to test the hypothesis. Thanks to the cumulative distribution of Theorem 1, we can compute the p-value of d s=1 η s D(X (s) , Y (s) ) under feature selection.

Theorem 1 (Lee et al., 2016) Suppose that z ∼ N (µ, Σ), and the feature selection event can be expressed as Az ≤ b for some matrix A and vector b, then for any given feature represented by η ∈ R n we have

F [V -(A,b),V + (A,b)] η µ,η Ση (η z) | Az ≤ b ∼ Unif(0, 1),

where

F [a,b] µ,σ 2 (x)

is the cumulative distribution function (CDF) of a truncated normal distribution with mean µ and variance σ 2 truncated at [a,b] . Given that α = A Ση η Ση , the lower and upper truncation points are given by

V -(A, b) = max j:αj <0 b j -(Az) j α j + η z, V + (A, b) = min j:αj >0 b j -(Az) j α j + η z.

Feature Selection with Discrepancy Measure: Assume we have an estimate of a discrepancy measure for each feature: z

= [ D(X (1) , Y (1) ), . . . , D(X (d) , Y (d) )] ∈ R d ∼ N (µ, Σ)

, where D(•, •) has large positive value when two distribution are different. We select the top-k features with largest discrepancy scores. We denote the index set of the selected k features by S, and that of the unselected k = d -k features by S. This feature selection event can be characterized by

D(X (s) , Y (s) ) ≥ D(X ( ) , Y ( ) )

, for all (s, ) ∈ S × S.

Note that we have in total k • k constraints.

The selection event can be rewritten as A s, z ≤ 0, for all (s, ) ∈ S × S, where 

A s, = [0 • • • 0 -1 s 0 • • • 0 1 0 • • • 0] and A s, ∈ R d is a row vector of A ∈ R (k• k)×d .

3.2. MAXIMUM MEAN DISCREPANCY (MMD):

We employ the maximum mean discrepancy (MMD) (Gretton et al., 2012) as divergence measure.

Let F be the unit ball in a reproducing kernel Hilbert space (RKHS) and K(x, x ) the corresponding positive definite kernel with

E x,x ∼p [ K(x, x)] < ∞, E x,x ∼q [ K(x, x)] < ∞, the squared population MMD is defined as MMD 2 [F, p, q] = E x,x [K(x, x )] -2E x,y [K(x, y)] + E y,y [K(y, y )],

where E x,y denotes the expectation over independent random variables x with distribution p and y with distribution q. It has been shown that if the kernel function et al., 2012) . In the following, we introduce two existing MMD estimators and then propose the incomplete U-statistics MMD estimator.

K(x, x ) is characteristic, then MMD[F, X, Y ] = 0 if and only if p = q (Gretton

(Complete) U-statistics estimator (Gretton et al., 2012) : We use the Gaussian kernel:

K(x, x ) = exp - x-x 2 2 2σ 2 x

, where σ x > 0 is the Gaussian width. When m = n, the complete U-statistics of MMD is defined as

MMD 2 u [F, X, Y ] = 1 m(m -1) i =j h(u i , u j ),

where h(u, u ) = K(x, x )+K(y, y )-K(x, y )-K(x , y) is the U-statistics kernel for MMD and u = [x y ] ∈ R 2d . However, since the complete U-statistics estimator of MMD is degenerated under p = q and does not follow normal distribution, this estimator cannot be used in PSI.

Algorithm 1 mmdInf (Feature Selection) H 0 :

EQUATION

sel )] ∈ R d . 2: Select k features (i.e., k = |S|) and compute A using z and Σ using Xc and Yc. Set b = 0.

3: Compute p-value as p = 1 -F [V -(A,b),V + (A,b)] 0,η Ση (η z).

Block estimator (Zaremba et al., 2013) : The block estimate of MMD is given by

MMD 2 b [F, X, Y ] = B 1 n n/B1 i=1 MMD 2 u [F, X i , Y i ],

where

X = [X 1 , . . . , X m/B1 ], Y = [Y 1 , . . . , Y n/B2 ],

and

X i ∈ R d×B1 and Y i ∈ R d×B2 are i-th partitioned data.

Here, we assume that the number of blocks m/B 1 = n/B 2 is an integer and

B = B 1 = B 2 .

A commonly used heuristic for the block size is to set B = √ n. This estimator asymptotically follows the normal distribution when B 1 and B 2 are finite and m and n go to infinity. The block estimator can be used for PSI, but the variance and normality depends on the partition of X and Y . Specifically, when the total number of samples is small, then a small block size would result in high variance, whereas larger block size tends to result in non-Gaussian response.

Incomplete U-statistics estimator: The described problems of the block-estimator motivated us to design a new MMD estimator that is normally distributed and has smaller variance. We therefore propose an MMD estimator based on the incomplete U-statistics (Blom, 1976; Janson, 1984; Lee, 1990) . The incomplete U-statistics estimator of MMD is given by

MMD 2 inc [F, X, Y ] = 1 (i,j)∈D h(u i , u j ),

where D is a subset of S n = {(i, j)} i =j , and = |D|. D can be fixed design or random design.

In particular, if we design D as D = {(1, 2), (3, 4), . . . , (n -1, n)} then the incomplete U-statistic corresponds to the linear-time MMD estimator (Gretton et al., 2012) .

3.3. COVARIANCE MATRIX ESTIMATION:

For PSI, we need to estimate the covariance matrix Σ from data.

Block estimator To estimate Σ, we first compute H b ∈ R d×(n/B) whose elements are the meansubtracted MMDs. Then, we regard H b as the data matrix and use a standard covariance estimator for estimating Σ b . Note that for small n, the number of blocks n/B tends to be smaller than the dimension d, and thus, the estimation accuracy of covariance matrix is low and affects the detection accuracy. To handle this issue, we employ the POET algorithm (Fan et al., 2013) .

Incomplete U-statistics estimator:

To estimate Σ, we first compute H inc ∈ R d× whose elements are the mean-subtracted U-stat kernels. Then, we regard H inc as the data matrix and can use a standard covariance estimator to compute Σ inc .

The estimation performance of the covariance estimation for the block estimate heavily depends on the block size B. That is, if B is large, we need to estimate Σ b from a small number of samples. On the other hand, in the incomplete U-statistics estimator, we can set = rn n/B, where r > 0 is a constant (See Section 4). Thus, in practice, the estimated covariance matrix tends to be more accurate than the block estimator counterpart. Figure 6 in the supplementary material shows the MSEs between the true and the estimated covariance matrices of the incomplete estimate and the block estimate, respectively. As can be seen, the error of the incomplete estimate is smaller than that of the block estimate, which is an advantage of the incomplete MMD over the block MMD in PSI.

3.4. ADDITIONAL APPLICATIONS (GANS ANALYSIS):

The proposed PSI framework is general and can be used for not only feature selection but also sample selection. In generative modeling, the generated distribution should match the data distribution; in other words, the discrepancy between the generated data and real data should be small. We can therefore apply the selective inference algorithm and use the significance value to evaluate the generation quality. In this paper we apply mmdInf to compare the performance of GANs. We first select the model whose generated samples has the smallest MMD score with the real data and then perform the hypothesis test.

Let x:

(s) i ∈ R p be a feature vector generated by s-th GAN model with random seed i and y j ∈ R p is a feature vector of an original image. Image features can be extracted by pre-trained Resnet (He et al., 2016) or auto-encoders. The hypothesis test can be written as (s) , Y ] = 0 | s generates samples closest to the real distribution, (s) , Y ] = 0| s generates samples closest to the real distribution.

• H 0 : MMD 2 inc [F, X

• H 1 : MMD 2 inc [F, X

Since we want to test the best generator that minimizes the discrepancy between generated and true samples (e.g., low MMD score), this sample selection event can be characterized by

MMD 2 inc [F, X (s) , Y ] ≤ MMD 2 inc [F, X ( ) , Y ], (s, ) ∈ S × S.

4. THEORETICAL ANALYSIS OF INCOMPLETE MMD:

We investigate the theoretical properties of the incomplete MMD estimator under the random design with replacement. For simplicity, we denote MMD

2 inc [F, X, Y ] = MMD 2 inc , MMD 2 [F, p, q] = MMD 2 , MMD 2 inc [F, X (s) , Y (s) ] = MMD 2,(s)

inc , and MMD 2 [F, p (s) , q (s) ] = MMD 2,(s) , respectively. We introduce the conditional expectations for the U-statistic kernel h(u, u ) as h 1 = E u [h(u, u )] and h 2 = h(u, u ). See the supplementary material for proof.

Theorem 2 Let c be the smallest integer such that h c = MMD 2 , and let n and tend to infinity such that γ = lim n, →∞ n -c , 0 ≤ γ ≤ ∞. For sampling with replacement, we have

     1 2 (MMD 2 inc -MMD 2 ) D -→ N (0, σ 2 ), if γ = 0. n c 2 (MMD 2 inc -MMD 2 ) D -→ V, if γ = ∞. 1 2 (MMD 2 inc -MMD 2 ) D -→ γ 1 2 V + T, if 0 < γ < ∞,

where V is the random variable of the limit distribution of n c 2 (MMD 2 u -MMD 2 ), T is the random variable of N (0, σ 2 ), σ 2 = V ar(h(u, u )), and T and V are independent.

Corollary 3 Assume lim n, →∞ n -2 = 0 and 0 < γ = lim n, →∞ n -1 < ∞. For sampling with replacement, incomplete U-statistics estimator of MMD is asymptotically normally distributed as

1 2 MMD 2 inc D -→ N (0, σ 2 ), if p = q. 1 2 (MMD 2 inc -MMD 2 ) D -→ N (0, σ 2 + γσ 2 u ), if p = q. where σ 2 = V ar(h(u, u )) and σ 2 u = 4(E u [(E u [h(u, u )]] -E u,u [h(u, u )]) 2 ).

Corollary 4 Assume lim n, →∞ n -1 = 0. For sampling with replacement, the incomplete Ustatistics estimator of MMD is asymptotically normally distributed as

1 2 (MMD 2 inc -MMD 2 ) D -→ N (0, σ 2 ).

Thus, in practice, by setting n 2 , the incomplete estimator is asymptotically normal and therefore can be applied in PSI. More specifically, we can set = rn n 2 , where r is a small constant. In practice, we found that r = 10 works well in general.

Figure 2 shows the empirical distribution under p = q and p = q for the complete estimator, the block estimator and the incomplete estimator. As can be observed, the empirical distribution of the incomplete estimator is normal for small sampling parameter r, and becomes similar to its complete counterpart if r is large; this is supported by Theorem 2 (γ = ∞).

Finally, the following theorem assure the joint normality of the z vector in Theorem 1. Theorem 5 Suppose 0 ≤ γ < ∞. Then, (d) converges to multivariate normal distribution.

1/2 MMD 2,(1) inc , . . . , MMD 2,(d) inc -MMD 2,(1) , . . . , MMD 2,

5. EXPERIMENTS:

We compared mmdInf with a naive testing baseline (mmd), which first selects features using MMD and estimates corresponding p-values with the same data of feature selection without adjustment for the selection event. For mmdInf, we used the three MMD estimators: the linear-time MMD (Gretton et al., 2012) , the block MMD (Zaremba et al., 2013) , and the incomplete MMD. We used 1/2 of data to calculate the covariance matrix of MMD and the rest to perform feature selection and inference. We fixed the number of selected features (prior to PSI) k to 30. In PSI the significance of each of the 30 selected features (from ranking MMD) is computed and features with p-value lower than the significance level α = 0.05 are selected as statistically significant features.

For block MMD, in each experiment we set the candidate of block size as B = {10, 20, 50} . For incomplete MMD, in each experiment the ratio between number of pairs (i, j) sampled to compute incomplete MMD score and sample size is fixed at r = n ∈ {0.5, 5, 10}. We reported the true positive rate (TPR) k k * where k is the number of true features selected by mmdInf or mmd and k * is the total number of true features in synthetic data. We further computed the false positive rate (FPR) k k-k * where k is the number of non-true features reported as positives. We ran all experiments over 200 times, and reported the average TPR and FPR.

5.1. SYNTHETIC EXPERIMENTS (PSI):

The number of features d is fixed to 50, and for each feature, data is randomly generated following a Gaussian distribution with set mean and variance. 10 out of the 50 features are set to be significantly different by shifting the distribution of one class away from the other (mean or variance). More specifically, we generate the synthetic data as

(a) Mean shift x ∼ N (0 50 , I 50 ), y ∼ N (µ, I), µ = [1 10 0 40 ] ∈ R 50 , (b) Variance shift x ∼ N (0 50 , I 50 ), y ∼ N (0, Σ), Σ = diag([1.51 10 1 40 ] ),

where N (µ, Σ) is a multivariate normal distribution with mean µ ∈ R d and covariance Σ ∈ R d×d , 1 p ∈ R d is a vector whose elements are all one, 0 d ∈ R d is a vector whose elements are all zero, and diag(a) ∈ R d×d is a diagonal matrix whose diagonal elements are a ∈ R d .

Published as a conference paper at ICLR 2019 As can be clearly seen, PSI successfully controls FPR with significance level α = 0.05 for all the three estimators, whereas the naive approach tends to have higher FPRs. Figures 3 shows the TPRs of the synthetic data. In both cases, the TPR of incomplete MMD converges to 1 significantly faster than the the other two empirical estimates.

B E G A N C R A M E R D C G A N D F M D R A G A N M I N I B A T C H S T D G A N W G A N

5.2. REAL-WORLD DATA (BENCHMARK):

We compared the proposed algorithm by using real-world datasets. Since it is difficult to decide what is a "true feature" in real-world data, we choose a few datasets for binary classification with small amount of features, and regard all the original features as true. We then concatenated random features to the true features (the total number of features d = 100). Table 1 shows TPRs and FPRs of mmdInf with different MMD estimators. It can be observed that the incomplete estimator significantly outperforms the other empirical estimates. Note that a higher TPR can be achieved with higher r, while the FPR is still controlled at 0.05 with r = 5. Additional comparison of different empirical estimates on real dataset (with varying sample size) can be found in Figures 8.

5.3. GANS ANALYSIS (SAMPLE SELECTION):

We also applied mmdInf for evaluating the generation quality of GANs. We trained BEGAN (Berthelot et al., 2017) , DCGAN (Radford et al., 2015) , STDGAN (Miyato et al., 2017) , Cramer GAN (Bellemare et al., 2017) , DFM (Warde-Farley & Bengio, 2016) , DRAGAN (Kodali et al., 2017) , Minibatch Discrimination GAN (Salimans et al., 2016) , and WGAN-GP (Gulrajani et al., 2017) , generated 5000 images (using Chainer GAN package1 with CIFAR10 datasets), and extracted 512 dimensional features by pre-trained Resnet18 (He et al., 2016) . For the true image sets, we subsampled 5000 images from CIFAR10 datasets and computed the 512 dimensional features using the same Resnet18. We then tested the difference between the generated images and the real images using mmdInf on extracted features (see Sec. 3.4).

We found that for all the members in the GAN family, the null hypothesis was rejected, i.e., the generated distribution and the real distribution are different. This result is consistent with the findings in (Sutherland et al., 2016) , which demonstrate that optimized MMD has perfect discriminative power in GANs evaluation. As sanity check, we evaluated mmdInf by constructing an "oracle" generative model that generates real images from CIFAR10. Next, we randomly selected 5000 images (a disjoint set from the oracle generative images) from CIFAR10 in each trial, and set the sampling ratio to r = 5. Figure 4 (a) shows the distribution of p-values computed by our algorithm. We can see that the p-values are distributed uniformly in the tests for the "oracle" generative model, which matches the theoretical result in Theorem 1. Thus the algorithm is able to detect the distribution difference and control the false positive rate. In other words, if the generated samples do not follow the original distribution, we can safely reject the null hypothesis with a given significance level α. shows the estimated MMD scores of each model. Based on the results, we could tell that DFM was the best model and DCGAN was the second best model to approximate the true distribution. However, the difference between various members is not obvious. Developing a validation pipeline based on mmdInf for GANs analysis would be one interesting line of future work.

6. CONCLUSION:

In this paper, we proposed a novel statistical testing framework mmdInf, which can find a set of statistically significant features that can discriminate two distributions. Through synthetic and realworld experiments, we demonstrated that mmdInf can successfully find important features and/or datasets. We also proposed a method for sample selection based on mmdInf and applied it in the evaluation of generative models.

Thus, in practice, by setting n 2 , the incomplete estimator is asymptotically normal and therefore can be applied in PSI. More specifically, we can set = rn n 2 , where r is a small constant. In practice, we found that r = 10 works well in general.

Theorem 5 Suppose 0 ≤ γ < ∞. Then, 2,(d) converges to multivariate normal distribution.

1/2 MMD 2,(1) inc , . . . , MMD 2,(d) inc -MMD 2,(1) , . . . , MMD

Proof: We use the fact that convergence in distribution of multivariate random variables results in convergence in distribution of univariate random variables. From Theorem29.4 in Billingsley (2008)(Cramér-Wold device), it is sufficient to prove that for any η s) ] converges to normal distribution derived in Corollary 3 when 0 ≤ γ < ∞, and from the continuous mapping theorem for g(x) = η x, we obtain the desired result.

EQUATION

ILLUSTRATIVE EXPERIMENTS:

Figure 5 shows the empirical distribution under p = q and p = q for the complete estimator, the block estimator and the incomplete estimator. As can be observed, the empirical distribution of the incomplete estimator is normal for small sampling parameter r, and becomes similar to its complete counterpart if r is large; this is supported by Theorem 2 (γ = ∞). Moreover, compared to the block estimator, the incomplete estimator tends to have a better trade-off between variance and normality.

Figure 6 shows the MSEs between the true and the estimated covariance matrices of the incomplete estimate and the block estimate, respectively. As can be seen, the error of the incomplete estimate is smaller than that of the block estimate, which is an advantage of the incomplete MMD over the block MMD in PSI. Moreover, for the incomplete U-statistics estimator, we have = rn samples to estimate the covariance matrix, while we have only n/B samples for the block estimate. That is, the number of samples that can be used for covariance estimation of the incomplete U-statistics estimator is rB times larger than that of the block estimator. For example, for n = 1000, we can use 5000 samples for the incomplete U-statistics, while we have only 50 samples for the block estimate. This is problematic for high-dimensional case (e.g., d = 100), and the estimated covariance matrix is not full rank. To alleviate this issue, we employ the POET algorithm for the block estimator.

Figure 7 (a) shows the Type II error comparison for two-sample test with one dimensional Gaussian mean shift data. The Type II error is computed when the Type I error is fixed at 0.05, and the incomplete MMD outperforms other estimators. Figure 7 (b) compares the computational time of the empirical estimates, and for small r the computational time of incomplete MMD is much less than that of the block MMD. Overall, the incomplete MMD has favorable properties in practice.

Figure 8 show the true positive rate (TPR) and the false positive rate (FPR) of the wine (white) dataset experiment. Here, we change the number of samples as 200, 400, . . . , 4000. As can be seen, both PSI algorithms can control the FPR with 2000 samples. Moreover, it is clear that the incomplete estimator can get better performance than the block estimator in both TPR and FPR with smaller number of samples.

EFFECTIVENESS OF PSI:

Here, we compared the PSI and non-PSI counterpart on benchmark data. As clearly seen, if we do not use the PSI, we cannot control FPR, while the proposed algorithm can successfully control FPR values. 

:

https://github.com/pfnet-research/chainer-gan-lib

Figures and Tables:

Figure 1: Comparison between Unbiased, Block, Linear, and Incomplete MMD estimation.

Under such construction, Az ≤ b can be satisfied by setting b = 0.

(a) MMD u . (b) MMD inc , r = 100. (c) MMD inc , r = 10. (d) MMD inc , r = 1.

Figure 2: Empirical distribution under p = q and p = q. (a) Complete U-statistics. (b)-(d): The incomplete MMD estimator with different sampling parameter r. For all plots, we fixed the number of samples as n = 200 and the dimensionality d = 1.

Figure 3: (a)(b): False positive rates at significant level α = 0.05 of the proposed incomplete estimator, block estimator and linear estimator with/without PSI. For incomplete MMD, we set = 10n. For block MMD, we set the block size B = √ n. The MMD without PSI computes the p-values without adjusting for the selection. (c)(d): True positive rate comparison of the following three empirical estimates for mmdInf.

Figure 4: (a) Histogram of p-values over 5000 runs. (b) Averaged incomplete MMD scores.

Figure3(a) and (b) show the FPRs of linear MMD, block MMD and incomplete MMD with or without PSI. As can be clearly seen, PSI successfully controls FPR with significance level α = 0.05 for all the three estimators, whereas the naive approach tends to have higher FPRs. Figures 3 shows the TPRs of the synthetic data. In both cases, the TPR of incomplete MMD converges to 1 significantly faster than the the other two empirical estimates.

Figure4(b) shows the estimated MMD scores of each model. Based on the results, we could tell that DFM was the best model and DCGAN was the second best model to approximate the true distribution. However, the difference between various members is not obvious. Developing a validation pipeline based on mmdInf for GANs analysis would be one interesting line of future work.

(a) MMD u . (b) MMD b , B = 100. (c) MMD b , B = 20. (d) MMD b , B = 5. (e) MMD inc , r = 100. (f) MMD inc , r = 50. (g) MMD inc , r = 10. (h) MMD inc , r = 1.

Figure 5: Empirical distribution under p = q and p = q. (a) Complete U-statistics. (b)-(d): The block MMD estimator with different block parameter B. (e)-(h): The incomplete MMD estimator with different sampling parameter r. For all plots, we fixed the number of samples as n = 200 and the dimensionality d = 1.

Post selection inference experimental results for real-world datasets. The average TPR and FPR over 200 trials are reported.

Table content: <table><tr><td/><td/><td/><td>Linear-Time</td><td/><td/><td colspan="2">Block</td><td/><td/><td/><td/><td colspan="2">Incomplete</td><td/><td/></tr><tr><td>Datasets</td><td>d</td><td>n</td><td/><td colspan="2">B = 10</td><td colspan="2">B = 20</td><td colspan="2">B = 50</td><td colspan="2">r = 0.5</td><td colspan="2">r = 5</td><td colspan="2">r = 10</td></tr><tr><td/><td/><td/><td>TPR FPR</td><td>TPR</td><td>FPR</td><td>TPR</td><td>FPR</td><td>TPR</td><td>FPR</td><td>TPR</td><td>FPR</td><td>TPR</td><td>FPR</td><td>TPR</td><td>FPR</td></tr><tr><td>Diabetis</td><td>8</td><td>768</td><td>0.05 0.06</td><td>0.20</td><td>0.04</td><td>0.40</td><td>0.12</td><td>0.46</td><td>0.18</td><td>0.13</td><td>0.06</td><td>0.52</td><td>0.07</td><td>0.65</td><td>0.10</td></tr><tr><td colspan="2">Wine (White) 11</td><td>4898</td><td>0.13 0.06</td><td>0.26</td><td>0.01</td><td>0.56</td><td>0.04</td><td>0.68</td><td>0.09</td><td>0.32</td><td>0.06</td><td>0.71</td><td>0.06</td><td>0.79</td><td>0.06</td></tr><tr><td>Wine (Red)</td><td>11</td><td>1599</td><td>0.08 0.06</td><td>0.26</td><td>0.04</td><td>0.37</td><td>0.09</td><td>0.49</td><td>0.15</td><td>0.15</td><td>0.06</td><td>0.51</td><td>0.07</td><td>0.61</td><td>0.07</td></tr><tr><td>Australia</td><td>7</td><td>690</td><td>0.08 0.06</td><td>0.36</td><td>0.06</td><td>0.47</td><td>0.11</td><td>0.65</td><td>0.21</td><td>0.16</td><td>0.05</td><td>0.66</td><td>0.07</td><td>0.79</td><td>0.09</td></tr><tr><td/><td>300</td><td/><td/><td/><td/><td/><td/><td/><td/><td/><td/><td/><td/><td/><td/></tr><tr><td/><td>200</td><td/><td/><td/><td/><td/><td/><td/><td/><td/><td/><td/><td/><td/><td/></tr><tr><td/><td>100</td><td/><td/><td/><td/><td/><td/><td/><td/><td/><td/><td/><td/><td/><td/></tr><tr><td/><td>0</td><td>0</td><td>0.5</td><td/><td/><td>1</td><td/><td/><td/><td/><td/><td/><td/><td/><td/></tr><tr><td/><td/><td/><td colspan="2">p-value</td><td/><td>(a)</td><td/><td/><td/><td/><td/><td/><td/><td/><td/></tr></table>

Peer review (iclr2019):
 The authors focus on the selection problem of k statistically significant features discriminating 2 probability distributions accessible via samples. They propose a non-parametric approach under the PSI (post selection inference) umbrella using MMD (maximum mean discrepancy) as a discrepancy measure between probability distributions. The idea is to apply (asymptotically) normal MMD estimators, rephrase the top-k selection problem as a linear constraint, and reduce the problem to Lee et al., 2016. The efficiency of the approach is illustrated on toy examples and in GAN (generative adversarial network) context. The technique complements the PSI-based independence testing approach recently proposed by Yamada et al., 2018. The submission is a well-organized, clearly written, nice contribution; it can be relevant to the machine learning community. Below I enlist a few suggestions to improve the manuscript: -Section 1: The notion of characteristic kernel (kernel when MMD is metric) has not been defined, but it was referred to. *Due to the mean embeddings in RKHS, all moment information is stored.*: This sentence is somewhat vague. -Section 1: *MMD can be computed in closed form*. This is rarely the case (except for e.g. Gaussian distributions with Gaussian or polynomial kernels). I assume that the authors wanted refer to the estimation of MMD. -Section 1: *K nearest neighbor approaches (Poczos & Schneider, 2011)*. The citation to this specific estimator can go under alpha-divergences. The Wasserstein metric could also be mentioned. -Section 3.1: k is used to denote the number of selected features and also the kernel used in MMD. I suggest using different notations. -Theorem 1: *.Phi is the CDF...*. There is no .Phi in the theorem. -Section 3.2: The existence of MMD (mean embedding) requires certain assumptions: E_{x.sim p}.sqrt{k(x,x)} < .infty, E_{x.sim q}.sqrt{k(x,x)} < .infty. -Section 3.2.: block estimator: *B_1 and B_2 are finite*. *fixed*? -Section 3.2.: MMD_{inc}: i) *S_{n,k}*: k looks superfluous. ii) *l*: it has not been introduced (cardinality of D). -Section 3.3: typo: *covraiance* (2x) -Section 3.3: Fan et al. 2013: The citation can go to .citep{}. -Theorem 2: i)*c* is left undefined. ii)Comma is missing before *where*. iii).xrightarrow{d} (Theorem 2, Corollary 3-4): Given that *d* also denotes dimension in the submission, I suggest using a different notation for convergence in distribution. -At the introduction of block-MMD the block size (B) was fixed, while in the experiments (e.g. Figure 3) it is growing with the sample size (B=.sqrt{n}). The assumption on B should be clearly stated. -Section 5.1: (b) mean shift: comma is missing before *where*. -References: i) Abbreviations and names in the titles should be capitalized (such as cramer, wasserstein, hilbert-schmidt, gan, nash). ii) Scholkopf should be Sch.{"o}lkopf (in the ALT 2005 work). iii) *Exact post-selection inference, with application to the lasso*: All the authors are listed; *et al.* is not needed.