Year: 2018

Abstract:
Implicit probabilistic models are models defined naturally in terms of a sampling procedure and often induces a likelihood function that cannot be expressed explicitly. We develop a simple method for estimating parameters in implicit models that does not require knowledge of the form of the likelihood function or any derived quantities, but can be shown to be equivalent to maximizing likelihood under some conditions. Our result holds in the non-asymptotic parametric setting, where both the capacity of the model and the number of data examples are finite. We also demonstrate encouraging experimental results.

heading: 1 INTRODUCTION
text: Generative modelling is a cornerstone of machine learning and has received increasing attention. Recent models like variational autoencoders (VAEs) (Kingma & Welling, 2013; Rezende et al., 2014) and generative adversarial nets (GANs) (Goodfellow et al., 2014; Gutmann et al., 2014), have delivered impressive advances in performance and generated a lot of excitement.
Generative models can be classified into two categories: prescribed models and implicit models (Diggle & Gratton, 1984; Mohamed & Lakshminarayanan, 2016). Prescribed models are defined by an explicit specification of the density, and so their unnormalized complete likelihood can be usually expressed in closed form. Examples include models whose complete likelihoods lie in the exponential family, such as mixture of Gaussians (Everitt, 1985), hidden Markov models (Baum & Petrie, 1966), Boltzmann machines (Hinton & Sejnowski, 1986). Because computing the normalization constant, also known as the partition function, is generally intractable, sampling from these models is challenging.
On the other hand, implicit models are defined most naturally in terms of a (simple) sampling procedure. Most models take the form of a deterministic parameterized transformation Tθ(·) of an analytic distribution, like an isotropic Gaussian. This can be naturally viewed as the distribution induced by the following sampling procedure:
1. Sample z ∼ N (0, I) 2. Return x := Tθ(z)
The transformation Tθ(·) often takes the form of a highly expressive function approximator, like a neural net. Examples include generative adversarial nets (GANs) (Goodfellow et al., 2014; Gutmann et al., 2014) and generative moment matching nets (GMMNs) (Li et al., 2015; Dziugaite et al., 2015). The marginal likelihood of such models can be characterized as follows:
pθ(x) = ∂ ∂x1 · · · ∂ ∂xd ∫ {z|∀i (Tθ(z))i≤xi } φ(z)dz
where φ(·) denotes the probability density function (PDF) of N (0, I). In general, attempting to reduce this to a closed-form expression is hopeless. Evaluating it numerically is also challenging, since the domain of integration could consist of an exponential number of disjoint regions and numerical differentiation is ill-conditioned.
These two categories of generative models are not mutually exclusive. Some models admit both an explicit specification of the density and a simple sampling procedure and so can be considered as both prescribed and implicit. Examples include variational autoencoders (Kingma & Welling,
2013; Rezende et al., 2014), their predecessors (MacKay, 1995; Bishop et al., 1998) and extensions (Burda et al., 2015), and directed/autoregressive models, e.g., (Neal, 1992; Bengio & Bengio, 2000; Larochelle & Murray, 2011; van den Oord et al., 2016).

heading: 1.1 CHALLENGES IN PARAMETER ESTIMATION
text: Maximum likelihood (Fisher, 1912; Edgeworth, 1908) is perhaps the standard method for estimating the parameters of a probabilistic model from observations. The maximum likelihood estimator (MLE) has a number of appealing properties: under mild regularity conditions, it is asymptotically consistent, efficient and normal. A long-standing challenge of training probabilistic models is the computational roadblocks of maximizing the log-likelihood function directly.
For prescribed models, maximizing likelihood directly requires computing the partition function, which is intractable for all but the simplest models. Many powerful techniques have been developed to attack this problem, including variational methods (Jordan et al., 1999), contrastive divergence (Hinton, 2002; Welling & Hinton, 2002), score matching (Hyvärinen, 2005) and pseudolikelihood maximization (Besag, 1975), among others.
For implicit models, the situation is even worse, as there is no term in the log-likelihood function that is in closed form; evaluating any term requires computing an intractable integral. As a result, maximizing likelihood in this setting seems hopelessly difficult. A variety of likelihood-free solutions have been proposed that in effect minimize a divergence measure between the data distribution and the model distribution. They come in two forms: those that minimize an f -divergence, and those that minimize an integral probability metric (Müller, 1997). In the former category are GANs, which are based on the idea of minimizing the distinguishability between data and samples (Tu, 2007; Gutmann & Hyvärinen, 2010). It has been shown that when given access to an infinitely powerful discriminator, the original GAN objective minimizes the Jensen-Shannon divergence, the − logD variant of the objective minimizes the reverse KL-divergence minus a bounded quantity (Arjovsky & Bottou, 2017), and later extensions (Nowozin et al., 2016) minimize arbitrary f -divergences. In the latter category are GMMNs which use maximum mean discrepancy (MMD) (Gretton et al., 2007) as the witness function.
In the case of GANs, despite the theoretical results, there are a number of challenges that arise in practice, such as mode dropping/collapse (Goodfellow et al., 2014; Arora & Zhang, 2017), vanishing gradients (Arjovsky & Bottou, 2017; Sinn & Rawat, 2017) and training instability (Goodfellow et al., 2014; Arora et al., 2017). A number of explanations have been proposed to explain these phenomena and point out that many theoretical results rely on three assumptions: the discriminator must have infinite modelling capacity (Goodfellow et al., 2014; Arora et al., 2017), the number of samples from the true data distribution must be infinite (Arora et al., 2017; Sinn & Rawat, 2017) and the gradient ascent-descent procedure (Arrow et al., 1958; Schmidhuber, 1992) can converge to a global pure-strategy Nash equilibrium (Goodfellow et al., 2014; Arora et al., 2017). When some of these assumptions do not hold, the theoretical guarantees do not necessarily apply. A number of ways have been proposed that alleviate some of these issues, e.g., (Zhao et al., 2016; Salimans et al., 2016; Donahue et al., 2016; Dumoulin et al., 2016; Arjovsky et al., 2017; Hjelm et al., 2017; Li et al., 2017; Zhu et al., 2017), but a way of solving all three issues simultaneously remains elusive.

heading: 1.2 OUR CONTRIBUTION
text: In this paper, we present an alternative method for estimating parameters in implicit models. Like the methods above, our method is likelihood-free, but can be shown to be equivalent to maximizing likelihood under some conditions. Our result holds when the capacity of the model is finite and the number of data examples is finite. The idea behind the method is simple: it finds the nearest sample to each data example and optimizes the model parameters to pull the sample towards it. The direction in which nearest neighbour search is performed is important: the proposed method ensures each data example has a similar sample, which contrasts with an alternative approach of pushing each sample to the nearest data example, which would ensure that each sample has a similar data example. The latter approach would permit all samples being similar to one data example. Such a scenario would be heavily penalized by the former approach.
The proposed method could sidestep the three issues mentioned above: mode collapse, vanishing gradients and training instability. Modes are not dropped because the loss ensures each data example has a sample nearby at optimality; gradients do not vanish because the gradient of the distance between a data example and its nearest sample does not become zero unless they coincide; training is stable because the estimator is the solution to a simple minimization problem. By leveraging recent advances in fast nearest neighbour search algorithms (Li & Malik, 2016; 2017), this approach is able to scale to large, high-dimensional datasets.

heading: 2 IMPLICIT MAXIMUM LIKELIHOOD ESTIMATOR

heading: 2.1 DEFINITION
text: We are given a set of n data examples x1, . . . ,xn and some unknown parameterized probability distribution Pθ with density pθ. We also have access to an oracle that allows us to draw independent and identically distributed (i.i.d.) samples from Pθ.
Let x̃θ1, . . . , x̃ θ m be i.i.d. samples from Pθ, where m ≥ n. For each data example xi, we define a random variable Rθi to be the distance between xi and the nearest sample. More precisely,
Rθi := min j∈[m] ∥∥x̃θj − xi∥∥22 where [m] denotes {1, . . . ,m}.
The implicit maximum likelihood estimator θ̂IMLE is defined as:
θ̂IMLE := arg min θ ERθ1,...,Rθn [ n∑ i=1 Rθi ]
= arg min θ Ex̃θ1,...,x̃θm [ n∑ i=1 min j∈[m] ∥∥x̃θj − xi∥∥22 ]

heading: 2.2 ALGORITHM
text: We outline the proposed parameter estimation procedure in Algorithm 1. In each outer iteration, we drawm i.i.d. samples from the current model Pθ. We then randomly select a batch of examples from the dataset and find the nearest sample from each data example. We then run a standard iterative optimization algorithm, like stochastic gradient descent (SGD), to minimize a sample-based version of the Implicit Maximum Likelihood Estimator (IMLE) objective.
Algorithm 1 Implicit maximum likelihood estimation (IMLE) procedure
Require: The dataset D = {xi}ni=1 and a sampling mechanism for the implicit model Pθ Initialize θ to a random vector for k = 1 to K do
Draw i.i.d. samples x̃θ1, . . . , x̃ θ m from Pθ Pick a random batch S ⊆ {1, . . . , n} σ(i)← arg minj ∥∥xi − x̃θj∥∥22 ∀i ∈ S for l = 1 to L do
Pick a random mini-batch S̃ ⊆ S θ ← θ − η∇θ ( n |S̃| ∑ i∈S̃ ∥∥∥xi − x̃θσ(i)∥∥∥2 2 ) end for
end for return θ
Because our algorithm needs to solve a nearest neighbour search problem in each outer iteration, the scalability of our method depends on our ability to find the nearest neighbours quickly. This was traditionally considered to be a hard problem, especially in high dimensions. However, this is no
longer the case, due to recent advances in nearest neighbour search algorithms (Li & Malik, 2016; 2017).
Note that the use of Euclidean distance is not a major limitation of the proposed approach. A variety of distance metrics are either exactly or approximately equivalent to Euclidean distance in some non-linear embedding space, in which case the theoretical guarantees are inherited from the Euclidean case. This encompasses popular distance metrics used in the literature, like the Euclidean distance between the activations of a neural net, which is often referred to as a perceptual similarity metric (Salimans et al., 2016; Dosovitskiy & Brox, 2016). The approach can be easily extended to use these metrics, though because this is the initial paper on this method, we focus on the vanilla setting of Euclidean distance in the natural representation of the data, e.g.: pixels, both for simplicity/clarity and for comparability to vanilla versions of other methods that do not use auxiliary sources of labelled data or leverage domain-specific prior knowledge. For distance metrics that cannot be embedded in Euclidean space, the analysis can be easily adapted with minor modifications as long as the volume of a ball under the metric has a simple dependence on its radius.

heading: 3 WHY MAXIMUM LIKELIHOOD
text: There has been debate (Huszár, 2015) over whether maximizing likelihood of the data is the appropriate objective for the purposes of learning generative models. Recall that maximizing likelihood is equivalent to minimizing DKL (pdata ‖pθ ), where pdata denotes the empirical data distribution and pθ denotes the model distribution. One proposed alternative is to minimize the reverse KLdivergence, DKL (pθ ‖pdata ), which is suggested (Huszár, 2015) to be better because it severely penalizes the model for generating an implausible sample, whereas the standard KL-divergence, DKL (pdata ‖pθ ), severely penalizes the model for assigning low density to a data example. As a result, when the model is underspecified, i.e. has less capacity than what’s necessary to fit all the modes of the data distribution, minimizing DKL (pθ ‖pdata ) leads to a narrow model distribution that concentrates around a few modes, whereas minimizingDKL (pdata ‖pθ ) leads to a broad model distribution that hedges between modes. The success of GANs in generating good samples is often attributed to the former phenomenon (Arjovsky & Bottou, 2017).
This argument, however, relies on the assumption that we have access to an infinite number of samples from the true data distribution. In practice, however, this assumption rarely holds: if we had access to the true data distribution, then there is usually no need to fit a generative model, since we can simply draw samples from the true data distribution. What happens when we only have the empirical data distribution? Recall that DKL (p ‖q ) is defined and finite only if p is absolutely continuous w.r.t. q, i.e.: q(x) = 0 implies p(x) = 0 for all x. In other words, DKL (p ‖q ) is defined and finite only if the support of p is contained in the support of q. Now, consider the difference between DKL (pdata ‖pθ ) and DKL (pθ ‖pdata ): minimizing the former, which is equivalent to maximizing likelihood, ensures that the support of the model distribution contains all data examples, whereas minimizing the latter ensures that the support of the model distribution is contained in the support of the empirical data distribution, which is just the set of data examples. In other words, maximum likelihood disallows mode dropping, whereas minimizing reverse KL-divergence forces the model to assign zero density to unseen data examples and effectively prohibits generalization. Furthermore, maximum likelihood discourages the model from assigning low density to any data example, since doing so would make the likelihood, which is the product of the densities at each of the data examples, small.
From the modelling perspective, because maximum likelihood is guaranteed to preserve all modes, it can make use of all available training data and can therefore be used to train high-capacity models that have a large number of parameters. In contrast, using an objective that permits mode dropping allows the model to pick and choose which data examples it wants to model. As a result, if the goal is to train a high-capacity model that can learn the underlying data distribution, we would not be able to do so using such an objective because we have no control over which modes the model chooses to drop. Put another way, we can think about the model’s performance along two axes: its ability to generate plausible samples (precision) and its ability to generate all modes of the data distribution (recall). A model that successfully learns the underlying distribution should score high along both axes. If mode dropping is allowed, then an improvement in precision may be achieved at the expense of lower recall and could represent a move to a different point on the same precision-recall
curve. As a result, since sample quality is an indicator of precision, improvement in sample quality in this setting may not mean an improvement in density estimation performance. On the other hand, if mode dropping is disallowed, since full recall is always guaranteed, an improvement in precision is achieved without sacrificing recall and so implies an upwards shift in the precision-recall curve. In this case, an improvement in sample quality does signify an improvement in density estimation performance, which may explain sample quality historically was an important way to evaluate the performance of generative models, most of which maximized likelihood. With the advent of generative models that permit mode dropping, however, sample quality is no longer a reliable indicator of density estimation performance, since good sample quality can be trivially achieved by dropping all but a few modes. In this setting, sample quality can be misleading, since a model with low recall on a lower precision-recall curve can achieve a better precision than a model with high recall on a higher precision-recall curve. Since it is hard to distinguish whether an improvement in sample quality is due to a move along the same precision-recall curve or a real shift in the curve, an objective that disallows mode dropping is critical tool that researchers can use to develop better models, since they can be sure that an apparent improvement in sample quality is due to a shift in the precision-recall curve.

heading: 4 ANALYSIS
text: Before formally stating the theoretical results, we first illustrate the intuition behind why the proposed estimator is equivalent to maximum likelihood estimator under some conditions. For simplicity, we will consider the special case where we only have a single data example x1 and a single sample x̃θ1. Consider the total density of Pθ inside a ball of radius of t centred at x1 as a function of t, a function that will be denoted as F̃ θ(t). If the density in the neighbourhood of x1 is high, then F̃ θ(t) would grow rapidly as t increases. If, on the other hand, the density in the neighbourhood of x1 is low, then F̃ θ(t) would grow slowly. So, maximizing likelihood is equivalent to making F̃ θ(t) grow as fast as possible. To this end, we can maximize the area under the function F̃ θ(t), or equivalently, minimize the area under the function 1 − F̃ θ(t). Observe that F̃ θ(t) can be interpreted as the cumulative distribution function (CDF) of the Euclidean distance between x1 and x̃θ1, which is a random variable because x̃θ1 is random and will be denoted as R̃
θ. Because R̃θ is nonnegative, recall that E [ R̃θ ] = ∫∞ 0 Pr ( R̃θ > t ) dt = ∫∞ 0 ( 1− F̃ θ(t) ) dt, which is exactly the
area under the function 1− F̃ θ(t). Therefore, we can maximize likelihood of a data example x1 by minimizing E [ R̃θ ] , or in other words, minimizing the expected distance between the data example
and a random sample. To extend this analysis to the case with multiple data examples, we show in the supplementary material that if the objective function is a summation, applying a monotonic transformation to each term and then reweighting appropriately preserves the optimizer under some conditions.
We now state the key theoretical result formally. Please refer to the supplementary material for the proof.
Theorem 1. Consider a set of observations x1, . . . ,xn, a parameterized family of distributions Pθ with probability density function (PDF) pθ(·) and a unique maximum likelihood solution θ∗. For any m ≥ 1, let x̃θ1, . . . , x̃θm ∼ Pθ be i.i.d. random variables and define r̃θ := ∥∥x̃θ1∥∥22, Rθ := minj∈[m]
∥∥x̃θj∥∥22 andRθi := minj∈[m] ∥∥x̃θj − xi∥∥22. Let F θ(·) be the cumulative distribution function (CDF) of r̃θ and Ψ(z) := minθ { E [ Rθ ] |pθ(0) = z } .

heading: If Pθ satisfies the following:
text: • pθ(x) is differentiable w.r.t. θ and continuous w.r.t. x everywhere.
• ∀θ,v, there exists θ′ such that pθ(x) = pθ′(x + v) ∀x.
• For any θ1, θ2, there exists θ0 such that F θ0(t) ≥ max { F θ1(t), F θ2(t) } ∀t ≥ 0 and
pθ0(0) = max {pθ1(0), pθ2(0)}.
• ∃τ > 0 such that ∀i ∈ [n] ∀θ /∈ Bθ∗(τ), pθ(xi) < pθ∗(xi), where Bθ∗(τ) denotes the ball centred at θ∗ of radius τ .
• Ψ(z) is differentiable everywhere.
• For all θ, if θ 6= θ∗, there exists j ∈ [d] such that〈 Ψ′(pθ(x1))pθ(x1) Ψ′(pθ∗ (x1))pθ∗ (x1) ...
Ψ′(pθ(xn))pθ(xn) Ψ′(pθ∗ (xn))pθ∗ (xn)
 ,  ∇θ (log pθ(x1))j... ∇θ (log pθ(xn))j 〉 6= 0. Then,
arg min θ n∑ i=1
E [ Rθi ]
Ψ′(pθ∗(xi))pθ∗(xi) = arg max θ n∑ i=1 log pθ(xi)
Furthermore, if pθ∗(x1) = · · · = pθ∗(xn), then,
arg min θ n∑ i=1 E [ Rθi ] = arg max θ n∑ i=1 log pθ(xi)
Now, we examine the restrictiveness of each condition. The first condition is satisfied by nearly all analytic distributions. The second condition is satisfied by nearly all distributions that have an unrestricted location parameter, since one can simply shift the location parameter by v. The third condition is satisfied by most distributions that have location and scale parameters, like a Gaussian distribution, since the scale can be made arbitrarily low and the location can be shifted so that the constraint on pθ(·) is satisfied. The fourth condition is satisfied by nearly all distributions, whose density eventually tends to zero as the distance from the optimal parameter setting tends to infinity. The fifth condition requires minθ { E [ Rθ ] |pθ(0) = z } to change smoothly as z changes. The final condition requires the two n-dimensional vectors, one of which can be chosen from a set of d vectors, to be not exactly orthogonal. As a result, this condition is usually satisfied when d is large, i.e. when the model is richly parameterized.
There is one remaining difficulty in applying this theorem, which is that the quantity 1/Ψ′(pθ∗(xi))pθ∗(xi), which appears as an coefficient on each term in the proposed objective, is typically not known. If we consider a new objective that ignores the coefficients, i.e. ∑n i=1 E [ Rθi ] ,
then minimizing this objective is equivalent to minimizing an upper bound on the ideal objective,∑n i=1 E [ Rθi ] /Ψ′(pθ∗(xi))pθ∗(xi). The tightness of this bound depends on the difference between the highest and lowest likelihood assigned to individual data points at the optimum, i.e. the maximum likelihood estimate of the parameters. Such a model should not assign high likelihoods to some points and low likelihoods to others as long as it has reasonable capacity, since doing so would make the overall likelihood, which is the product of the likelihoods of individual data points, low. Therefore, the upper bound is usually reasonably tight.

heading: 5 EXPERIMENTS
text: We trained generative models using the proposed method on three standard benchmark datasets, MNIST, the Toronto Faces Dataset (TFD) and CIFAR-10. All models take the form of feedforward neural nets with isotropic Gaussian noise as input.
For MNIST, the architecture consists of two fully connected hidden layers with 1200 units each followed by a fully connected output layer with 784 units. ReLU activations were used for hidden layers and sigmoids were used for the output layer. For TFD, the architecture is wider and consists of two fully connected hidden layers with 8000 units each followed by a fully connected output layer with 2304 units. For both MNIST and TFD, the dimensionality of the noise vector is 100.
For CIFAR-10, we used a simple convolutional architecture with 1000-dimensional Gaussian noise as input. The architecture consists of five convolutional layers with 512 output channels and a kernel size of 5 that all produce 4 × 4 feature maps, followed by a bilinear upsampling layer that doubles the width and height of the feature maps. There is a batch normalization layer followed by leaky ReLU activations with slope −0.2 after each convolutional layer. This design is then repeated for each subsequent level of resolution, namely 8 × 8, 16× 16 and 32× 32, so that we have 20 convolutional layers, each with output 512 channels. We then add a final output layer with three output channels on top, followed by sigmoid activations. We note that this architecture has more capacity than typical architectures used in other methods, like (Radford et al., 2015). This is because our method aims to capture all modes of the data distribution and therefore needs more modelling capacity than methods that are permitted to drop modes.
Evaluation for implicit generative models in general remains an open problem. Various intrinsic and extrinsic evaluation metrics have been proposed, all of which have limitations. Extrinsic evaluation metrics measure performance indirectly via a downstream task (Salimans et al., 2016). Unfortunately, dependence on the downstream task could introduce bias and may not capture de-
sirable properties of the generative model that do not affect performance on the task. Intrinsic evaluation metrics measure performance without relying on external models or data. Popular examples include estimated log-likelihood (Bengio et al., 2014; Wu et al., 2016) and visual assessment of sample quality. While recent literature has focused more on the latter and less on the former, it should be noted that they evaluate different properties – sample quality reflects precision, i.e.: how accurate the model samples are compared to the ground truth, whereas estimated log-likelihood focuses on recall, i.e.: how much of the diversity in the data distribution the model captures. Consequently, both are important metrics; one is not a replacement for the other. As pointed out by (Theis et al., 2015), “qualitative as well as quantitative analyses based on model samples can be misleading about a model’s density estimation performance, as well as the probabilistic model’s performance in applications other than image synthesis.” Two models that achieve different levels of precision may simply be at different points on the same precision-recall curve, and therefore may not be directly comparable. Models that achieve the same level of recall, on the other hand, may be directly compared. So, for methods that maximize likelihood, which are guaranteed to preserve all modes and achieve full recall, both sample quality and estimated log-likelihood capture precision. Because most generative models traditionally maximized likelihood or a lower bound on the likelihood, the only property that differed across models was precision, which may explain why sample quality has historically been seen as an important indicator of performance. However, in heterogenous experimental settings with different models optimized for various objectives, sample quality does not necessarily reflect how well a model learns the underlying data distribution. Therefore, under these settings, both precision and recall need to be measured. While there is not yet a reliable way to measure recall (given the known issues of estimated log-likelihoods in high dimensions), this does not mean that sample quality can be a valid substitute for estimated log-likelihoods, as it cannot detect the lack of diversity of samples. A secondary issue that is more easily solvable is that samples presented in papers are sometimes cherry-picked; as a result, they capture the maximum sample quality, but not necessarily the mean sample quality.
To mitigate these problems to some extent, we avoid cherry-picking and visualize randomly chosen samples, which are shown in Figure 1. We also report the estimated log-likelihood in Table 1. As mentioned above, both evaluation criteria have biases/deficiencies, so performing well on either of these metrics does not necessarily indicate good density estimation performance. However, not performing badly on either metric can provide some comfort that the model is simultaneously able to achieve reasonable precision and recall.
As shown in Figure 1, despite its simplicity, the proposed method is able to generate reasonably good samples for MNIST, TFD and CIFAR-10. While it is commonly believed that minimizing reverse KL-divergence is necessary to produce good samples and maximizing likelihood necessarily leads to poor samples (Grover et al., 2017), the results suggest that this is not necessarily the case. Even
though Euclidean distance was used in the objective, the samples do not appear to be desaturated or overly blurry. Samples also seem fairly diverse. This is supported by the estimated log-likelihood results in Table 1. Because the model achieved a high score on that metric on both MNIST and TFD, this suggests that the model did not suffer from significant mode dropping.
In Figure 4 in the supplementary material, we show samples and their nearest neighbours in the training set. Each sample is quite different from its nearest neighbour in the training set, suggesting that the model has not overfitted to examples in the training set.
Next, we visualize the learned manifold by walking along a geodesic on the manifold between pairs of samples. More concretely, we generate five samples, arrange them in arbitrary order, perform linear interpolation in latent variable space between adjacent pairs of samples, and generate an image from the interpolated latent variable. As shown in Figure 3, the images along the path of interpolation appear visually plausible and do not have noisy artifacts. In addition, the transition from one image to the next appears smooth, including for CIFAR-10, which contrasts with findings in the literature that suggest the transition between two natural images tends to be abrupt. This indicates that the support of the model distribution has not collapsed to a set of isolated points and that the proposed method is able to learn the geometry of the data manifold, even though it does not learn a distance metric explicitly.
Finally, we illustrate the evolution of samples as training progresses in Figure 2. As shown, the samples are initially blurry and become sharper over time. Importantly, sample quality consistently improves over time, which demonstrates the stability of training.
While our sample quality may not be state-of-the-art, it is important to remember that these results are obtained under the setting of full recall. So, this does not necessarily mean that our method models the underlying data distribution less accurately than other methods that achieve better sample quality, as some of them may drop modes and therefore achieve less than full recall. As previously mentioned, this does not suggest a fundamental tradeoff between precision and recall that cannot be overcome – on the contrary, our method provides researchers with a way of designing models that can improve the precision-recall curve without needing to worry that the observed improvements are due to a movement along the curve. With refinements to the model, it is possible to move the curve upwards and obtain better sample quality at any level of recall as a consequence. This is left for future work; as this is the initial paper on this approach, its value stems from the foundation it lays for a new research direction upon which subsequent work can be built, as opposed to the current results themselves. For this paper, we made a deliberate decision to keep the model simple, since non-essential practically motivated enhancements are less grounded in theory, may obfuscate the key underlying idea and could impart the impression that they are critical to making the approach work in practice. The fact that our method is able to generate more plausible samples on CIFAR-10 than other methods at similar stages of development, such as the initial versions of GAN (Goodfellow et al., 2014) and PixelRNN (van den Oord et al., 2016), despite the minimal sophistication of our method and architecture, shows the promise of the approach. Later iterations of other methods incorporate additional supervision in the form of pretrained weights and/or make task-specific modifications to the architecture and training procedure, which were critical to achieving state-of-the-art sample quality. We do believe the question of how the architecture should be refined in the context of our method to take advantage of task-specific insights is an important one, and is an area ripe for future exploration.

heading: 6 DISCUSSION
text: In this section, we consider and address some possible concerns about our method.

heading: 6.1 DOES MAXIMIZING LIKELIHOOD NECESSARILY LEAD TO POOR SAMPLE QUALITY?
text: It has been suggested (Huszár, 2015) that maximizing likelihood leads to poor sample quality because when the model is underspecified, it will try to cover all modes of the empirical data distribution and therefore assign high density to regions with few data examples. There is also empirical evidence (Grover et al., 2017) for a negative correlation between sample quality and log likelihood, suggesting an inherent trade-off between maximizing likelihood and achieving good sample quality. A popular solution is to minimize reverse KL-divergence instead, which trades off recall for pre-
cision. This is an imperfect solution, as the ultimate goal is to model all the modes and generate high-quality samples.
Note that this apparent trade-off exists that the model capacity is assumed to be fixed. We argue that a more promising approach would be to increase the capacity of the model, so that it is less underspecified. As the model capacity increases, avoiding mode dropping becomes more important, because otherwise there will not be enough training data to fit the larger number of parameters to. This is precisely a setting appropriate for maximum likelihood. As a result, it is possible that a combination of increasing the model capacity and maximum likelihood training can achieve good precision and recall simultaneously.

heading: 6.2 WOULD MINIMIZING DISTANCE TO THE NEAREST SAMPLES CAUSE OVERFITTING?
text: When the model has infinite capacity, minimizing distance from data examples to their nearest samples will lead to a model distribution that memorizes data examples. The same is true if we maximize likelihood. Likewise, minimizing any divergence measure will lead to memorization of data examples, since the minimum divergence is zero and by definition, this can only happen if the model distribution is the same as the empirical data distribution, whose support is confined to the set of data examples. This implies that whenever we have a finite number of data examples, any method that learns a model with infinite capacity will memorize the data examples and will hence overfit.
To get around this, most methods learn a parametric model with finite capacity. In the parametric setting, the minimum divergence is not necessarily zero; the same is true for the minimum distance from data examples to their nearest samples. Therefore, the optimum of these objective functions is not necessarily a model distribution that memorizes data examples, and so overfitting will not necessarily occur.

heading: 6.3 DOES DISJOINT SUPPORT BREAK MAXIMUM LIKELIHOOD?
text: Arjovsky et al. (2017) observes that the data distribution and the model distribution are supported on low-dimensional manifolds and so they are unlikely to have a non-negligible intersection. They point out DKL (pdata ‖pθ ) would be infinite in this case, or equivalently, the likelihood would be zero. While this does not invalidate the theoretical soundness of maximum likelihood, since the maximum of a non-negative function that is zero almost everywhere is still well-defined, it does cause a lot of practical issues for gradient-based learning, as the gradient is zero almost everywhere. This is believed to be one reason that models like variational autoencoders (Kingma & Welling, 2013; Rezende et al., 2014) use a Gaussian distribution with high variance for the conditional likelihood/observation model rather than a distribution close to the Dirac delta, so that the support of the model distribution is broadened to cover all the data examples (Arjovsky et al., 2017).
This issue does not affect our method, as our loss function is different from the log-likelihood function, even though their optima are the same (under some conditions). As the result, the gradients of our loss function are different from those of log-likelihood. When the supports of the data distribution and the model distribution do not overlap, each data example is likely far away from its nearest sample and so the gradient is large. Moreover, the farther the data examples are from the samples, the larger the gradient gets. Therefore, even when the gradient of log-likelihood can be tractably computed, there may be situations when the proposed method would work better than maximizing likelihood directly.

heading: 7 CONCLUSION
text: We presented a simple and versatile method for parameter estimation when the form of the likelihood is unknown. The method works by drawing samples from the model, finding the nearest sample to every data example and adjusting the parameters of the model so that it is closer to the data example. We showed that performing this procedure is equivalent to maximizing likelihood under some conditions. The proposed method can capture the full diversity of the data and avoids common issues like mode collapse, vanishing gradients and training instability. The method combined with vanilla model architectures is able to achieve encouraging results on MNIST, TFD and CIFAR-10.

heading: 8 SUPPLEMENTARY MATERIAL
text: Before proving the main result, we first prove the following intermediate results: Lemma 1. Let Ω ⊆ Rd and V ⊆ R. For i ∈ [N ], let fi : Ω → V be differentiable on Ω and Φ : V → R be differentiable on V and strictly increasing. Assume arg minθ∈Ω ∑N i=1 fi(θ)
exists and is unique. Let θ∗ := arg minθ∈Ω ∑N i=1 fi(θ) and wi := 1/Φ ′(fi(θ ∗)). If the following conditions hold:
• There is a bounded set S ⊆ Ω such that bd(S) ⊆ Ω, θ∗ ∈ S and ∀fi, ∀θ ∈ Ω\S, fi(θ) > fi(θ ∗), where bd(S) denotes the boundary of S. • For all θ ∈ Ω, if θ 6= θ∗, there exists j ∈ [d] such that〈 w1Φ ′(f1(θ))
... wnΦ ′(fn(θ))
 ,  ∂f1/∂θj(θ)...
∂fn/∂θj(θ) 〉 6= 0. Then arg minx∈Ω ∑N i=1 wiΦ(fi(θ)) exists and is unique. Furthermore,
arg minθ∈Ω ∑N i=1 wiΦ(fi(θ)) = arg minθ∈Ω ∑N i=1 fi(θ).
Proof. Let S ⊆ Ω be the bounded set such that bd(S) ⊆ Ω, θ∗ ∈ S and ∀fi, ∀θ ∈ Ω \ S, fi(θ) > fi(θ
∗). Consider the closure of S := S ∪ bd(S), denoted as S̄. Because S ⊆ Ω and bd(S) ⊆ Ω, S̄ ⊆ Ω. Since S is bounded, S̄ is bounded. Because S̄ ⊆ Ω ⊆ Rd and is closed and bounded, it is compact.
Consider the function ∑N i=1 wiΦ(fi(·)). By the differentiability of fi’s and Φ, ∑N i=1 wiΦ(fi(·)) is differentiable on Ω and hence continuous on Ω. By the compactness of S̄ and the continuity of ∑N i=1 wiΦ(fi(·)) on S̄ ⊆ Ω, Extreme Value Theorem applies, which im-
plies that minθ∈S̄ ∑N i=1 wiΦ(fi(θ)) exists. Let θ̃ ∈ S̄ be such that ∑N i=1 wiΦ(fi(θ̃)) =
minθ∈S̄ ∑N i=1 wiΦ(fi(θ)).

heading: By definition of S, ∀fi, ∀θ ∈ Ω \ S, fi(θ) > fi(θ∗), implying that Φ(fi(θ)) > Φ(fi(θ∗))
text: since Φ is strictly increasing. Because Φ′(·) > 0, wi > 0 and so ∑N i=1 wiΦ(fi(θ)) >∑N
i=1 wiΦ(fi(θ ∗)) ∀θ ∈ Ω \ S. At the same time, since θ∗ ∈ S ⊂ S̄, by definition of θ̃,∑N i=1 wiΦ(fi(θ̃)) ≤ ∑N i=1 wiΦ(fi(θ ∗)). Combining these two facts yields ∑N i=1 wiΦ(fi(θ̃)) ≤∑N
i=1 wiΦ(fi(θ ∗)) < ∑N i=1 wiΦ(fi(θ)) ∀θ ∈ Ω \ S. Since the inequality is strict, this implies that
θ̃ /∈ Ω \ S, and so θ̃ ∈ S̄ \ (Ω \ S) ⊆ Ω \ (Ω \ S) = S.
In addition, because θ̃ is the minimizer of ∑N i=1 wiΦ(fi(·)) on S̄, ∑N i=1 wiΦ(fi(θ̃)) ≤∑N
i=1 wiΦ(fi(θ)) ∀θ ∈ S̄. So, ∑N i=1 wiΦ(fi(θ̃)) ≤ ∑N i=1 wiΦ(fi(θ)) ∀θ ∈ S̄ ∪
(Ω \ S) ⊇ S ∪ (Ω \ S) = Ω. Hence, θ̃ is a minimizer of ∑N i=1 wiΦ(fi(·)) on Ω, and so
minθ∈Ω ∑N i=1 wiΦ(fi(θ)) exists. Because ∑N i=1 wiΦ(fi(·)) is differentiable on Ω, θ̃ must be a
critical point of ∑N i=1 wiΦ(fi(·)) on Ω.
On the other hand, since Φ is differentiable on V and fi(θ) ∈ V for all θ ∈ Ω, Φ′(fi(θ)) exists for all θ ∈ Ω. So,
∇ ( N∑ i=1 wiΦ(fi(θ)) ) = N∑ i=1 wi∇ (Φ(fi(θ)))
= N∑ i=1 wiΦ ′(fi(θ))∇fi(θ)
= N∑ i=1 Φ′(fi(θ)) Φ′(fi(θ∗)) ∇fi(θ)
At θ = θ∗,
∇ ( N∑ i=1 wiΦ(fi(θ ∗)) ) = N∑ i=1 Φ′(fi(θ ∗)) Φ′(fi(θ∗)) ∇fi(θ∗)
= N∑ i=1 ∇fi(θ∗)
Since each fi is differentiable on Ω, ∑N i=1 fi is differentiable on Ω. Combining this with the fact
that θ∗ is the minimizer of ∑N i=1 fi on Ω, it follows that ∇ (∑N i=1 fi(θ ∗) ) = ∑N i=1∇fi(θ∗) = 0.
Hence, ∇ (∑N i=1 wiΦ(fi(θ ∗)) ) = 0 and so θ∗ is a critical point of ∑N i=1 wiΦ(fi(·)).
Because ∀θ ∈ Ω, if θ 6= θ∗, ∃j ∈ [d] such that 〈 w1Φ ′(f1(θ))
... wnΦ ′(fn(θ))
 ,  ∂f1/∂θj(θ)...
∂fn/∂θj(θ) 〉 6= 0, ∑N i=1 wiΦ ′(fi(θ))∇fi(θ) = ∇ (∑N i=1 wiΦ(fi(θ)) ) 6= 0 for any θ 6= θ∗ ∈ Ω. Therefore, θ∗
is the only critical point of ∑N i=1 wiΦ(fi(·)) on Ω. Since θ̃ is a critical point on Ω, we can con-
clude that θ∗ = θ̃, and so θ∗ is a minimizer of ∑N i=1 wiΦ(fi(·)) on Ω. Since any other minimizer must be a critical point and θ∗ is the only critical point, θ∗ is the unique minimizer. So, arg minθ∈Ω ∑N i=1 fi(θ) = θ ∗ = arg minθ∈Ω ∑N i=1 wiΦ(fi(θ)).
Lemma 2. Let P be a distribution on Rd whose density p(·) is continuous at a point x0 ∈ Rd and x ∼ P be a random variable. Let r̃ := ‖x− x0‖2, κ := πd/2/Γ ( d 2 + 1 ) , where Γ(·) denotes the gamma function 1, and r := κr̃d. Let G(·) denote the cumulative distribution function (CDF) of r and ∂+G(·) denote the one-sided derivative of G from the right. Then, ∂+G(0) = p(x0).
Proof. By definition of ∂+G(·),
∂+G(0) = lim h→0+ G(h)−G(0) h = lim h→0+ G(h) h
= lim h→0+ Pr (r ≤ h) h = lim h→0+
Pr ( r̃ ≤ d √ h/κ ) h
1The constant κ is the the ratio of the volume of a d-dimensional ball of radius r̃ to a d-dimensional cube of side length r̃.
If we define h̃ := d √ h/κ, the above can be re-written as:
∂+G(0) = lim h̃→0+
Pr ( r̃ ≤ h̃ ) κh̃d = lim h̃→0+ ∫ Bx0 (h̃) p(u)du κh̃d
We want to show that limh̃→0+ (∫
Bx0 (h̃) p(u)du
) /κh̃d = p(x0). In other words, we want to show
∀ > 0 ∃δ > 0 such that ∀h̃ ∈ (0, δ), ∣∣∣∣ ∫Bx0 (h̃) p(u)duκh̃d − p(x0) ∣∣∣∣ < . Let > 0 be arbitrary.
Since p(·) is continuous at x0, by definition, ∀̃ > 0 ∃δ̃ > 0 such that ∀u ∈ Bx0(δ̃), |p(u)− p(x0)| < ̃. Let δ̃ > 0 be such that ∀u ∈ Bx0(δ̃), p(x0) − < p(u) < p(x0) + . We choose δ = δ̃.
Let 0 < h̃ < δ be arbitrary. Since p(x0)− < p(u) < p(x0)+ ∀u ∈ Bx0(δ̃) = Bx0(δ) ⊃ Bx0(h̃),∫ Bx0 (h̃) p(u)du < ∫ Bx0 (h̃) (p(x0) + ) du
= (p(x0) + ) ∫ Bx0 (h̃) du
Observe that ∫ Bx0 (h̃) du is the volume of a d-dimensional ball of radius h̃, so ∫ Bx0 (h̃) du = κh̃d.
Thus, ∫ Bx0 (h̃) p(u)du < κh̃d (p(x0) + ), implying that (∫ Bx0 (h̃) p(u)du ) /κh̃d < p(x0) + . By
similar reasoning, we conclude that (∫
Bx0 (h̃) p(u)du
) /κh̃d > p(x0)− .
Hence,
∣∣∣∣∣ ∫ Bx0 (h̃) p(u)du κh̃d − p(x0) ∣∣∣∣∣ < ∀h̃ ∈ (0, δ) Therefore,
∂+G(0) = lim h̃→0+
∫ Bx0 (h̃) p(u)du
κh̃d = p(x0)
Lemma 3. Let Pθ be a parameterized family of distributions on Rd with parameter θ and probability density function (PDF) pθ(·) that is continuous at a point xi. Consider a random variable x̃θ1 ∼ Pθ and define r̃θi := ∥∥x̃θ1 − xi∥∥22, whose cumulative distribution function (CDF) is denoted by F θi (·). Assume Pθ has the following property: for any θ1, θ2, there exists θ0 such that F θ0i (t) ≥ max { F θ1i (t), F θ2 i (t) } ∀t ≥ 0 and pθ0(xi) = max {pθ1(xi), pθ2(xi)}. For any m ≥ 1,
let x̃θ1, . . . , x̃ θ m ∼ Pθ be i.i.d. random variables and define Rθi := minj∈[m] ∥∥x̃θj − xi∥∥22. Then the function Ψi : z 7→ minθ { E [ Rθi ] |pθ(xi) = z } is strictly decreasing.
Proof. Let rθi := κ ( r̃θi )d/2 = κ ∥∥x̃θ1 − xi∥∥d2 be a random variable and let Gθi (·) be the CDF of rθi . Since Rθi is nonnegative,
E [ Rθi ] = ∫ ∞ 0 Pr ( Rθi > t ) dt
= ∫ ∞ 0 ( Pr (∥∥x̃θ1 − xi∥∥22 > t))m dt
= ∫ ∞ 0 ( Pr ( κ ∥∥x̃θ1 − xi∥∥d2 > κtd/2))m dt
= ∫ ∞ 0 ( Pr ( rθi > κt d/2 ))m dt
= ∫ ∞ 0 ( 1−Gθi ( κtd/2 ))m dt
Also, by Lemma 2, pθ(xi) = ∂+Gθi (0). Using these facts, we can rewrite minθ { E [ Rθi ] |pθ(xi) = z } as minθ {∫∞ 0 ( 1−Gθi ( κtd/2 ))m dt ∣∣∂+Gθi (0) = z}. By definition
of Ψi, minθ {∫∞
0
( 1−Gθi ( κtd/2 ))m dt ∣∣∂+Gθi (0) = z} exists for all z. Let φi(z) be a value of
θ that attains the minimum. Define G∗i (y, z) := G φi(z) i (y). By definition, ∂+ ∂yG ∗ i (0, z) = z, where ∂+ ∂yG ∗ i (y, z) denotes the one-sided partial derivative from the right w.r.t. y. Also, since G ∗ i (·, z) is the CDF of a distribution of a non-negative random variable, G∗i (0, z) = 0.
By definition of ∂+∂yG ∗ i (0, z), ∀ > 0 ∃δ > 0 such that ∀h ∈ (0, δ), ∣∣∣G∗i (h,z)−G∗i (0,z)h − z∣∣∣ < . Let z′ > z. Let δ > 0 be such that ∀h ∈ (0, δ),
∣∣∣G∗i (h,z)−G∗i (0,z)h − z∣∣∣ < z′−z2 and δ′ > 0 be such that ∀h ∈ (0, δ′),
∣∣∣G∗i (h,z′)−G∗i (0,z′)h − z′∣∣∣ < z′−z2 . Consider h ∈ (0,min(δ, δ′)). Then, G ∗ i (h,z)−G ∗ i (0,z)
h = G∗i (h,z) h < z + z′−z 2 = z+z′
2 and G∗i (h,z ′)−G∗i (0,z ′)
h = G∗i (h,z ′) h > z ′ − z ′−z 2 = z+z′
2 . So, G∗i (h, z)
h < z + z′ 2 < G∗i (h, z ′) h
Multiplying by h on both sides, we conclude that G∗i (h, z) < G ∗ i (h, z ′) ∀h ∈ (0,min(δ, δ′)). Let α := d √ min(δ, δ′)/κ. We can break ∫∞
0
( 1−G∗i ( κtd/2, z ))m dt into two terms:∫ ∞
0
( 1−G∗i ( κtd/2, z ))m dt
= ∫ α 0 ( 1−G∗i ( κtd/2, z ))m dt+ ∫ ∞ α ( 1−G∗i ( κtd/2, z ))m dt
We can also do the same for ∫∞
0
( 1−G∗i ( κtd/2, z′ ))m dt.
Because G∗i (h, z) < G ∗ i (h, z ′) ∀h ∈ (0,min(δ, δ′)), G∗i (κtd/2, z) < G∗i (κtd/2, z′) ∀t ∈ (0, α). It follows that 1 − G∗i (κtd/2, z) > 1 − G∗i (κtd/2, z′) and ( 1−G∗i ( κtd/2, z ))m >(
1−G∗i (κtd/2, z′) )m ∀t ∈ (0, α). So, ∫ α
0
( 1−G∗i ( κtd/2, z ))m dt >∫ α
0
( 1−G∗i ( κtd/2, z′ ))m dt.
We now consider the second term. First, observe that F θi (t) = Pr (∥∥x̃θ1 − xi∥∥22 ≤ t) =
Pr ( κ ∥∥x̃θ1 − xi∥∥d2 ≤ κtd/2) = Gθi (κtd/2) for all t ≥ 0. So, by the property of Pθ, for
any θ1, θ2, there exists θ0 such that Gθ0i (κt d/2) = F θ0i (t) ≥ max { F θ1i (t), F θ2 i (t) } =
max { Gθ1i (κt d/2), Gθ2i (κt d/2) } ∀t ≥ 0 and ∂+Gθ0i (0) = pθ0(xi) = max {pθ1(xi), pθ2(xi)} =
max { ∂+G θ1 i (0), ∂+G θ2 i (0) } .
Take θ1 = φi(z) and θ2 = φi(z′). Let θ0 be such that Gθ0i (κt d/2) ≥ max { Gθ1i (κt d/2), Gθ2i (κt d/2) } ∀t ≥ 0 and ∂+Gθ0i (0) = max { ∂+G θ1 i (0), ∂+G θ2 i (0) } . By
definition of φi(·), ∂+Gθ1i (0) = z and ∂+G θ2 i (0) = z ′. So, ∂+Gθ0i (0) = max {z, z′} = z′. Since Gθ0i (κt d/2) ≥ Gθ2i (κtd/2) ∀t ≥ 0, 1 − G θ0 i ( κtd/2 ) ≤ 1 − Gθ2i ( κtd/2 ) ∀t ≥ 0
and so ∫∞
0
( 1−Gθ0i ( κtd/2 ))m dt ≤ ∫∞ 0 ( 1−Gθ2i ( κtd/2 ))m dt. On the other hand, because
θ2 = φi(z ′) minimizes ∫∞ 0 ( 1−Gθi ( κtd/2 ))m dt among all θ’s such that ∂+Gθi (0) = z ′ and ∂+G θ0 i (0) = z ′, ∫∞
0
( 1−Gθ2i ( κtd/2 ))m dt ≤ ∫∞ 0 ( 1−Gθ0i ( κtd/2 ))m dt. We can therefore
conclude that ∫∞
0
( 1−Gθ0i ( κtd/2 ))m dt = ∫∞ 0 ( 1−Gθ2i ( κtd/2 ))m dt. Since 1−Gθ0i ( κtd/2 ) ≤
1 − Gθ2i ( κtd/2 ) ∀t ≥ 0, the only situation where this can happen is when Gθ0i ( κtd/2 ) =
Gθ2i ( κtd/2 ) ∀t ≥ 0.
By definition of G∗i , G ∗ i
( κtd/2, z ) = G
φi(z) i (κt d/2) = Gθ1i (κt d/2) and G∗i
( κtd/2, z′ ) =
G φi(z ′) i (κt d/2) = Gθ2i (κt d/2) = Gθ0i ( κtd/2 ) . By definition of θ0, Gθ0i ( κtd/2 ) ≥ Gθ1i (κtd/2) ∀t ≥
0. So, G∗i ( κtd/2, z′ ) = Gθ2i (κt d/2) ≥ Gθ1i (κtd/2) = G∗i ( κtd/2, z ) ∀t ≥ 0. Hence,∫∞
α
( 1−G∗i ( κtd/2, z′ ))m dt ≤ ∫∞ α ( 1−G∗i ( κtd/2, z ))m dt.
Combining with the previous result that ∫ α
0
( 1−G∗i ( κtd/2, z′ ))m dt <∫ α
0
( 1−G∗i ( κtd/2, z ))m dt, it follows that:∫ ∞
0
( 1−G∗i ( κtd/2, z′ ))m dt
= ∫ α 0 ( 1−G∗i ( κtd/2, z′ ))m dt+ ∫ ∞ α ( 1−G∗i ( κtd/2, z′ ))m dt
< ∫ α 0 ( 1−G∗i ( κtd/2, z ))m dt+ ∫ ∞ α ( 1−G∗i ( κtd/2, z′ ))m dt
≤ ∫ α 0 ( 1−G∗i ( κtd/2, z ))m dt+ ∫ ∞ α ( 1−G∗i ( κtd/2, z ))m dt
= ∫ ∞ 0 ( 1−G∗i ( κtd/2, z ))m dt
By definition, ∫ ∞ 0 ( 1−G∗i ( κtd/2, z ))m dt
= ∫ ∞ 0 ( 1−Gφi(z)i (κt d/2) )m dt
= min θ {∫ ∞ 0 ( 1−Gθi ( κtd/2 ))m dt ∣∣∂+Gθi (0) = z}
= min θ
{ E [ Rθi ] |pθ(xi) = z } = Ψi(z)
Similarly, ∫∞
0
( 1−G∗i ( κtd/2, z′ ))m dt = Ψi(z
′). We can therefore conclude that Ψi(z′) < Ψi(z) whenever z′ > z.
We now prove the main result. Theorem 1. Consider a set of observations x1, . . . ,xn, a parameterized family of distributions Pθ with probability density function (PDF) pθ(·) and a unique maximum likelihood solution θ∗. For any m ≥ 1, let x̃θ1, . . . , x̃θm ∼ Pθ be i.i.d. random variables and define r̃θ :=
∥∥x̃θ1∥∥22, Rθ := minj∈[m]
∥∥x̃θj∥∥22 andRθi := minj∈[m] ∥∥x̃θj − xi∥∥22. Let F θ(·) be the cumulative distribution function (CDF) of r̃θ and Ψ(z) := minθ { E [ Rθ ] |pθ(0) = z } .
If Pθ satisfies the following:
• pθ(x) is differentiable w.r.t. θ and continuous w.r.t. x everywhere.
• ∀θ,v, there exists θ′ such that pθ(x) = pθ′(x + v) ∀x. • For any θ1, θ2, there exists θ0 such that F θ0(t) ≥ max { F θ1(t), F θ2(t) } ∀t ≥ 0 and
pθ0(0) = max {pθ1(0), pθ2(0)}.
• ∃τ > 0 such that ∀i ∈ [n] ∀θ /∈ Bθ∗(τ), pθ(xi) < pθ∗(xi), where Bθ∗(τ) denotes the ball centred at θ∗ of radius τ .
• Ψ(z) is differentiable everywhere.
• For all θ, if θ 6= θ∗, there exists j ∈ [d] such that〈 Ψ′(pθ(x1))pθ(x1) Ψ′(pθ∗ (x1))pθ∗ (x1) ...
Ψ′(pθ(xn))pθ(xn) Ψ′(pθ∗ (xn))pθ∗ (xn)
 ,  ∇θ (log pθ(x1))j... ∇θ (log pθ(xn))j 〉 6= 0. Then,
arg min θ n∑ i=1
E [ Rθi ]
Ψ′(pθ∗(xi))pθ∗(xi) = arg max θ n∑ i=1 log pθ(xi)
Furthermore, if pθ∗(x1) = · · · = pθ∗(xn), then,
arg min θ n∑ i=1 E [ Rθi ] = arg max θ n∑ i=1 log pθ(xi)
Proof. Pick an arbitrary i ∈ [n]. We first prove a few basic facts. By the second property of Pθ, ∀θ ∃θ′ such that pθ(u) = pθ′(u − xi) ∀u. In particular, pθ(xi) = pθ′(xi − xi) = pθ′(0). Let F θi be as defined in Lemma 3.
F θi (t) = Pr ( r̃θi ≤ t ) = Pr (∥∥x̃θ1 − xi∥∥2 ≤ √t) =
∫ Bxi ( √ t) pθ(u)du = ∫ Bxi ( √ t) pθ′(u− xi)du
= ∫ B0( √ t) pθ′(u)du = Pr ( r̃θ ′ ≤ t ) = F θ ′ (t)
Similarly, ∀θ′ ∃θ such that pθ′(u) = pθ(u + xi) ∀u. In particular, pθ′(0) = pθ(0 + xi) = pθ(xi).
F θ ′ (t) = Pr ( r̃θ ′ ≤ t ) = ∫ B0( √ t) pθ′(u)du
= ∫ B0( √ t) pθ(u + xi)du = ∫ Bxi ( √ t) pθ(u)du
=Pr (∥∥x̃θ1 − xi∥∥2 ≤ √t) = Pr (r̃θi ≤ t) = F θi (t)
Let θ1, θ2 be arbitrary. The facts above imply that there exist θ′1 and θ ′ 2 such that F θ1 i (t) = F θ′1(t), F θ2i (t) = F θ′2(t), pθ1(xi) = pθ′1(0) and pθ2(xi) = pθ′2(0).
By the third property of Pθ, let θ′0 be such that F θ′0(t) ≥ max { F θ ′ 1(t), F θ ′ 2(t) } ∀t ≥ 0 and
pθ′0(0) = max { pθ′1(0), pθ′2(0) } . By the facts above, it follows that there exists θ0 such that F θ ′ 0(t) = F θ0i (t) and pθ′0(0) = pθ0(xi).
So, we can conclude that for any θ1, θ2, there exists θ0 such that F θ0i (t) ≥ max { F θ1i (t), F θ2 i (t) } ∀t ≥ 0 and pθ0(xi) = max {pθ1(xi), pθ2(xi)}.
By Lemma 3, Ψi(z) = minθ { E [ Rθi ] |pθ(xi) = z } is strictly decreasing.
Consider any θ. By the facts above, there exists θ′ such that pθ(xi) = pθ′(0) and F θi (t) = F θ′(t) ∀t. Therefore,
E [ Rθi ] = ∫ ∞ 0 Pr ( Rθi > t ) dt
= ∫ ∞ 0 ( Pr (∥∥x̃θ1 − xi∥∥22 > t))m dt
= ∫ ∞ 0 ( 1− F θi (t) )m dt
= ∫ ∞ 0 ( 1− F θ ′ (t) )m dt
= ∫ ∞ 0 Pr ( Rθ ′ > t ) dt
=E [ Rθ ′ ]
So, ∀z
Ψi(z) = min θ
{ E [ Rθi ] |pθ(xi) = z } = min
θ′
{ E [ Rθ ′ ] |pθ′(0) = z } =Ψ(z)
Because Ψi(·) is strictly decreasing, Ψ(·) is also strictly decreasing. We would like to apply Lemma 1, with fi(θ) = − log pθ(xi) ∀i ∈ [n] and Φ(y) = Ψ(exp(−y)). By the first property of Pθ, pθ(·) is differentiable w.r.t. θ and so fi(θ) is differentiable for all i. By the fifth property of Pθ, Ψ(·) is differentiable and so Φ(·) is differentiable. Since y 7→ exp(−y) is strictly decreasing and Ψ(·) is strictly decreasing, Φ(·) is strictly increasing. Since there is a unique maximum likelihood solution θ∗, minθ ∑n i=1 fi(θ) = maxθ ∑n i=1 log pθ(xi) exists and has a unique minimizer. By the fourth property of Pθ, the first condition of Lemma 1 is satisfied. By the sixth property of Pθ, the second condition of Lemma 1 is satisfied. Since all conditions are satisfied, we apply Lemma 1 and conclude that
min θ n∑ i=1 wiΦ(fi(θ)) = min θ n∑ i=1 wiΨ(pθ(xi))
= min θ n∑ i=1 wiΨi(pθ(xi))
= min θ n∑ i=1
E [ Rθi ]
Ψ′(pθ∗(xi))pθ∗(xi)
exists and has a unique minimizer. Furthermore,
arg min θ n∑ i=1
E [ Rθi ]
Ψ′(pθ∗(xi))pθ∗(xi) = arg min θ n∑ i=1 − log pθ(xi)
= arg max θ n∑ i=1 log pθ(xi)
If pθ(x1) = · · · pθ(xn), then w1 = · · · = wn, and so arg minθ ∑n i=1 wiE [ Rθi ] =
arg minθ ∑n i=1 E [ Rθi ] = arg maxθ ∑n i=1 log pθ(xi).

Review summaries:
The paper proposes a new algorithm for implicit maximum likelihood estimation based on a fast Nearest Neighbor search. The algorithm can be used to implicitly maximize the likelihood of models for which the former quantity is not intractable but for which sampling is easy which is typically the case for implicit models.  The paper shows that under some conditions the optimal solution of the algorithm corresponds to the MLE solution and provides some experimental evidence that the method leads to a higher likelihood. However, The paper lacks clarity and the experiments are not really convincing. Here are some remarks:

This paper proposes a nearest-neighbor-based algorithm for implicit maximum likelihood.  Samples are produced by the generator network and then a nearest neighbors algorithm is run to match the samples with their nearest data point.  The generator is then updated using the Euclidean distance between samples and neighbors as the optimization objective.  Six conditions are then provided, and if they are met, then the authors show that this method is performing maximum likelihood on the implied density.  Experiments report Parzen window density estimates, samples from the model, and latent-space interpolations for MNIST, Toronto Faces, and CIFAR-10.