Title: CLASSIFICATION LOGIT TWO-SAMPLE TESTING BY NEURAL NETWORKS

Abstract: The recent success of generative adversarial networks and variational learning suggests training a classifier network may work well in addressing the classical two-sample problem. Network-based tests have the computational advantage that the algorithm scales to large samples. This paper proposes to use the difference of the logit of a trained neural network classifier evaluated on the two finite samples as the test statistic. Theoretically, we prove the testing power to differentiate two smooth densities given that the network is sufficiently parametrized, by comparing the learned logit function to the log ratio of the densities, the latter maximizing the population training objective. When the two densities lie on or near to lowdimensional manifolds embedded in possibly high-dimensional space, the needed network complexity is reduced to only depending on the intrinsic manifold geometry. In experiments, the method demonstrates better performance than previous network-based tests which use the classification accuracy as the test statistic, and compares favorably to certain kernel maximum mean discrepancy (MMD) tests on synthetic and hand-written digits datasets.

1. INTRODUCTION:

The two-sample test problem concerns the comparison of two unknown distributions p and q from finitely observed data samples (Lehmann & Romano, 2006 ). As a central problem in statistics, it is widely encountered in general data analysis of biomedical data, audio and imaging data, etc. (Borgwardt et al., 2006; Chwialkowski et al., 2015; Jitkrittum et al., 2016; Lopez-Paz & Oquab, 2016; Cheng et al., 2017) , and particularly, in the machine learning application of training and evaluating generative models (Li et al., 2015; 2017; Goodfellow et al., 2014; Arjovsky et al., 2017; Nowozin et al., 2016; Lloyd & Ghahramani, 2015; Sutherland et al., 2016; Chwialkowski et al., 2016; Liu et al., 2016; Jitkrittum et al., 2017) . Many existing tests are based on certain estimators of a distance or divergence between p and q. Important examples include Maximum Mean Discrepancy (MMD), especially kernel-based MMD (Anderson et al., 1994; Gretton et al., 2012a) and distance of Reproducing Kernel Hilbert Space Mean Embedding (Chwialkowski et al., 2015; Jitkrittum et al., 2016) , divergence based methods which may involve non-parametric estimation of density difference or density ratio (Sriperumbudur et al., 2009; Wornowizki & Fried, 2016; Sugiyama et al., 2013; Kanamori et al., 2012) . While these methods have been intensively studied and theoretically wellunderstood, the application is often restricted to data of small dimensionality and/or small sample size, or certain specific classes of densities p and q, due to model and computational limitations.

The powerful expressiveness of neural networks and the recent progress in optimization of deep networks suggest the natural idea of using a network for the two-sample problem, as revisited in (Lopez-Paz & Oquab, 2016) . In the training of generative adversarial networks (GAN) (Goodfellow et al., 2014; Arjovsky et al., 2017; Nowozin et al., 2016) , in each iteration a discriminative network (D-net) is trained to distinguish the model density q produced by a generative network from the data density p which is only accessible via observed samples, that is, a two-sample problem. Strictly speaking, the task of D-net is a goodness-of-fit problem as the model density q is analytically given (Chwialkowski et al., 2016; Liu et al., 2016; Jitkrittum et al., 2017) . Since batch sampling is commonly used in training GAN and other generative networks, the ability of a trained network, such as the D-net in GAN, to detect the difference between two densities from finite samples is crucial for such applications. We review these connections in more detail in Section 1.1.

The current paper studies the method of training a network for the two-sample problem, where the test statistic is the log ratio of the class probabilities averaged over samples, which can be computed once a classifier network is trained. Our contributions include: (1) We introduce a network-based two-sample test statistic based on classification logit, and the algorithm inherits the scalable computational efficiency of neural networks; (2) Theoretical guarantee of testing power is proved for smooth densities p and q in R D , and the needed network complexity is reduced to be intrinsic when p and q lie on or near to low-dimensional manifolds embedded in possibly high dimensional ambient space. (3) Numerical experiments show that the proposed test compares favorably to kernel MMD tests and earlier neural network classifier test based on classification accuracy, on both synthetic manifold data and hand-written digits datasets.

The proposed statistic belongs to a general class of f -divergence between p and q, which means that the method may extend to a broad class of classification networks as suggested by f -GAN (Nowozin et al., 2016) beyond maximizing softmax loss. Since KL divergence (corresponding to softmax loss) is a prototypical case of f -divergence, we focus on softmax classifier network in this paper.

1.1. RELATED WORKS:

Classification and two-sample testing. The relations between two-sample testing, divergence estimation and binary classification has been pointed out in earlier statistical literature Friedman (2004) ; Sriperumbudur et al. (2009); Reid & Williamson (2011) . (Ramdas et al., 2016) studied Fisher LDA classifier used for testing mean shift of Gaussian distributions. Discriminative approach has also been used to detect and correct covariant shifts (Bickel et al., 2007; 2009) . Training classifier provides an estimator of density ratio, as has been pointed out in (Menon & Ong, 2016) and in the formulation of learning generative models (Mohamed & Lakshminarayanan, 2016) . While distribution divergence estimation has been studied and used for two-sample problems (Kanamori et al., 2011; 2012; Sugiyama et al., 2013; Wornowizki & Fried, 2016) , the use of neural network as a divergence estimator for two-sample testing was less investigated. In terms of theoretical guarantee of test power, the analysis in (Lopez-Paz & Oquab, 2016) assumes a non-zero population test statistic under H 1 which is not specified, along with other approximations. Theoretical analysis of neural network two-sample testing power remains limited.

MMD and kernel MMD tests. MMD encloses a wide class of two-sample statistics such as Kolmogorov-Smirnov statistic, Wasserstein metric, and general integral probability metrics. Particularly, kernel-based MMD (Anderson et al., 1994; Gretton et al., 2012a) has been widely applied due to its non-parametric form, and recently in training moment matching networks (MMN) (Li et al., 2015; 2017) and evaluating generative models (Sutherland et al., 2016) . To optimize kernel parameters, (Gretton et al., 2012a) considers the selection of kernel bandwidth from data, (Jitkrittum et al., 2016) studies the optimization of reference locations in the Mean Embedding test, (Gretton et al., 2012b) optimizes the kernel via a convex combination of multiple kernels, (Cheng et al., 2017) choses anisotropic kernels. The combination of neural network feature learning and kernel MMD has been studied in (Li et al., 2017) , where the training is typically more costly than that of a classifier network. Compared to kernel methods, neural networks are algorithmically more efficient and scalable, and the current paper investigates if it also has advantage in testing power.

Relation to goodness-of-fit test and GAN. The goodness-of-fit problem differs from the twosample problem in that one of the two densities is analytically accessible. Using the explicit formula of q, methods based on kernel Stein discrepancy have been developed in (Chwialkowski et al., 2016; Liu et al., 2016; Jitkrittum et al., 2017) and applied to generative model evaluation. However, the computation of the score function ∇ log q may be difficult for certain generative models, including many generative networks. Meanwhile, in many generative models including MMN and GAN the goodness-of-fit is evaluated by batch sampling, i.e. the two-sample setting: Kernel MMD is used in MMN, and GAN, Wasserstein GAN (Arjovsky et al., 2017) and f-GAN (Nowozin et al., 2016) estimate density divergence by a trained network (the D-net). Since the success of GAN training relies on the discriminative power of the D-net, the efficiency of using a neural network for the two-sample test is important for the training and evaluating of such models.

2. LOG-RATIO TEST BY NETWORK CLASSIFIER:



2.1. TWO-SAMPLE PROBLEM:

Formally, the two-sample problem asks to test the null hypothesis H 0 : p = q given datasets X = {x i } n X i=1 and Y = {y j } n Y j=1 where x i ∼ p i.i.d., y j ∼ q i.i.d., and X is independent from Y . The test method is usually based upon a statistic T = T (X, Y ), which is computed from the two datasets, and a test threshold τ , and the null hypothesis H 0 is rejected if T > τ . To control the false discovery under the null, the threshold τ is usually set to the smallest value s.t. Pr[ T > τ |H 0 ] ≤ α, where α ∈ (0, 1) is a pre-specified number called the significance level of the test (typically α = 0.05). Algorithm-wise, τ is given either by theory (the probabilistic distribution of T under H 0 ) or computed from data (Higgins, 2003; Gretton et al., 2012a) .

2.2. TEST STATISTIC AND DENSITY LOG RATIO ESTIMATION:

The proposed test statistic is computed in the following way: given X and Y as above, without loss of generality suppose n = n X + n Y is even integer. Same as in (Lopez-Paz & Oquab, 2016) , we split the dataset z) e u θ (z) +e v θ (z) , where u θ and v θ are activations in the last hidden layer of the network, θ denoting the network parametrization. We define f θ := u θ -v θ , which is the logit, and let the test statistic be

D = {(x i , 0)} n X i=1 ∪ {(y j , 0)} n Y j=1 = {(z i , l i )} n i=1 , l i ∈ {0,

e u θ (z) +e v θ (z) , Pr[l = 1|z] = e v θ (

T = 1 |X te | x∈Xte f θ (x) - 1 |Y te | y∈Yte f θ (y).

(1)

which can be written as T = f θ (x)(p te (x) -qte (x))dx where pte and qte stand for the empirical density of X te and Y te respectively.

As has been pointed out in (Menon & Ong, 2016) , the training of the classifier can be interpreted as minimizing a Bregman divergence between the estimated logit f θ and the true log ratio f * = log p q . Thus the proposed statistic T can be viewed as estimating the symmetric KL divergence (p -q) log p q = KL(p||q) + KL(q||p), which is an f -divergence with f (u) = (u -1) log u (Kanamori et al., 2011; Nowozin et al., 2016) . The testing power of (1) will be theoretically analyzed in Section 4, particularly when p and q lie on or near to low-dimensional manifolds.

2.3. ALGORITHM IN PRACTICE:

Threshold τ . In practice, the test threshold τ can be computed by a permutation test (Higgins, 2003) : randomly permute the |D te | many labels l i on the test set, and recompute the test statistics for m perm times, typically a few tens. This gives an empirical distribution of T under the null hypothesis where both densities equal a mixture of the original p and q. Then τ is set to be the (1-α)-quantile of the empirical distribution so as to control the type-I error to be at most α. Density difference indicator. Many two-sample methods also provide an indication of where q differs from p, which is often of more application interest, e.g., via the witness function in kernel MMD (Gretton et al., 2012a) . For the proposed test, such a differential indicator is provided by the logit f θ of the trained classifier network, which can be viewed as an estimator of log density ratio f * as discussed above. Following kernel MMD, we call f θ (and f * ) the empirical (and population) witness function of the proposed classification logit test.

Network training and computational complexity. Our training of network is conducted via Adam (Kingma & Ba, 2014), the convergence of which has been analyzed in many places such as (Reddi et al., 2019) . We use fixed learning rate over a fixed number of epochs, and it is entirely possible that our training procedure is over simplified and better usage of stochastic gradient descent method as studied in (Bottou, 2010; Zeiler, 2012; Sutskever et al., 2013; Shamir & Zhang, 2013; Du et al., 2018) may lead to improved performance. Given n data samples, evaluating the network output on each sample takes a fixed amount of flops, and thus computing the test statistic takes O(n) operations. The permutation test to determine τ adds negligible cost since f θ has been evaluated at each test sample, and permuting the class labels only reorders these computed values. The training phase is certainly more expensive, though theoretically the overall complexity is O(n) assuming that training is terminated after a fixed number of epochs. Note that the computation can be conducted by batch sampling so the algorithm scales to large sample size and also to multiple sample problems.

3. A ONE-DIMENSIONAL EXAMPLE:

Setting-up. We compare the proposed test based on log ratio (net-logit) to (1) gaussian kernel MMD (gmmd) and ( 2) test based on classification accuracy (net-acc) (Lopez-Paz & Oquab, 2016), on 1D example where the distributions are Test power. The table in Fig. 1 lists the power for the three methods, where net-logit gives significantly better average power ∼ 80%, and the power of net-acc and gmmd are similar, both ∼ 20%. The variation of the power is much larger for the two net-based tests though (c.f. Fig. A.2). We note that such large variation is due to the instability of network training, possibly due to small training size, and is a limitation of the current net-based methods.

EQUATION

One may observe that the above comparison to kernel MMD is not fair: First, kernel MMD with median-distance σ does not use the training set, thus it would be a more fair comparison if gmmd can use all the data samples without training-test splitting. Second, the median setting of σ may not be optimal and can be improved by existing methods in literature. We thus repeat the experiment of gmmd with median-distance σ which uses all the 400 samples ("gmmd+"), and also test over a range of values of σ which are {2 -3 , 2 -2 , • • • 2 3 }, and select the best test power ("gmmd++"). Other experimental setting being the same as above. The results are reported in Table A .1, where gmmd+ achieves a test power of 47% and gmmd++ a power of 57%, remaining inferior to net-logit, while both with small variation (std 2) and thus are more stable than net-based tests. Results with other values of δ and sample sizes are reported in Sec. 5.1, Fig. 2 .

Witness function.:

The three tests all provide witness functions to indicate where q differs from p:

The population witness functions for gmmd is w σ (x) := g σ (x -y)(p(y) -q(y))dy, g σ (z) := e -|z| 2 /(2σ 2 ) , and its empirical counterpart is by replacing p and q with the empirical densities of X te and Y te respectively in the integral. Recall that the population and empirical witness function for net-logit test are f * (x) = log p(x) q(x) and f θ respectively. For net-acc, when |X te | = |Y te |, it is equivalent to using the sign (taking value of ±1) of f θ instead of f θ in computing the test statistic in (1), as shown in (A.1). Thus we call Sign(f θ (x)) the empirical witness function for the net-acc test, and Sign(f * (x)) its population version.

The population and empirical witness functions (in one test run) are plotted in Fig. 1 . Comparing to gmmd, the witness function of net-logit, i.e., the log density ratio, weighs larger at the differential region which is at the tail of the density p. This is also reflected in the empirical witness functions. Taking the sign of f θ as done in classification accuracy test introduces discontinuity of at the decision boundary neat x = 2, which leads to comparatively larger variance in view of the mean of the statistic under H 1 . This intuitively explains why the net-logit test is more powerful here, and a quantitative comparison of mean vs. standard deviation of the three tests is given in Appendix A.2.

4. ANALYSIS OF TEST POWER:

In this section we prove the power of the network logit test assuming that f θ is identified by minimizing the population classification loss. The proof is based on network approximation analysis, and a key question is the needed network complexity, particularly how it scales with the data dimensionality. We first analyze the general case of smooth densities p and q in R D , and then consider the important case where p and q lie on or near a smooth low-dimensional manifold, where we reduce the needed network complexity to depend on the intrinsic manifold geometry only.

4.1. IDENTIFICATION OF f θ BY POPULATION LOSS:

Training with the population loss of the classification network can be expressed as (samples from X and Y are of same number)

EQUATION

where F Θ denotes the class of functions that can be expressed as the difference of the outputs in the last hidden layer of the classification network. A direct verification shows that f * = arg max f L[f ] = log p q (see e.g. (Goodfellow et al., 2014) where it is proved in terms of D = e f 1+e f ), which characterizes the solution of (3) when the function class is arbitrarily large or something large enough to contains f * . Then L[f * ] = 1 2 p log 2p p+q + q log 2q p+q = JSD(p, q), which is the Jensen-Shannon divergence. For certain classification network of finite complexity, f * may not be contained in F Θ , and instead the optimization (3) finds some

EQUATION

We analyze the exact optimization of the population loss only. Since we do not address training from finite samples in this section, we abuse the notation of n (which used to denote

|D tr | + |D te |) to be n = |X te | = |Y te |.

Then after f = f θ is identified by ( 4), the test statistic T n = T can be computed as in (1), and we denote its expectation as

T [f θ ] = f θ (p -q). 4.2 TESTING POWER FOR DISTINGUISHING GENERAL p AND q IN R D

We assume that p and q are smooth densities in R D compactly supported on a bounded region Ω, and without loss of generality Ω = {x ∈ R D , |x| < 1}. The analysis proceeds in the following three steps, all proofs in Appendix B:

Step 1. Use neural network approximation result to construct a network function f con ∈ F Θ that uniformly approximates f * , which bounds

|L[f con ] -L[f * ]| proportionally. This implies that L[f con ] > C for some C ≈ L[f * ] = JSD(p, q) > 0 (c.f. Proposition B.1). The relation L[f θ ] ≥ L[f con ] then gives that L[f θ ] > C > 0.

Step 2. The lower bound of L[f θ ] serves to show that T [f θ ] > 4C under H 1 via a relation between the two (c.f. Lemma B.2).

Step 3. We compute and bound the variance of T n to be at order O(n -1 ), under both H 0 and H 1 (c.f. Proposition B.4). This together with the proved O(1) mean (bias) of T n under H 1 will prove the test power to be asymptotically 1.

Combining these steps gives the following theorem, where we set r = 2 in Proposition B.1: Theorem 4.1. Let the densities p = e u , q = e v be C 2 and supported on Ω, the unit ball in R D , and

u, v ∈ C 2 (Ω), f * = u-v.

If p = q such that JSD(p, q) > 0, then for any small 0 < ε 1 < JSD(p, q), there is a neural network architecture Θ with O(ε

-D/2 1

) many trainable parameters, where the constant depends on the regularity of the second derivative of f * , such that f θ identified by (3) satisfies that

(1) ET n = T [f θ ] > 4C, where C = JSD(p, q) -ε 1 > 0.

(2) For all n, Var(

T n ) = 1 n (Var x∼p (f θ (x)) + Var x∼q (f θ (x))) ≤ 2 n B 2 Θ , B

Θ is a constant depending on the network function family Θ, as defined in (A.3).

(3

) Under H 0 , √ nT n → N (0, σ 2 H0 ) in distribution; Under H 1 , √ n(T n -T [f θ ]) → N (0, σ 2 H1 ) in distribution; σ 2 H0 = 2Var x∼p (f θ (x)), σ 2 H1 = Var x∼p (f θ (x)) + Var x∼q (f θ (x)), both ≤ 2B 2 Θ .

We discuss the extension to less regular f * after Proposition B.1, including the case where p and q nearly non-overlap. In particular, the more relevant situation for two-sample problem is where p and q weakly differ. The above Theorem directly gives the asymptotic test consistency:

Corollary 4.2. Notations and settings as in Theorem 4.1, and for given p, q, the network architecture Θ has been fixed to satisfy C = JSD(p, q) -ε 1 > 0. Let 0 < α < 1 be the two-sample test level, typically α = 0.05, and the test threshold be

τ n = σ H 0 √ n Ψ -1 (α), where Ψ(x) := ∞ x 1 √ 2π e -y 2 /2 dy, then, as n → ∞, Pr[T n > τ n |H 0 ] → α and Pr[T n > τ n |H 1 ] → 1 -Ψ √ n T [f θ ] -τ n σ H1 ≥ 1 -Ψ √ n(4C -τ n ) σ H1

which ≥ 1 -c 0 e -c1n for positive constants c 0 and c 1 .

The analysis reveals that the critical regime for two-sample test is when the divergence between p and q is ∼ O(n -1/2 ). Actually, one may obtain a lower bound of testing power based on the bound of ET n and Var(T n ) in Theorem 4.1 and Chebyshev inequality to control the large deviation. This will prove a finite-sample testing power which is positive and approaching 1 as long as JSD(p, q) exceeds a constant multiple of n -1/2 , where the constant depends on the choice of Θ which guarantees both small ε 1 and B Θ . We omit the details here.

4.3. REDUCING NETWORK COMPLEXITY TO BE INTRINSIC:

The above analysis is for general p and q in R D , and when D is large the needed network complexity grows exponentially. We now show that when p and q are on/near manifold M as in many applications, the network complexity can be reduced to be intrinsic, i.e., depending on the manifold only, and replacing D by the intrinsic dimension d in Theorem 4.1.

On-manifold densities. When both p and q are degenerate and support on the manifold M, all the integrals above, L[f ] and T [f ], are carried out on the manifold only. Specifically, assume that M is compact smooth manifold without boundary, and p and q are smooth on M with respect to the Riemannian geometry, then the log ratio f * is also smooth on M. This allows to apply the manifold function approximation result in (Shaham et al., 2018) to obtain

f con -f * L ∞ (M) < ε 1 with needed network complexity of O( -d/2 1

) (c.f. Theorem B.6). The rest of the proof remains the same by replacing all the integrals in Ω ⊂ R D to be on M.

Near-manifold densities. The analysis extends when p and q decay sufficiently fast away from the manifold, as proved in Proposition B.5. This is because the manifold is locally near Euclidean and there exists differentiable one-to-one mapping between the manifold chart and the local tangent plane on local neighborhoods in R D of radius δ > 0, which is an absolute constant determined by the manifold M. This allows to approximate integrals L[f con ] and L[f * ] in R D by their counterparts on M, using the off-manifold decay of the densities, and the two integrals on M are close due to uniform approximation of f * on the manifold. Proofs in Appendix B.

5. EXPERIMENTS:

This section conducts numerical experiments of the proposed two-sample test and compares with alternatives, on synthetic 1D and manifold densities and evaluating hand-written digits generating models. Codes to produce all experiments will be publicly online.

5.1. SYNTHETIC DATA:

1D normal density departure. The following three examples all have p = N (0, 1), and Eg.1. Mean shift, q = N (δ, 1); Eg.2. Dilation of variance, q = N (0, (1 + δ) 2 ). Eg.3. Mixture with bump at tail, q as in (2). We examine the tests: (1) net-acc, (2) net-logit, (3) gmmd which sets σ to be median distance, (4) gmmd-ad which selects σ by maximizing kernel MMD discrepancy on the training set, and (5) gmmd+, (6) gmmd++ as described in Sec. 3. ( 1)-(4) only use the test set when measuring the power, while (5)-( 6) access both the training and test sets. More details about experimental setting-up in Appendix A.

The test powers of all the methods are plotted in Fig. 2 for the three examples. For Eg.1 and Eg.2, gmmd+, gmmd++ are performing consistently better than the other four which only access the test data set, particularly in Eg.1. Eg.3 has been discussed in Sec. 3, and net-logit gives stronger power than gmmd+, gmmd++ when n all > 200, where net-acc remains inferior to the two. Note that gmmd-ad does worse than gmmd, that is, the median-distance choice of kernel bandwidth σ works better than the adaptive choice here. Among the four methods (1)-( 4), the performances on Eg.1 are comparable, and net-logit gives better power on Eg.2 and Eg. 3. This is especially the case of Eg. 3, where net-logit shows the most significant advantage. Densities on a 2D manifold. The example consists of p and q which lie on the sphere S 2 , a 2dimensional manifold embedded in R 3 . A realization of samples X and Y is shown in the left of Fig. 3 . Fig. 3 plots the test power of the 5 methods over increasing density departure δ and sample size. It can been seen that net-logit gives the fastest growth of power as δ increases and the strongest average power for all n all , but the variation can be large if the power is not close to 1. Unlike some of the 1D cases, gmmd with median σ does not do better than the trivial power for all cases (blue solid), even with access to the full data samples gmmd+ the power is only 0.2 with the largest n all (blue dash). The adaptive choice by maximizing MMD discrepancy on training set improves the power significantly (gmmd-ad, green solid), but does not do as well as net-acc, which again performs inferior to net-logit. The optimal choice of σ (gmmd++, green dash) achieves better power than net-acc at n all = 200 and comparable performance with larger n all . net-logit performs better than gmmd++ and the advantage is more evident when n all > 200. This indicates that larger sample size can be particularly in favor of network-based tests, which rely on the search in the network parameter space optimized on a separated training set.

5.2. GENERATED VS AUTHENTIC MNIST DATA:

As a real-world data example, we study the task of distinguishing "faked" MNIST samples produced by a pre-trained generative network from authentic ones. The MNIST dataset consists of gray-scale hand-written digits of size 28 × 28 falling into 10 classes, which is relatively simple and thus is viewed to lie near to low-dimensional manifolds in the ambient space of R 784 . The classifier network used in net-logit test is a convolutional neural network (CNN) with 2 convolutional layers. More details about the generative and classification networks in Appendix A.

We compare (1) net-acc (2) net-logit (3) gmmd (4) gmmd-ad on two samples X and Y , half of D = X ∪ Y used for training. X consists of authentic MNIST samples, and Y of a mixture of authentic and faked ones, i.e. p = p data and q = (1 -δ)p data + δp model , δ ∈ [0, 1]. The test power is evaluated on 400 test runs and the training is repeated for 20 replicas. The results for increasing δ and sample size n all = |D| up to 500 is shown in Fig. 5 , where net-logit gives the strongest power throughout all cases, and the two network-based tests significantly outperforms the other two when n all ≥ 300. The adaptive choice of kernel bandwidth also improves the power over the mediandistance choice, shown by the better power of gmmd-ad than gmmd. The standard deviation of the net-acc and net-logit power is less than that of gmmd-ad power when n all = 300 and δ ≥ 0.4, when the former two give near to 1 power. We also observe that the training of the CNN classifier in this experiment is more stable than that of the previous fully-connected network on low-dimensional synthetic data, as revealed in the training error evolution plots, c.f. shows the best discriminative power, net-acc gives comparable performance starting n all = 300, while gmmd and gmmd-ad gives trivial power up to n all = 500, c.f. 

6. CONCLUSION:

The paper proposes to use estimated log ratio to compute a two-sample statistic once a classification network has been trained on a split training set. The proposed statistic empirically demonstrates stronger testing power than previously studied neural network classifier tests based on classification accuracy. It also compares favorable to gaussian kernel-based MMD in certain settings, especially for higher dimensional data, including distinguishing generated MNIST digits from authentic ones.

The proposed test has more advantage with large samples, due to that larger training set makes the training more stable, as well as its linear computational complexity and scalable algorithm. Theoretically, we prove the power of the proposed test when the network is sufficiently parametrized, and reduce the needed network complexity to be intrinsic when p and q lie on or sufficiently near to low-dimensional manifolds in possibly high-dimensional space,

The analysis in this paper gives a positive result towards justifying the power of two-sample tests based on training a classification network. However, the proof is based on a network approximation analysis and optimization of population loss, which means that the derived testing power only applied when the training perfectly identifies the global optimizer, a situation not necessarily happening in practice. In experiments we observe that the performance of the network-based tests has larger variance than traditional methods like kernel MMD, due to the instability of the training, particularly with small training size. Thus, more understanding of the network optimization, which is not addressed in the current paper, is needed so as to better understand network classification twosample tests and to develop better methods. At last, we have not systematically explored the effects of different network architecture on the performance, which surely has an influence. Training via SGD with momentum 0.9 produces similar result. The result is qualitatively the same when the number of hidden units varies from 16 to 1024. We have not investigated the optimal choice of network architecture hyperparameters for the two-sample problem.

Adaptive choice of σ in gmmd-ad. In the training phase, the algorithm computes the gaussian kernel MMD discrepancy

TMMD (X, Y ) = 1 |X| 2 x,x ∈X k σ (x, x ) + 1 |Y | 2 y,y ∈Y k σ (y, y ) - 2 |X||Y | x∈X, y∈Y k σ (x, y) on the training set X = X tr , Y = Y tr , for a range of values of the kernel bandwidth σ, i.e. σ = {2 -3 , • • • , 2 3 }. k σ (x, y) = exp{-|x-y| 2

2σ 2 } is the isotropic gaussian kernel. A plot of MMD discrepancy as a function of varying σ is given in Fig. A.1. The σ which maximizes the MMD discrepancy on the training set is then chosen to compute the test statistic on the test set. The method is equivalent to the training process in (Li et al., 2017) with only one trainable parameter which is the kernel bandwidth σ. The MMD test statistic also takes the form as TMMD in (Gretton et al., 2012a; Li et al., 2017) .

A.2 MORE DETAILS ABOUT RESULTS IN SEC. 3:

Test power. The way of computing the test power is empirical and has randomness: for kernel mmd the variation is due to the finite number of runs (n run times), and for network based tests there is extra variation due to the stochastic optimization of the network. Thus we use experiment replicas to recored the variations of the test power.

The empirical distribution of the test power of the three methods over training replicas is given in A.2, corresponding to the experiment in Fig. 1 . The plots of the two net-based methods indicate large variation of the power given by each trained network, that is, the "quality" of the trained net to discriminate the two densities varies. This instability is due to limited training samples as well as the randomness in the optimization algorithm. We observe decreased power variation with larger training set, where the stochastic optimization converges to solutions of lower error and the resulted net gives better two-sample test power. However, the two-sample problem itself is expected to be easier with larger n too. 

| = |Y te | = 1 2 |D te |, Sign(z) = 1 if z ≥ 0 and -1 if z < 0, Tnet-acc = 1 2   1 |X te | x∈Xte 1 {f θ (x)≥0} + 1 |Y te | y∈Yte 1 {f θ (y)<0}   = 1 2   1 |X te | x∈Xte 1 2 (1 + Sign(f θ (x))) + 1 |Y te | y∈Yte 1 2 (1 -Sign(f θ (x)))   = 1 2 + 1 4   1 |X te | x∈Xte Sign(f θ (x)) - 1 |Y te | y∈Yte Sign(f θ (x))   . (A.1)

Quantitative comparison of mean and std of test statistics. Let w be the population witness function of the three methods respectively, and define

Mean := E x∼p,Y ∼q (w(X) -w(Y )), Std := Var x∼p (w(X)) + Var Y ∼q (w(Y )).

For tests using T = w(p -q) as the statistic, such as in net-logit and net-acc, by independence of the samples the mean and variance of T are Mean and Std/ √ n respectively. For kernel MMD, the actually test statistic is computed via quadratic sums, however the mean remains the same, and Std/ √ n will be a lower bound of the standard deviation of the MMD statistic (Serfling, 1981) . Strictly speaking, the test statistic is computed from empirical rather than the population witness function. For net-logit and net-acc, considering population witness function is as if the training is able to identify the exact optimizer which lies inside the representable function family of the network, an idealized scenario. With this idealization, the relation of Mean and Std to the test power can be made rigorous making use of the asymptotic normality of the test statistic (as independent sums), as done in Section 4. The conclusion gives that the larger the Mean, and the smaller the Std, the more powerful the test is going to be. For kernel MMD, these two quantities similarly indicate the testing power, see e.g. (Serfling, 1981; Cheng et al., 2017) . Thus we will use Mean and Std for all three methods for comparison.

To remove the scaling equivalence of test statistics (a test statistic multiplied by a positive constant gives the same test power), we will use the ratio of Mean and Std as an indicator of test power. For the 1D example in Section 3, due to the explicit formula of p and q the values of Mean and Std can be analytically computed, which are shown in Training and testing split is half-and-half in all cases.

Densities p and q on the 2D manifold. The construction of x i ∼ p and y j ∼ q are as below:

x i = T (u i ), y j = T (v j ) where T : R 2 → R 3 is a smooth mapping from unit square to the spherical surface given by

T (x 1 , x 2 ) = 1 R x 1 , x 2 , R 2 -x 2 1 -x 2 2 , R = 1.5,

and u i , v j are i.i.d. copies of random variables u and v in R 2 distributed as

EQUATION

2 ) by angle δ, where the 4 random variables are all independent.

A.4 EXPERIMENT ON MNIST DATA IN SEC. 5.2

The pre-trained generated model is based on a convolutional auto-encoder: c5x5x1x16 -re -ap 2x2 -c5x5x16x32 -re -ap 2x2 -fc128 -re -fc10 -re ← code space R 10 -fc128 -re -ct 5x5x128x32 -re -ct5x5x32x16 (upsample 2x2) -re -ct5x5x16x1 (upsample 2x2) -Euclidean loss where "c" stands for convolutional layer, "ct" for transposed convolutional layers, "re" for Relu activation, and "ap" for average pooling. The auto-encoder is trained on 50000 MNIST dataset for 20 epochs using Adam with learning rate decreasing from 10 -3 to 10 -6 and batch size 100. The sampling of generative model is conducted by adding a small isotropic gaussian noise ("giggering") to the 10-dimensional codes of authentic MNIST digits computed by an encoder, and then mapping through the decoder to R 784 .

We also prepare another generative model by removing the bottleneck layer in the above autoencoder architecture and retrain the model, which gives smaller reconstruction error and a higherdimensional code space of R 128 . The generative model is conducted in the same way by sampling in the code space using gaussian noise of smaller variance per coordinate. This produces faked images that are closer to the authentic ones in Euclidean distance in R 784 , however less explore the "manifold" of p data . The test power of the four methods is shown in Suppose that p and q are non-vanishing inside Ω, and let p = e u and q = e v for smooth potential functions u and v, then f * = u -v is also smooth. We first assume that f * is has properly bounded derivatives, then standard network approximation theory, e.g. that in (Yarotsky, 2017) , guarantees that one can approximate f * by a network output function f con of a multi-layer fully-connected network with

f con -f * L ∞ (Ω) = O(N -r/D

) where N is the number of parameters in the network, and the constant depends on the regularity of r-th derivative of f * . This leads to the following guaranteed positive lower bound of L[f con ]: Proposition B.1. Let the densities p = e u , q = e v be C r and supported on Ω, the unit ball in R D , and u, v ∈ C r (Ω), f * = u -v. If p = q such that JSD(p, q) > 0, then for any 0 < ε 1 < JSD(p, q), are both finite, Ω is the unit sphere in R D , where densities p and q are supported on. Then for any f ∈ F Θ , Var x∼p (f (x)), Var x∼q (f (x)) ≤ min{2B

EQUATION

Θ } 2 := B 2 Θ .

(A.

3)

The boundedness of B

Θ and B

(1)

Θ relies on the regularization of the network function family F Θ . For a given Θ, B (0) Θ and B

(1) Θ may be very large compared to f L ∞ (Ω) and Lip(f ) of the trained f θ , which leads to a loose upper bound of the variance. In practice, regularization techniques may lead to smaller values of |f | and Lip(f ) while at a price of smaller L[f ], revealing the trade-off between stability of the test statistic and the sensitivity to detect differential density departure at large samples. We do not further pursue this problem in the current paper.

B.1.4 PROOFS:

Proof of Proposition B.1. Under the assumptions, f * can by approximated by f con ∈ F Θ such that

f con -f * L ∞ (Ω) ≤ ε 1 .

Writing f con as f , we then have that

|L[f ] -L[f * ]| ≤ 1 2 Ω p log 2e f 1 + e f -log 2e f * 1 + e f * + Ω q log 2 1 + e f -log 2 1 + e f * ≤ 1 2 Ω p(x)|f (x) -f * (x)| + Ω q|f (x) -f * (x)| , ≤ ε 1 ,

where the second inequality uses that log e ξ 1+e ξ and log(1 + e ξ ) are Lip 1. This implies the claim.

Proof of Lemma B.2. By definition,

2L:

[f ] = p log 2 1 + e -f + q log 2 1 + e f = -p log 1 + e -f 2 -q log 1 + e f 2 ≤ -p log e -f 2 -q log e f 2

(for any real number ξ,

1 + e ξ 2 ≥ e ξ 2 ) = f 2 (p -q) = 1 2 T [f ].

Proof of Lemma B.3.

1 2 T [f ] -2L[f ] = p f 2 -q f 2 -p log 2e f 1 + e f -q log 2 1 + e f = p log

1 + e f 2e f /2 + q log 1 + e f 2e f /2 = (p + q) log e -f /2 + e f /2 2 , and by that e x +e -x 2 ≤ e x 2 /2 ,

1 2 T [f ] -2L[f ] ≤ (p + q) log e f 2 /2 = f 2 2 (p + q).

This immediately gives the following lemma Lemma B.9. For any ξ : Ω → R which is Lip continuous, and ξ| M = 0, then for any p ∈ P σ , where pi is an atlas dependent projection of the density p to U i , the explicit formula to be given below, then 

:



Figures and Tables:

1}, into two halves, i.e. D = D tr ∪ D te , |D tr | = |D te | = n 2 , D te = X te ∪ Y te and similarly for the training set. A binary classification neural network is trained on D tr using softmax loss, which gives estimated class probabilities as Pr[l = 0|z] = e u θ (z)

Figure 1: Plots: Top-left: Two densities p and q as in (2); Right three columns: The test statistic T on |Xte| = |Yte| = 100 samples i.e. under H1 (red cross) and the histogram of T under 1000 permutation tests i.e. H0 (blue curve). The population witness function (black curve) and the empirical one evaluated on test data (red cross for Xte, blue crosses for Yte) of the three methods are shown in the bottom row. Table:The mean, standard deviation ("std") and median of the test power computed over nrep = 20 replicas, as described in Sec. 3.the number δ ∈ [0, 1] controls the difference of the densities. We set δ = 0.08, and p and q are illustrated in Fig.1. For both network tests, 200 training and 200 testing samples are used, and same samples in X and Y . Testing power is computed by the frequency of rejecting H 0 in n run = 400 test runs, and n rep = 20 replicas of whole experiment are used to compute mean and standard deviation of the estimated power. The kernel bandwidth σ in gmmd is set to be the median of the pairwise distances among all samples(Gretton et al., 2012a). gmmd can make use of the training set for a more fair comparison, see below. More experimental details in Appendix A.

Figure 2: Three examples of 1D data in Sec. 5.1. Test power of: gmmd (blue), gmmd-ad (green), net-acc (pink), net-logit (red), error bar standing for the standard deviation of the estimated power over 20 replicas, and gmmd+ (blue dash), gmmd++ (green dash). n all = |X| + |Y | including half-half training-testing split.

Figure 3: Test power of the different tests on data on sphere in R 3 in Sec. 5.1. Markers same as in Fig. 2.

Figure 4: Two-sample problem of differentiating p, the density of authentic MNIST digits, and q which contains a δ = 0.4 fraction of digits "faked" by a generative model. |X| = |Y | = 500. The gmmd-ad and net-logit tests use half as training set, and test on the other |Dte| = 500 samples. Left and middle: the most likely fake digits identified by the empirical witness functions of the two tests, red box indicates authentic digits incorrectly identified. Right: The test statistic T (H1) and the histogram of its value under 1000 permutation tests (H0).

Figure 5: Test power of gmmd (blue) gmmd-ad (green) net-acc (pink) net-logit (red) on differentiating authentic vs synthesized MNIST digits produced by a generative model, where sample X has all authentic ones, and δ stands for the fraction of synthesized ones in Y , n all = |X| + |Y | including half-half training-testing split.

Fig. A.4.Setting n all = 1000, δ = 0.4, the results of gmmd-ad and net-logit in one test run is shown in Fig.4. Based on the n all = 500 plot in Fig.5, both tests shall have non-trivial power, and that of net-logit shall be close to 1. In this test, both methods correctly rejects H 0 , yet the net-logit statistic deviates from the distribution of T |H 0 more significantly, indicating stronger power (shown in the histogram plots). To compare the detecting ability of the empirical witness function ŵ of gmmd-ad and net-logit, for each method, we sort the 250 samples in Y te (among which 100 are faked ones) in ascending order of the value of ŵ and select the first 100 samples. These are samples which the model views as most likely to be faked ones. The success rate of identifying faked samples is ∼ 60 by gmmd-ad ŵ, and ∼ 90 by net-logit ŵ. The first 48 most likely faked digits identified with both witness functions are plotted in Fig.4, where gmmd-ad ŵ incorrectly includes 5 authentic samples, and none by net-logit ŵ.

Figure A.1: Left: MMD discrepancy on trained set used by gmmd-ad to select kernel bandwidth σ. Middle and right: training of the classification network used in net-based tests. On the synthetic 1D dataset.

Figure A.2: Histogram of estimated test power from 400 test runs of gmmd, net-acc and net-logit over 20 replicas of training (no training for gmmd), on the example in Sec. 3, Fig. 1.

2: The values of Mean, Std, and their ratio of the three tests, where p and q are as in the 1D example in Sec. 3, c.f. Fig.1.

Figure A.3: Plots of the population witness functions normalized to have Std = 1, of the three tests in Sec. 3, c.f. Fig. 1, Table A.2.

Figure A.4: Same plot as Fig. 5 with another pre-trained generative model which produces faked images that are closer to the authentic ones.

Figure A.5: Same plot as Fig. A.1 on MNIST data.

Fig. A.4. The classification network used in net-logit is the following CNN c5x5x1x16 -re -ap 2x2 -c5x5x16x32 -re -ap 2x2 -fc128 -re -fc2 -softmax loss where dropout is used between the last 2 fully-connected layers. The classification CNN is trained for 100 epochs using Adam with learning rate 10 -3 and batch size 100. A typical plot of evolution of training loss and training error is given in Fig. A.5. The procedure of adaptive selection of σ by gmmd-ad is same as in Sec. 5.1, where the bandwidth search range is σ = {2 -1 , • • • , 2 6 }. B PROOFS AND DETAILS OF THE ANALYSIS IN SECTION 4 B.1 PROOFS AND DETAILS FOR THEOREM 4.1 B.1.1 CONSTRUCTION OF f CON BY NETWORK

|ξ(x)|p(x)dx < Lip(ξ)c 3 σ.Proof of Lemma B.9. By compactness and smoothness of M, for any x ∈ Ω, there exists ψ(x) ∈ M s.t. ψ(x) -x 2 = d(x, M). thus |ξ(x)| = |ξ(x) -ξ(ψ(x))| ≤ Lip(ξ) x -ψ(x) 2 = Lip(ξ)d(x, M). Then Ω |ξ(x)|p(x)dx ≤ Ω Lip(ξ)d(x, M)p(x)dx < Lip(ξ) • c 3 σ,where the last < is by Lemma B.8.The following Lemma fulfills element (1) in the proof. Lemma B.10. Let g : Ω → R has finte Lip(g) and g L ∞ (Ω) , p ∈ P σ , and definep(x) = K i=1 η i (x)p i (x),

p(x)d M (x) ≤ K( g L ∞ (Ω) L η,M + Lip(g)(1 + β M ))c 3 σ. Proof of Lemma B.10. Let H i := φ i (U i ) = φ i (N i ) for each i = 1, • • • , K, η i (x)p(x)dx (p supported on ∪ i N i , c.f. comment after Theorem B.6) ≈ i Ni g(ψ i • φ i (x))η i (x)p(x)dx (error 2) = i Hi (g • η i )(ψ i (u)) η i (z)p i (z)d M (z) (definition of pi ) z)p i (z) d M (z) = M g(z)p(z)d M (z),where d M (z) stands for the Reimannian volume measure on M. Thus it suffices to bound the error in (error 1) and (error 2) and show that the sum ≤ the right hand side of the claim in the Lemma.Bound of (error 1):Ω p(x)g(x)dx -Ω p(x)g(x) K i=1 ηi (x) ≤ g L ∞ (Ω) Ω p(x) 1 -K i=1 ηi (x) dx,and the Lip constant of the function ξ := (1 -K i=1 ηi ) is upper bounded by K i=1 Lip(η i ) ≤ KL η,M by Lemma B.7. Also ξ vanishes on M. Applying Lemma B.9 gives that (error 1) ≤ g L ∞ (Ω) KL η,M c 3 σ.

1: The mean, standard deviation ("std") and median of the test power of the various methods computed from nrun = 400 test runs over nrep = 20 replicas on the 1D example in Sec. 3. The gmmd, net-acc, net-logit tests are computed on |Xte| = |Yte| = 100 samples, where net-acc and net-logit train a classification network on another training set of size |Xtr| = |Ytr| = 100. gmmd only uses the test set and sets the kernel bandwidth σ to be the median distance. gmmd+ and gmmd++ accesses both the training and test sets, where gmmd+ uses the median distance as σ, and gmmd++ reports the best power over varying range of choices of σ, as described in Sec. 3. The results of gmmd, net-acc, net-logit are also reported in Fig.1.

Table content: <table><tr><td/><td colspan="5">gmmd gmmd+ gmmd++ net-acc net-logit</td></tr><tr><td>mean</td><td>19.14</td><td>46.63</td><td>57.29</td><td>19.98</td><td>78.09</td></tr><tr><td>std</td><td>1.95</td><td>2.49</td><td>1.578</td><td>10.43</td><td>20.56</td></tr><tr><td colspan="2">median 19.63</td><td>47.13</td><td>57.38</td><td>17.63</td><td>84.13</td></tr></table>

1 gives the full table of test power including that of the methods gmmd+ and gmmd++. Equivalent form of net-acc test. Here we show that the net-acc test studied in (Lopez-Paz & Oquab, 2016) is equivalent to using Sign(f θ ) instead of f θ in (1) when n X = n Y , up to multiplying and adding constants. Specifically, by the definition of test statistic in (Lopez-Paz & Oquab, 2016), and recall that |X te

Table content: <table/>

Table A.2. The net-logit gives the largest ratio in this comparison. The normalized witness functions are plotted in Fig. A.3, where a constant

Table content: <table><tr><td/><td>Mean</td><td>Std</td><td>Mean/Std</td></tr><tr><td>gmmd</td><td colspan="2">0.0087 0.0421</td><td>0.2075</td></tr><tr><td colspan="3">net-acc 0.1579 0.6087</td><td>0.2594</td></tr><tr><td colspan="3">net-logit 0.2445 0.9011</td><td>0.2714</td></tr></table>

Peer review (iclr2020):
The paper proposes a Neural Networks based algorithm for the two-sample test problem, which is one of the most fundamental problem in statistics and machine learning. Given independent samples from two distributions p and q, the goal is to tell whether p = q or p = q. The test is usually based on some test statistic T and a threshold T_0 where Pr(T > T_0 | p = q) < a (a is the tolerance for the false alarm rate and usually set to be 0.05 for statistical significance). The decision rule is simply testing whether T > T_0 holds. In this paper, the authors propose a test statistic based on the difference of the sums of output logit s(the difference between the output neurons of the last layer) of a trained neural network classifier on both sample sets. The training of the network is done using half of the samples and the threshold is set using a permutation test on the training set. The final decision is based on the outputs of the network on the other half of the data (test set). The paper shows that the test is consistent asymptotically if the network provides a good approximation for the JL divergence between the two distributions, which connects recent advances in network approximation theory. In the experiment part, the paper considers synthetic experiments using 1-d and 2-d Gaussian distribution and real-data experiments using the MNIST dataset to compare the proposed method with other methods including net-acc (using the difference in the output probabilities) and variance kinds of Gaussian kernel MMD methods. The paper presents some nice ideas, but overall I find the theory part not deep and fail to provide enough insights. In the training phase of the network, not only the number of the parameters in the network matters, but the number of training samples is also crucial for the network to approximate the optimal function. It would be nice to have some discussions on this. I also find the experiment part not convincing. In the experiment part, I think the fairest comparison to Gaussian kernel MMD is to use gmmd+ and gmmd++. For the 1-d Gaussian case (figure 2), the experiments show that the proposed test is not always better than these two in some distribution cases. In the 2-d case, the experiment is only done for one distribution. For real dataset (MNIST, figure 5), results for gmmd+ and gmmd++ are not presented. I would expect more thorough experiments on different distributions to make a more convincing statement. So overall I would not recommend this paper for admission. Minor comments: Section 2.2, line 3: y_i*s should have label 1.