Title: AE-OT: A NEW GENERATIVE MODEL BASED ON EX-TENDED SEMI-DISCRETE OPTIMAL TRANSPORT

Abstract: Current generative models like generative adversarial networks (GANs) and variational autoencoders (VAEs) have attracted huge attention due to its capability to generate visual realistic images. However, most of the existing models suffer from the mode collapse or mode mixture problems. In this work, we give a theoretic explanation of the both problems by Figalli's regularity theory of optimal transportation maps. Basically, the generator compute the transportation maps between the white noise distributions and the data distributions, which are in general discontinuous. However, deep neural networks (DNNs) can only represent continuous maps. This intrinsic conflict induces mode collapse and mode mixture. In order to tackle the both problems, we explicitly separate the manifold embedding and the optimal transportation; the first part is carried out using an autoencoder (AE) to map the images onto the latent space; the second part is accomplished using a GPU-based convex optimization to find the discontinuous transportation maps. Composing the extended optimal transport (OT) map and the decoder, we can finally generate new images from the white noise. This AE-OT model avoids representing discontinuous maps by DNNs, therefore effectively prevents mode collapse and mode mixture.

1. INTRODUCTION:

Generative adversarial networks (GANs) (Goodfellow et al. (2014) ) and variational autoencoders (VAEs) (Kingma & Welling (2013) ) emerge as the dominant approaches for unconditional image generation. When trained on appropriate datasets, they are able to produce realistic and visual appealing samples. GAN methods train an unconditional generator that regresses real images from random noise and a discriminator that measures the difference between generated samples and real images. Despite GANs' advantages, they have critical drawbacks. 1) Training of GANs are tricky and sensitive to hyperparameters. 2) GANs suffer from mode collapse, in which the generator only learns to generate few modes of data distribution while missing others, although samples from the missing modes occur throughout the training data (see e.g. Goodfellow (2016) ). While for the VAEs, the encoder is used to map the data distribution to a Gaussian latent distribution, which is then mapped back to the data distribution by the decoder. While standard VAEs tend to capture all modes, they often generate ambiguous images on multimodal real data distributions. We propose that these phenomena relates deeply with the singularities of distribution transport maps.

Manifold Distribution Hypothesis:

In deep learning, the manifold distribution hypothesis is well accepted, which assumes the distribution of a specific class of natural data is concentrated on a low dimensional manifold embedded in the high dimensional data space Tenenbaum et al. (2000) . Therefore, GANs and VAEs implicitly aim to accomplish two major tasks: 1) manifold embedding: to find the encoding/decoding maps between the data manifold embedded in the image space and the latent space; 2) probability distribution transport: to transport a given white noise distribution to the data distribution, either in the latent or in the image space. Top row shows real data distributions, the bottom row gives the noise distributions. On top, each cluster represents a mode, the spurious generated samples are red crosses (mode mixture); at the bottom, red dotted lines are the singularity set, red crosses are mapped to be spurious samples by DNNs. (a) If the support of the target distribution is convex, DNN (f 1 ) is able to approximate the transport map f1 well. When the support of the target distributions are concave, there are two situations: (b) single mode and (c) multi modes. In (b), DNN, represented by f 2 , cannot approximate the transport map f2 well and generates some spurious samples. f3 gives the transport map of multi-mode, when approximating it with continuous DNNs, either mode collapse f 31 or mode mixture f 32 will happen.

Distribution Transformation:

The generator of GAN model and the decoder of VAE model are trained to compute a transport map that transforms a known continuous distribution (e.g. Gaussian white noise) to the real data distribution. Namely, the transport map pushes forward the given noise distribution to a generated distribution to approximate the real data distribution, the similarity between the two distributions determines the generalization ability of the generator Ben-David et al. (2010) .

Discontinuity and Mode Collapse/Mixture It is a common practice among GAN/VAE models that the generators/decoders are expressed by deep neural networks, which can only represent continuous mappings. Unfortunately, as pointed out by works Nagarajan & Kolter (2017); Khayatkhoei et al. (2018) ; Xiao et al. (2018) , the transport maps may be discontinuous when there are multiple modes in the data distribution. This intrinsic conflict can cause mode collapse or mode mixture. The later means that the generated samples are mixtures of multiple modes and look spurious or ambiguous.1 

Furthermore, even when the real data distribution has a single mode, ambiguous data (e.g. a human face image with mismatched eye colors) can still present. This can be explained by Brenier's polar factorization theorem Brenier (1991b; 1987; 1991a ) and Figalli's regularity theorem Figalli (2010); Chen & Figalli (2017) (Thm. 5 in Appendix B), which asserts that if the support of the target distribution is not convex, then there will be singularity sets on the support of the source distribution, such that the transport map is discontinuous on these sets. This shows the intrinsic difficulties of conventional GANs/VAEs cannot be eliminated, as shown in Fig. 1 .

Conquering Mode Collapse/Mixture However, according to Brenier (1987; 1991a) theorem, the optimal transport map can be represented as the gradient map of the Brenier potential. At the regular points, the Brenier potential is differentiable, its gradient map (the transport map) is continuous; at the singularities, the Brenier potential is continuous but not differentiable, and its gradient map is discontinuous. Conventional GANs and VAEs model the gradient map directly and encounter the trouble of discontinuity. In contrast, we propose to model the globally continuous Brenier potential to avoid mode collapse/mixture.

More specifically, our proposed AE-OT model separates the manifold embedding step and the probability distribution transformation step, the former is carried out by an autoencoder (AE), the latter is accomplished by a convex optimization framework (OT). The OT step computes the Brenier potential explicitly and is able to locate the singularity set (the discontinuous points of the gradient map) based on Figalli's theory. Our experimental results demonstrate that the proposed method can not only cover all of the modes, but also avoid generating spurious samples (mode mixture).

Contributions (i) From theoretical aspect, this work gives a thorough explanation of mode collapse and mode mixture by the regularity theory of optimal transportation developed by Figalli ( 2018Fields medalist) and the reasons why standard GANs/VAEs cannot solve this problem perfectly.

(ii) From practical aspect, this work proposes a novel model called AE-OT, which first encodes the data manifold into the latent space, then compute the Brenier potential to represent the optimal transportation map in the latent space. The Figalli's singularity set can be located efficiently and avoided when generating new samples. In this way, our model eliminates mode collapse and mode mixture successfully. (iii) The algorithm for finding the Brenier potential and the optimal transportation map can be accelerated with GPU based convex optimization algorithm. The method converges to the unique global optimum with bounded error estimate. (iv) Our experiment results demonstrate the efficiency and efficacy of the proposed method.

2. RELATED WORK:

Optimal Transport Optimal transport (OT) has been successfully applied in the areas like economy, optics and machine learning, as surveyed in Solomon (2018) and Peyré & Cuturi (2018) . In Gu et al. (2016) , the intrinsic connection between Brenier theory in OT and Alexandroff theory in convex geometry was established, and applied for deep learning in Lei et al. (2017) by an convex optimization. Figalli and the collaborators Figalli (2010); Chen & Figalli (2017) proposed that when the support of the data distribution is non-convex, the transport map will be discontinuous. Besides, various non-adversarial methods has also been proposed recently. GLO Bojanowski et al. (2017) adopts the "encoder-less autoencoder" method to generate new images with a non-adversarial loss function. IMLE Li & Malik (2018) proposed an ICP related generative model training approach.

Generative models:

Later GLANN Hoshen & Malik (2019) combines the advantages of GLO and IMLE, an embedding from image space to latent space was first found using GLO and then a transformation between an arbitrary distribution and latent code was computed using IMLE.

Mitigating Mode Collapsing Recently, Nagarajan & Kolter (2017); Khayatkhoei et al. (2018) ; Xiao et al. (2018) also realize the training difficulties of GANs come from the approximation of discontinuous functions with continuous DNNs. By the gradient-based regularization, GDGAN Nagarajan & Kolter (2017) do relieve the mode collapse phenomenon of GANs, but mode mixture still exists. Khayatkhoei et al. (2018) proposes to use multiple GANs to overcome the mode collapse. Xiao et al. (2018) proposed to embed the images into a latent space according to Bourgain's theorem, and train the generator by sampling a Gaussian mixture distribution in the latent space instead of a unimodal Gaussian. The recently introduced normalized diversification by Liu et al. (2018) can also help overcome mode collapse successfully. However, all of them cannot solve the mode mixture well.

All these works Nagarajan & Kolter (2017); Khayatkhoei et al. (2018) ; Xiao et al. (2018) explain that if the target data distribution has multiple modes, the transport map is discontinuous, but DNNs can only represent continuous mappings, the intrinsic conflict causes mode collapse.

3. COMPUTATIONAL ALGORITHMS:

Overview of AE-OT Model Our AE-OT model is summarized in Fig. 2 , it has two major components: i) (AE) An autoencoder is trained to encode (f θ ) the data manifold from the image space X to the latent space Z, and map the data distribution to the latent code distribution; then the decoder g ξ decodes the latent code back to the data manifold. ii) (OT) This module computes the optimal transportation map T from the noise distribution to the latent code distribution. First, the Brenier potential is found by a convex optimization process according to Gu et al. (2016) , whose gradient is the semi-discrete optimal transport map, where the target is the discrete set of latent codes of training samples; then the transport map is piece-wise linearly extended to a global continuous map T , where the image domain becomes a simplicial complex obtained by triangulating the above latent codes. Finally, the singularity set in the source domain is located and avoided when generating new samples. As a result, given a random noise x, we can get the generated image by g ξ • T (x).

Semi-Discrete OT Map Suppose the source measure µ (Gaussian or uniform distribution) is absolutely continuous defined on a convex domain

Ω ⊂ R d , the target domain is a discrete set, Y = {y 1 , y 2 , • • • , y n }, y i ∈ R d , the target measure is a Dirac measure, ν = n i=1 ν i δ(y -y i ), i = 1, 2, .

. . , n, with the equal total mass as the source measure, µ(Ω)

= n i=1 ν i . Under a semi- discrete transport map T : Ω → Y , a cell decomposition is induced Ω = n i=1 W i , such that every x in each cell W i is mapped to the target y i , T : x ∈ W i → y i . The map T is measure preserving, denoted as T # µ = ν, if the µ-volume of each cell W i equals to the ν-measure of the image T (W i ) = y i , µ(W i ) = ν i .

The cost function is given by c : Ω × Y → R, where c(x, y) represent the cost for transporting a unit mass from x to y. The total cost of T is given by Ω c(x, T (x))dµ(x) = n i=1 Wi c(x, y i )dµ(x). Semi-discrete optimal transport map is the measure-preserving map that minimizes the total cost, T * := arg min T # µ=ν Ω c(x, T (x))dµ(x).

When the cost function is the L 2 distance c(x, y) = 1/2 x -y 2 , Brenier's theorem claims that the semi-discrete OT map is given by the gradient map of a piece-wise (PL) convex function, the so-called Brenier potential

u h : Ω → R, u h (x) := max n i=1 {π h,i (x)}, where π h,i (x) = x, y i + h i is the hyperplane corresponding to y i ∈ Y . As shown in Fig. 3(a), the projection of the graph of u h decomposes Ω into cells W i (h), each cell W i (h) is the projection of the supporting plane π h,i (x).

The height vector h is the unique optimizer of the following convex energy under the condition that

EQUATION

where w i (η) is the µ-volume of W i (η). The convex energy E(h) can be optimized simply by gradient descend method with

∇E(h) = (w i (h) -ν i ) T .

The key is to compute the µ-volume w i (h) of each cell W i (h), which can be estimated using conventional Monte Carlo method. We draw N random samples from µ distribution,

{x j } ∼ i.i.d. µ, ∀j ∈ J , the estimated µ-volume of each cell is ŵi (h) = #{j ∈ J | x j ∈ W i (h)}/N . Given x j , we can find W i in which x j ∈ W i by i = arg max i { x j , y i + h i }, i = 1, 2, . . . , n. When N is large enough, ŵi (h) converges to w i (h).

Then the gradient of the energy is approximated as

∇E ≈ ( ŵi (h) -ν i ) T .

Once the gradient is estimated, we can use Adam algorithm Kingma & Ba (2015) to minimize the energy. Sampling of x is independent of each other and finding the cell that x is located only involves matrix multiplication and sorting. Hence the Monte Carlo method has a natural parallel computation implementation on GPUs.

The approximation error is proportional to the inverse of the square root the amount of Monte Carlo samples. Asymptotically the number of Monte Carlo samples increases exponentially with respect to the dimension d (see e.g. Weed & Bach (2017) ). This brings huge computational burdens. To find a good balance between precision and speed, we adaptively adjust the number of random samples. In practice, we apply the following strategy: if the energy E(h) ceases decreasing for a number of consecutive steps, we double the amount of Monte Carlo samples. The algorithmic details of semi-discrete OT map are summarized in Alg. 1. Generate N uniformly distributed samples {xj } N j=1 .

6:

Calculate ∇h = ( ŵi(h) -νi) T .

7:

∇h = ∇h -mean(∇h).

8::

Update h by Adam algorithm with β1 = 0.9, β2 = 0.5.

9::

if E(h) has not decreased for s steps then 10:

N ← N × 2.

11:

end if 12: until Converge 13: OT map T (•) ← ∇(maxi •, yi + hi).

Algorithm 2 Generate latent code 1: Input: Optimal transport map T (•), number of samples to generate n, angle threshold θ. 2: Output: Generated latent code P . 

6::

Compute dihedral angles θi k between πi 0 and πi k .

7::

Select θi k with θi k ≤ θ, result in îk = 0, 1, . . . , d1.

8: if ∀k, θi k > θ then Abandon x 9: else Generate latent code T (x) = d 1 k=0 λ k T (ĉ îk ) with λ k = d -1 (x, ĉî k )/ d 1 j=0 d -1 (x, ĉî j ).

10:

end if 11: until Generate n new latent code

Piece-wise Linear Extension The semi-discrete OT map ∇u h : Ω → Y maps all x ∈ Ω to the latent codes of training samples {y i }'s and won't generate new samples. Therefore, we extend the semi-discrete OT map T = ∇u h to a piecewise linear (PL) mapping T as follows. The projection of u h in the source domain induces a cell decomposition of Ω, of which each cell is of µ-volume ν i and is mapped to the corresponding y i . By representing the cells by their µ-mass centers as c i := Wi(h) xdµ(x), we can get the point-wise map t : c i → y i . The Poincaré of the cell decomposition induces a triangulation of the centers . Because all of the y i s are used to construct the simplicial complex T (Y ) in the support of the target distribution, we can guarantee that no mode is lost.

C = {c i }: if W i ∩ W j = ∅, then c i is connected with c j to form an edge [c i , c j ]. Similarly, if W i0 ∩ W i1 • • • ∩ W i k = ∅, then there is a k-dimensional simplex [c i0 , c i1 , . . . , c i k ].

In practice, the µ-mass center c i is approximated by the mean value of all the Monte-Carlo samples inside W i (h), ĉi = xj ∈Wi x j /#{x j ∈ W i }, where x j ∼ µ. The connectivity information T (C) is too complicated to construct and to store in high dimensional space, thus T (C) is not explicitly built. Instead, we find the simplex σ ∈ T (C) containing x as follows: given a random point x ∈ Ω, evaluate and sort its Euclidean distances to the centers d(x, ĉi ), i = 1, 2, . . . , n in the ascending order. Suppose the first d + 1 items are {d(x, ĉi0 ), d(x, ĉi1 ), . . . , d(x, ĉi d )}, then σ is formed by {ĉ i k }. The bary-centric coordinates λi k are estimated as λi k = d -1 (x, ĉi k )/ d k=0 d -1 (x, ĉi k ). However, this may generate some spurious samples. To overcome it, we need further to detect the singular set.

Singular Set Detection According to Figalli's theory Figalli (2010) ; Chen & Figalli (2017) , if there are multiple modes or the support of the target distribution is concave, there will be singular sets Σ ⊂ Ω, where the Brenier potential is continuous but not differentiable, making its gradient map, i.e. the transport map, discontinuous.

As shown in Fig. 3 (a), the source distribution is uniformly defined on Ω, and the target empirical distribution has two modes. There is one ridge (the red line) on the Brenier potential u h , whose projection is the singular set Σ (the red line in Ω). Ω \ Σ consists of two connected components, each of them is mapped onto a single mode. Σ consists of codimension 1 facets of cells. If W i (h) ∩ W j (h) ⊂ Σ, then the dihedral angle between two supporting planes π h,i and π h,j of u h is prominently large. Therefore, on the graph of Brenier potential, we pick the pairs of facets whose dihedral angles are larger than a given threshold, the projection of their intersection gives a co-dimension 1 cell in the singular set Σ. During the generation process, if a random sample x is around Σ, it will be mapped by T to the gaps among the modes. When generating new latent codes, we just abandon such samples. And this help our method prevent the mode mixture phenomenon.

Given the extended OT map T (x), some of the polyhedrons transverse the singular set (the red lines of Fig . 3(b) ), which means that different vertices of the polyhedron belongs to different mode. If the sample x falls into such a polyhedron (the dotted red triangle), we just abandon it (as shown in Fig. 3 (b), the red dot x2 is just abandoned). Specifically, given x, we can detect if it belongs to the singular set by checking the angles θ i k between π i0 and π i k , k = 1, 2, . . . , d as θ i k = y i0 , y i k / y i0 • y i k . If all of the angles θ i k is larger than a threshold θ, we say x belongs to the singular set and just abandon it. Or we just select a subset {π i k } with θ i k ≤ θ, denoted as {π îk , k = 0, 1, . . . , d 1 }. Then we can compute λ k = d -1 (x, ĉî k )/ d1 j=0 d -1 (x, ĉî j ) and T (x) = d1 k=0 λ k T (ĉ îk ). Intuitively, T (•) smooths the discrete function T (•) in regions where latent codes are dense and keep the discontinuity of T (•) where latent codes are very sparse. In this way we avoid generating spurious latent code and thus improve the generation quality. The algorithm to generate new code is shown in Alg. 2 and the effect of threshold filtering is further investigated in Appendix 4.1.

4. EXPERIMENTS:

In order to validate that the proposed method can solve the mode collapse/mixture problems and generate controllable high quality images, several experiments are conducted.

The first experiment explores the influence of the angle thresholds for the singularity detection on synthetic dataset.

The second experiment focuses on toy sets, so that the complexity of the tasks can be manually controlled and the mode and quality of the generated samples can be accurately computed. Lin et al. (2018) did a large-scale comparison with previous methods that explicitly proposed to mitigate mode collapse and thus established a baseline for comparison. For consistent evaluation, we set up our experiment on the same benchmark dataset as theirs, and make the comparison.

In the last experiment, we run the proposed method mainly on 4 public datasets, MNIST LeCun & Cortes ( 2010 2018) did in their papers. Besides, the architecture of the decoder is the same as Lucic et al. (2018) , in which the authors did a large-scale study to evaluated the best performance of 8 different generative models including various GAN models and VAE, and the encoder is set to be the mirror of decoder.

4.1. SINGLE PARAMETER SELECTIVE INTERPOLATION:

On synthetic datasets, effects of angle threshold filtering can be visually inspected. As illustrated in Fig. 4 , number of mode is a monotonically increasing function with respect to angle threshold θ. Quality of generated samples is effected directly by choosing different θ. Generally, small θ encourages interpolation in between closely related real samples while too large θ will result in interpolation between samples from different modes, which might in turn lower generation quality. On synthetic datasets, where modes are isotropic and different modes are clearly separable, an ideal θ that captures all modes while avoids generating low quality samples can be chosen within a relatively wide band. For real world datasets of unknown modes, an ideal θ needs to be hand tuned as the separability of different modes depends largely on input data pattern and quality of the embedding map. 

4.2. MITIGATION OF MODE COLLAPSE AND MODE MIXTURE IN SYNTHETIC DATASET:

Since synthetic dataset consists of explicit distributions and known modes, mode collapse and the quality of the generated sample can be accurately measured. We choose the same synthetic datasets and metrics as in Lin et al. (2018) . Specifically, we use 2D-ring and 2D-grid for test sets and Number of modes, Percentage of high-quality samples, reverse Kullback-Leibler (KL) divergence as evaluation metrics. Number of modes counts the amount of modes captured by samples produced a generative model. Percentage of high-quality samples measures the proportion of samples that generated within three standard deviations of the nearest mode. reverse KL divergence measures how well generated samples balance among all modes regarding the real distribution. In Lin et al. (2018) , the authors evaluated GAN, ALI, MD and PacGAN on synthetic sets with above three metrics. Each networks are trained under the same generator architecture with a total of approximated 400K training parameters. For AE-OT test, since the source space and the target space are both 2-dimensional, there is no need to train any autoencoder. A two dimensional extended OT map is directly computed. Our results are included in table 1, and benchmarks of previous methods are copied from Lin et al. (2018) and Xiao et al. (2018) . Generally speaking, the samples generated by the proposed method can not only cover all of the modes, the quality of them is also better than others.

Besides, we experiment on stack MNIST dataset and CelebA dataset to further illustrate the performance of the proposed method, and the results are shown in Section C.1 and C.2 in the Appendix.

4.3. QUANTITATIVE COMPARISON WITH FID:

FID is computed by: (1) extract the visual-meaningful features of both the generated and real images through the inception network, (2) fit the features in both the generated and real feature spaces with Gaussian distribution, and then (3) compute the distance between the two Gaussian distributions with the following formula

F ID = µ r -µ g 2 2 + T r(Σ r + Σ g -2(Σ r Σ g ) 1 2

), where µ r , µ g mean the means of the real and generated features, Σ r , Σ g represent the variances of both distributions.

We report our results in Tab. (2), in which the compared data comes from Lucic et al. (2018 )Hoshen & Malik (2019) . In general, the proposed model achieves better than or comparable scores to other state-of-the-art generative models. Theoretically, the FID scores of our proposed generative models should be close to that of the pre-trained autoencoders, and this is also validated in our experiments.

The autoencoder architecture we use here cannot find a good encoding for the CelebA dataset due to the limited capacity. But the FID score of the generation model is still approach to the autoencoder.

In order to verify that with appropriate capacity of autoencoder, the proposed model works. We use the generator of DCGAN Radford et al. (2016) as the decoder of the autoencoder, then the reported FID score is 28.6, outperforming other models. Also, some of the generated images are displayed in Fig. 7 of the Appendix.

We also display the generating results for the four dataset in Fig. 5 . It includes the original images, the best generating results of Lucic et al. (2018) , including various GANs and VAE, the results of Hoshen & Malik (2019) and ours, row by row.

5. CONCLUSION:

This work gives a theoretic explanation for mode collapse/mixture by Brenier's theory and Figalli's regularity theory of optimal transport maps. When the target measure has concave support, the OT map is discontinuous on the signular sets. But DNNs can only represent continuous functions, this conflict causes the both problems. In order to solve this problem, the AE-OT model is proposed by separating manifold embedding and measure transformation. The former step is computed using an autoencoder, the latter is carried out using the extended semi-discrete OT map based on GPUs. The model is tested thoroughly and extensively by both synthetic and real data sets. The experimental results validates the discontinuity of the OT maps and demonstrate the advantages comparing to the state-of-the-arts. Theorem 3 (Kitagawa-Mérigot-Thibert Kitagawa et al. (2019) ). Assume the cost function is quadratic distance, µ has convex support and also that (i) The probability density of µ is C 0,α (Ω) for α in (0; 1]. (ii) µ has positive Poincaré-Wirtinger constant. Then the Newton algorithm for semi-discrete optimal transport converges globally with linear rate and locally with rate 1 + α.

Though in our method the gradient descend method is applied, the above theorem also ensures its convergence because of the convexity of the energy function E(h) we adopted.

B FIGALLI'S THEORY:

In this section, we show the fact that even for the case of single mode, the transport map may still be discontinuous, which will cause the instability of the training process of GANs. The arguments are mainly based on the regularity theory of transport maps developed by Figalli Chen & Figalli (2017) ; Figalli (2010) and so on.

According to Brenier's Theorem Brenier (1987; 1991a) , any transport map can be decomposed into a measure preserving map and a solution to the Monge-Ampére equation, which is the optimal transport map under the L 2 cost function. Therefore, the continuity of the transport map can be reduced to the regularity (smoothness) of the solution to the Monge-Ampére equation. When the support of the target measure is convex and the density functions are smooth, Caffarelli showed the map is differentiable; otherwise if the target domain is not convex, Figalli showed the map is discontinuous, and gave precise description of the singularity set. In this section, we briefly introduce Figallis' theory, and conduct an experiment using CelebA data set to show the existence of the singularity set, hence demonstrate the fact that the transport maps computed in GANs are discontinuous.

B.1 CONVEX DOMAINS -CAFFARELLI THEOREM:

Let Ω and Λ are two bounded open sets in R n , and let f : R n → R and g : R n → R be two positive functions such that f = 0 in R n \ Ω, g = 0 in R n \ Λ, and

Ω f = Λ g = 1.

According to Brenier's Theorem Brenier (1987; 1991a) , there exists a globally Lipschitz convex function ϕ : R n → R such that ∇ϕ # f = g and and ∇ϕ(x) ∈ Λ for L 2 -a.e. x ∈ R n . We say ϕ weakly solves the Monge-Ampére equation

EQUATION

together with the boundary condition ∇ϕ(R n ) ⊂ Λ. ϕ is called the Briener potential.

As shown by Caffarelli [9] , if Λ is convex, then ϕ is strictly convex, and it solves the Monge-Ampére equation 6. The regularity theory has been estabilished (see Caffarelli (1990a; b; 1991) ), such as

1. if λ ≤ f, g ≤ 1/λ for some λ > 0, then ϕ ∈ C 1,α(λ) loc (Ω). 2. if f ∈ C k,α loc (Ω) and g ∈ C k,α loc (Λ), then ϕ ∈ C k+2,α loc (Ω), (k ≥ 0, α ∈ (0, 1)).

B.2 NON-CONVEX DOMAINS -FIGALLI THEOREM:

However, if Λ is not convex, the regularity of the Brenier potential can not be guaranteed. Figalli gave examples in Figalli (2010) , such that 1. Ω is convex, Λ is simply connected, but non-convex; 2. the density functions f and g are smooth, f ∈ C ∞ (Ω) and g ∈ C ∞ (Λ); 3. the Brenier potential ϕ ∈ C 1 (Ω), the transport map ∇ϕ is not continuous.

In this scenario, the transport map can not be learned using DNNs, and training process is unstable or the GAN model generates unrealistic samples.

Figalli's construction Let ϕ : R n → R be a convex function. Its subdifferential at a point x is defined by ∂ϕ(x) := {y ∈ R n |ϕ(z) ≥ ϕ(x) + y • (z -x), ∀z ∈ R n }. ϕ is differentiable at a point x if and only if ∂ϕ(x) is a singleton. Figalli decomposes the set of non-differentiability points according to the dimension of the singular set:

EQUATION

For any k = 0, . . . , n, the set Σ k (ϕ) is (n -k)-rectifiable. The set of reachable subgradients at x as

∇ * ϕ := lim k→+∞ ∇ϕ(x k )|x k ∈ Σ 0 , x k → x .

It is known that the convex hull of ∇ * ϕ(x), coincides with ∂ϕ(x).

Theorem 4 (Figalli) . Assume that there exists λ > 0 such that λ 

≤ f ≤ 1/λ in Ω, λ ≤ g ≤ 1/λ in Λ,

(Ω \ Ω ) = L 2 (Λ \ Λ ) = 0, such that ϕ ∈ C 1,α (Ω )

, ∇ϕ is a bi-Hölder homeomorphism between Ω and Λ , and ϕ is an Alexandrov solution of 6 inside Ω . In particular, Caffarelli's regularity theory for strictly convex Alexandrov solutions of the Monge-Ampére equations applies to ϕ inside Ω .

Figalli studies the singular set of ϕ in Ω, i.e. the set of points x ∈ Ω where ϕ is not differentiable, denoted as Sing. Figalli shows the singularity set has the following characterization,

Sing = {x ∈ Ω|∂ϕ(x) ∩ Λ = ∅, ∇ * ϕ(x) ⊂ ∂Λ, ∂ϕ(x) ⊂ Λ} .

it can be decomposed into connected components Sing := ∪ i S i . For planar case, Theorem 5 (Figalli Singularity Set) . The number of connected components of Sing is at most countable. Moreover:

1. either S i coincides with an isolated point {xi} for some x i ∈ Ω, and in this case the boundary of ∂(xi) is entirely contained inside ∂Λ (so that ∂ϕ(xi) completely fills a hole in Λ);

2. or S i can be written as a disjoint union as follows:

S i = j γ ij ,

where γ ij : I ij → Sing are embedded Lipschitz curves parameterized by arc-lengh, I ij is an interval.

B.3 ELEMENTARY EXPERIMENTS:

We have designed a numerical experiment to verify Figalli's theorems in low dimensional case.

Fig. 6 shows another computational result, which demonstrates the singularity structure in Figalli's theorem. The source domain Ω is the unit disk, the target domain Λ is with complicated geometry. The singularity set of the optimal transport map satisfies the description of Figalli's theorem 4,

Σ 1 = 3 i=0 γ k , Σ 2 = 1 j=0 x j .

∂ϕ(x 0 ) fills the hole on Λ. For any interior point p ∈ γ 1 , ∂ϕ(p) is a line segment connecting two points on the boundary of Λ. We test diversity of generated samples from our AE-OT method on stack MNIST dataset that consists of 128,000 samples in 1,000 modes with each sample stacking three handwritten digit images from MNIST dataset LeCun et al. (1998) . Number of modes counts the amount of modes captured by samples produced a generative model. The reverse KL divergence is computed by first assign each samples to their nearest mode, and compute the KL divergence between histogram of sample count on each mode and the histogram of real data. We choose angle threshold θ = 0.5 for AE-OT method. Details of network architectures are listed in supplementary materials. Experiments results are summarized in table 5, which show our method achieves best performance in terms of modes captured and reverse KL divergence on stacked MNIST dataset.

C ADDITIONAL EXPERIMENTS:

C.2 CELEBA EXPERIMENT we evaluate our method on CelebA dataset by measuring collision probability in a batch of 1024 generated images of size 64-by-64. If a pair of identical images appear, a collision is declared, and thus higher collision probability means lower generation diversity. The same metric has been used in Lin et al. (2018) for evaluation of PacGAN. To make a consistent comparison, we design our autoencoder network with encoder having the same architecture as in previous work and decoder having a mirrored architecture of encoder. Angle threshold θ is chosen to be 0.7 for AE-OT test. Results are listed in table 6, with corresponding images can be downloaded here. Results have shown that our method achieves best result in terms of probability of collision. Autoencoder network structures can be found at table 7 and 8. Given any two images in the dataset, we can find the images between them by linear interpolation in the noise space because the one to one correspondence between µ mass centers in the noise space and the images in the dataset is provided by the proposed algorithm. For other generation models, though the interpolation can be done successfully in the noise space, they cannot find the correspondence from the noise space and the image space. The results of the linear interpolation are shown in Fig. 8 . 

:

For example, a generator generates obscure digits mixing 0 and 8 but neither 0 nor 8 on the MNIST dataset.

Figures and Tables:

Figure1: Mode collapse/mixture caused by the discontinuity of the transport map. Top row shows real data distributions, the bottom row gives the noise distributions. On top, each cluster represents a mode, the spurious generated samples are red crosses (mode mixture); at the bottom, red dotted lines are the singularity set, red crosses are mapped to be spurious samples by DNNs. (a) If the support of the target distribution is convex, DNN (f 1 ) is able to approximate the transport map f1 well. When the support of the target distributions are concave, there are two situations: (b) single mode and (c) multi modes. In (b), DNN, represented by f 2 , cannot approximate the transport map f2 well and generates some spurious samples. f3 gives the transport map of multi-mode, when approximating it with continuous DNNs, either mode collapse f 31 or mode mixture f 32 will happen.

Figure2: AE-OT model. AE: f θ and g ξ represent the encoding and decoding maps respectively, where θ and ξ are their corresponding network parameters. In the latent space Z, the latent codes are clustered into three different modes, represented as marks with different shapes (i.e. disks, squares and circles). OT: The singular set between different modes is plotted with dashed lines. Finally, the generator of our model, which generates realistic images from random noise samples, is the composition of the extended OT map T and the decoding map g ξ .

Illustration of the proposed algorithm with two modes in Y . (a) Brenier potential and the corresponding power diagram. Each cell W i is mapped to y i , which is the slope of π h,i . The red line in Ω gives the singular set. (b) The extended semi-discrete map. By computing the weighted center of each cell, and then triangulating the centers according to the power diagram, we get the PL map T (x) from Ω to Y . If the samples x 1 is not in the triangles transverse the singular set, we map it to the corresponding T (x 1 ) in Y .

Semi-discrete OT Map 1: Input: Latent codes Y = {yi} i∈I , empirical latent code distribution ν = 1 |I| i∈I δy i , number of Monte Carlo samples N , positive integer s. 2: Output: Optimal transport map T (•). 3: Initialize h = (h1, h2, . . . , h |I| ) ← (0, 0, . . . , 0). 4: repeat 5:

3: Compute ĉi by Monte Carlo method. 4: repeat 5: Sample x ∼ µ, Find the smallest d + 1 vertex around x as {d(x, ĉi 0 ), d(x, ĉi 1 ), . . . , d(x, ĉi d )}.

All these simplices form a triangulation of C (a simplicial complex), denoted as T (C) (the green triangles in the left of Fig. 3(b)). We can triangulate Y in the same way to obtain the triangulation T (Y ) (the green triangles in the right of Fig. 3(b)). Once a random sample x is drawn from the distribution µ, we can find the simplex σ in T (C) containing x. Assume the simplex σ has d + 1 vertices {c i0 , c i1 , . . . , c i d }, the bary-centric coordinates of x in σ is defined as x = d k=0 λ k c i k , and d k=0 λ k = 1 with all λ k non-negative. Then the generated latent code of x under this piece-wise linear map is given by T (x) = d k=0 λ k y i k (In Fig. 3(b), the green dot x1 is mapped to be T (x1))

), MNIST-FANSION Han Xiao & Vollgraf (2017), CIFAR-10 Krizhevsky (2009) and CelebA Zhang et al. (2018), just like the authors of Hoshen & Malik (2019) Sajjadi et al. (2018) Lucic et al. (

Figure 4: Effect of increasing angle threshold θ. (a) shows target distribution. (b) and (c) shows AE-OT results when θ is too small (as in (b)) or too large (as in (c)). (d) shows a proper choice of θ that precisely captured and generalized all modes.

Figure 5: The comparison between the proposed method and the state-of-the-art on MNIST LeCun & Cortes (2010), Fashion MNIST Han Xiao & Vollgraf (2017), Cifar10 Krizhevsky (2009) and CelebA Zhang et al. (2018). The first row shows some of the real images in each dataset. The second row corresponds to the best results of Lucic et al. (2018); The third row gives the results of Hoshen & Malik (2019); then we display our generating results in the last row.

Figure 6: Figalli's example: Singularity structure of an optimal transport map, Fig.3.2 in Figalli (2010), Lei et al. (2020).

Generative model is one of the main tasks in the machine learning field.

Table content: <table><tr><td>One of</td></tr><tr><td>the most used image generation methods is Variational Autoencoders (VAEs), where the decoders</td></tr><tr><td>approximate real data distributions from a Gaussian distribution in a variational approach (e.g Kingma</td></tr><tr><td>&amp; Welling (2013) and Rezende et al. (2014)). Later, Adversarial Autoencoders (AAEs) Makhzani</td></tr><tr><td>et al. (2015) and Wasserstein Autoencoders (WAEs) Tolstikhin et al. (2018) are proposed following</td></tr><tr><td>the similar scheme. Although VAEs are relatively simple to train, images they generate look blurry.</td></tr><tr><td>Generative Adversarial Networks (GANs) Goodfellow et al. (2014) proposed by Goodfellow et.al can</td></tr><tr><td>produce better quality images. While being a powerful tool in generating realistically looking samples,</td></tr><tr><td>GANs suffer from the mode collapse problem. To solve this problem, a huge number of methods,</td></tr><tr><td>including changing the loss function (e.g. Wasserstein GAN Arjovsky et al. (2017)), regularizing the</td></tr><tr><td>discriminators to be Lipschitz (clipping Arjovsky et al. (2017), gradient regularization Gulrajani et al.</td></tr><tr><td>(2017), Mescheder et al. (2018) or spectral normalization Miyato et al. (2018)), were proposed.</td></tr></table>

Experiments on synthetic datasets. Under standard benchmark settings, AE-OT achieves best performances over an average of 10 independent experiment results in terms of modes captured, probability of high quality samples and reverse KL divergence. The mean values and standard deviations of the experiment results are reported here.

Table content: <table><tr><td/><td/><td colspan="2">2D-ring</td><td/><td/><td/><td/><td>2D-grid</td></tr><tr><td/><td>Modes</td><td colspan="2">high quality</td><td/><td/><td/><td>Modes</td><td>high quality</td></tr><tr><td/><td>(Max 8)</td><td colspan="2">samples</td><td colspan="2">reverse KL</td><td colspan="2">(Max 25)</td><td>samples</td><td>reverse KL</td></tr><tr><td>GAN</td><td>6.3±0.5</td><td colspan="6">98.2±0.2% 0.45±0.09 17.3±0.8</td><td colspan="2">94.8±0.7% 0.70±0.07</td></tr><tr><td>ALI</td><td>6.6±0.3</td><td colspan="6">97.6±0.4% 0.36±0.04 24.1±0.4</td><td colspan="2">95.7±0.6% 0.14±0.03</td></tr><tr><td>MD</td><td>4.3±0.8</td><td colspan="6">36.6±8.8% 1.93±0.11 23.8±0.5</td><td colspan="2">79.9±3.2% 0.18±0.03</td></tr><tr><td>PacGAN2</td><td>7.9±0.1</td><td colspan="6">95.6±2.0% 0.07±0.03 23.8±0.7</td><td colspan="2">91.3±0.8% 0.13±0.04</td></tr><tr><td>PacGAN3</td><td>7.8±0.1</td><td colspan="6">97.7±0.3% 0.10±0.02 24.6±0.4</td><td colspan="2">94.2±0.4% 0.06±0.02</td></tr><tr><td>PacGAN4</td><td>7.8±0.1</td><td colspan="6">95.9±1.4% 0.07±0.02 24.8±0.2</td><td colspan="2">93.6±0.6% 0.04±0.01</td></tr><tr><td>BourGAN</td><td>8.0±0.0</td><td colspan="6">99.8±2.9% 4e-4±2e-4 24.9±0.3</td><td colspan="2">95.9±0.2% 0.01±0.02</td></tr><tr><td>AE-OT</td><td>8.0±0.0</td><td colspan="6">99.6±0.3% 0.004±0.001 25.0±0.0</td><td colspan="2">99.8±0.2% 0.007±0.002</td></tr><tr><td/><td/><td colspan="6">Table 2: Quantitative comparison with FID</td><td/></tr><tr><td/><td/><td colspan="2">Adversarial</td><td/><td/><td/><td colspan="2">Non-Adversarial</td><td>Reference</td></tr><tr><td>Dataset</td><td>NS GAN</td><td>LSGAN</td><td colspan="2">WGAN</td><td colspan="2">BEGAN</td><td>VAE</td><td>GLANN</td><td>AE</td><td>Ours</td></tr><tr><td>MNIST</td><td>6.8±0.5</td><td>7.8±0.6</td><td colspan="2">6.7±0.4</td><td colspan="3">13.1±1.0 23.8±0.6</td><td>8.6±0.1</td><td>5.5</td><td>6.2±0.2</td></tr><tr><td>Fansion</td><td>26.5±1.6</td><td>30.7±2.2</td><td colspan="6">21.5±1.6 22.9±0.9 58.7±1.2 13.0±0.1</td><td>4.7</td><td>10.1±0.3</td></tr><tr><td colspan="10">CIFAR-10 58.5±1.9 87.1±47.5 55.2±2.3 71.4±1.6 65.4±0.2 46.5±0.2 28.2</td><td>38.3±0.5</td></tr><tr><td>CelebA</td><td>55.0±3.3</td><td>53.9±2.8</td><td colspan="7">41.3±2.0 38.9±0.9 85.7±3.8 46.3±0.1 67.5 68.4± 0.5</td></tr></table>

and that ∂Ω and ∂Λ are continuous. Then ϕ is strictly convex inside Ω. Moreover there exist two open sets Ω ⊂ Ω and Λ ⊂ Λ, with L 2

Table content: <table/>

Encoder architecture for stack MNIST layer number of outputs kernel size stride BN activation

Table content: <table><tr><td colspan="2">Input x ∼ P data 28*28*3</td><td/><td/><td/></tr><tr><td>Convolution</td><td>14*14*16</td><td>5*5</td><td>2</td><td>LeakyReLU</td></tr><tr><td>Convolution</td><td>7*7*32</td><td>5*5</td><td>2</td><td>Yes LeakyReLU</td></tr><tr><td>Convolution</td><td>4*4*64</td><td>5*5</td><td>2</td><td>Yes LeakyReLU</td></tr><tr><td>Convolution</td><td>2*2*128</td><td>5*5</td><td>2</td><td>Yes LeakyReLU</td></tr><tr><td colspan="2">Fully connected 100</td><td/><td/><td/></tr><tr><td colspan="5">architecture same as the consistent generator architectures in GANs, and encoder having mirrored</td></tr><tr><td>architecture.</td><td/><td/><td/><td/></tr></table>

Decoder architecture for stack MNIST layer number of outputs kernel size stride BN activation

Table content: <table><tr><td>Input z ∼ P latent</td><td>100</td><td/><td/></tr><tr><td>Fully connected</td><td>2*2*128</td><td/><td/><td>Yes ReLU</td></tr><tr><td colspan="2">Transposed Convolution 4*4*64</td><td>5*5</td><td>2</td><td>Yes ReLU</td></tr><tr><td colspan="2">Transposed Convolution 7*7*32</td><td>5*5</td><td>2</td><td>Yes ReLU</td></tr><tr><td colspan="2">Transposed Convolution 14*14*16</td><td>5*5</td><td>2</td><td>Yes ReLU</td></tr><tr><td colspan="2">Transposed Convolution 28*28*3</td><td>5*5</td><td>2</td><td>Tanh</td></tr></table>

Experiments on stacked MNIST. Results have shown that our method achieves best results in terms of mode captured and reverse KL divergence. (*) In WGAN, PacWGAN and AE-OT experiments, number of feature maps in each network layer is a quarter of those in other experiments.

Table content: <table><tr><td>Stacked MNIST</td><td/></tr><tr><td>Modes</td><td>KL</td></tr></table>

Probability of identical images in a batch of 1024 generated images from DCGAN, PacGAN2 and AE-OT. Results have shown that our method achives best result in terms of collision probability on CelebA dataset.

Table content: <table><tr><td>Discriminator size</td><td colspan="2">Probability of collision</td><td/></tr><tr><td>(Decoder size)</td><td colspan="3">DCGAN PacDCGAN2 AE-OT</td></tr><tr><td>273K</td><td>1</td><td>0.33</td><td>0</td></tr><tr><td>4×273K</td><td>0.42</td><td>0</td><td>0</td></tr><tr><td>16×273K</td><td>0.86</td><td>0</td><td>0</td></tr><tr><td>25×273K</td><td>0.65</td><td>0.17</td><td>0</td></tr></table>

Encoder architecture in CelebA experiment layer number of outputs kernel size stride BN activation

Table content: <table><tr><td colspan="2">Input x ∼ P data 64*64*3</td><td/><td/><td/></tr><tr><td>Convolution</td><td>32*32*dim_f</td><td>4*4</td><td>2</td><td>LeakyReLU</td></tr><tr><td>Convolution</td><td>16*16*dim_f*2</td><td>4*4</td><td>2</td><td>Yes LeakyReLU</td></tr><tr><td>Convolution</td><td>8*8*dim_f*4</td><td>4*4</td><td>2</td><td>Yes LeakyReLU</td></tr><tr><td>Convolution</td><td>4*4*dim_f*8</td><td>4*4</td><td>2</td><td>Yes LeakyReLU</td></tr><tr><td>Convolution</td><td>100</td><td>4*4</td><td>1</td><td/></tr></table>

Decoder architecture in CelebA experiment layer number of outputs kernel size stride BN activation

Table content: <table><tr><td>Input z ∼ P latent</td><td>100</td><td/><td/></tr><tr><td colspan="2">Transposed Convolution 4*4*dim_f*8</td><td/><td/></tr><tr><td colspan="2">Transposed Convolution 8*8*dim_f*4</td><td>4*4</td><td>2</td><td>Yes ReLU</td></tr><tr><td colspan="2">Transposed Convolution 16*16*dim_f*2</td><td>4*4</td><td>2</td><td>Yes ReLU</td></tr><tr><td colspan="2">Transposed Convolution 32*32*dim_f</td><td>4*4</td><td>2</td><td>Yes ReLU</td></tr><tr><td colspan="2">Transposed Convolution 64*64*3</td><td>4*4</td><td>2</td><td>Tanh</td></tr></table>

Peer review (iclr2020):
This paper deals with an important problem of mode collapse and mode mixture. In order to tackle the both problems, the paper proposes to separate the manifold embedding and the optimal transportation problems; the first part being carried out using an autoencoder to map the images onto the latent space and the second part is accomplished using a GPU-based convex optimization to find the discontinuous transportation maps. I have some doubts about moving from the *semi-discrete OT map* to the piece-wise linear extension. The illustration in Fig. 3, and implicit in all the explanation charts is the fact that discontinuity can be found by a linear separation. This seems to be an extremely simplifying assumption, which leads to not so great visual results from the paper. Although the numerical results seems promising, I feel that fewer images, but larger in size, and analysis of mode collapse phenomenon in real images would have been much better. Singular set detection seems to be the most tricky part in this paper, which should have been explained further. The Simplex projection assumption, renders this part not that tricky, but that is where I feel the biggest doubt about this paper lies. The authors themselves mention the need for a high quality auto encoder model to encode celebA dataset, which has been improved upon by numerous other papers, the claims seems not too strong. Also, the method does not have any adversarial training and hence, it studies the GAN idea from only fixing the generator point of view.