GENERATIVE LEARNING WITH EULER PARTICLE TRANSPORT
Year: 2020

Abstract:
We propose an Euler particle transport (EPT) approach for generative learning. The proposed approach is motivated by the problem of finding the optimal transport map from a reference distribution to a target distribution characterized by the Monge-Ampere equation. Interpreting the infinitesimal linearization of the MongeAmpere equation from the perspective of gradient flows in measure spaces leads to a stochastic McKean-Vlasov equation. We use the forward Euler method to solve this equation. The resulting forward Euler map pushes forward a reference distribution to the target. This map is the composition of a sequence of simple residual maps, which are computationally stable and easy to train. The key task in training is the estimation of the density ratios or differences that determine the residual maps. We estimate the density ratios (differences) based on the Bregman divergence with a gradient penalty using deep density-ratio (difference) fitting. We show that the proposed density-ratio (difference) estimators do not suffer from the “curse of dimensionality” if data is supported on a lower-dimensional manifold. Numerical experiments with multi-mode synthetic datasets and comparisons with the existing methods on real benchmark datasets support our theoretical results and demonstrate the effectiveness of the proposed method.

heading: 1 INTRODUCTION
text: The ability to efficiently sample from complex distributions plays a key role in a variety of prediction and inference tasks in machine learning and statistics (Salakhutdinov, 2015). The long-standing methodology for learning an underlying distribution relies on an explicit statistical data model, which can be difficult to specify in many applications such as image analysis, computer vision and natural language processing. In contrast, implicit generative models do not assume a specific form of the data distribution, but rather learn a nonlinear map to transform a reference distribution to the target distribution. This modeling approach has been shown to achieve impressive performance in many machine learning tasks (Reed et al., 2016; Zhu et al., 2017). Generative adversarial networks (GAN) (Goodfellow et al., 2014), variational auto-encoders (VAE) (Kingma & Welling, 2014) and flow-based methods (Rezende & Mohamed, 2015) are important representatives of implicit generative models.
In this paper, we propose an Euler particle transport (EPT) approach for learning a generative model by integrating ideas from optimal transport, numerical ODE, density-ratio estimation and deep neural networks. We formulate the problem of generative learning as that of finding a nonlinear transform that pushes forward a reference to the target based on the quadratic Wasserstein distance. Since it is challenging to solve the resulting Monge-Ampère equation, we consider the continuity equation derived from the linearization of the Monge-Ampère equation, which is a gradient flows converging to the target distribution. We solve the Mckean-Vlasov equation associated with the gradient flow using the forward Euler method. The resulting EPT that pushes forward the reference distribution to the target distribution is a composition of a sequence of simple residual maps, which are computationally stable and easy to train. The residual maps are completely determined by the density ratios between the distributions at the current iterations and the target distribution. We estimate density ratios based on the Bregman divergence with a gradient regularizer using deep density-ratio fitting.
We establish bounds on the approximation errors due to linearization of the Monge-Ampère equation, Euler discretization of the Mckean-Vlasov equation, and deep density-ratio estimation. Our result on
the error rate for the proposed density-ratio estimators improves the minimax rate of nonparametric estimation via exploring the low-dimensional structure of the data and circumvents the “curse of dimensionality”. Experimental results on multi-mode synthetic data and comparisons with stateof-the-art GANs on benchmark data support our theoretical findings and demonstrate that EPT is computationally more stable and easier to train than GANs. Using simple ReLU ResNets without batch normalization and spectral normalization, we obtained results that are better than or comparable with those using GANs trained with such tricks.

heading: 2 EULER PARTICLE TRANSPORT
text: Let X ∈ Rm be a random vector with distribution ν, and let Z be a random vector with distribution µ. We assume that µ has a known and simple form. Our goal is to construct a transformation T such that T#µ = ν, where T#µ denotes the push-forward distribution of µ by T , that is, the distribution of T (Z). Then we can sample from ν by first generating a Z ∼ µ and calculate T (Z). In practice, ν is unknown and only a random sample {Xi}ni=1 i.i.d. ν is available. We must construct T based on the sample.
There may exist multiple transports T with T#µ = ν. The optimal transport is the one that minimizes the quadratic Wasserstein distance between µ and ν defined by
W2(µ, ν) = { inf γ∈Γ(µ,ν)
E(Z,X)∼γ [‖Z −X‖22]} 1 2 , (1)
where Γ(µ, ν) denotes the set of couplings of (µ, ν) (Villani, 2008; Ambrosio et al., 2008). Suppose that µ and ν have densities q and p with respect to the Lesbeque measure, respectively. Then the optimal transport map T such that T#µ = ν is characterized by the Monge-Ampère equation (Brenier, 1991; McCann, 1995; Santambrogio, 2015). Specifically, the minimization problem in (1) admits a unique solution γ = (1, T )#µ with T = ∇Ψ, µ-a.e., where 1 is the identity map and ∇Ψ is the gradient of the potential function Ψ : Rm → R. This function is convex and satisfies the Monge-Ampère equation
det(∇2Ψ(z)) = q(z) p(∇Ψ(z)) , z ∈ Rm. (2)
Therefore, to find the optimal transport T , it suffices to solve (2) for Ψ. However, it is challenging to solve this degenerate elliptic equation due to its highly nonlinear nature.
Below we describe the proposed EPT method for obtaining an approximate solution of the MongeAmpère equation (2). It consists of the following steps: (a) linearizing (2) via residual maps, (b) determining the velocity fields governing the stochastic McKean-Vlasov equation resulting from the linearization, (c) calculating the forward Euler particle transport map and, (d) training the EPT map by estimating the velocity fields from data. Since velocity fields are completely determined by density ratios, this step amounts to nonparametric density ratio estimation. We also provide bounds on the errors due to linearization, discretization and estimation. Mathematical details and proofs are given in the appendix.
Linearization via residual map A basic approach to addressing the difficulty due to nonlinearity is linearization. We use a linearization method based on the residual map
Tt,Φt = ∇Ψ = 1+ t∇Φt, t ≥ 0, (3)
where Φt : Rm → R1 is a function to be chosen such that the law of Tt,Φt(Z) approaches ν as t increases (Villani, 2008). We give the specific form of Φt below, see Theorem B.1 in the appendix for details.
This linearization scheme leads to the stochastic process Xt : Rm → Rm satisfying the McKeanVlasov equation
d dt Xt(x) = vt(Xt(x)), t ≥ 0, with X0 ∼ µ, µ- a.e. x ∈ Rm, (4)
where vt is the velocity vector field of Xt. In addition, we have vt = ∇Φt. Thus vt also determines the residual map (3). The details of the derivation are given in Theorems B.2 and B.1. in the appendix. Therefore, estimating the residual map (3) is equivalent to estimating vt.
The movement of Xt along t is completely governed by vt, given the initial value. We choose a vt to decrease the discrepancy between the distribution of Xt, say µt, at time t and the target ν with respect to a properly chosen measure. An equivalent formulation of (4) is through the gradient flow {µt}t≥0 with {vt}t≥0 as its velocity fields, see Proposition B.1 in the appendix. Computationally it is more convenient to work with (4).
Determining velocity field The basic intuition is that we should move in the direction that decreases the differences between µt and the target ν. We use an energy functional L[µt] to measure such differences. An important energy functional L[µt] is the f -divergence (Ali & Silvey, 1966),
L[µt] = Df (µt‖ν) = ∫ Rm p(x)f ( qt(x) p(x) ) dx, (5)
where qt is the density of µt, p is the density of ν and f : R+ → R is assumed to be a twicedifferentiable convex function with f(1) = 0. We choose Φt such that L[µt] is minimized. We show in Theorem B.1 in the appendix that Φt(x) = −f ′(rt(x)) and vt(x) = ∇Φt(x). Therefore,
vt(x) = −f ′′(rt(x))∇rt(x), where rt(x) = qt(x)
p(x) , x ∈ Rm.
For example, if we use the χ2-divergence with f(c) = (c−1)2/2, then vt(x) = ∇rt(x) is simply the gradient of the density ratio. Other types of velocity fields can be obtained by using different energy functionals such as the Lebesgue norm of the density difference, i.e., L[µt] = ∫ Rm |qt(x)−p(x)|
2dx, see Section B.2 for details.
The forward Euler method Numerically, we need to discretize the McKean-Vlasov equation (4). Let s > 0 be a small step size. We use the forward Euler method defined iteratively by:
Tk = 1+ svk, (6) Xk+1 = Tk(Xk), (7) µk+1 = (Tk)#µk, (8)
where X0 ∼ µ, µ0 = µ, vk is the velocity field at the kth step, k = 0, 1, ...,K for some large K. The particle process {Xk}k≥0 is a discretized version of the continuous process {Xt}t≥0 in (4). The final transport map is the composition of a sequence of simple residual maps T0, T1 . . . , TK , i.e., T = TK ◦ TK−1 · · · ◦ T0. This updating scheme is based on the forward Euler method for solving equation (4). This is the reason we refer to the proposed method as Euler particle transport (EPT).
Training EPT When the target ν is unknown and only a random sample is available, it is natural to learn ν by first estimating the discrete velocity fields vk at the sample level and then plugging the estimator of vk in (6). For example, if we use the f -divergence as the energy functional, estimating vk(x) = −f ′′(rk(x))∇rk(x) boils down to estimating the density ratios rk(x) = qk(x)/p(x) dynamically at each iteration k. Nonparametric density-ratio estimation using Bregman divergences and gradient regularizer are discussed in Section 4 below. Let v̂k be the estimated velocity fields at the kth iteration. The kth estimated residual map is T̂k = 1+ sv̂k. Finally, the trained map is
T̂ = T̂K ◦ T̂K−1 ◦ · · · ◦ T̂0. (9)
Theoretical guarantees We establish the following bound on the approximation error due to the linearization of the Monge-Ampère equation under appropriate conditions:
W2(µt, ν) = O(e−λt), (10) for some λ > 0, see Proposition B.1 in the appendix. Therefore, µt converges to ν exponentially fast as t→∞. For an integer K ≥ 1 and a small s > 0, let {µst : t ∈ [ks, (k + 1)s), k = 0, . . . ,K} be a piecewise constant interpolation between µks and µ(k+1)s, k = 0, 1, . . . ,K. Under the assumption that the velocity fields vt are Lipschitz continuous with respect to (x, µt), the discretization error of µst can be bounded in a finite time interval [0, T ) as follows:
sup t∈[0,T )
W2(µt, µst ) = O(s). (11)
The proof of (11) is given in Proposition B.2 in the appendix. The error bounds (10) and (11) imply that the distributions of the particles Xk generated by the EPT map defined in (7) with a small s and a sufficiently large k converges to the target ν at the rate of discretization size s.
When training the EPT map, we use the deep neural networks to estimate the density ratios (density differences) with samples. In Theorem 4.1, we provide an estimation error bound that improves the minimax rate of deep nonparametric estimation via exploring the low-dimensional structure of data and circumvents the “curse of dimensionality.” Thus this result is of independent interest in nonparametric estimation using deep neural networks.

heading: 3 IMPLEMENTATION
text: We now described how to implement EPT and train the optimal transport T with an i.i.d. sample {Xi}ni=1 ⊂ Rm from an unknown target distribution ν. The EPT map is trained via the forward Euler iteration (6)-(8) with a small step size s > 0. The resulting map is a composition of a sequence of residual maps, i.e., TK ◦ TK−1 ◦ ... ◦ T0 for a large K. As implied by Theorem 4.1 in Section 4, each Tk, k = 0, ...,K can be estimated with high accuracy by T̂k = 1 + sv̂k, where v̂k(x) = −f ′′(R̂φ(x))∇R̂φ(x). Here R̂φ is the density-ratio estimator defined in (14) below based on {Yi}ni=1 ∼ qk and the data {Xi}ni=1 ∼ p. Therefore, according to the EPT map (9), the particles

heading: T̂ (Ỹi) ≡ T̂K ◦ T̂K−1 ◦ ... ◦ T̂0(Ỹi), i = 1, . . . , n
text: serve as samples drawn from the target distribution ν, where particles {Ỹi}ni=1 ⊂ Rm are sampled from a simple reference distribution µ.
In many applications, high-dimensional complex data such as images, texts and natural languages, tend to have low-dimensional latent features. To learn generative models with latent low-dimensional structures, it is beneficial to have the option of first sampling particles {Zi}ni=1 from a low-dimensional reference distribution µ̃ ∈ P2(R`) with ` d. Then we apply T̂ to particles Ỹi = Gθ(Zi), i = 1, ..., n, where we introduce another deep neural network Gθ : R` → Rm with parameter θ. We can estimate Gθ via fitting the pairs {(Zi, Ỹi)}ni=1. We describe the EPT algorithm below.
• Outer loop for modeling low dimensional latent structure (optional) – Sample {Zi}ni=1 ⊂ R` from a low-dimensional reference distribution µ̃ and let Ỹi = Gθ(Zi), i = 1, 2, . . . , n.
– Inner loop for finding the push-forward map ∗ If there are no outer loops, sample Ỹi ∼ µ, i = 1, . . . , n. ∗ Get v̂(x) = −f ′′(R̂φ(x))∇R̂φ(x) via solving (14) below with Yi = Ỹi. Set T̂ = 1 + sv̂ with a small step size s.
∗ Update the particles Ỹi = T̂ (Ỹi), i = 1, . . . , n. – End inner loop – If there are outer loops, update the parameter θ of Gθ(·) via solving
minθ ∑n i=1 ‖Gθ(Zi)− Ỹi‖22/n.
• End outer loop

heading: 4 DEEP DENSITY-RATIO AND DENSITY-DIFFERENCE FITTING
text: The evaluation of velocity fields depends on the dynamic estimation of a discrepancy between the push-forward distribution qt and the target distribution p. Density-ratio and density-difference fitting with the Bregman score provides a unified framework for such discrepancy estimation without estimating each density separately (Gneiting & Raftery, 2007; Dawid, 2007; Sugiyama et al., 2012a;b; Kanamori & Sugiyama, 2014).
Let r(x) = q(x)/p(x) be the density ratio between a given density q(x) and the target p(x). Let g : R → R be a differentiable and strictly convex function. The separable Bregman score with the base probability density p for measuring the discrepancy between r and a measurable function R : Rm → R1 is
B(r,R) = EX∼p[g′(R(X))R(X)− g(R(X))]− EX∼q[g′(R(X))].
Here we focus on the widely used least-squares density-ratio (LSDR) fitting with g(c) = (c− 1)2 as a working example, i.e.,
BLSDR(r,R) = EX∼p[R(X)2]− 2EX∼q[R(X)] + 1. (12)
For other choice of g, such as g(c) = c log c− (c+ 1) log(c+ 1) corresponding to estimating r via the logistic regression (LR), and the scenario of density difference fitting will be presented in detail in Section B.3.1.
Gradient regularizer The distributions of real data may have a low-dimensional structure with their support concentrated on a low-dimensional manifold, which may cause the f -divergence to be ill-posed due to non-overlapping supports. To exploit such underlying low-dimensional structures and avoid ill-posedness, we derive a simple weighted gradient regularizer 12Ep[g
′′(R)‖∇R‖22], motivated by recent works on smoothing via noise injection (Sønderby et al., 2017; Arjovsky & Bottou, 2017). This serves as a regularizer for deep density-ratio fitting. For example, with g(c) = (c − 1)2, the resulting gradient regularizer is Ep[‖∇R‖22], (13) which recovers the well-known squared Sobolev semi-norm in nonparametric statistics. Gradient regularization stabilizes and improves the long time performance of EPT. The detailed derivation is presented in Section B.3.2.
LSDR estimation with gradient regularizer Let {Xi}ni=1 and {Yi}ni=1 be two collections of i.i.d data from densities p(x) and q(x), respectively. Let H ≡ HD,W,S,B be the set of ReLU neural networks Rφ with parameter φ, depth D, widthW , size S, and ‖Rφ‖∞ ≤ B. We combine the least squares loss (12) with the gradient regularizer (13) as our objective function. The resulting gradient regularized LSDR estimator of r = p/q is given by
R̂φ ∈ arg min Rφ∈H
1
n n∑ i=1 [Rφ(Xi) 2 − 2Rφ(Yi)] + α 1 n n∑ i=1 ‖∇Rφ(Xi)‖22, (14)
where α ≥ 0 is a regularization parameter. Estimation error bound We first show that the density ratio r is identifiable through the objective function by proving that, at the population level, we can recover the density ratio r via minimizing
BαLSDR(R) = BLSDR(r,R) + αEp[‖∇R‖22] + C,
where BLSDR is defined in (12) and C = EX∼q[r2(X)]− 1. Lemma 4.1. For any α ≥ 0, we have r ∈ arg minRBαLSDR(R). In addition, BαLSDR(R) ≥ 0 for any R with EX∼pR2(X) <∞, and BαLSDR(R) = 0 iff R(x) = r(x) = 1 (q, p)-a.e. x ∈ Rm.
This identifiabiity result shows that the target density ratio is the unique minimizer of the population version of the empirical criterion in (14). This provides a the basis for establishing the convergence result of deep nonparametric density-ratio estimation.
Next we bound the nonparametric estimation error ‖R̂φ − r‖L2(ν) under the assumptions that the support of ν is concentrated on a compact low-dimensional manifold and r is Lipsichiz continuous. Let M ⊆ [−c, c]m be a Riemannian manifold (Lee, 2010) with dimension m, condition number 1/τ , volume V , geodesic covering regularityR, and m M = O (m ln(mVR/τ)) m. Denote M = {x ∈ [−c, c]m : inf{‖x− y‖2 : y ∈M} ≤ } , ∈ (0, 1). Theorem 4.1. Assume supp(r) = M and r(x) satisfies |r(x)| ≤ B for a finite constant B > 0 and is Lipschitz continuous with Lipschitz constant L. Suppose the topological parameter ofHD,W,S,B in (14) with α = 0 satisfies D = O(log n),W = O(n M 2(2+M) / log n), S = O(n M−2 M+2 / log4 n), and B = 2B. Then,
E{Xi,Yi}ni=1 [‖R̂φ − r‖ 2 L2(ν)] ≤ C(B 2 + cLmM)n−2/(2+M),
where C is a universal constant.
The error bound established in Theorem 4.1 for the nonparametric deep density-ratio fitting is new. This result is of independent interest for nonparametric estimation with deep neural networks.
The above derived rate O(n− 2
2+M lnm ) is faster than the optimal rate of convergence for nonparametric estimation of a Lipschitz target in Rm, where the optimal rate is O(n− 2 2+m ) (Stone, 1982; Schmidt-Hieber, 2020) as long as the intrinsic dimensionM of the data is much smaller than the ambient dimension m. Therefore, the proposed density-ratio estimators circumvent the “curse of dimensionality” if data is supported on a lower-dimensional manifold.

heading: 5 RELATED WORK
text: We discuss connections between EPT and the existing related works. The existing generative models, such as VAEs, GANs and flow-based methods, parameterize a transform map with a neural network, say G, that solves
min G D(G#µ, ν), (15)
where D(·, ·) is an integral probability discrepancy. The original GAN (Goodfellow et al., 2014), f -GAN (Nowozin et al., 2016) and WGAN (Arjovsky et al., 2017) solve the dual form of (15) by parameterizing the dual variable using another neural network with D as the JS-divergence, the f -divergence and the 1-Wasserstein distance, respectively. Based on the fact that the 1-Wasserstein distance can be evaluated from samples via linear programming (Sriperumbudur et al., 2012), Liu et al. (2018) and Genevay et al. (2018) proposed training the primal form of WGAN via a two-stage method that solves the linear programm. SWGAN (Deshpande et al., 2018) and MMDGAN (Li et al., 2017; Binkowski et al., 2018) use the sliced quadratic Wasserstein distance and the maximum mean discrepancy (MMD) as D, respectively.
Vanilla VAE (Kingma & Welling, 2014) approximately solves the primal form of (15) with the KL-divergence loss under the framework of variational inference. Several authors have proposed methods that use optimal transport losses, such as various forms of Wasserstein distances between the distribution of learned latent codes and the prior distribution as the regularizer in VAE to improve performance. These methods include WAE (Tolstikhin et al., 2018), Sliced WAE (Kolouri et al., 2019) and Sinkhorn AE (Patrini et al., 2019).
Discrete time flow-based methods minimize (15) with the KL divergence loss (Rezende & Mohamed, 2015; Dinh et al., 2015; 2017; Kingma et al., 2016; Papamakarios et al., 2017; Kingma & Dhariwal, 2018). Grathwohl et al. (2019) proposed an ODE flow approach for fast training in such methods using the adjoint equation (Chen et al., 2018b). By introducing the optimal transport tools into maximum likelihood training, Chen et al. (2018a) and Zhang et al. (2018) considered continuous time flow. Chen et al. (2018a) proposed a gradient flow in measure spaces in the framework of variational inference and then discretized it with the implicit movement minimizing scheme (De Giorgi, 1993; Jordan et al., 1998). Zhang et al. (2018) considered gradient flows in measure spaces with time invariant velocity fields. CFGGAN (Johnson & Zhang, 2018) derived from the perspective of optimization in the functional space is a special form of EPT with L[·] taken as the KL divergence. SW flow (Liutkus et al., 2019) and MMD flow (Arbel et al., 2019) are gradient flows in measure spaces. MMD flow can be recovered from EPT by first choosing L[·] as the Lebesgue norm and then projecting the corresponding velocity vector fields onto reproducing kernel Hilbert spaces, please see Appendix B.4 for a proof. However, neither SW flow nor MMD flow can model hidden low-dimensional structure with the particle sampling procedure.
SVGD in (Liu, 2017) and the proposed EPT are both particle methods based on gradient flow in measure spaces. However, the SVGD samples from an un-normalized density, while EPT focuses on generative leaning, i.e., learning the distribution from samples. At the population level, projecting the velocity fields of EPT with KL divergence onto reproducing kernel Hilbert Spaces will recover the velocity fields of SVGD. The proof is given in Appendix B.5. Score-based methods in (Song & Ermon, 2019; 2020; Ho et al., 2020) are also particle methods based on unadjusted Langevin flow and deep score estimators. At the population level, the velocity fields of these score-based methods are random since they have a Brownian motion term, while the velocity fields of EPT are deterministic. At the sample level, these score-based methods need to learn a vector-valued deep score function. while in EPT we need to estimate the density ratios which are scalar functions.

heading: 6 EXPERIMENTS
text: The implementation details on numerical settings, network structures, SGD optimizers and hyperparameters are given in the appendix. All experiments are performed using NVIDIA Tesla K80 GPUs. The PyTorch code of EPT is available at https://github.com/anonymous/EPT.
2D Examples. We use EPT to learn 2D distributions adapted from Grathwohl et al. (2019) with multiple modes and density ridges. The first row in Figure 1 shows kernel density estimation (KDE) plots of 50k samples from target distributions including (from left to right) 8Gaussians, pinwheel, moons, checkerboard, 2spirals, and circles.
The second and third rows show the KDE plots of learned samples via EPT with f -divergence/ Lebesgue norm (left six of the second/ third row), and surface plots of estimated density ratio/ difference after 20k iterations of EPT with f -divergence/ Lebesgue norm (right six of the second/ third row), respectively. Clearly, the generated samples via EPT are nearly indistinguishable from those of the target samples and the estimated density-ratio/ difference functions are approximately equal to 1/0, indicating the learnt distribution matches the target well.
We further visualize the transport maps learned with 5squares and large4gaussians from 4squares and small4gaussians, respectively. We use 200 particles connected with grey lines to manifest the learned transport maps. As shown in the left two figures in Figures 2, the central squares of 5squares were learned better with the gradient penalty, which is consistent with the result of the estimated density-ratio right two figures in Figure 2. For large4gaussians, the learned transport map exhibited some optimality under quadratic Wasserstein distance due to the obvious correspondence between the samples in left two figures in Figure 2.
4 4
Results on Benchmark Image Data. We show the performance of applying EPT to benchmark data MNIST (LeCun et al., 1998), CIFAR10 (Krizhevsky & Hinton, 2009) and CelebA (Liu et al., 2015) using ReLU ResNets without batch normalization and spectral normalization. The particle evolutions on MNIST and CIFAR10 without using outer loop are shown in Figure 3. Clearly, EPT can transport samples from a multivariate normal distribution into a target distribution.
We further compare EPT using the outer loop with the generative models including WGAN, SNGAN and MMDGAN. We considered different f -divergences, including Pearson’s χ2, KL, JS and logD (Gao et al., 2019)) and different deep density-ratio fitting methods (LSDR and LR). Table 1 shows FID (Heusel et al., 2017) evaluated with five bootstrap sampling of EPT with four divergences on CIFAR10. We can see that EPT using ReLU ResNets without batch normalization and spectral normalization attains (usually better) comparable FID scores with the state-of-the-art generative models. Comparisons of the real samples and learned samples on MNIST, CIFAR10 and CelebA are shown in Figure 4, where high-fidelity learned samples are comparable to real samples visually.

heading: 7 CONCLUSION
text: EPT is a new approach for generative learning via training a transport map that pushes forward a reference to the target. This approach uses the forward Euler method for solving the McKean-Vlasov equation, which results from linearizing the Monge-Ampère equation that characterizes the optimal transport map. The EPT map is a composition of a sequence of simple residual maps. The key task in training is the estimation of density ratios that completely determine the residual maps. We estimate density ratios based on the Bregman divergence with gradient penalty using deep density-ratio fitting. We establish bounds on the approximation errors due to linearization, discretization, and density-ratio estimation. Our results provide strong theoretical guarantees for the proposed method and ensure that the EPT map converges fast to the target. We also show that the proposed density-ratio (difference) estimators do not suffer from the “curse of dimensionality” if data is supported on a lower-dimensional manifold. This is an interesting result in itself since density-ratio estimation is an important basic problem in machine learning and statistics. Because EPT is easy to train, computationally stable, and
enjoys strong theoretical guarantees, we expect it to be a useful addition to the methods for generating learning.
The proposed EPT method is motivated from the Monge-Ampère equation that characterizes the optimal transport map. However, while the EPT map pushes forward a reference distribution to the target, it is not an estimate of the optimal transport map itself. How to consistently estimate the Monge-Ampére optimal map is a challenging and open problem.

heading: APPENDIX
text: In the appendix, we provide the implementation details on numerical settings, network structures, SGD optimizers, and hyper-parameters in the paper. We show the numerical convergence of EPT with simulated datasets and compare the learning and inference of EPT with other generative models. We give detailed theoretical background and proofs of the results mentioned in the paper. We also provide proofs MMD flow and SVGD can be derived from EPT by choosing appropriate f -divergences.

heading: A APPENDIX: NUMERICAL EXPERIMENTS
text: A.1 IMPLEMENTATION DETAILS, NETWORK STRUCTURES, HYPER-PARAMETERS
We provide the details of two versions of the EPT algorithm, EPTv1 in Algorithm 1 and EPTv2 in Algorithm 2 below. In Algorithm 1, we describe the algorithm without outer loops. In Algorithm 2, we describe the algorithm with a latent structure and outer loops.
A.1.1 2D EXAMPLES
Experiments on 2D examples in our work were performed with deep LSDR fitting and the Pearson χ2 divergence. We use the EPTv1 (Algorithm 1) without outer loops. In inner loops, only a multilayer perceptron (MLP) was utilized for dynamic estimation of the density ratio between the model distribution qk and the target distribution p. The network structure and hyper-parameters in EPT and deep LSDR fitting were shared in all 2D experiments. We adopt EPT to push particles from a predrawn pool consisting of 50k i.i.d. Gaussian particles to evolve in 20k steps. We used RMSProp with the learning rate 0.0005 and the batch size 1k as the SGD optimizer. The details are given in Table A1 and Table A2. We note that s is the step size, n is the number of particles, α is the penalty coefficient, and T is the mini-batch gradient descent times of deep LSDR fitting or deep logistic regression in each inner loop hereinafter.
Table A1: MLP for deep LSDR fitting.
Layer Details Output size
1 Linear, ReLU 64
2 Linear, ReLU 64
3 Linear, ReLU 64
4 Linear 1
Table A2: Hyper-parameters in EPT on 2D examples.
Parameter s n α T
Value 0.005 50k 0 or 0.5 5

heading: A.1.2 REAL IMAGE DATA
text: Datasets. We evaluated EPT on three benchmark datasets including two small datasets MNIST, CIFAR10 and one large dataset CelebA from GAN literature. MNIST contains a training set of 60k examples and a test set of 10k examples as 28× 28 bilevel images which were resized to 32× 32 resolution. There are a training set of 50k examples and a test set of 10k examples as 32× 32 color images in CIFAR10. We randomly divided the 200k celebrity images in CelebA into two sets for training and test according to the ratio 9:1. We also pre-processed CelebA images by first taking a 160× 160 central crop and then resizing to the 64× 64 resolution. Only the training sets are used to train our models.
Algorithm 1: EPTv1: Euler particle transport Input: K ∈ N∗, s > 0, α > 0 // maximum loop count, step size, regularization coeficient
Xi ∼ ν, Ỹ 0i ∼ µ, i = 1, 2, · · · , n // real samples, initial particles k ← 0 while k < K do
R̂kφ ∈ arg minRφ 1n ∑n i=1[Rφ(Xi)
2 + α‖∇Rφ(Xi)‖22 − 2Rφ(Ỹ ki )] via SGD // determine the density ratio v̂k(x) = −f ′′(R̂kφ(x))∇R̂kφ(x) // approximate the velocity field T̂ k = 1+ sv̂k // define the forward Euler map Ỹ k+1i = T̂ k(Ỹ ki ), i = 1, 2, · · · , n // update particles k ← k + 1 end Output: Ỹ ki ∼ µ̃k, i = 1, 2, · · · , n // transported particles
Evaluation metrics. Fréchet Inception Distance (FID) (Heusel et al., 2017) computes the Wasserstein distance W2 with summary statistics (mean µ and variance Σ) of real samples xs and generated samples gs in the feature space of the Inception-v3 model (Szegedy et al., 2016), i.e., FID = ‖µx−µg‖22 +Tr(Σx+Σg−2(ΣxΣg) 1 2 ). Here, FID is reported with the TensorFlow implementation and lower FID is better.
Network architectures and hyper-parameter settings. We employed the ResNet architectures used by Gao et al. (2019) in our EPT algorithm. Especially, the batch normalization (Ioffe & Szegedy, 2015) and the spectral normalization (Miyato et al., 2018) of networks were omitted for EPT-LSDRχ2. To train neural networks, we set SGD optimizers as RMSProp with the learning rate 0.0001 and the batch size 100. Inputs {Zi}ni=1 in EPTv2 (Algorithm 2) were vectors generated from a 128-dimensional standard normal distribution on all three datasets. Hyper-parameters are listed in Table A3 where IL expresses the number of inner loops in each outer loop. Even without outer loops, EPTv1 (Algorithm 1) can generate images on MNIST and CIFAR10 as well by making use of a large set of particles. Table A4 shows the hyper-parameters.
Table A3: Hyper-parameters in EPT with outer loops on real image datasets.
Parameter ` s n α T IL
Value 128 0.5 1k 0 1 20
Table A4: Hyper-parameters in EPT without outer loops on real image datasets.
Parameter s n α T
Value 0.5 4k 0 5

heading: A.2 NUMERICAL CONVERGENCE
text: We illustrate the convergence property of the learning dynamics of EPTv1 on synthetic datasets pinwheel, checkerboard and 2spirals. As shown in Figure 5, on the three test datasets, the dynamics of both the estimated LSDR fitting losses in (14) with α = 0 and the estimated value of the gradient norms EX∼qk [‖∇Rφ(X)‖2] demonstrate the estimated LSDR loss converges to the theoretical value −1.

heading: A.3 LEARNING AND INFERENCE
text: The learning process of EPT performs particle evolution via solving the McKean-Vlasov equation using forward Euler iterations. The iterations rely on the estimation of the density ratios (difference) between the pushforward distributions and the target distribution. To make the inference of EPTv1 more amendable, we propose EPTv2 based on EPTv1. EPTv2 takes advantage of a neural network to fit the pushforward map. The inference of EPTv2 is fast since the pushforward map is parameterized as a neural network and only forward propagation is involved. These aspects distinguish EPTv2 from score-based generative models Song & Ermon (2019; 2020) which simulate Langevin dynamics to generate samples.

heading: B APPENDIX: THEORETICAL BACKGROUNDS AND PROOFS

heading: B.1 THEORETICAL BACKGROUNDS FOR SECTION 2
text: For convenience, we first give the following notation to be used in this section. Let P2(Rm) denote the space of Borel probability measures on Rm with finite second moments, and let Pa2 (Rm) denote the subset of P2(Rm) in which measures are absolutely continuous with respect to the Lebesgue measure (all distributions are assumed to satisfy this assumption hereinafter). TanµP2(Rm) denotes the tangent space to P2(Rm) at µ. Let ACloc(R+,P2(Rm)) := {µt : I → P2(Rm) is absolutely continuous, |µ′t| ∈ L2(I), I ⊂ R+}. Liploc(Rm) denotes the set of functions that are Lipschitz continuous on any compact set of Rm. For any ` ∈ [1,∞], we use L`(µ,Rm) (L`loc(µ,Rm)) to denote the L` space of µ-measurable functions on Rm (on any compact set of Rm). With 1, det and tr, we refer to the identity map, the determinant and the trace. We use ∇, ∇2 and ∆ to denote the gradient or Jacobian operator, the Hessian operator and the Laplace operator, respectively.
We are now ready to describe the proposed method in a mathematically rigorous fashion and provide theoretical guarantees. Let X ∼ q, X̃ = Tt,Φ(X), and denote the distribution of X̃ as q̃. With a small t, the map Tt,φ is invertible according to the implicit function theorem, and we have the change of variables formula
det(∇2Ψ)(x) = |det(∇Tt,Φ)(x)| = q(x)
q̃(x̃) , (B-1)
where x̃ = Tt,Φ(x). (B-2)
Using the fact ddt ∣∣ t=0 det(A + tB) = det(A)tr ( A−1B ) ∀A,B ∈ Rm×m with A invertible, and applying the first order Taylor expansion to (B-1), we have
log q̃(x̃)− log q(x) = −t∆Φ(x) + o(t). (B-3)
Let t→ 0 in (B-2) and (B-3), we obtain a random process {xt} and its law qt satisfying dxt dt = ∇Φ(xt), with x0 ∼ q, (B-4)
d ln qt(xt)
dt = −∆Φ(xt), with q0 = q. (B-5)
Equations (B-4) and (B-5) resulting from linearization of the Monge-Ampère equation (2) can be interpreted as gradient flows in measure spaces (Ambrosio et al., 2008). And thanks to this connection, we can resort to solving a continuity equation characterized by a type of McKean-Vlasov equation, an ODE system that is easier to handle.
B.2 GRADIENT FLOWS IN Pa2 (Rm)
For µ ∈ Pa2 (Rm) with density q, let L[µ] = ∫ Rm F (q(x))dx : Pa2 (Rm)→ R+ ∪ {0} (B-6)
be an energy functional satisfying ν ∈ arg minL[·], where F (·) : R+ → R1 is a twice-differentiable convex function. Among the widely used metrics on Pa2 (Rm) in implicit generative learning, the following two are important examples of L[·] : (1) f -divergence given in (5) (Ali & Silvey, 1966); (2) Lebesgue norm of density difference:
‖µ− ν‖2L2(Rm) = ∫ Rm |q(x)− p(x)|2dx. (B-7)
Definition. We call {µt}t∈R+ ⊂ ACloc(R+,P2(Rm)) a gradient flow of the functional L[·], if {µt}t∈R+ ⊂ Pa2 (Rm) a.e., t ∈ R+ and the velocity vector field vt ∈ TanµtP2(Rm) satisfies vt ∈ −∂L[µt] a.e. t ∈ R+, where ∂L[·] is the subdifferential of L[·]. The gradient flow {µt}t∈R+ of L[·] enjoys the following nice properties.
Proposition B.1. (i) The following continuity equation holds in the sense of distributions.
∂ ∂t µt = −∇ · (µtvt) in R+ × Rm with µ0 = µ, (B-8)
(ii) Energy decay along the gradient flow: ddtL[µt] = −‖vt‖ 2 L2(µt,Rm) a.e. t ∈ R +. In addition, W2(µt, ν) = O(exp−λt), if L[µ] is λ-geodetically convex with λ > 0 1.
(iii) Conversely, if {µt}t is the solution of continuity equation (B-8) in (i) with vt(x) specified by (B-9) in (ii), then {µt}t is a gradient flow of L[·].
Remark B.1. In part (ii) of Proposition B.1, for general f -divergences, we assume the functional L to be λ-geodesically convex for the convergence of µt to the target ν in the quadratic Wasserstein distance. However, for the KL divergence, the convergence can be guaranteed if ν satisfies the logSobolev inequality(Otto & Villani, 2000). In addition, the distributions that are strongly log-concave outside a bounded region, but not necessarily log-concave inside the region satisfy the log-Sobolev inequality, see, for example, Holley & Stroock (1987). Here the functional L can even be nonconvex, an example includes the densities with double-well potential.
Remark B.2. Equation (8.48) in Proposition 8.4.6 of and Ambrosio et al. (2008) shows the connection (locally) of the velocity vt of the gradient flow µt and the optimal transport along µt, i.e., let T µt+h µt be the optimal transport from µt to µt+h, then T µt+h µt = I + hvt + o(h) in Lp. So locally, I + hvt approximates the optimal transport map from µt to µt+h on [t, t+ h] for a small h.
Proof. (i) The continuity equation (B-8) follows from the definition of the gradient flow directly, see, page 281 in (Ambrosio et al., 2008). (ii) The first equality follows from the chain rule and integration by part, see, Theorem 24.2 of Villani (2008). The second one on linear convergence follows from Theorem 24.7 of Villani (2008), where the assumption on λ in equation (24.6) is equivalent to the λ-geodetically convex assumption here. (iii) Similar to (i) see, page 281 in Ambrosio et al. (2008).
Theorem B.1. (i) Representation of the velocity fields: if the density qt of µt is differentiable, then
vt(x) = −∇F ′(qt(x)) µt-a.e. x ∈ Rm. (B-9)
(ii) If we let Φ be time-dependent in (B-4)–(B-5), i.e., Φt, then the linearized Monge-Ampère equations (B-4)–(B-5) are the same as the continuity equation (B-8) by taking Φt(x) = −F ′(qt(x)).
Proof. (i) Recall L[µ] is a functional on Pa2 (Rm). By the classical results in calculus of variation (Gelfand & Fomin, 2000),
∂L[q] ∂q (x) = d dt L[q + tg] |t=0= F ′(q(x)),
where ∂L[q]∂q denotes the first order of variation of L[·] at q, and q, g are the densities of µ and an arbitrary ξ ∈ Pa2 (Rm), respectively. Let
LF (z) = zF ′(z)− F (z) : R1 → R1.
Some algebra shows, ∇LF (q(x)) = q(x)∇F ′(q(x)).
Then, it follows from Theorem 10.4.6 in (Ambrosio et al., 2008) that
∇F ′(q(x)) = ∂oL(µ), 1 We say that L is λ-geodetically convex if there exists a constant λ > 0 such that for every µ1, µ2 ∈
Pa2 (Rm), there exists a constant speed geodestic γ : [0, 1]→ Pa2 (Rm) such that γ0 = µ1, γ1 = µ2 and
L(γs) ≤ (1− s)L(µ1) + sL(µ2)− λ
2 s(1− s)d(µ1, µ2), ∀s ∈ [0, 1],
where d is a metric defined on Pa2 (Rm) such as the quadratic Wasserstein distance.
where, ∂oL(µ) denotes the one in ∂L(µ) with minimum length. The above display and the definition of gradient flow implies the representation of the velocity fields vt.
(ii) The time dependent form of (B-4)-(B-5) reads dxt dt = ∇Φt(xt), with x0 ∼ q,
d ln qt(xt)
dt = −∆Φt(xt), with q0 = q.
By chain rule and substituting the first equation into the second one, we have 1
qt ( dqt dt + dqt dxt dxt dt ) = 1 qt ( dqt dt +∇qt∇Φt(xt))
= −∆Φt(xt), which implies,
dqt dt = −qt∆Φt(xt)−∇qt∇Φt(xt) = −∇ · (qt∇Φt).
By (B-9), the above display coincides with the continuity equation (B-8) with vt = ∇Φt = −∇F ′(qt(x)).
Theorem B.1 and Proposition B.1 imply that {µt}t, the solution of the continuity equation (B-8) with vt(x) = −∇F ′(qt(x)), converges rapidly to the target distribution ν. Furthermore, the continuity equation has the following representation under mild regularity conditions on the velocity fields. Theorem B.2. Assume ‖vt‖L1(µt,Rm) ∈ L1loc(R+) and vt(·) ∈ Liploc(Rm) with upper bound Bt and Lipschitz constant Lt such that (Bt + Lt) ∈ L1loc(R+). Then the solution of the continuity equation (B-8) can be represented as µt = (Xt)#µ, where Xt(x) : R+ × Rm → Rm satisfies the McKean-Vlasov equation (4).
Proof. The Lipschitz assumption of vt implies the existence and uniqueness of the McKean-Vlasov equation (4) according to the classical results in ODE (Arnold, 2012). By the uniqueness of the continuity equation, see Proposition 8.1.7 in Ambrosio et al. (2008), it is sufficient to show that µt = (Xt)#µ satisfies the continuity equation (B-8) in a weak sense. This can be done by the standard test function and smoothing approximation arguments, see, Theorem 4.4 in Santambrogio (2015) for details.
As shown in Lemma B.1 below, the velocity fields associated with the f -divergence (5) and the Lebesgue norm (B-7) are determined by density ratio and density difference respectively. Lemma B.1. The velocity fields vt satisfy
vt(x) = { −f ′′(rt(x))∇rt(x), L[µ] = Df (µ‖ν), where rt(x) = qt(x)p(x) , −2∇dt(x), L[µ] = ‖µ− ν‖2L2(Rm), where dt(x) = qt(x)− p(x).
Proof. By definition,
F (qt(x)) = { p(x)f( qt(x)p(x) ), L[µ] = Df (µ‖ν), (qt(x)− p(x))2, L[µ] = ‖µ− ν‖2L2(Rm).
Direct calculation shows
F ′(qt(x)) = { f ′( qt(x)p(x) ), L[µ] = Df (µ‖ν), 2(qt(x)− p(x)), L[µ] = ‖µ− ν‖2L2(Rm).
Then, the desired result follows from the above display and (B-9).
Several methods have been developed to estimate density ratio and density difference in the literature. Examples include probabilistic classification approaches, moment matching and direct density-ratio (difference) fitting, see Sugiyama et al. (2012a;b); Kanamori & Sugiyama (2014); Mohamed & Lakshminarayanan (2016) and the references therein.
Proposition B.2. Suppose that the velocity fields vt are Lipschitz continuous with respect to (x, µt), that is, there exists a finite constant Lv > 0 such that
‖vt(x)− vt̃(x̃)‖ ≤ Lv[‖x− x̃‖+W2(µt, µt̃)], t, t̃ > 0 and x, x̃ ∈ Rm. (B-10)
Then for any finite T > 0, the bound (11) on the discretization error holds:
sup t∈[0,T ]
W2(µt, µst ) = O(s).
Remark B.3. If we take f(x) = (x− 1)2/2 in Lemma B.1, then the velocity fields vt(x) = ∇rt(x), where rt(x) = qt(x)/p(x). In the proof of Theorem B.1, part (ii), it is shown that qt satisfies dqt/dt = −∇ · (qt∇Φt). Thus for this simple f -divergence function, the verification of the Lipschitz condition (B-10) amounts to verifying that∇rt(x) is Lipschitz in the sense of (B-10).
Proof. Without loss of generality let K = Ts > 1 be an integer. Recall {µ s t t ∈ [ks, (k + 1)s) is the piecewise constant interpolation between µk and µk+1 defined as
µst = (T k,s t )#µk,
where, T k,st = 1+ (t− ks)vk,
µk is defined in (16)-(18) with vk = vks, i.e., the continuous velocity in (B-9) at time ks, k = 0, ..,K − 1, µ0 = µ. Under assumption (B-10) we can first show in a way similar to the proof of Lemma 10 in Arbel et al. (2019) that
W2(µks, µk) = O(s). (B-11)
Let Γ be the optimal coupling between µk and µks, and (X,Y ) ∼ Γ. Let Xt = T k,st (X) and Yt be the solution of (4) with X0 = Y and t ∈ [ks, (k + 1)s). Then
Xt ∼ µst , Yt ∼ µt and
Yt = Y + ∫ t ks vt̃(Yt̃)dt̃.
It follows that
W22 (µt, µks) ≤ E[‖Yt − Y ‖22] (B-12) = E[‖ ∫ t ks vt̃(Yt̃)dt̃‖22]
≤ E[( ∫ t ks ‖vt̃(Yt̃)‖2dt̃)2] ≤ O(s2).
where, the first inequality follows from the definition ofW2, and the last equality follows from the the uniform bounded assumption of vt. Similarly,
W22 (µk, µst ) ≤ E[‖X −Xt‖22] = E[‖(t− ks)vk(X)‖22] ≤ O(s2). (B-13)
Then,
W2(µt, µst ) ≤ W2(µt, µks) +W2(µks, µk) +W2(µk, µst ) ≤ O(s),
where the first inequality follows from the triangle inequality, see for example Lemma 5.3 in Santambrogio (2015), and the second one follows from (B-11)-(B-13).

heading: B.3 DERIVATION AND PROOFS OF THE RESULTS IN SECTION 4.

heading: B.3.1 BREGMAN SCORE FOR DENSITY RATIO/DIFFERENCE
text: The separable Bregman score with the base probability measure p to measure the discrepancy between a measurable function R : Rm → R1 and the density ratio r is
Bratio(r,R) = EX∼p[g′(R(X))(R(X)− r(X))− g(R(X))] = EX∼p[g′(R(X))R(X)− g(R(X))]− EX∼q[g′(R(X))].
It can be verified that Bratio(r,R) ≥ Bratio(r, r), where the equality holds iff R = r. For deep density-difference fitting, a neural network D : Rm → R1 is utilized to estimate the density-difference d(x) = q(x)− p(x) between a given density q and the target p. The separable Bregman score with the base probability measure w to measure the discrepancy between D and d can be derived similarly,
Bdiff(d,D) = EX∼p[w(X)g′(D(X))]− EX∼q[w(X)g′(D(X))] + EX∼w[g′(D(X))D(X)− g(D(X))].
Here, we focus on the widely used least-squares density-ratio (LSDR) fitting with g(c) = (c− 1)2 as a working example for estimating the density ratio r. The LSDR loss function is
BLSDR(r,R) = EX∼p[R(X)2]− 2EX∼q[R(X)] + 1.

heading: B.3.2 GRADIENT PENALTY
text: We consider a noise convolution form of Bratio(r,R) with Gaussian noise ∼ N (0, αI), Bαratio(r,R) = EX∼pE [g′(R(X + ))R(X + )− g(R(X + ))]− EX∼qE [g′(R(X + ))].
Taylor expansion applied to R gives
E [R(x+ )] = R(x) + α
2 ∆R(x) +O(α2).
Using equations (13)-(17) in Roth et al. (2017), we get
Bαratio(r,R) ≈ Bratio(r,R) + α
2 Ep[g′′(R)‖∇R‖22],
i.e., 12Ep[g ′′(R)‖∇R‖22] serves as a regularizer for deep density-ratio fitting when g is twice differentiable.

heading: B.3.3 PROOFS IN SECTION 4
text: Below we prove Theorem 4.1 in Section 4. Lemma B.2. For given densities p(x) and q(x), let r(x) = q(x)/p(x) with C = EX∼q[r(X)]−1 < ∞. For any α ≥ 0, define a nonnegative functional
BαLSDR(R) = BLSDR(r,R) + αEp[‖∇R‖22] + C, R is measurable. Then, r ∈ arg minRB0LSDR(R). Moreover, Bα(R) = 0 iff R(x) = r(x) = 1, (q, p)-a.e. x ∈ Rm.
Proof. By definition, it is easy to check B0LSDR(R) = Bratio(r,R)−Bratio(r, r),
where Bratio(r,R) is the Bregman score with the base probability measure p between R and r. Then r ∈ arg minmeasureableRB0LSDR(R) follow from the fact Bratio(r,R) ≥ Bratio(r, r) and the equality holds iff R = r. Since Bα(R) = B0LSDR(R) + αEp[‖∇R‖22] ≥ 0, Then,
Bα(R) = 0
iff B0LSDR(R) = 0 and Ep[‖∇R‖22] = 0, which is further equivalent to R = r = constant (q, p)-a.e. ,
and the constant = 1 since r is a density ratio.
Theorem B.3. Assume supp(r) = M and r(x) is Lipschitz continuous with the bound B and the Lipschitz constant L. Suppose the topological parameter ofHD,W,S,B in (20) with α = 0 satisfies D = O(log n),W = O(n M 2(2+M) / log n), S = O(n M−2 M+2 / log4 n), and B = 2B. Then,
E{Xi,Yi}n1 [‖R̂φ − r‖ 2 L2(ν)] ≤ C(B 2 + cLmM)n−2/(2+M),
where C is a universal constant.
Proof. We use B(R) to denote B0LSDR − C for simplicity, i.e.,
B(R) = EX∼p[R(X)2]− 2EX∼q[R(X)]. (B-14)
Rewrite (20) with α = 0 as
R̂φ ∈ arg min Rφ∈HD,W,S,B B̂(Rφ) = n∑ i=1 1 n (Rφ(Xi) 2 − 2Rφ(Yi)). (B-15)
By Lemma B.2 and Fermat’s rule (Clarke, 1990), we know 0 ∈ ∂B(r). Then, ∀R direct calculation yields,
‖R− r‖2L2(ν) = B(R)−B(r)− 〈∂B(r), R− r〉 = B(R)−B(r). (B-16)
∀R̄φ ∈ HD,W,S,B we have,
‖R̂φ − r‖2L2(ν) = B(R̂φ)−B(r) (B-17)
= B(R̂φ)− B̂(R̂φ) + B̂(R̂φ)− B̂(R̄φ) + B̂(R̄φ)−B(R̄φ) + B(R̄φ)−B(r) ≤ 2 sup
R∈HD,W,S,B |B(R)− B̂(R)|+ ‖R̄φ − r‖2L2(ν),
where the inequality uses the definition of R̂φ, R̄φ and (B-16). We prove the theorem by upper bounding the expected value of the right hand side term in (B-17). To this end, we need the following auxiliary results (B-18)-(B-20).
E{Zi}ni [sup R |B(R)− B̂(R)|] ≤ 4C1(2B + 1)G(H), (B-18)
where
G(H) = E{Zi, i}ni
[ sup
R∈HD,W,S,B | 1 n n∑ i=1 iR(Zi)| ] is the Gaussian complexity ofHD,W,S,B (Bartlett & Mendelson, 2002).
Proof of (B-18). Let g(c) = c2 − c, z = (x,y) ∈ Rm × Rm,
R̃(z) = (g ◦R)(z) = R2(x)−R(y).
Denote Z = (X,Y ), Zi = (Xi, Yi), i = 1, ..., n with X,Xi i.i.d. ∼ p, Y, Yi i.i.d. ∼ q. Let Z̃i be an i.i.d. copy of Zi, and σi( i) be i.i.d. Rademacher random (standard normal) variables that are independent of Zi and Z̃i. Then,
B(R) = EZ [R̃(Z)] = 1
n EZ̃i [R̃(Z̃i)],
and
B̂(R) = 1
n n∑ i=1 R̃(Zi).
Denote
R(H) = 1 n E{Zi,σi}ni [ sup R∈HD,W,S,B | n∑ i=1 σiR(Zi)|]
as the Rademacher complexity ofHD,W,S,B (Bartlett & Mendelson, 2002). Then,
E{Zi}ni [sup R |B(R)− B̂(R)|] = 1 n E{Zi}ni [sup R | n∑ i=1 (EZ̃i [R̃(Z̃i)]− R̃(Zi))|]
≤ 1 n E{Zi,Z̃i}ni [supR |R̃(Z̃i)− R̃(Zi)|]
= 1
n E{Zi,Z̃i,σi}ni [supR | n∑ i=1 σi(R̃(Z̃i)− R̃(Zi))|]
≤ 1 n E{Zi,σi}ni [sup R | n∑ i=1 σiR̃(Zi)|] + 1 n E{Z̃i,σi}ni [supR | n∑ i=1 σiR̃(Z̃i)|]
= 2R(g ◦ H) ≤ 4(2B + 1)R(H) ≤ 4C1(2B + 1)G(H),
where, the first inequality follows from the Jensen’s inequality, and the second equality holds since the distribution of σi(R̃(Z̃i)− R̃(Zi)) and R̃(Z̃i)− R̃(Zi) are the same, and the last equality holds since the distribution of the two terms are the same, and last two inequality follows from the Lipschitz contraction property where the Lipschitz constant of g onHD,W,S,B is bounded by 2B + 1 and the relationship between the Gaussian complexity and the Rademacher complexity, see for Theorem 12 and Lemma 4 in Bartlett & Mendelson (2002), respectively.
G(H) ≤ C2B √
n
DS logS log
n DS logS exp(− log2 n DS logS ). (B-19)
Proof of (B-19). SinceH is negation closed,
G(H) = E{Zi, i}ni [ sup R∈HD,W,S,B
1
n n∑ i=1 iR(Zi)]
= EZi [E i [ sup R∈HD,W,S,B
1
n n∑ i=1 iR(Zi)]|{Zi}ni=1].
Conditioning on {Zi}ni=1, ∀R, R̃ ∈ HD,W,S,B it easy to check
V i [ 1
n n∑ i=1 i(R(Zi)− R̃(Zi))] = dH2 (R, R̃)√ n ,
where, dH2 (R, R̃) = 1√ n √∑n i=1(R(Zi)− R̃(Zi))2. Observing the diameter ofHD,W,S,B under dH2 is at most B, we have
G(H) ≤ C3√ n E{Zi}ni=1 [ ∫ B 0 √ logN (H, dH2 , δ)dδ]
≤ C3√ n E{Zi}ni=1 [ ∫ B 0 √ logN (H, dH∞, δ)dδ]
≤ C3√ n ∫ B 0 √ VCH log 6Bn δVCH dδ, ≤ C4B( n
VCH )1/2 log(
n VCH ) exp(− log2( n VCH ))
≤ C2B √
n
DS logS log
n DS logS exp(− log2 n DS logS )
where, the first inequality follows from the chaining Theorem 8.1.3 in Vershynin (2018), and the second inequality holds due to dH2 ≤ dH∞, and in the third inequality we used the relationship between
the matric entropy and the VC-dimension of the ReLU networks HD,W,S,B (Anthony & Bartlett, 2009), i.e.,
logN (H, dH∞, δ) ≤ VCH log 6Bn δVCH ,
and the fourth inequality follows by some calculation, and the last inequality holds due to the upper bound of VC-dimension for the ReLU networkHD,W,S,B satisfying
VCH ≤ C5DS logS,
see Bartlett et al. (2019). For any two integer M,N , there exists a R̄φ ∈ HD,W,S,B with width W = max{8MN1/M + 4M, 12N + 14} and depth D = 9M + 12, and B = 2B, such that
‖r − R̄φ‖2L2(ν) ≤ C6cLmM(NM) −4/M. (B-20)
Proof of (B-20). We use Lemma 4.1, Theorem 4.3, 4.4 and following the proof of Theorem 1.3 in Shen et al. (2019). Let A be the random orthoprojector in Theorem 4.4, then it is to check A(M ) ⊂ A([−c, c]m) ⊂ [−c √ m, √ mc]M. Let r̃ be an extension of the restriction of r on M , which is defined similarly as g̃ on page 30 in Shen et al. (2019). Since we assume the target r is Lipschitz continuous with the bound B and the Lipschitz constant L, let small enough, then by Theorem 4.3, there exist a ReLU network R̃φ ∈ HD,W,S,B with width
W = max{8MN1/M + 4M, 12N + 14},
and depth D = 9M + 12,
and B = 2B, such that
‖r̃ − R̃φ‖L∞(M \N ) ≤ 80cL √ mM(NM)−2/m,
and ‖R̃φ‖L∞(M ) ≤ B + 3Lc √ mM,
where, N is a ν− negligible set with ν(N ) can be arbitrary small. Define R̄φ = R̃φ ◦A. Then, following the proof after equation (4.8) in Theorem 1.3 of Shen et al. (2019), we get our (B-20) and
‖R̄φ‖L∞(M \N ) ≤ 2B, ‖R̄φ‖L∞(N ) ≤ 2B + 3cL √ mM.
Let DS logS < n, combing the results (B-17) - (B-20), we have
E{Xi,Yi}n1 [‖R̂φ − r‖ 2 L2(ν)]
≤ 8C1(2B + 1)G(H) + C6cLmM(NM)−4/M ≤ 8C1(2B + 1)C2B √ DS logS
n log
n
DS logS + C6cLmM(NM)−4/M
≤ C(B2 + cLmM)n−2/(2+M),
where, last inequality holds since we choose M = log n, N = n M 2(2+M) / log n, S = n M−2 M+2 / log4 n, i.e., D = 9 log n+ 12,W = 12n M 2(2+M) / log n+ 14.

heading: B.4 THE RELATIONSHIP BETWEEN EPT AND MMD FLOW
text: Here we show that MMD flow can be considered a special case of EPT.
Proof. Let H be a reproducing kernel Hilbert space with characteristic kernel K(x, z). Recall in MMD flow,
L[µ] = 1 2 ‖µ− ν‖2mmd,
and ∂L[µ] ∂µ (x) =
∫ K(x, z)dµ(z)− ∫ K(x, z)dν(z),
and the vector fields
vmmdt = −∇ ∂L[µ] ∂µt
= ∫ ∇xK(x, z)dν(z)− ∫ ∇xK(x, z)dµt(z)
= ∫ ∇xK(x, z)p(z)dz − ∫ ∇xK(x, z)qt(z)dz
By Lemma B.1, the vector fields corresponding the Lebesgue norm 12‖µ−ν‖ 2 L2(Rm) = 1 2 ∫ Rm |q(x)− p(x)|2dx are defined as vt = ∇p(x)−∇qt(x).
Next, we will show the vector fields vmmdt is exactly by projecting the vector fields vt on to the reproducing kernel Hilbert spaceHm = H⊗m. By the definition of reproducing kernel we have,
p(x) = 〈p(·),K(x, ·)〉H = ∫ K(x, z)p(z)dz,
and qt(x) = 〈qt(·),K(x, ·)〉H = ∫ K(x, z)qt(z)dz.
Hence,
vt(x) = ∇p(x)−∇qt(x)
= ∫ ∇xK(x, z)(p(z)− qt(z))dz
=vmmdt (x).
This completes the proof.

heading: B.5 PROOF OF THE RELATION BETWEEN EPT AND SVGD
text: Proof. Let f(u) = u log u in (5). With this f the velocity fields vt = −f ′′(rt)∇rt = −∇rt(x)rt(x) Let g in a Stein class associated with qt.
〈vt,g〉H(qt) =− ∫
g(x)T ∇rt(x) rt(x) qt(x)dx
=− ∫ g(x)T∇ log rt(x)qt(x)dx
=− EX∼qt(x)[g(x) T∇ log qt(X) + g(x)T∇ log p(X)]
=− EX∼qt(x)[g(x) T∇ log qt(X) +∇ · g(x)]
+ EX∼qt(x)[g(x) T∇ log p(X) +∇ · g(x)]
=− EX∼qt(x)[Tqtg] + EX∼qt(x)[Tpg] =EX∼qt(x)[Tpg],
where the last equality is obtained by restricting g in a Stein class associated with qt, i.e., EX∼qt(x)Tqtg = 0. This is exactly the velocity fields of SVGD (Liu, 2017).

Review summaries:
This paper considers generative learning by discretizing a Wasserstein gradient with Euler methods. More precisely, some samples of a target distribution are given and the goal is to pushforward some samples of an initial distribution to the target distribution. The proposed method is obtained by  minimizing the f-divergence between the initial distribution and the target distribution, but considering the Wasserstein gradient flow of the f-divergence w.r.t. the target distribution (= the objective function). This Wasserstein gradient flow is discretized via Euler method to obtain the proposed algorithm. This Euler method involves the Wasserstein gradient of the objective, which is intractable. The authors describe a statistical methodology to compute this Wasserstein gradient based on the samples of the target distribution and samples from the current distribution. They also prove a bound for the estimated Wasserstein gradient wrt the true Wasserstein gradient. Finally, the paper presents relevant numerical experiments.

This paper tackles generative modeling (sampling, in particular) via finding the push forward functions T (equivalently, the velocity fields v) that iteratively moves particles from a reference distribution toward the target data distribution. The velocity fields are solved by minimizing the f-divergence between the particle density at iteration k and the target data density, which is shown to be in the form of gradient of density ratio. Based on this intuition, the training stage becomes estimating the density ratio via neural networks, for each iteration k=1,...,K. However, estimating the density ratio can be quite difficult when two densities have little overlapping support. Thus, the author proposes to add gradient regularizer to the density ratio estimating function. The experiment on real-world computer vision benchmarks demonstrate reasonable sampling quality, and the FID score on CIFAR-10 is comparable to some GAN baselines in the generative modeling literature.