Title: STEIN BRIDGING: ENABLING MUTUAL REINFORCE-MENT BETWEEN EXPLICIT AND IMPLICIT GENERA-TIVE MODELS

Abstract: Deep generative models are generally categorized into explicit models and implicit models. The former defines an explicit density form, whose normalizing constant is often unknown; while the latter, including generative adversarial networks (GANs), generates samples without explicitly defining a density function. In spite of substantial recent advances demonstrating the power of the two classes of generative models in many applications, both of them, when used alone, suffer from respective limitations and drawbacks. To mitigate these issues, we propose Stein Bridging, a novel joint training framework that connects an explicit (un-normalized) density estimator and an implicit sample generator with Stein discrepancy. We show that the Stein Bridge induces new regularization schemes for both explicit and implicit models. Convergence analysis and extensive experiments demonstrate that the Stein Bridging i) improves the stability and sample quality of the GAN training, and ii) facilitates the density estimator to seek more modes in data and alleviate the mode-collapse issue. Additionally, we discuss several applications of Stein Bridging and useful tricks in practical implementation used in our experiments.

1. INTRODUCTION:

Deep generative model, as a powerful unsupervised framework for learning the distribution of highdimensional multi-modal data, has been extensively studied in recent literature. Typically, there are two types of generative models (Goodfellow et al., 2014) . Explicit models define an explicit (unnormalized) density function, while implicit models learn to sample from the distribution without explicitly defining a density function.

Explicit models have wide applications in undirected graphical models (LeCun et al., 2006; Salakhutdinov & Hinton, 2009; Hinton et al., 2006; Ngiam et al., 2011) , random graph theory (Robins et al., 2007) , energy-based reinforcement learning (Haarnoja et al., 2017) , etc. However, the unknown normalizing constant makes the model hard to train and sample from, and the explicit models might not be able to capture the complex structure of true samples while maintaining tractability. In contrast, implicit models are more flexible in training and easy to sample from, and in pariticular, generative adverarial networks (GANs) have shown great power in learning representations of images, natural languages, graphs, etc. (Goodfellow et al., 2014; Radford et al., 2016; Arjovsky et al., 2017; Brock et al., 2019) . Nevertheless, due to the minimax game between generator and discriminator/critic in GANs, the training process often suffers from instability, and produces undesirable samples often associated with missing modes in data or generating extra modes out of data. Therefore, it motivates us to consider jointly learning of two models that can presumably compensate and reinforce each other in the training process. (More discussions about related works are in Appendix A.) Furthermore, there are many of situations where we do need both an explicit (unnormalized) density and a flexible implicit sampler. For sample evaluation, it is not enough to merely distinguish samples between real and faked one, and one may also expect to provide fine-grained evaluation on generated samples, where the energy values given by the explicit models can be a good metric (Dai et al., 2017) . Another situation is outlier detection. Implicit models often leverage all true samples (possibly mixed with corrupted samples) as true examples for training. To make up for the issue, explicit models could help to detect out-of-distribution samples via the estimated densities (Zhai et al., 2016) . Also, when given insufficient observed samples, explicit models may fail to capture an accurate distribution, in which case implicit model may help with data augmentation and facilitate training for density estimation. These situations inspire us to combine both of the worlds so as to take advantage of two models in an effective way.

In this work, we aim at jointly learning explicit and implicit generative models. In our framework, an explicit energy model is used to estimate the unnormalized densities of true samples via minimizing a Stein discrepancy; in the meantime, an implicit generator model is exploited to minimize the Wasserstein metric (or Jensen-Shannon divergence) between distributions of true and generated samples. On top of these, another Stein discrepancy, acting as a bridge between implicit generated samples and explicit estimated densities, is introduced and pushes the two models to achieve a consensus. We show that the Stein bridge allows the two generative models to reinforce each other by imposing new regularizations on both models, which help the generator to output high-quality samples and facilitate the energy model to avoid mode-collapse. Moreover, we show that the joint training helps to stabilize GAN training via a convergence analysis. Extensive experiments on various tasks verify our theoretical findings as well as demonstrate the superiority of proposed methods compared with existing deep energy models and GAN-based models.

2. BACKGROUND:

In this section, we briefly provide some technical background used in our model.

Energy Model. The energy model assigns each data x ∈ R d with a scalar energy value E φ (x), where E φ (•) is called the energy function and is parameterized by φ. The model is expected to assign low energy to true samples according to a Gibbs distribution p φ (x) = exp{-E φ (x)}/Z φ , where Z φ is a normalizing constant dependent on φ. The normalizing term Z φ is often hard to compute, making the training intractable, and various methods are proposed to detour such term (see Appendix A).

Stein Discrepancy. Stein discrepancy (Gorham & Mackey, 2015; Liu et al., 2016; Chwialkowski et al., 2016; Oates et al., 2017) is a measure of closeness between two probability distributions that does not require the knowledge for the normalizing constant of one of the compared distributions. Let P and Q be two probability distributions on X ⊂ R d , and assume Q has a (unnormalized) density q. The Stein discrepancy S(P, Q) is defined as S(P, Q) := sup

EQUATION

where F is often chosen to be a Stein class (see, e.g., Definition 2.1 in Liu et al. (2016) ), f : R d → R d is a vector-valued function called Stein critic and φ is an operation that transforms a d × d matrix into a scalar value. One common choice1 of φ is trace operation when d = d. If F is a unit ball in some reproducing kernel Hilbert space (RKHS) with a positive definite kernel k, it induces Kernel Stein Discrepancy (KSD). We provide more details in Appendix B.

Wasserstein Metric. Wasserstein metric is suitable for measuring distances between two distributions with non-overlapping supports (Arjovsky et al., 2017) . The Wasserstein-1 metric between distributions P and Q is defined as W(P, Q) := min γ E (x,y)∼γ [ xy ], where the minimization is over all joint distributions with marginals P and Q. By Kantorovich-Rubinstein duality, it has a dual representation W(P, Q) := max

EQUATION

where the maximization is over all 1-Lipschitz continuous functions.

Sobolev space and Sobolev dual norm. Use L 2 to denote the canonical Hilbert space on R d equipped with an inner product u, v L 2 := R d uvdx. The Sobolev space H 1 is defined as the closure of C ∞ 0 , the set of smooth functions on R d with compact support, with respect to the norm

u H 1 := R d (u 2 + ∇u 2 2 )dx 1/2 . For v ∈ L 2 , its Sobolev dual norm v H -1 is defined by (Evans & Society, 2010) v H -1 := sup u∈H 1 v, u L 2 : R d ∇u 2 2 dx ≤ 1, R d u(x)dx = 0 .

The constraint R d u(x)dx = 0 is necessary in ensuring the finiteness of the supremum. • H -1 can be viewed as a measure of smoothness, which measures the similarity (in terms of largest L 2 -norm) between v and a subset of smooth functions in H 1 .

3. PROPOSED MODEL:

In this section, we formulate our model, Stein Bridging, and highlight its regularization effects.

3.1. MODEL FORMULATION:

We denote by P real the underlying real distribution from which the data {x} are sampled. We simultaneously learn two generative models -one explicit and one implicit -that represent estimates of P real . The explicit generative model has an explicit probability density P E proportional to exp(-E(x)), where E is referred to as an energy function. The implicit generative model transforms an easy-to-sample random noise z with distribution P 0 via a generator G to a generated sample x = G(z) with distribution P G . We use the Stein discrepancy as a measure of closeness between the explicit unnormalized density P E and the real distribution P real , and use the Wasserstein metric as a measure of closeness between the implicit distribution P G and P real .

To jointly learn the two generative models P G and P E , arguably the most straightforward approach is to minimize the sum of the Stein discrepancy and the Wasserstein metric:

min E,G W(P real , P G ) + λS(P real , P E ),

where λ ≥ 0 is a weight coefficient. However, this approach appears no different than learning the two generative models separately. To better train the model, we incorporate the objective another term S(P G , P E ) -called the Stein bridge -that measures the closeness between the explicit unnormalized density P E and the implicit distribution P G :

EQUATION

where λ 1 , λ 2 ≥ 0 are weight coefficients. Although the Stein bridge might seem redundant mathematically, we show that it helps regularize the models in Section 3.2.

The Wasserstein term in (3) is implemented using its equivalent dual representation (2). The two Stein terms in (3) can be implemented using (1) with either a Stein critic parameterized by a neural network, or the Kernel Stein Discrepancy. To reduce the computational cost, the two Stein critics share their parameters, namely, kernels or neural networks. A scheme of our framework is presented in Fig. 1 . We also discuss some related works that attempt to combine both of the worlds (such as energy-based GAN, contrastive learning and cooperative learning) in Appendix A.3, and highlight the difference between our method and theirs in terms of the objective in Table 1 .

Remark. In general, we can also choose other statistical distances in (3) to measure closeness between probability distributions. For example, the Wasserstein metric W(P real , P G ) can be replaced by other common choices for implicit generative models, such as Jensen-Shannon divergence used in the original GAN paper (Goodfellow et al., 2014) . If the normalizing constant of P E is known or easy to calculate, one can replace the Stein discrepancy by the Kullback-Leibler divergence, which is equivalent to the maximum likelihood estimation. We present details for model specifications in various forms and training algorithm in Appendix E.2.

3.2. REGULARIZATION EFFECTS BY VIRTUE OF THE STEIN BRIDGE:

The intuitive motivation of the Stein bridge term in (3) is to push the two models to achieve a consensus. In this subsection, we theoretically show that the Stein bridge allows the two models to reinforce each other by imposing regularizations on the critics.

3.2.1. KERNEL SOBOLEV DUAL NORM REGULARIZATION ON THE WASSERSTEIN CRITIC:

We show the regularization effect of the Stein bridge on the Wasserstein critic. Define the kernel Sobolev dual norm as which can be viewed as a kernel generalization of the Sobolev dual norm defined in Section 2, which reduces to the Sobolev dual norm when k(x, x ) = I(x = x ) and P being the Lebesgue measure.

D H -1 (P;k) := sup u∈C ∞ 0 D, u L 2 (P) : E x,x ∼P [∇u(x) k(x, x )∇u(x )] ≤ 1, E P [h] = 0 .

We have the following result.

Theorem 1. Assume that {P G } G exhausts all continuous probability distributions and S is chosen as kernel Stein discrepancy. Then problem (3) is equivalent to

min E max D E y∼P E [D(y)] -E x∼P real [D(x)] - 1 4λ 2 D H -1 (P E ;k) .

According to Section 2, the regularization term would penalize the non-smoothness of the Wasserstein critic D, which is in the same spirit of gradient-based penalty (e.g., Gulrajani et al. (2017) ; Roth et al. (2017) ), but with a new way to encouraging smoothness.

Another way to interpret the Sobolev dual norm penalty is by observing that if k(x, x ) = I(x = x ) and and E P E [D] = 0 (Villani, 2008) , then

D H -1 (P E ;k) = lim →0 W 2 ((1 + D)P E , P E ) ,

where W 2 denotes the 2-Wasserstein metric. Therefore, the regularization ensures that D would not change suddenly on the high-density region of P E , and the explicit model reinforces the learning of the Wasserstein critic. 

3.2.2. LIPSCHITZ REGULARIZATION:



4. CONVERGENCE ANALYSIS:

In Section 3.2, we justify Stein Bridging by showing the regularization effects. In this section, we further show that it could help to stabilize GAN training with local convergence guarantee. To this end, we first compare the behaviors of WGAN, likelihood-and entropy-regularized WGAN, and our Stein Bridging under SGD via an easy to comprehend toy example. Then we give a formal result that interprets why the introduction of (unnormalized) density estimator could stablize GAN training and help for convergence.

4.1. ANALYSIS OF A LINEAR SYSTEM:

The training for minimax game in GAN is difficult. When using traditional gradient methods, the training would suffer from some oscillatory behaviors (Goodfellow (2017) ; Liang & Stokes (2019) ).

In order to better understand the optimization behaviors, we first study a one-dimension linear system that provides some insights on this problem. 

EQUATION

Unfortunately, such simple objective cannot guarantee convergence by traditional gradient methods like SGD with alternate updating3 :

θ k+1 = θ k + ηψ k ,, ψ k+1 = ψ k + η(1 -θ k+1

). Such optimization would suffer from an oscillatory behavior, i.e., the updated parameters go around the optimum point ([ψ * , θ * ] = [0, 1]) forming a circle without converging to the centrality, which is shown in Fig. 2 (a). A recent study in Liang & Stokes (2019) theoretically show that such oscillation is due to the interaction term in (4).

One solution to the instability of GAN training is to add (likelihood) regularization, which has been widely studied by recent literatures (Warde-Farley & Bengio (2017) ; Li & Turner (2018) ). With regularization term, the objective changes into

min θ max |ψ|≤1 ψE[x] -ψθE[z] -λE[log µ(θz)],

where µ(•) denotes the likelihood function and λ is a hyperparameter. A recent study (Tao et al. (2019) ) proves that when λ < 0 (likelihood-regularization), the extra term is equivalent to maximizing sample evidence, helping to stabilize GAN training; when λ > 0 (entropy-regularization), the extra term maximizes sample entropy, which encourages diversity of generator. Here we consider a Gaussian likelihood function for generated sample x , µ(x ) = exp(-1 2 (x -b)2 ) which is up to a constant, and then the objective becomes (see Appendix D.1 for details): The above system would converge with λ < 0 and diverge with λ > 0 in gradient-based optimization, shown in Fig. 2(a) . Another issue of likelihood-regularization is that the extra term changes the optimum point and makes the model converge to a biased distribution, as proved by Tao et al. (2019) . In this case, one can verify that the optimum point becomes [ψ * , θ * ] = [-λ, 1], resulting a bias. To avoid this issue, Tao et al. (2019) proposes to temporally decrease |λ| through training. However, such method would also be stuck in oscillation when |λ| gets close to zero as is shown in Fig. 2(a) .

EQUATION

1 = 1, 2 = 1 1 = 1, 2 = 5 1 = 1, 2 = 0.2 1 = 5, 2 = 1 1 = 0.2, 2 = 1 (b)

Finally, let us consider our proposed model. We also simplify the density estimator as a basic energy model p φ (x) = exp(-1 2 x 2 -φx) whose score function is ∇ x log p φ (x) = -x -φ. Then if we specify the two Stein discrepancies in (3) as KSD, we have the objective (see Appendix D.1 for details),

EQUATION

)

Interestingly, one can verify that for ∀λ 1 , λ 2 , the optimum point remains the same [ψ * , θ * , φ * ] = [0, 1, -1]. Then we show that the optimization can guarantee convergence to [ψ * , θ * , φ * ]. Proposition 1. Using alternate SGD for (6) geometrically decreases the square norm

EQUATION

In Fig. 2 (a), we can see that Stein Bridging achieves a good convergence to the right optimum.

Compared with (4), the objective (6) adds a new bilinear term φ • θ, which acts like a connection between the two generator and estimator, and two other quadratic terms, which help to push the values to decrease through training. The added terms and the original terms in ( 6) cooperate to guarantee convergence to a unique optimum. (More discussions in Appendix D.1).

We further generalize the analysis to multi-dimensional bilinear system 2019)). For any bilinear system, with the added term

F (ψ, θ) = θ Aψ -b θ - c ψ

H(φ, θ) = 1 2 (θ + φ) B(θ + φ) where B = (AA ) 1 2

to the objective, we can prove that i) the optimum point remains the same as the original system (Proposition 2) and ii) using alternate SGD algorithm for the new objective can guarantee convergence (Theorem 4). The results are given in Appendix D.3.

4.2. LOCAL CONVERGENCE FOR A GENERAL MODEL:

To study the convergence for Stein Bridging, we proceed to consider a general optimization objective

min θ max ψ min φ L(θ, ψ, φ),

where L(θ, ψ, φ) = F (θ, ψ)+H(θ, φ), and

ω f = [θ, ψ] and ω h = [θ, φ] (θ is a shared parameter set). Use ω * = [θ * , ψ * , φ * ] to denote the optimum point of L and ω * f = [θ * , ψ * ], ω * h = [θ * , φ * ] represent the optimum points of F and H respectively. Define Ω f = Ω θ × Ω ψ and Ω h = Ω θ × Ω φ ,

where Ω θ , Ω ψ , Ω φ denote constraint sets for θ, ψ, φ respectively. Function H is µ-strongly convex, and F is µ-strongly convex for θ and µ-strongly concave for ψ (see Appendix E.4 for definition of strongly convex condition). Here we define h

(ω h ) = ∇ θ H + ∇ φ H, f (ω f ) = ∇ θ F -∇ ψ F

, and then we have the following theorem. Theorem 3. If F is µ-strongly convex-concave and H is µ-strongly convex, we can leverage the alternate SGD algorithm, i.e.

EQUATION

EQUATION

where ω

t+1/2 h = [θ t , ψ t ], ω t+1 h = [θ t+1/2 , ψ t+1 ], ω t+1/2 f = [θ t+1/2 , φ t ], ω t+1 f = [θ t+1 , φ t+1 ], and P Ω (ω) = arg min ω ∈Ω ω -ω 2

2 denotes the projection mapping to Ω. Then we can achieve the convergence by using 1 2µ < η < 1 µ .

Theorem 3 shows that Stein Bridging could converge to at least a local optimum. Due to the unknown and intricate landscape of deep neural networks, the global optimization and convergence analysis for GAN has remained as an unexplored problem. Despite the fact that strong convexity assumption cannot be guaranteed with deep neural networks, the optimization could converge to a stable point once there exists a local region that satisfies the strongly convex conditions. In the experiments, we will empirically compare the training stability of each method on various datasets to validate our theoretical discussions. 

5. EXPERIMENTS:

In this section, we conduct experiments to verify the effectiveness of proposed method from multifaceted views. First, we select three tasks with different evaluation metrics in Section 5.1, 5.2 and 5.3. Then we further discuss some applications of joint training as well as some useful tricks in Section 5.4, 5.5 and 5.6. The codes will be released later.

We consider two synthetic datasets with mixtures of Gaussian distributions: Two-Circle and Two-Spiral. The first one is composed of 24 Gaussian mixtures that lie in two circles. Such dataset is extended from the 8-Gaussian-mixture scenario which is widely used in previous GAN papers and is more difficult, so that we can use it to test the quality of generated samples and mode coverage of learned energy. The second synthetic dataset consists of 100 Gaussian mixtures whose centers are densely arranged on two centrally symmetrical spiral-shaped curves. This dataset can be used to examine the power of generative model on complicated data distributions. The ground-truth distributions and samples are shown in Fig. 3 (a) and Fig. 4 (a). Furthermore, we also apply the methods to MNIST and CIFAR datasets which require the model to deal with high-dimensional data. In each dataset, we use observed samples as input of the model and leverage them to train the generators and the estimators. The details for each dataset are reported in Appendix E.1.

In our experiments, we also replace the Wasserstein metric in (3) by JS divergence. To well distinguish different specifications, we term the model Joint-W if using Wasserstein metric and Joint-JS if using JS divergence in this section. We consider several competitors. First, for implicit generative models, we consider valina GAN, WGAN-GP (Gulrajani et al. (2017) ), likelihood-regularized GAN/WGAN-GP (short as GAN+LR/WGAN+LR), entropy-regularized GAN/WGAN-GP (short as GAN+ER/WGAN+ER) and a recently proposed variational annealing regularization (Tao et al. (2019) ) for GAN (short as GAN+VA/WGAN+VA) to compare the quality of generated samples. We employ the denoising auto-encoder to estimate the gradient for regularization penalty, which is proposed by Alain & Bengio (2014) and utilized by Tao et al. (2019) . Second, for explicit models, we consider Deep Energy Model (DEM) which is optimized based on Stein discrepancy, and energybased GAN (EGAN) (Dai et al. (2017) ). Besides, we also compare with Deep Directed Generative (DGM) Model (Kim & Bengio (2017) ) which adopts contrastive divergence to unite sample generator and energy estimator. See Appendix A for brief introduction of these methods and Appendix E.3 for implementation details for each method.

5.1. SAMPLE QUALITY OF IMPLICIT MODEL:

Calibrating explicit (unnormalized) density model with implicit generator is expected to improve the quality of generated samples. In Fig. 3 and Fig. 4 erated samples given by GAN, WGAN-GP and DGM (the worst one in this case) locating between two Gaussian components, and the boundary for each component is not distinguishable. Since the ground-truth densities of regions between two components are very low, such generated samples possess low-quality, which depicts that these models capture the combinations of two dominated features (i.e., modes) in the data but such combination does not make sense in practice. By contrast, Joint-JS and Joint-W could alleviate such issue, reduce the low-quality samples and produce more distinguishable boundaries for components. In Two-Spiral, similarly, the generated samples given by GAN and WGAN-GP form a circle instead of two spirals while the samples of DGM 'link' two spirals. Joint-JS manages to focus more on true high densities compared to GAN and Joint-W provides the best results. To quantitatively measure the sample quality, we adopt two metrics: Maximum Mean Discrepancy (MMD) and High-quality Sample Rate (HSR). The detailed definitions are given in Appendix E.4 and we report the results in Table 5 .

We visualize the generated digits/images on MNIST/CIFAR-10 datasets in Fig. 9 and Fig. 10 and use Inception Score and conditional entropy of predicted classes (CEPC) to measure the sample quality (See Appendix E.4 for details). As shown in Table 2 , Joint-W (resp. Joint-JS) is superior than WGAN-GP (resp. DCGAN), regularized WGAN (resp. DCGAN) and DGM. The CEPC characterizes how well the picture can be distinguished by a pre-trained classifier, i.e., the quality of picture, so the results depict that proposed method could give higher-quality generated pictures. We also rank the generated digits (and true digits) on MNIST w.r.t the densities given by the energy model in Fig. 11 , Fig. 12 and Fig. 13 . As depicted in the figures, the digits with high densities (or low densities) given by Joint-JS possess enough diversity (the thickness, the inclination angles as well as the shapes of digits diverses). By constrast, all the digits with high densities given by DGM tend to be thin and digits with low densities are very thick. Also, as for EGAN, digits with high (or low) densities appear to have the same inclination angle (for high densities, '1' keeps straight and '9' 'leans' to the left while for low densities, just the opposite). Such phenomenon indicates that DGM and EGAN tend to allocate high (or low) densities to data with certain modes and would miss some modes that possibly possess high densities in ground-truth distributions. Fortunately, our method overcomes the issue and manages to capture complicated distributions.

5.2. DENSITY ESTIMATION:



5.3. ENHANCING THE STABILITY OF GAN:

Our discussions and analysis show that joint training helps to stabilize GAN training. In Fig. 5 we present the learning curves of Joint-W (resp. Joint-JS) compared with WGAN (resp. GAN or DCGAN) and its regularization-based variants on different datasets. One can clearly see from the curves that joint training could reduce the variance of metric values especially during the second half of training. Furthermore, we visualize the generated pictures given by the same noise z in adjacent epoches in Fig. 6 . The results show that Joint-W outputs more stable generation in adjacent epoches while the generated samples given by WGAN-GP and WGAN+VA exhibit an obvious variation.

Especially, some digits generated by WGAN-GP and WGAN+VA change from one class to another. Such phenomenon is quite similar to the oscillatory behavior with non-convergence in optimization that we discuss in Section 4.1.

Another issue discussed in Section 4.1 is the bias of model distribution for regularized GAN methods. To quantify this evaluation, we calculate l 1 and l 2 distances between the means of 50000 generated digits (resp. images) and 50000 true digits (resp. images) in MNIST (reps. CIFAR-10).

The results are shown in Table 3 . The smaller distances given by Joint-W indicate that it converges to a better local optimum with smaller bias from the original data distribution. Also, in Table 6 (resp. Table 7 ), we report the distances for digits (resp. images) in each class on MNIST (resp. CIFAR).

5.4. DETECTING OUT-OF-DISTRIBUTION SAMPLES:

The explicit model estimates densities for each sample and one of its applications is to detect outliers in the input data. Here, we adopt CIFAR-10 to measure the ability of our estimator to distinguish the in-distribution samples and (true/false) out-of-distribution samples. We consider four situations and in each case, we consider the test images of CIFAR-10 as positive set (expected to allocate high densities) and construct a negative set (expected to allocate low densities). We let the model output densities for images in two sets, rank them according to the densities and plot the ROC curve for false-positive rate v.s. true-positive rate in Fig. 8 . In the first case, we flip each image in the positive set as negative set. Note that such flipped images are not out-of-distribution samples, so the model is expected to allocate high densities to them, i.e., the ROC curve should be close to a straight line from (0, 0) to (1, 1). The results show that Joint-W, EGAN and DEM give the exact results while DGM assigns all flipped images with lower densities, which means that it fails to capture the semantics in images. In the following three cases, we i) generate random noise, ii) average two images with different CIFAR classes, and iii) adopt Lsun Bedroom dataset as the negative set, respectively. In these situations, the model is expected to distinguish the images in two sets. The results in Fig. 8 show that DGM provides the best results while the performance of Joint-W is quite close to DGM and much better than DEM and EGAN. By contrast, the AUC of Joint-W exhibits a small decline when the sample size is more than 500 and suffers from an obvious decline when it is less than 300. Such phenomenon demonstrates lower sensitivity of joint training to observed sample size.

5.6. WHEN TO START JOINT LEARNING:

In our experiment, we also observe an interesting phenomenon: the performance achieved at convergence would be better if we start joint training after some iterations with independent training for the generator and the estimator. In other words, at the beginning, we could set λ 2 = 0 (or some very small values) in (3) and after some iterations set it as a normal level. We report the inception scores on MNIST with different numbers of iterations for independent training in Fig. 7 (c) where we can see that the score firstly goes up and then goes down when we increase iterations for independent training. Such phenomenon is quite similar to the 'warm up' trick used for training deep networks where one can use small learning rates at iterations in the begining and amplify its value for further training. One intuitive reason behind this phenomenon is that at the beginning, both the generator and estimator are weak and if we minimize the discrepancy between them at this point, they would possibly constrain each other and get limited in some bad local optima. When they become strong enough after some training iterations, uniting them through joint training would help them compensate and reinforce each other as our discussions.

6. CONCLUSIONS:

In this paper, we aim at uniting the training for implicit generative model (represented by GAN) and explicit generative model (represented by a deep energy-based model). Besides two loss terms for GAN and energy-based model, we introduce the third loss characterized via Stein discrepancy between the generator in GAN and the energy-based model. Theoretically, we show that joint training could i) help to stablize GAN training and facilitate its convergence, and ii) enforcing dual regularization effects on both models and help to escape from local optima in optimization. We also conduct extensive experiments with different tasks and application senarios to verify our theoretical findings as well as demonstrate the superiority of our method compared with various GAN models and deep energy-based models.

A LITERATURE REVIEWS:

We discuss some of related literatures and shed lights on the relationship between our work with others.

A 2019)). Apart from DEM, there are also some other forms of deep explicit models based on restricted Boltzmann machines like deep belief networks (Hinton et al. (2006) ) and deep Boltzmann machines (Salakhutdinov & Hinton (2009) ).

The normalized constant under the energy function requires an intractable integral over all possible instances, which makes the model hard to learn via Maximum Likelihood Estimation (MLE). To solve this issue, some works propose to approximate the constant by MCMC methods (Geman & Geman (1984) ; Neal (2011)). However, MCMC requires an inner-loop samples in each training, which induces high computational costs. Another solution is to optimize an alternate surrogate loss function. For example, contrastive divergence (CD) (Liu & Wang (2017) ) is proposed to measure how much KL divergence can be improved by running a small numbers of Markov chain steps towards the intractable likelihood, while score matching (SM) (Hyvärinen (2005) ) detours the constant by minimizing the distance for gradients of log-likelihoods. Moreover, the intractable normalized constant makes it hard to sample from. To obtain an accurate samples from unnormalized densities, many studies propose to approximate the generation by diffusion-based processes, like generative flow (Nguyen et al. (2017) ) and variational gradient descent (Liu & Wang (2016) ). Also, a recent work (Hu et al. (2018) ) leverages Stein discrepancy to design a neural sampler from unnormalized densities. The fundamental disadvantage of explicit model is that the energy-based learning is difficult to accurately capture the distribution of true samples due to the low manifold of real-world instances (Liu & Wang (2017) ).

A.2 IMPLICIT GENERATIVE MODELS:

Implicit generative models focus on a generation mapping from random noises to generated samples. Such mapping function is often called as generator and possesses better flexibility compared with explicit models. Two typical implicit models are Variational Auto-Encoder (VAE) (Kingma & Welling ( 2014)) and Generative Adversarial Networks (GAN) (Goodfellow et al. (2014) ). VAE introduces a latent variable and attempts to maximize the variational lower bound for likelihood of joint distribution of latent variable and observable variable, while GAN targets an adversarial game between the generator and a discriminator (or critic in WGAN) that aims at discriminating the generated and true samples. In this paper, we focus on GAN and its variants (e.g., WGAN (Arjovsky et al. ( 2017)), WGAN-GP (Gulrajani et al. (2017) ), DCGAN (Radford et al. (2016) ), etc.) as the implicit generative model and we leave the discussions on VAE as future work.

Two important issues concerning GAN and its variants are instability of training and local optima. The typical local optima for GAN can be divided into two categories: mode-collapse (the model fails to capture all the modes in data) and mode-redundance (the model generates modes that do not exist in data). Recently there are many attempts to solve these issues from various perspectives. One perspective is from regularization. Two typical regularization methods are likelihood-based and entropy-based regularization with the prominent examples Warde-Farley & Bengio (2017) and Li & Turner (2018) that respectively leverage denoising feature matching and implicit gradient approximation to enforce the regularization constraints. The likelihood and entropy regularizations could respectively help the generator to focus on data distribution and encourage more diverse samples, and a recent work (Tao et al. (2019) ) uses Langevin dynamics to indicate that i) the entropy and likelihood regularizations are equivalent and share an opposite relationship in mathematics, and ii) both regularizations would make the model converge to a surrogate point with a bias from original data distribution. Then Tao et al. (2019) proposes a variational annealing strategy to empirically unite two regularizations and tackle the biased distributions.

To deal with the instability issue, there are also some recent literatures from optimization perspectives and proposes different algorithms to address the non-convergence of minimax game optimization (for instance, Gemp & Mahadevan (2018) ; Liang & Stokes (2019) ; Gidel et al. ( 2019)). Moreover, the disadvantage of implicit models is the lack of explicit densities over instances, which disables the black-box generator to characterize the distributions behind data.

A 2018)) leverages the samples of generator as the MCMC initialization for energy-based model. The synthesized samples produced from finite-step MCMC are closer to the energy model and the generator is optimized to make the finite-step MCMC revise its initial samples. Also, a recent work Du et al. (2018) proposes to regard the explicit model as a teacher net who guides the training of implicit generator as a student net to produce samples that could overcome the mode-collapse issue.

The main drawback of cooperative training is that they indirectly optimize the discrepancy between the generator and data distribution via the energy model as a 'mediator', which leads to a fact that once the energy model gets stuck in a local optimum (e.g., mode-collapse or mode-redundance) the training for the generator would be affected. In other words, the training for two models would constrain rather than exactly compensate each other. In Table 1 , we do a high-level comparison among the above-mentioned generative models w.r.t the objectives. Different from other methods, our model considers three discrepancies simultaneously as a triangle to jointly train the generator and the estimator, enabling them to compensate and reinforce each other.

B BACKGROUND FOR STEIN DISCREPANCY:

Assume q(x) to be a continuously differentiable density supported on X ⊂ R d and f : R d → R d a smooth vector function. Define A q [f (x)] = ∇ x log q(x)f (x) + ∇ x f (x) as a Stein operator. If f is a Stein class (satisfying some mild boundary conditions) then we have the following Stein identity property:

E x∼q [A q [f (x)]] = E x∼q [∇ x log q(x)f (x) + ∇ x f (x)] = 0.

Such property induces the Stein discrepancy between distributions P : p(x) and Q : q(x), x ∈ X :

EQUATION

)

where f is what we call Stein critic that exploits over function space F and if F is large enough then S(Q, P) = 0 if and only if Q = P. Note that in (1), we do not need the normalized constant for p(x) which enables Stein discrepancy to deal with unnormalized density.

If F is a unit ball in a Reproducing Kernel Hilbert Space (RKHS) with a positive definite kernel function k(•, •), then the supremum in (1) would have a close form (see Liu et al. (2016) ; Chwialkowski et al. (2016) ; Oates et al. (2017) for more details):

EQUATION

where where the minimization with respect to E is over all energy functions, the minimization with respect to P is over all probability distributions with continuous density, and the maximization is over all 1-Lipschitz continuous functions. From the definition of kernel Stein discrepancy

u p (x, x ) = ∇ x log p(x) k(x, x )∇ x log p(x ) + ∇ x log p(x) ∇ x k(x, x ) + ∇ x k(x, x ) ∇ x log p(x ) + tr(∇ x,x k(x, x )).

S(P, P E ) = E x,x ∈P [(∇ x log P(x) -∇ x log P E (x)) k(x, x )(∇ x log P(x ) -∇ x log P E (x ))],

S(P, P E ) is infinite if P is not absolutely continuous with respect to P E . Hence, to minimize the objective, it suffices to consider those P's that are absolutely continuous with respect to P E . We set h(x) := dP/dP E (x) -1, where dP/dP E is the Radon-Nikodym derivative. It follows that the problem becomes

min E,h max D E P E [D] + E P E [hD] -E P real [D] + λ 1 S(P real , P E ) + λ 2 • E x,x ∼P [∇ x log(1 + h(x)) k(x, x )∇ x log(1 + h(x ))] ,

where the minimization with respect to h is over all L 1 (P E ) functions with P E -expectation zero.

Fixing E, we claim that we can swap min h and max D . Indeed, without loss of generality, we can restrict D to be such that D(x 0 ) = 0 for some element x 0 , as a constant shift does not change the value of

E P E [(1 + h)D] -E P real [D].

The set of Lipschitz functions that vanish at x 0 is a Banach space, and the subset of 1-Lipschtiz functions is compact Weaver (1999) . Moreover, L 1 (P E ) is also a Banach space. The above verifies the condition of Sion's minimax theorem, and thus the claim is proved.

Swapping min h and max D , we consider

min E max D min h:E P E [h]=0 {E P E [hD] + λ 2 • E x,x ∼P [∇ x log(1 + h(x)) k(x, x )∇ x log(1 + h(x ))]} = min E max D min h:E P E [h]=0 E P E [hD] + λ 2 • E x,x ∼P ∇ x h(x) 1 + h(x) k(x, x ) ∇ x h(x ) 1 + h(x ) = min E max D min h:E P E [h]=0 E P E [hD] + λ 2 • E x,x ∼P E ∇ x h(x) k(x, x )∇ x h(x ) ,

where the second equality follows from the chain rule of the derivative. For the inner minimum, we have that min

h:E P E [h]=0 E P E [hD] + λ 2 • E x,x ∼P E ∇ x h(x) k(x, x )∇ x h(x ) = min r≥0 min h:E P E [h]=0 E P E [hD] + λ 2 r 2 : E x,x ∼P E ∇ x h(x) k(x, x )∇ x h(x ) ≤ r 2 = min r≥0 min h:E P E [h]=0 rE P E [hD] + λ 2 r 2 : E x,x ∼P E ∇ x h(x) k(x, x )∇ x h(x ) ≤ 1 = min r≥0 λ 2 r 2 -r D H -1 (P E ;k) = - 1 4λ 2 D H -1 (P E ;k) ,

where the first equality follows the introduction of the auxiliary variable r, the second equality follows a change of variable from h to rh; and the second equality the third equality follows from the definition of the kernel Sobolev dual norm. Plugging in the previous equation yields the ideal result.

C.2 PROOF FOR THEOREM 2

Proof. Essentially, the result is a consequence of distributionally robust optimization with Wasserstein metric (Gao & Kleywegt, 2016; Gao et al., 2017) . Here we provide a simplified version for completeness.

Using the definition of Stein discrepancy, we rewrite the problem as min

E,P max f {λ 1 S(P real , P E ) + λ 2 E y∼P [A P E f (y)] + W(P real , P)}

Using the compactness of the Stein class and similar argument to the proof of Theorem 1, we can swap the minimization over P and the maximization over f . Fixing f , consider

min P {λ 2 E y∼P [A P E f (y)] + W(P real , P)}.

Using the definition of Wasserstein metric

W(P real , P) = min γ E (x,y)∼γ [ x -y ],

where the minimization is over all joint distributions of (x, y) with x-marginal P real and y-marginal P, we write the problem above as

min P,γ {E (x,y)∼γ [λ 2 A P E f (y) + ||x -y||]},

where γ has marginals P real and P. Since P is unconstrained, the above problem is further equivalent to min

γ {E (x,y)∼γ [λ 2 A P E f (y)] + ||x -y||]},

where the minimization is over all joint distributions of (x, y) with x-marginal being P real . Using the law of total expectation, the problem above is equivalent to

min {γx}x E x∼P real [E y∼γx [λ 2 A P E f (y) + ||x -y|| | x]] = E x∼P real min γx {E y∼γx [λ 2 A P E f (y) + ||x -y|| | x]} = E x∼P real min y∈X {λ 2 A P E f (y) + ||x -y||}

where the minimization in the first equation is over γ x , all conditional distributions of y given x, and the exchanging of min and E follows from the interchangebability principle (Shapiro et al., 2009) ; the second equality holds because the infimum can be restricted to the set of point masses.

We finally have the original problem is equivalent to

min E max f λ 1 S(P real , P E ) + E x∼P real min y∈X {λ 2 A P E f (y) + ||x -y||} .

Hence the proof is completed using the definition of Moreau-Yosida regularization.

D PROOFS AND MORE DISCUSSIONS IN SECTION 4 D.1 DETAILS FOR ONE-DIMENSIONAL CASE:

For the analysis of 1-dim regularized WGAN in section 3.1.1, we assume a Gaussian likelihood function for sample x, µ(x) = exp(-1 2 (x -b) 2 ) which is up to a constant. Its parameter can be estimated by b = E[x]. Then for generated sample x = θz, we have

E(log µ(θz)) = -1 2 E[z 2 ]θ 2 + E[z]E[x]θ -1 2 E[x] 2 . Like the case in WGAN, we consider E[x] = E[z] = 1. Assume Var[z] = 1 and we have E[z 2 ] = 1 + E[z].

Hence, for the analysis on likelihood-(and entropy-) regularized WGAN, we can study the following system: min

θ max ψ ψ -ψ • θ -λ(θ 2 -θ).

When λ = 1, the above objective degrades to (4); when λ < 0 (likelihood-regularization), the the gradient of regularization term pushes θ to shrink, which helps for convergence; when λ > 0 (entropy-regularization), the added term forms an amplifiying strength on θ and leads to divergence.

We proceed to consider 1-dim case of Stein Bridging with energy model

p φ (x) = exp(-1 2 x 2 -φx). If using KSD with kernel k(x 1 , x 2 ) = I(x 1 = x 2 ), then S(P real , P E ) = E x1,x2 [(∇ x1 log p φ (x 1 ) -∇ x1 log µ(x 1 ))k(x 1 , x 2 )(∇ x2 log p φ (x 2 ) -∇ x2 log µ(x 2 ))] = E x [(∇ x log p φ (x) -∇ x log µ(x)) 2 ] = (φ + E[x]) 2 . Similarly, one can obtain S(P G , P E ) = (φ + θE[z]) 2 .

Therefore we arrive at the objective in ( 6)

EQUATION

)

Interestingly, the added terms λ1 2 (1 + φ) 2 + λ2 2 (θ + φ) 2 in ( 6) and the original terms ψ -ψ • θ in WGAN play both necessary roles to guarantee the convergence to the unique optimum points

[ψ * , θ * , φ * ] = [0, 1, -1].

If we remove the critic and optimize θ and φ with the remaining loss terms, we would find that the training would converge but not necessarily to [ψ * , θ * ] = [0, 1] (since the optimum points are not unique in this case). On the other hand, if we remove the estimator, the system degrades to (4) and would not converge to the unique optimum point [ψ * , θ * ] = [0, 1]. If we consider both of the world and optimize three terms together, the training would converge to a unique global optimum

[ψ * , θ * , φ * ] = [0, 1, -1].

D.2 PROOF FOR PROPOSITION 1:

Proof. Instead of directly studying the optimization for (6), we first prove the following problem will converge to the unique optimum,

EQUATION

)

Applying alternate SGD we have the following iterations:

ψ t+1 = ψ t + η * θ t , φ t+1 = φ t -η * (θ t + 2φ t ) = (1 -2η)φ t -ηθ t , θ t+1 = θ t -η(ψ t+1 + φ t+1 + θ t ) = -η(1 -2η)φ t + (1 -η)θ t -ηψ t .

Then we obtain the relationship between adjacent iterations:

ψ t+1 φ t+1 θ t+1 = 1 0 η 0 1 -2η -η -η -η(1 -2η) 1 -η • ψ t φ t θ t = M • ψ t φ t θ t

We further calculate the eigenvalues for matrix M and have the following equations (assume the eigenvalue as λ):

(λ -1) 3 + 3η(λ -1) 2 + 2η 2 (1 + η)(λ -1) + 2η 3 = 0.

One can verify that the solutions to the above equation satisfy |λ|

< (1 -η + η 2 )(1 + η -η 2 ).

Then we have the following relationship

ψ t+1 φ t+1 θ t+1 2 2 = [ψ t φ t θ t ] • M M • ψ t φ t θ t 2 2 ≤ λ 2 m • ψ t φ t θ t 2 2

where λ m denotes the eigenvalue with the maximum absolute value of matrix M . Hence, we have

ψ 2 t+1 + φ 2 t+1 + θ 2 t+1 ≤ (1 -η + η 2 )(1 + η -η 2 )[ψ 2 t + φ 2 t + θ 2 t ].

We proceed to replace ψ, φ and θ in (13) by ψ , φ and θ respectively and conduct a change of variable: let θ = 1 -θ and φ = -1 -φ. Then we get the conclusion in the proposition.

D.3 GENERALIZATION TO BILINEAR SYSTEMS:

Our analysis in the one-dimension case inspires us that we can add affiliated variable to modify the objective and stabilize the training for general bilinear system. The bilinear system is of wide interest for researchers focusing on stability of GAN training (Goodfellow (2017) ; Liang & Stokes (2019) ; Gidel et al. (2019) ; Gemp & Mahadevan (2018) ). The general bilinear function can be written as

EQUATION

where ψ, θ are both r-dimensional vectors and the objective is min θ max ψ F (ψ, θ) which can be seen as a basic form of various GAN objectives. Unfortunately, if we directly use simultaneous (resp. alternate) SGD to optimize such objectives, one can obtain divergence (resp. fluctuation).

To solve the issue, some recent papers propose several optimization algorithms, like extrapolation from the past (Gidel et al. ( 2019)), crossing the curl (Gemp & Mahadevan (2018) ) and consensus optimization (Liang & Stokes (2019) ). Also, Liang & Stokes (2019) shows that it is the interaction term which generates non-zero values for ∇ θψ F and ∇ ψθ F that leads to such instability of training. Different from previous works that focused on algorithmic perspective, we propose to add new affiliated variables which modify the objective function and allow the SGD algorithm to achieve convergence without changing the optimum points.

Based on the minimax objective of ( 14) we add affiliated r-dimensional variable φ (corresponding to the estimator in our model) the original system and tackle the following problem:

EQUATION

where

H(φ, θ) = 1 2 (θ + φ) B(θ + φ), B = (AA )

1 2 and α is a non-negative constant. Theoretically, the new problem keeps the optimum points of ( 14 Proof. The condition tells us that ∇ θ F (ψ * , θ) = 0 and ∇ ψ F (ψ, θ * ) = 0. Then we derive the gradients for L(ψ, φ, θ),

EQUATION

EQUATION

EQUATION

Combining ( 17) and ( 18) we get φ * = -θ * . Hence, the optimum point of ( 15) is [ψ * , θ * , φ * ] where φ * = -θ * .

The advantage of the new problem is that it can be solved by SGD algorithm and guarantees convergence theoretically. We formulate the results in the following theorem. 

EQUATION

we can achieve convergence to [ψ * , θ * , φ * ] where φ * = -θ * with at least linear rate of (1

-η 1 + η 2 2 )(1 + η 2 -η 2 1 )

where η 1 = ησ min , η 2 = ησ max and σ min (resp. σ max ) denotes the maximum (resp. minimum) singular value of matrix A.

To prove Theorem 3, we can prove a more general argument.

Lemma 1. If we consider any first-order optimization method on (15), i.e.,

ψ t+1 ∈ ψ 0 + span(L(ψ 0 , φ, θ), • • • , F (ψ t , φ, θ)), ∀t ∈ N, φ t+1 ∈ ψ 0 + span(L(ψ, φ 0 , θ), • • • , L(ψ, φ t , θ)), ∀t ∈ N, θ t+1 ∈ ψ 0 + span(L(ψ, φ, θ 0 ), • • • , L(ψ, φ, θ t )), ∀t ∈ N, Then we have ψ t = V (ψ t -ψ * ), φ t = U (φ t -φ * ), θ t = U (θ t -θ * ),

where U and V are the singular vectors decomposed by matrix A using SVD decomposition, i.e., A = UDV and the triple ([

ψ t ] i , [ φ t ] i , [ θ t ] i ) 1≤i≤r

follows the update rule with step size σ i η as the same optimization method on a unidimensional problem

EQUATION

)

with step size η, where σ i denotes the i-th singular value on the diagonal of D.

Proof. The proof is extended from the proof of Lemma 3 in Gidel et al. (2019) . The general class of first-order optimization methods derive the following updations:

ψ t+1 = ψ 0 + t+1 s=0 ρ st (A θ s -c) = ψ 0 + t+1 s=0 ρ st A (θ s -θ * ), φ t+1 = φ 0 + 1 2 t+1 s=0 δ st (B + B )(θ s + φ s ), θ t+1 = θ 0 + t+1 s=0 µ st [A(ψ s -ψ * ) + 1 2 (B + B )(θ s + φ s )],

where ρ st , δ st , µ st ∈ R depend on specific optimization method (for example, in SGD, ρ tt = δ tt = µ tt remain as a non-zero constant for ∀t and other coefficients are zero).

Using SVD A = UDV and the fact θ * = -φ * , B = (UDD U ) = D, we have

V (ψ t+1 -ψ * ) = V (ψ 0 -ψ * ) + t+1 s=0 ρ st D U (θ s -θ * ) U (φ t+1 -φ * ) = U (φ 0 -φ * ) + t+1 s=0 δ st U D(θ s -θ * ) + U D(φ s -φ * ), U (θ t+1 -θ * ) = U (θ 0 -θ * ) + t+1 s=0 ρ st [DV (ψ s -ψ * ) + U D(θ s -θ * ) + U D(φ s -φ * )],

and equivalently,

ψ t+1 = ψ 0 + t+1 s=0 ρ st D θ t , φ t = φ 0 + t+1 s=0 δ st D( θ t + φ t ), θ t+1 = θ 0 + t+1 s=0 ρ st D( ψ t + θ t + φ t ).

Note that D is a rectangular matrix with non-zero elements on a diagonal block of size r. Hence, the above r-dimensional problem can be reduced to r unidimensional problems:

[ ψ t+1 ] i = [ ψ 0 ] i + t+1 s=0 ρ st σ i [ θ t ] i , [ φ t ] i = [ φ 0 ] i + t+1 s=0 δ st σ i ([ θ t ] i + [ φ t ] i ), [ θ t+1 ] i = [ θ 0 ] i + t+1 s=0 ρ st σ i ([ ψ t ] i + [ θ t ] i + [ φ t ] i ).

The above iterations can be conducted independently in each dimension where the optimization in i-th dimension follows the same updating rule with step size σ i η as problem in (20).

Furthermore, since problem (20) can achieve convergence with a linear rate of (1-η+η 2 )(1+η-η 2 ) using alternate SGD (the proof is similar to that of (( 13))), the multi-dimensional problem in (15) can achieve convergence by SGD with at least a rate of (1-η 1 +η 2 2 )(1+η 2 -η 2 1 ) where η 1 = ησ max , η 2 = ησ min and σ max (resp. σ min ) denotes the maximum (resp. minimum) singular value of matrix A. We conclude the proof for Theorem 4.

Theorem 4 suggests that the added term H(φ, θ) with affiliated variables φ could help the SGD algorithm achieve convergence to the the same optimum points as directly optimizing F (ψ, θ). Our method is related to consensus optimization algorithm (Liang & Stokes (2019) ) which adds a regularization term ∇ θ F (ψ, θ) + ∇ ψ F (ψ, θ) to ( 14) resulting extra quadratic terms for θ and ψ. The disadvantage of such method is the requirement of Hessian matrix of F (ψ, θ) which is computational expensive for high-dimensional data. By contrast, our solution only requires the first-order derivatives.

D.4 STRONGLY CONVEXITY:

In section 3.1.2, we assume H(θ, φ) as a µ-strongly convex function which indicates that it satisfies the conditions:

(∇ θ H(θ, •) -∇ θ H(θ , •)) (θ -θ ) ≥ µ θ -θ 2 2 , ∀θ, θ ∈ Ω θ , (∇ φ H(•, φ) -∇ φ H(•, φ )) (φ -φ ) ≥ µ φ -φ 2 2 , ∀φ, φ ∈ Ω φ .

Bedises, F (θ, ψ) is µ-strongly convex for θ and µ-strongly concave for ψ so it satisfies:

(∇ θ F (θ, •) -∇ θ F (θ , •)) (θ -θ ) ≥ µ θ -θ 2 2 , ∀θ, θ ∈ Ω θ , (∇ ψ F (•, ψ ) -∇ ψ F (•, ψ)) (ψ -ψ ) ≥ µ ψ -ψ 2 2 , ∀ψ, ψ ∈ Ω ψ .

In section 3.1.2, we also define h(ω h ) = ∇ θ H + ∇ φ H and f (ω f ) = ∇ θ F -∇ ψ F , so the above condition can be written in a more compact form,

(h(ω h ) -h(ω h )) (ω h -ω h ) ≥ µ ω h -ω h 2 2 , ∀ω h , ω h ∈ Ω h , (f (ω f ) -f (ω f )) (ω f -ω f ) ≥ µ ω f -ω f 2 2 , ∀ω f , ω f ∈ Ω f .

D.5 PROOF FOR THEOREM 2:

The proof relies on two lemmas, Lemma 2. For any ω ∈ Ω and ω + = P Ω (ω + u), then we have

ω + -ω 2 2 ≤ u (ω + -ω).

Proof. Since ω + is a projection of ω + u on a convex set Ω, we have

EQUATION

Rearranging the above inequality one can easily get the lemma.

Lemma 3. If function Φ(ω) is µ-strongly convex, we have µ ω -ω * 2 2 ≤ ∇Φ(ω) (ω -ω * ). Similarly, if Φ(ω) is µ-strongly concave, we have µ(ω -ω * ) ≤ ∇ -Φ(ω) (ω -ω * ).

Proof. By optimality of ω * , we have

∇Φ(ω * ) (ω -ω * ) ≥ 0. Since Φ is µ-convex, we can further derive µ ω -ω * 2 2 ≤ ∇Φ(ω * ) (ω -ω * ) + µ ω -ω * 2 2 ≤ ∇Φ(ω) (ω -ω * ).

Proof. (Proof for Theorem 3) We apply Lemma 2 to (9) with (ω, u, ω

+ ) = (ω t+1/2 f , -ηf (ω t+1/2 f ), ω t+1

f ) and we have

EQUATION

Then we have

EQUATION

According to Lemma 3, we have

EQUATION

)

Plug ( 24) into ( 23) and we get

EQUATION

)

The above inequality is equivalent to

EQUATION

Similarly, one can obtain

EQUATION

i.e.,

EQUATION

) Combining ( 26) and (28) we have

ω t+1 -ω * 2 2 = θ t+1 -θ * 2 2 + ψ t+1 -ψ * 2 2 + φ t+1 -φ * 2 2 ≤ (1 -2ηµ) θ t -θ * 2 2 + (2 -2ηµ) ψ t+1 -ψ * 2 2 + (2 -2ηµ) φ t+1 -φ * 2 2 ≤ (2 -2ηµ) ω t -ω * 2 2 . (29) Hence, if 1 2µ < η < 1 µ we have 0 < 2 -2ηµ < 1 and ω t -ω * 2 2 ≤ (2 -2ηµ) t ω 0 -ω * 2 2 E DETAILS FOR EXPERIMENT SETUP E.1 SYNTHETIC DATASETS

We provide the details for two synthetic datasets. The Two-Circle dataset consists of 24 Gaussian mixtures where 8 of them are located in an inner circle with radius r 1 = 4 and 16 of them lie in an outer circle with radius r 2 = 8. For each Gaussian component, the covariance matrix is 0.2 0 0 0.2 = σ 1 I and the mean value is [r 1 cos t, r 1 sin t], where t = 2π•k 8 , k = 1, • • • , 8, for the inner circle, and [r 2 cos t, r 2 sin t], where t = 2π•k 16 , k = 1, • • • , 16 for the outer circle. We sample N 1 = 2000 points as true observed samples for model training. In section 5.5, we consider noised data scenario. In this case, we randomly add n noise points sampled from Gaussian distribution N (0, σ 0 I) where σ 0 = 2 to the original true samples. Here we set n = [40, 100, 160, 300, 400, 600, 800, 1000] .

The Two-Spiral dataset contains 100 Gaussian mixtures whose centers locate on two spiral-shaped curves. For each Gaussian component, the covariance matrix is 0.5 0 0 0.5 = σ 2 I and the mean value is [-c 1 cos c 1 , c 1 sin c 1 ], where c 1 = 2π 3 + linspace(0, 0.5, 50) • 2π, for one spiral, and [c 2 cos c 2 , -c 2 sin c 2 ], where c 2 = 2π 3 + linspace(0, 0.5, 50) • 2π for another spiral. We sample N 2 = 5000 points as true observed samples. In section 5.5, we consider insufficient data scenario. In this case, the sample size N 2 is reduced to [100, 200, 300, 500, 700, 1000, 2000] . 

D 1 D 2 D 3 Objective W S S min θ min φ max ψ max π E x∼P data [d ψ (x)] -E z∼p0 [d ψ (G θ (z))] +λ 1 E x∼P data [A p φ [f π (x)]] + λ 2 E z∼ 0 [A p φ [f π (G θ (z))]]

where D 1 , D 2 , D 3 denote three general discrepancy measures for distributions. As stated in our remark, D 1 can be specified as arbitrary discrepancy measures for implicit generative models. Here we also use JS divergence, the objective for valina GAN. To well distinguish them, we call the model using Wasserstein distance (resp. JS divergence) as Joint-W (resp. Joint-JS) in our experiments. On the other hand, the two Stein discrepancies in (3) can be specified by KSD (as defined by S k in (11)) or general Stein discrepancy with an extra critic (as defined by S in (1)). Hence, the two specifications for D 1 and the two for D 2 (D 3 ) compose four different combinations in total, and we organize the objectives in each case in Table 4 .

In our experiments, we use KSD with RBF kernels for D 2 and D 3 in Joint-W and Joint-JS on two synthetic datasets. For MNIST with conditional training (given the digit class as model input), we also use KSD with RBF kernels. For MNIST and CIFAR with unconditional training (the class is not given as known information), we find that KSD cannot provide desirable results so we adopt general Stein discrepancy for higher model capacity.

The objectives in Table 4 appear to be comutationally expensive. In the worst case (using general Stein discrepancy), there are two minimax operations where one is from GAN or WGAN and one is from Stein discrepancy estimation. To guarantee training efficiency, we alternatively update the generator, estimator, Wasserstein critic and Stein critic over the parameters θ, φ, ψ and π respectively. Specifically, in one iteration, we optimize the generator over θ and the estimator over φ with one step respectively, and then optimize the Wasserstein critic over ψ with n d steps and the Stein critic over π with n c steps. Such training approach guarantees the same time complexity order of proposed method as that of GAN or WGAN, and the training time for our model can be bounded within constant times the time for training GAN model. In our experiment, we set n d = n c = 5 and empirically find that our model Stein Bridging would be two times slower than WGAN on average. We present the training algorithm for Stein Bridging in Algorithm 1.

E.3 IMPLEMENTATION DETAILS:

We give the information of network architectures and hyper-parameter settings for our model as well as each competitor in our experiments. '2', '3', '4', '5', '6', '7', '8' and '9' stand for 'airplane', 'automobile', 'bird', 'cat', 'deer', 'dog', 'frog', 'horse', 'ship' and 'truck' respectively.) (a) Randomly sampled over all digits (b) Randomly sampled over digits with top 50% densities 

:

One can also use other forms for φ, like matrix norm when d = d(Liu et al. (2016)).

Our theoretical discussions focus on WGAN, and we also compare with original GAN in the experiments.

Here, we adopt the most widely used alternate updating strategy. The simultaneous updating, i.e., θ k+1 = θ k + ηψ k and ψ k+1 = ψ k + η(1 -θ k ), would diverge in this case.

Figures and Tables:

Figure 1: Model framework for Stein Bridging which jointly train an implicit sample generator and an explicit (unnormalized) density estimator via a Stein bridge.

Figure 2: (a) Numerical iterations for SGD training of WGAN, likelihood-regularized WGAN (WGAN+LR), variational annealing for WGAN+LR (WGAN+LR+VA), entropy-regularized WGAN (WGAN+ER) and our Stein Bridging. (b) Stein Bridging with different λ 1 and λ 2 .

Figure 3: (a) True samples and (b)∼(f) generated samples produced by the generators of different methods on Two-Circle (upper line) and Two-Spiral (bottom line) datasets.

Figure 4: (a) True densities and (b)∼(f) estimated densities given by the estimators of different methods on Two-Circle (upper line) and Two-Spiral (bottom line) datasets.

Figure 5: Learning curves of Joint-W (resp. Joint-JS) compared with WGAN (resp. GAN or DC-GAN) and its regularization-based variants.

Figure 6: Generated digits (resp. images) given by the same noise z in adjacent training epoches on MNIST (reps. CIFAR) dataset.

Figure 8: ROC curves for evaluation of outlier detection on CIFAR-10.

5.5 ADDRESSING DATA INSUFFICIENCY AND NOISY DATAWe proceed to test the model performance in some extreme situations where the observed samples are mixed with noises or the observed samples are quite insufficient. The results are presented in Fig.7(a) where we add different ratios of random noise to the true samples in Two-Circle dataset and Fig. 7(b) where we only sample insufficient data for training in Two-Spiral dataset. The details are in Appendix E.1. The noise in data impacts the performance of WGAN and Joint-W, but comparatively, the performance decline for Joint-W is less insignificant than WGAN, which indicates better robustness of joint training w.r.t noised data. In Fig. 7(b), when the sample size decreases from 2000 to 100, the AUC value of DEM declines dramatically, showing its dependency on sufficient training samples.

) unchanged. Let L(ψ, φ, θ) = F (ψ, θ) + αG(φ, θ) Proposition 2. Assume the optimum point of min θ max ψ (ψ, θ) are [ψ * , θ * ],then the optimum points of (15) would be [ψ * , θ * , φ * ] where φ * = -θ * .

φ, θ) using alternate SGD algorithm, i.e.,

S k S k min θ min φ max ψ E x∼P data [d ψ (x)] -E z∼p0 [d ψ (G θ (z))] +λ 1 E x,x ∼P data [u p φ (x, x )] + λ 2 E z,z ∼p0 [u p φ (G θ (z), G θ (z ))] J S S S min θ min φ max ψ max π E x∼Pr [log(d ψ (x))] + E z∼p0 [log(1 -d ψ (G θ (z)))] +λ 1 E x∼P data [A p φ [f π (x)]] + λ 2 E z∼ 0 [A p φ [f π (G θ (z))]] J S S k S k min θ min φ max ψ E x∼Pr [log(d ψ (x))] + E z∼p0 [log(1 -d ψ (G θ (z)))] +λ 1 E x,x ∼P data [u p φ (x, x )] + λ 2 E z,z ∼p0 [u p φ (G θ (z), G θ (z ))]

Figure 9: Generated digits given by Joint-W on MNIST.

Figure 10: Generated images given by Joint-W on CIFAR.

Figure11: The generated digits (and real digits) with the highest densities and the lowest densities given by Joint-W.

Figure12: The generated digits (and real digits) with the highest densities and the lowest densities given by DGM.

Figure13: The generated digits (and real digits) with the highest densities and the lowest densities given by EGAN.

Comparison of objectives between different generative models, where D

Table content: <table/>

ON THE STEIN CRITICWe next investigate how the Stein bridge helps to regularize the Stein critic. Recall that the two Stein terms in (3) share the same Stein critic. We have the following result.Theorem 2 shows that the Stein bridge, together with the Wasserstein metric W(P real , P G ), plays as a smoothness regularization on the Stein critic f via Moreau-Yosida regularization, which smoothens the Stein operator A P E f and further encourages the energy model to seek more modes in data instead of focusing on some dominated modes, thus helps alleviate the mode-collapse issue. To the best of our knowledge, this suggests a novel regularization scheme for Stein-based GAN.

Table content: <table><tr><td colspan="3">Theorem 2. Assume {P G } G exhausts all continuous probability distributions, and the Stein class</td></tr><tr><td colspan="3">defining the Stein discrepancy is compact (in some linear topological space). Then problem (3) is</td></tr><tr><td>equivalent to</td><td/><td/></tr><tr><td>min E</td><td>max f</td><td>{λ 1 S(P real , P E ) + λ 2 E x∼P real [M λ2A P E f (x)]},</td></tr><tr><td colspan="3">where M λ2A P E f (•) denotes the (generalized) Moreau-Yosida regularization of the function A P E f with parameter λ 2 , i.e., M λ2A P E f (x) = min y∈X {A P E f (y) + 1 λ2 ||x -y||}.</td></tr></table>

Inception scores (higher is better) and conditional entropies (short as CEPC and lower is better) on MNIST and CIFAR-10. † We directly use the best result reported in their paper.

Table content: <table><tr><td colspan="3">MNIST (Conditional)</td><td colspan="3">MNIST (Unconditional)</td><td colspan="3">CIFAR-10 (Unconditional)</td></tr><tr><td>Method</td><td colspan="2">Score CEPC</td><td>Method</td><td colspan="2">Score CEPC</td><td>Method</td><td colspan="2">Score CEPC</td></tr><tr><td>DCGAN</td><td>8.43</td><td>0.168</td><td>WGAN-GP</td><td>7.71</td><td>0.256</td><td>WGAN-GP</td><td>6.80</td><td>0.153</td></tr><tr><td>DCGAN+LR</td><td>8.40</td><td>0.171</td><td>WGAN+LR</td><td>7.82</td><td>0.243</td><td>WGAN+LR</td><td>6.89</td><td>0.154</td></tr><tr><td>DCGAN+ER</td><td>8.33</td><td>0.179</td><td>WGAN+ER</td><td>7.75</td><td>0.252</td><td>WGAN+ER</td><td>6.99</td><td>0.156</td></tr><tr><td>DCGAN+VA</td><td>8.40</td><td>0.172</td><td>WGAN+VA</td><td>7.74</td><td>0.254</td><td>WGAN+VA</td><td colspan="2">6.95 † 0.154</td></tr><tr><td>DGM</td><td>8.15</td><td>0.201</td><td>DGM</td><td>6.87</td><td>0.372</td><td>DGM</td><td>4.79</td><td>0.146</td></tr><tr><td>Joint-JS(ours)</td><td>8.53</td><td>0.156</td><td>Joint-W(ours)</td><td>7.90</td><td>0.231</td><td>Joint-W(ours)</td><td>7.11</td><td>0.151</td></tr></table>

OF EXPLICIT MODELAnother advantage of joint learning is that the generator could help the density estimator to capture more accurate distribution. As shown in Fig3, both Joint-JS and Joint-W manage to capture all Gaussian components while other methods miss some of modes. In Fig4, Joint-JS and Joint-W exactly fit the ground-truth distribution. By contrast, DEM misses one spiral while EGAN degrades to a uniform-like distribution. DGM manages to fit two spirals but allocate high densities to regions that have low densities in the groung-truth distribution. To quantitatively measure the performance, we introduce three evaluation metrics: KL & JS divergence between the ground-truth and estimated densities and Area Under the Curve (AUC) for false-positive rate v.

Table content: <table/>

Distances between means of generated digits (resp. images) and ground-truth digits (resp. images) on MNIST (resp. CIFAR-10).

Table content: <table/>

.3 ATTEMPTS TO COMBINE BOTH OF THE WORLDS Recently, there are several studies that attempt to combine explicit and implicit generative models from different ways. For instance,Zhao et al. (2017) proposes energy-based GAN that leverages energy model as discriminator to distinguish the generated and true samples. The similar idea is also used byKim & Bengio (2017) andDai et al. (2017) which let the discriminator estimate a scaler energy value for each sample. Such discriminator is optimized to give high energy to generated samples and low energy to true samples while the generator aims at generating samples with low energy. The fundamental difference is that Zhao et al. (2017) andDai et al. (2017) both aim at minimizing the discrepancy between distributions of generated and true samples while the motivation ofKim & Bengio (2017) is to minimize the KL divergence between estimated densities and true samples.Kim & Bengio (2017) adopts contrastive divergence (CD) to link MLE for energy model over true data with the adversarial training of energy-based GAN. However, both CD-based method and energy-based GAN have limited power for both generator and discriminator. Firstly, if the generated samples resemble true samples, then the gradients for discriminator given by true and generated samples are just the opposite and will counteract each other, and the training will stop before the discriminitor captures accurate data distribution. Second, since the objective boils down to minimizing the KL divergence (forKim & Bengio (2017)) or Wasserstein distance (forDai et al.

Table content: <table/>

P real [D] + λ 1 S(P real , P E ) + λ 2 S(P, P E )},

Table content: <table><tr><td/><td>The (11) gives the Kernel Stein Discrepancy</td></tr><tr><td>(KSD).</td><td/></tr><tr><td colspan="2">C PROOFS OF RESULTS IN SECTION 3.2</td></tr><tr><td colspan="2">C.1 PROOF OF THEOREM 1</td></tr><tr><td colspan="2">Proof. Using Kantorovich's duality, we rewrite the problem as</td></tr><tr><td>min E,P</td><td>max</td></tr></table>

Objectives for different specifications of D 1 (P real , P G ), D 2 (P real , P E ) and D 3 (P G , P E ). We specify D 1 as Wasserstein distance or JS divergence in our paper and for D 2 and D 3 we consider the general Stein discrepancy or kernel Stein discrepancy. Here we use W, J S to denote Wasserstein distance and JS divergence respectively, and S, S k to represent general Stein discrepancy and kernel Stein discrepancy respectively. We omit the gradient penalty term for Wasserstein distance here but use it in experiments.E.2 MODEL SPECIFICATIONS AND TRAINING ALGORITHMIn different tasks, we consider different model specifications in order to meet the demand of capacify as well as test the effectiveness under various settings. Our proposed framework (3) adopts Wasserstein distance for the first term and two Stein discrepancies for the second and the third terms. We can write (3) as a more general form min θ,φ D 1 (P real , P G ) + λ 1 D 2 (P real , P E ) + λ 2 D 3 (P G , P E ),

Table content: <table/>

l1 and l2 distances between means of true images and generated images in each class on CIFAR. (Class '0', '1',

Table content: <table/>

Peer review (iclr2020):
This paper proposes a training objective that combines three terms: * A Stein discrepancy for learning a energy model with intractable normalizing constant * A Wasserstein GAN objective for learning an implicit neural sampler. * A Stein discrepancy for minimizing the distance between distributions defined by the energy model and the GAN. The third term is called *Stein bridging* by the authors. It seems pretty difficult to motivate such a bridging term because from the first glance this term does not add anything to learning the two models from data. So I*m wondering how the authors motivate themselves to study this modification. This is my main concern about the paper if this term appears simply because that energy model has an unnormalized density while from GANs we can sample, and Stein discrepancy is best applicable to such pair of distributions. Apart from the concern on motivations, I tried to follow the arguments in Section 3, as the bridging term is justified as regularization to both models. However, I think the proof of Theorem 1 is incorrect: * I don*t think it makes sense from 
abla log (1 + h(x)) to (1 + h(x))
abla h(x), even if the taylor expansion suggested by the authors is applied. * In the next step, *Consider a further approximation*, this approximation basically sets 1 + h(x) to 1, if h(x) is approximately 0, then P=P_E.. Minor: * IN proof of Theorem 1, the expectation should be always over x,x*~P_E instead of x~P, right?