Year: 2019

Abstract:
Modelling statistical relationships beyond the conditional mean is crucial in many settings. Conditional density estimation (CDE) aims to learn the full conditional probability density from data. Though highly expressive, neural network based CDE models can suffer from severe over-fitting when trained with the maximum likelihood objective. Due to the inherent structure of such models, classical regularization approaches in the parameter space are rendered ineffective. To address this issue, we develop a model-agnostic noise regularization method for CDE that adds random perturbations to the data during training. We demonstrate that the proposed approach corresponds to a smoothness regularization and prove its asymptotic consistency. In our experiments, noise regularization significantly and consistently outperforms other regularization methods across seven data sets and three CDE models. The effectiveness of noise regularization makes neural network based CDE the preferable method over previous nonand semi-parametric approaches, even when training data is scarce.

heading: 1 INTRODUCTION
text: While regression analysis aims to describe the conditional mean E[y|x] of a response y given inputs x, many problems such as risk management and planning under uncertainty require gaining insight about deviations from the mean and their associated likelihood. The stochastic dependency of y on x can be captured by modeling the conditional probability density p(y|x). Inferring such a density function from a set of empirical observations {(xn, yn)}Nn=1 is typically referred to as conditional density estimation (CDE) and is the focus of this paper.
In the recent machine learning literature, there has been a resurgence of interest in high-capacity density models based on neural networks (Dinh et al., 2017; Ambrogioni et al., 2017; Kingma & Dhariwal, 2018). Since this line of work mainly focuses on the modelling of images based on large scale data sets, over-fitting and noisy observations are of minor concern in this context. In contrast, we are interested in CDE in settings where data may be scarce and noisy. When combined with maximum likelihood estimation, the flexibility of such high-capacity models results in over-fitting and poor generalization. While regression typically assumes Gaussian conditional noise, CDE uses expressive distribution families to model deviations from the conditional mean. Hence, the overfitting problem tends to be even more severe in CDE than in regression.
Classical regularization of the neural network weights such as weight decay (Pratt & Hanson, 1989) has been shown to be effective for regression and classification. However, in the context of CDE, the output of the neural network merely controls the parameters of a density model such as a Gaussian Mixture or Normalizing Flow. This makes the standard regularization methods in the parameter space less effective and harder to analyze.
Aiming to address this issue, we propose and analyze noise regularization, a method well-studied in the context of regression and classification, for the purpose of conditional density estimation. In that, the paper attempts to close a gap in previous research. By adding small random perturbations to the data during training, the conditional density estimate is smoothed and tends to generalize better. In fact, we show that adding noise during maximum likelihood estimation is equivalent to penalizing the second derivatives of the conditional log-probability. Visually, the respective regularization term punishes very curved or even spiky density estimators in favor of smoother variants, which proves to be a favorable inductive bias in many applications. Moreover, under some regularity conditions, we show that the proposed regularization scheme is asymptotically consistent, converging to the unbiased maximum likelihood estimator. This does not only support the soundness of the proposed
method but also endows us with useful insight in how to set the regularization intensity relative to the data dimensionality and training set size.
Overall, the proposed noise regularization scheme is easy to implement and agnostic to the parameterization of the CDE model. We empirically demonstrate its effectiveness on three different neural network based models. The experimental results show that noise regularization outperforms other regularization methods significantly and consistently across various data sets. Finally, we demonstrate that, when properly regularized, neural network based CDE is able to improve upon state-of-the art non-parametric estimators, even when only 400 training observations are available.

heading: 2 BACKGROUND
text: Density Estimation. Let X be a random variable with probability density function (PDF) p(x) defined over the domain X ⊆ Rdx . Given a collection D = {x1, ..., xn} of observations sampled from p(x), the goal is to find a good estimate f̂(x) of the true density function p. In parametric estimation, the PDF f̂ is assumed to belong to a parametric family F = {f̂θ(·)|θ ∈ Θ} where the density function is described by a finite dimensional parameter θ ∈ Θ. The standard method for estimating θ is maximum likelihood estimation, wherein θ is chosen so that the likelihood of the data D is maximized. This is equivalent to minimizing the Kullback-Leibler divergence between the empirical data distribution pD(x) = 1n ∑n i=1 δ(||x − xi||) (i.e., mixture of point masses in the observations xi) and the parametric distribution f̂θ:
θ∗ = arg max θ∈Θ n∑ i=1 log f̂θ(xi) = arg min θ∈Θ DKL(pD||f̂θ) (1)
From a geometric perspective, (1) can be viewed as an orthogonal projection of pD(x) onto F w.r.t. the reverse KL-divergence. Hence, (1) is also commonly referred to as an M-projection (Murphy, 2012; Nielsen, 2018). In contrast, non-parametric density estimators make implicit smoothness assumptions through a kernel function. The most popular non-parametric method, kernel density estimation (KDE), places a symmetric density function K(z), the so-called kernel, on each training data point xn (Rosenblatt, 1956; Parzen, 1962). The resulting density estimate reads as q̂(x) =
1 nhd ∑n i=1K ( x−xi h ) . One popular choice of K(·) is a Gaussian K(z) = (2π)− d2 exp (− 12z
2). Beyond the appropriate choice of K(·), a central challenge is the selection of the bandwidth parameter h which controls the smoothness of the estimated PDF (Li & Racine, 2007).
Conditional Density Estimation (CDE). Let (X,Y ) be a pair of random variables with respective domains X ⊆ Rdx and Y ⊆ Rdy and realizations x and y. Let p(y|x) = p(x, y)/p(x) denote the conditional probability density of y given x. Typically, Y is referred to as a dependent variable (explained variable) and X as conditional (explanatory) variable. Given a dataset of observations D = {(xn, yn)}Nn=1 drawn from the joint distribution (xn, yn) ∼ p(x, y), the aim of conditional density estimation (CDE) is to find an estimate f̂(y|x) of the true conditional density p(y|x). In the context of CDE, the KL-divergence objective is expressed as expectation over p(x):
Ex∼p(x) [ DKL(p(y|x)||f̂(y|x)) ] = E(x,y)∼p(x,y) [ log p(y|x)− log f̂(y|x) ] (2)
Corresponding to (1), we refer to the minimization of (2) w.r.t. θ as conditional M-projection. Given a dataset D drawn i.i.d. from p(x, y), the conditional MLE following from (2) can be stated as
θ∗ = arg min θ − n∑ i=1 log f̂θ(yi|xi) (3)

heading: 3 RELATED WORK
text: The first part of this section discusses relevant work in the field of CDE, focusing on high-capacity models that make little prior assumptions. The second part relates our approach to previous regularization and data augmentation methods.
Non-parametric CDE. A vast body of literature in statistics and econometrics studies nonparametric kernel density estimators (KDE) (Rosenblatt, 1956; Parzen, 1962) and the associated bandwidth selection problem, which concerns choosing the appropriate amount of smoothing (Silverman,
1982; Hall et al., 1992; Cao et al., 1994). To estimate conditional probabilities, previous work proposes to estimate both the joint and marginal probability separately with KDE and then computing the conditional probability as their ratio (Hyndman et al., 1996; Li & Racine, 2007). Other approaches combine non-parametric elements with parametric elements (Tresp, 2001; Sugiyama & Takeuchi, 2010; Dutordoir et al., 2018). Despite their theoretical appeal, non-parametric density estimators suffer from poor generalization in regions where data is sparse (e.g., tail regions), causing rapid performance deterioration as the data dimensionality increases (Scott & Wand, 1991).
CDE based on neural networks. Most work in machine learning focuses on flexible parametric function approximators for CDE. In our experiments, we use the work of Bishop (1994) and Ambrogioni et al. (2017), who propose to use a neural network to control the parameters of a mixture density model. A recent trend in machine learning are latent density models such as cGANs (Mirza & Osindero, 2014) and cVAEs (Sohn et al., 2015). Although such methods have been shown successful for estimating distributions of images, the probability density function (PDF) of such models is intractable. More promising in this sense are normalizing flows (Rezende & Mohamed, 2015; Dinh et al., 2017; Trippe & Turner, 2018), since they provide the PDF in tractable form. We employ a neural network controlling the parameters of a normalizing flow as our third CDE model to showcase the empirical efficacy of our regularization approach.
Regularization. Since neural network based CDE models suffer from severe over-fitting when trained with the MLE objective, they require proper regularization. Classical regularization of the parameters such as weight decay (Pratt & Hanson, 1989; Krogh & Hertz, 1992; Nowlan & Hinton, 1992), l1/l2-penalties (Mackay, 1992; Ng, 2004) and Bayesian priors (Murray & Edwards, 1993; Hinton & Van Camp, 1993) have been shown to work well in the regression and classification setting. However, in the context of CDE, it is less clear what kind of inductive bias such a regularization imposes on the density estimate. In contrast, our regularization approach is agnostic w.r.t. parametrization and is shown to penalize strong variations of the log-density function.
Regularization methods such as dropout are closely related to ensemble methods (Srivastava et al., 2014). Thus, they are orthogonal to our work and can be freely combined with noise regularization.
Adding noise during training. Adding noise during training is a common scheme that has been proposed in various forms. This includes noise on the neural network weights or activations (Wan et al., 2013; Srivastava et al., 2014; Gal & Uk, 2016) and additive noise on the gradients for scalable MCMC posterior inference (Welling & Teh, 2011; Chen et al., 2014). While this line of work corresponds to noise in the parameter space, other research suggests to augment the training data through random and/or adversarial transformations of the data (Sietsma & Dow, 1991; Burges & Schölkopf, 1996; Goodfellow et al., 2015; Yuan et al., 2017). Our approach transforms the training observations by adding small random perturbations. While this form of regularization has been studied in the context of regression and classification problems (Holmstrom & Koistinen, 1992a; Webb, 1994; Bishop, 1995; Natarajan et al., 2013; Maaten et al., 2013), this paper focuses on the regularization of CDE. In particular, we build on top of the results of Webb (1994) showing that training with noise corresponds to a penalty on strong variations of the log-density and extend previous consistency results for regression of Holmstrom & Koistinen (1992a) to the more general setting of CDE. To our best knowledge, this is also the first paper to evaluate the empirical efficacy of noise regularization for density estimation.

heading: 4 NOISE REGULARIZATION
text: When considering expressive families of conditional densities, standard maximum likelihood estimation of the model parameters θ is ill suited. As can be observed in Figure 1, simply minimizing the negative log-likelihood of the data leads to severe over-fitting and poor generalization beyond the training data. Hence, it is necessary to impose additional inductive bias, for instance, in the form of regularization. Unlike in regression or classification, the form of inductive bias imposed by popular regularization techniques such as weight decay (Krogh & Hertz, 1991; Kukačka et al., 2017) is less clear in the CDE setting, where the neural network weights often only indirectly control the probability density through a unconditional density model, e.g., a Gaussian Mixture.
We propose to add noise perturbations to the data points during the optimization of the log-likelihood objective. This can be understood as replacing the original data points (xi, yi) by random variables x̃i = xi + ξx and ỹi = yi + ξy where the perturbation vectors are sampled from noise distribu-
tions Kx(ξx) and Ky(ξy) respectively. Further, we choose the noise to be zero centered as well as identically and independently distributed among the data dimensions, with standard deviation h:
Eξ∼K(ξ) [ξ] = 0 and Eξ∼K(ξ) [ ξξ> ] = h2I (4)
This can be seen as data augmentation, where “synthetic” data is generated by randomly perturbing the original data. Since the supply of noise vectors is technically unlimited, an arbitrary large augmented data set can be generated by repetitively sampling data points from D, and adding a random perturbation vector to the respective data point. This procedure is formalized in Algorithm 1.
For notational brevity, we set Z := X × Y , z := (x>, y>)> and denote f̂θ(z) := f̂θ(y|x). The presented noise regularization approach is agnostic to whether we are concerned with unconditional or conditional MLE. Thus, the generic notation also allows us to generalize the results to both settings (derived in the remainder of the paper).
Algorithm 1 (Conditional) MLE with Noise Regularization - Generic Procedure
Require: D = {z1, ..., zn}, noise intensity h Require: number of perturbed samples r,
1: for j = 1 to r do 2: Select i ∈ {1, ..., n} with equal prob. 3: Draw perturbation ξ ∼ K 4: Set z̃j = zi + hξ 5: return arg minθ∈Θ− ∑r j=1 log f̂θ(z̃j)
Algorithm 2 (Conditional) MLE with Noise Regularization - Mini-Batch Gradient Descent
Require: D = {z1, ..., zn}, noise intensity h Require: learning rate α, mini-batch size m
1: Initialize θ 2: while θ not converged do 3: Sample minibatch {z1, ..., zm} ⊂ D 4: for j = 1 to m do 5: Draw perturbation ξ ∼ K 6: Set z̃j = zj + hξ 7: θ ← θ + α∇θ ∑m j=1 log f̂θ(z̃j) 8: return optimized parameter θ
When considering highly flexible parametric families such as Mixture Density Networks (MDNs) (Bishop, 1994), the maximum likelihood solution in line 5 of Algorithm 1 is no longer tractable. In such case, one typically resorts to numerical optimization techniques such as mini-batch gradient descent and variations thereof. In this context, the generic procedure in Algorithm 1 can be transformed into a simple extensions of mini-batch gradient descent on the MLE objective (see Algorithm 2). Specifically, each mini-batch is perturbed with i.i.d. noise before computing the MLE objective function (forward pass) and the respective gradients (backward pass).

heading: 4.1 VARIABLE NOISE AS SMOOTHNESS REGULARIZATION
text: Intuitively, the previously presented variable noise can be interpreted as “smearing” the data points during the maximum likelihood estimation. This alleviates the jaggedness of the density estimate arising from an un-regularized maximum likelihood objective in flexible density classes. We will now give this intuition a formal foundation, by mathematically analyzing the effect of the noise perturbations.
Before discussing the particular effects of randomly perturbing the data during conditional maximum likelihood estimation, we first analyze noise regularization in a more general case. Let l(D) be a loss function over a set of data points D = {z1, ..., zn}, which can be partitioned into a sum of losses l(D) = ∑n i=1 l(zi), corresponding to each data point zi: The expected loss l(zi+ξ), resulting from adding random perturbations, can be approximated by a second order Taylor expansion around
zi. Using the assumption about ξ in (4), the expected loss an be written as
Eξ∼K(ξ) [l(zi + ξ)] = l(zi) + 1
2 Eξ∼K(ξ)
[ ξ>H(i)ξ ] +O(ξ3) ≈ l(zi) + h2
2 tr(H(i)) (5)
where l(zi) is the loss without noise and H(i) = ∂ 2l ∂z2 (z) ∣∣ zi
the Hessian of l w.r.t z, evaluated at zi. Assuming that the noise ξ is small in its magnitude, O(ξ3) is negligible. This effect has been observed earlier by Webb (1994) and Bishop (1994). See Appendix A for derivations.
When concerned with maximum likelihood estimation of a conditional density f̂θ(y|x), the loss function coincides with the negative conditional log-likelihood l(yi, xi) = − log f̂θ(yi|xi). Let the standard deviation of the additive data noise ξx, ξy be hx and hy respectively. Maximum likelihood estimation (MLE) with data noise is equivalent to minimizing the loss
l(D) ≈ − n∑ i=1 log f̂θ(yi|xi)− h2y 2 n∑ i=1 dy∑ j=1 ∂2 log f̂θ(y|x) ∂y(j)∂y(j) ∣∣∣∣x=xi y=yi − h 2 x 2 n∑ i=1 dx∑ j=1 ∂2 log f̂θ(y|x) ∂x(j)∂x(j) ∣∣∣∣x=xi y=yi
(6) In that, the first term corresponds to the standard MLE objective, while the other two terms constitute a form of smoothness regularization. The second term in (6) penalizes large negative second derivatives of the conditional log density estimate log f̂θ(y|x) w.r.t. y. As the MLE objective pushes the density estimate towards high densities and strong concavity in the data points yi, the regularization term counteracts this tendency to over-fit and overall smoothes the fitted distribution. The third term penalizes large negative second derivatives w.r.t. the conditional variable x, thereby regularizing the sensitivity of the density estimate to changes in the conditional variable. The intensity of the noise regularization can be controlled through the variance (h2x and h 2 y) of the random perturbations.
Figure 1 illustrates the effect of the introduced noise regularization scheme on MDN estimates. Plain maximum likelihood estimation (left) leads to strong over-fitting, resulting in a spiky distribution that generalizes poorly beyond the training data. In contrast, training with noise regularization (center and right) results in smoother density estimates that are closer to the true conditional density.

heading: 4.2 CONSISTENCY OF NOISE REGULARIZATION
text: We now establish asymptotic consistency results for the proposed noise regularization. In particular, we show that, under some regularity conditions, concerning integrability and decay of the noise regularization, the solution of Algorithm 1 converges to the asymptotic MLE solution.
Let f̂θ(z) : Rdz × Θ → (0,∞) a continuous function of z and θ. Moreover, we assume that the parameter space Θ is compact. In the classical MLE setting, the idealized loss, corresponding to a (conditional) M-projection of the true data distribution onto the parametric family, reads as
l(θ) = −Ep(z) [ log f̂θ(z) ] (7)
As we typically just have a finite number of samples from p(z), the respective empirical estimate l̂n(θ) = − 1n ∑n i=1 log f̂θ(zi), zi
i.i.d∼ p(z) is used as training objective. Note that we now define the loss as function of θ, and, for fixed θ, treat ln(θ) as a random variable. Under some regularity conditions, one can invoke the uniform law of large numbers to show consistency of the empirical ML objective in the sense that supθ∈Θ |l̂n(θ)− l(θ)|
a.s.−−→ 0 (see Appendix B for details). In case of the presented noise regularization scheme, the maximum likelihood estimation is performed using on the augmented data {z̃j} rather than the original data {zi}. For our analysis, we view Algorithm 1 from a slightly different angle. In fact, the data augmentation procedure of uniformly selecting a data point from {z1, ..., zn} and perturbing it with a noise vector drawn from K can be viewed as drawing i.i.d. samples from a kernel density estimate q̂(h)n (z) = 1 n ∑n i=1 1 hdz K ( z−zi h ) . Hence, maximum likelihood estimation with variable noise can be understood as
1. forming a kernel density estimate q̂(h)n of the training data
2. followed by a (conditional) M-projection of q̂(h)n onto the parametric family.
In that, step 2 aims to find the θ∗ that minimizes the following objective: l(h)n (θ) = −Eq̂(h)n (z) [ log f̂θ(z) ] (8)
Since (8) is generally intractable, r samples are drawn from the kernel density estimate, forming the following Monte Carlo approximation of (8) which corresponds to the loss in line 5 Algorithm 1:
l̂(h)n,r(θ) = − 1
r r∑ j=1 log f̂θ(z̃j) , z̃j ∼ q̂(h)n (z) (9)
We are concerned with the consistency of the training procedure in Algorithm 1, similar to the classical MLE consistency result discussed above. Hence, we need to show that supθ∈Θ ∣∣∣l̂(h)n,r(θ)− l(θ)∣∣∣ a.s.−−→ 0 as n, r →∞. We begin our argument by decomposing the problem into easier sub-problems. In particular, the triangle inequality is used to obtain the following upper bound:
sup θ∈Θ ∣∣∣l̂(h)n,r(θ)− l(θ)∣∣∣ ≤ sup θ∈Θ ∣∣∣l̂(h)n,r(θ)− l(h)n (θ)∣∣∣+ sup θ∈Θ ∣∣∣l(h)n (θ)− l(θ)∣∣∣ (10) Note that l̂(h)n,r(θ) is based on samples from the kernel density estimate, which are obtained by adding random noise vectors ξ ∼ K(·) to our original training data. Since we can sample an unlimited amount of such random noise vectors, r can be chosen arbitrarily high. This allows us to make supθ∈Θ |l̂ (h) n,r(θ) − l(h)n (θ)| arbitrary small by the uniform law of large numbers. In order to make supθ∈Θ |l (h) n (θ)− l(θ)| small in the limit n→∞, the sequence of bandwidth parameters hn needs to be chosen appropriately. Such results can then be combined using a union bound argument. In the following we outline the steps leading us to the desired results. In that, the proof methodology is similar to Holmstrom & Koistinen (1992b). While they show consistency results for regression with a quadratic loss function, our proof deals with generic and inherently unbounded log-likelihood objectives and thus holds for a much more general class of learning problems. The full proofs can be found in the Appendix.
Initially, we have to make asymptotic integrability assumptions that ensure that the expectations in l(h)n (θ) and l(θ) are well-behaved in the limit (see Appendix C for details). Given respective integrability, we are able to obtain the following proposition.
Proposition 1 Suppose the regularity conditions (28) and (29) are satisfied, and that lim n→∞ hn = 0, lim n→∞ n(hn) d =∞ (11)
Then, lim n→∞ sup θ∈Θ ∣∣∣l(h)n (θ)− l(θ)∣∣∣ = 0 (12) almost surely.
In (11) we find conditions on the asymptotic behavior of the smoothing sequence (hn). These conditions also give us valuable guidance on how to properly choose the noise intensity in line 4 of Algorithm 1 (see Section 4.3 for discussion). The result in (12) demonstrates that, under the discussed conditions, replacing the empirical data distribution with a kernel density estimate still results in an asymptotically consistent maximum likelihood objective. However, as previously discussed, l (h) n (θ) is intractable and, thus, replaced by its sample estimate l̂ (h) n,r. Since we can draw an arbitrary amount of samples from q̂(h)n , we can approximate l (h) n (θ) with arbitrary precision. Given a fixed
data setD of size n > n0, this means that limr→∞ supθ∈Θ ∣∣∣l̂(h)n,r(θ)− l(h)n (θ)∣∣∣ = 0 almost surely, by
(29) and the uniform law of large numbers. Since our original goal was to also show consistency for n→∞, this result is combined with Proposition 1, obtaining the following consistency theorem.
Theorem 1 Suppose the regularity conditions (28) and (29) are satisfied, hn fulfills (11) and Θ is compact. Then,
lim n→∞ lim r→∞ sup θ∈Θ ∣∣∣l̂(h)n,r(θ)− l(θ)∣∣∣ = 0 (13) almost surely.
In that, lim used to denote the limit superior (“lim sup”) of a sequence.
Training a (conditional) density model with noise regularization means minimizing l̂(h)n,r(θ) w.r.t. θ. As result of this optimization, one obtains a parameter vector θ̂(h)n,r, which we hope is close to the minimizing parameter θ̄ of the ideal objective function l(θ). In the following, we establish asymptotic consistency results, similar to Theorem 1, in the parameter space. Therefore we first have to formalize the concept of closeness and optimality in the parameter space. Since a minimizing parameter θ̄ of l(θ) may not be unique, we define Θ∗ = {θ∗ | l(θ∗) ≤ l(θ) ∀θ ∈ Θ} as the set of global minimizers of l(θ), and d(θ,Θ∗) = minθ∗∈Θ∗{||θ− θ∗||2} as the distance of an arbitrary parameter θ to Θ∗. Based on these definitions, it can be shown that Algorithm 1 is asymptotically consistent in a sense that the minimizer of θ̂(h)n,r converges almost surely to the set of optimal parameters Θ∗.
Theorem 2 Suppose the regularity conditions (28) and (29) are satisfied, hn fulfills (11) and Θ is compact. For r > 0 and n > n0, let θ̂ (h) n,r ∈ Θ be a global minimizer of the empirical objective l̂(h)n,r. Then lim n→∞ lim r→∞ d(θ̂(h)n,r,Θ ∗) = 0 (14)
almost surely.
Note that Theorem 2 considers global optimizers, but equivalently holds for compact neighborhoods of a local minimum θ∗ (see discussion in Appendix C).

heading: 4.3 CHOOSING THE NOISE INTENSITY
text: After discussing the properties of noise regularization, we are interested in how to properly choose the noise intensity h, for different training data sets. Ideally, we would like to choose h so that |l(h)n (θ) − l(θ)| is minimized, which is practically not feasible since l(θ) is intractable. Inequality (30) gives as an upper bound on this quantity, suggesting to minimize l1 distance between the kernel density estimate q(h)n and the data distribution p(z). This is in turn a well-studied problem in the kernel density estimation literature (see e.g., Devroye & Luc (1987)). Unfortunately, general solutions of this problem require knowing p(z) which is not the case in practice. Under the assumption that p(z) and the kernel function K are Gaussian, the optimal bandwidth can be derived as h = 1.06σ̂n− 1 4+d (Silverman, 1986). In that, σ̂ denotes the estimated standard deviation of the data, n the number of data points and d the dimensionality of Z . This formula is widely known as the rule of thumb and often used as a heuristic for choosing h.
In addition, the conditions in (11) give us further intuition. The first condition tells us that hn needs to decay towards zero as n becomes large. This reflects the general theme in machine learning that the more data is available, the less inductive bias / regularization should be imposed. The second condition suggests that the bandwidth decay must happen at a rate slower than n− 1 d . For instance, the rule of thumb fulfills these two criteria and thus constitutes a useful guideline for selecting h. However, for highly non-Gaussian data distributions, the respective hn may decay too slowly and a faster decay rate such as n− 1 1+d may be appropriate.

heading: 5 EXPERIMENTS
text: This section provides a detailed experimental analysis of the proposed method, aiming to empirically validate the theoretical arguments outlined previously and investigating the practical efficacy of our regularization approach. In all experiments we use Gaussian pertubations of the data, i.e., K(ξ) = N (0, I). Since one of the key features of our noise regularization scheme is that it is agnostic to the choice of model, we evaluate its performance on three different neural network based CDE models: Mixture Density Networks (MDN) (Bishop, 1994), Kernel Mixture Networks (KMN) (Ambrogioni et al., 2017) and Normalizing Flows Networks (NFN) (Rezende & Mohamed, 2015; Trippe & Turner, 2018). In our experiments, we consider both simulated as well as real-world data sets. In particular, we simulate data from a 4-dimensional Gaussian Mixture (dx = 2, dy = 2) and a Skew-Normal distribution whose parameters are functionally dependent on x (dx = 1, dy = 1). In terms of real-world data, we use the following three data sources:
Euro Stoxx: Daily stock-market returns of the Euro Stoxx 50 index conditioned on various stock return factors relevant in finance (dx = 14, dy = 1).
NYC Taxi: Drop-off locations of Manhattan taxi trips conditioned on the pickup location, weekday and time (dx = 6, dy = 2).
UCI datasets: Standard data sets from the UCI machine learning repository (Dua & Graff, 2017). In particular, Boston Housing (dx = 13, dy = 1), Concrete (dx = 8, dy = 1), Energy (dx = 9, dy = 1).
The reported scores are test log-likelihoods, averaged over at least 5 random seeds alongside the respective standard deviation. For further details regarding the data sets and simulated data, we refer to Appendix E. The experiment data and code is available at TODO

heading: 5.1 NOISE INTENSITY SCHEDULES
text: We complement the discussion in 4.3 with an empirical investigation of different schedules of hn. In particular, we compare a) the rule of thumb hn ∝ n− 1 4+d b) a square root decay schedule hn ∝ n− 1
1+d c) a constant bandwidth hn = const. ∈ (0,∞) and d) no noise regularization, i.e. hn = 0. Figure 2 plots the respective test log-likelihoods against an increasing training set size n for the two simulated densities Gaussian Mixture and Skew Normal.
First, we observe that bandwidth rates which conform with the decay conditions seem to converge in performance to the non-regularized maximum likelihood estimator (red) as n becomes large. This reflects the theoretical result of Theorem 1. Second, a fixed bandwidth across n (green), violating (11), imposes asymptotic bias and thus saturates in performance vastly before its counterparts. Third, as hypothesized, the relatively slow decay of hn through the rule of thumb works better for data distributions that have larger similarities to a Gaussian, i.e., in our case the Skew Normal distribution. In contrast, the highly non-Gaussian data from the Gaussian Mixture requires faster decay rates like the square root decay schedule. Most importantly, noise regularization substantially improves the estimator’s performance when only little training data is available.

heading: 5.2 REGULARIZATION COMPARISON
text: We now investigate how the proposed noise regularization scheme compares to classical regularization techniques. In particular, we consider an l1 and l2-penalty on the neural network weights as regularization term, the weight decay technique of Loshchilov & Hutter (2019)1, as well a Bayesian neural network (Neal, 2012) trained with variational inference using a Gaussian prior and posterior (Blei et al., 2017).
First, we study the performance of the regularization techniques on our two simulation benchmarks. Figure 3 depicts the respective test log-likelihood across different training set sizes. For each regularization method, the regularization hyper-parameter has been optimized via grid search.
As one would expect, the importance of regularization, i.e., performance difference to un-regularized model, decreases as the amount of training data becomes larger. The noise regularization scheme
1Note that an l2 regularizer and weight decay are not equivalent since we use the adaptive learning rate technique Adam. See Loshchilov & Hutter (2019) for details.
yields similar performance across the different CDE models while the other regularizers vary greatly in their performance depending on the different models. This reflects the fact that noise regularization is agnostic to the parameterization of the CDE model while regularizers in the parameter space are dependent on the internal structure of the model. Most importantly, noise regularization performs well across all models and sample sizes. In the great majority of configurations it outperforms the other methods. Especially when little training data is available, noise regularization ensures a moderate test error while the other approaches mostly fail to do so.
Next, we consider real world data sets. Since now the amount of data we can use for hyper-parameter selection, training and evaluation is limited, we use 5-fold cross-validation to select the parameters for each regularization method. The test log-likelihoods, reported in Table 1, are averages over 3 different train/test splits and 5 seeds each for initializing the neural networks. The held out test set amounts to 20% of the overall data sets. Consistent with the results of the simulation study, noise regularization outperforms the other methods across the great majority of data sets and CDE models.

heading: 5.3 CONDITIONAL DENSITY ESTIMATOR BENCHMARK STUDY
text: We benchmark neural network based density estimators against state-of-the art CDE approaches. While neural networks are the obvious choice when a large amount of training data is available, we pose the questions how such estimators compete against well-established non-parametric methods in small data regimes. In particular, we compare to the three following CDE methods:
Conditional Kernel Density Estimation (CKDE). Non-parametric method that forms a KDE of both p(x, y) and p(x) to compute its estimate as p̂(y|x) := p̂(x, y)/p̂(x) (Li & Racine, 2007).
-Neighborhood kernel density estimation (NKDE). Non-parametric method that considers only a local subset of training points to form a density estimate.
Least-Squares Conditional Density Estimation (LSCDE). Semi-parametric estimator that computes the conditional density as linear combination of fixed kernels (Sugiyama & Takeuchi, 2010).
For the kernel density estimation based methods CKDE and NKDE, we perform bandwidth selection via the rule of thumb (R.O.T) (Silverman, 1982; Sheather & Jones, 1991) and via maximum likelihood leave-one-out cross-validation (CV-ML) (Rudemo, 1982; Hall et al., 1992). In case of LSCDE, MDN, KMN and NFN, the respective hyper-parameters are selected via 5-fold cross-validation grid search on the training set. Note that, in contrast to Section 5.2 which focuses on regularization parameters, the grid search here extends to more hyper-parameters. The respective test log-likelihood scores are listed in Table 2. For the majority of data sets, all three neural network based methods outperform all of the non- and semi-parametric methods. Perhaps surprisingly, it can be seen that, when properly regularized, neural network based CDE works well even when training data is scarce, such as in case of the Boston Housing data set.

heading: 6 CONCLUSION
text: This paper addresses conditional density estimation with high-capacity models. In particular, we propose to add small random perturbations to the data during training. We demonstrate that the resulting noise regularization method corresponds to a smoothness regularization and prove its asymptotic consistency. The experimental results underline the effectiveness of the proposed method, demonstrating that it consistently outperforms other regularization methods across various conditional density models and data sets. This makes neural network based CDE the preferable method, even when only little training data is available. While we assess the estimator performance in terms of the test log-likelihood, an interesting question for future research is whether the noise regularization also improves the respective uncertainty estimates for downstream tasks such as safe control and decision making.

heading: A DERIVATION SMOOTHNESS REGULARIZATION
text: Let l(D) be a loss function over a set of data points D = {z1, ..., zN}, which can be partitioned into a sum of losses corresponding to each data point xn:
lD(D) = n∑ i=1 l(zi) (15)
Also, let each zi be perturbed by a random noise vector ξ ∼ K(ξ) with zero mean and i.i.d. elements, i.e. Eξ∼K(ξ) [ξ] = 0 and Eξ∼K(ξ) [ ξnξ > j ] = h2I (16) The resulting loss l(zi + ξ) can be approximated by a second order Taylor expansion around zi
l(zi + ξ) = l(zi) + ξ >∇zl(z) ∣∣ zi + 1 2 ξ>∇2zl(z) ∣∣ zi ξ +O(ξ3) (17)
Assuming that the noise ξ is small in its magnitude, O(ξ3) may be neglected. The expected loss under K(ξ) follows directly from (17):
Eξ∼K(ξ) [l(zi + ξ)] = l(zi) + Eξ∼K(ξ) [ ξ>∇xl(z) ∣∣ zi ] + 1 2 Eξ∼K(ξ) [ ξ>∇2xl(z) ∣∣ zi ξ ]
(18)
Using the assumption about ξ in (16) we can simplify (18) as follows:
Eξ∼K(ξ) [l(zi + ξ)] = l(zi) + Eξ∼K(ξ) [ξ] >∇zl(z) ∣∣ zi + 1 2 Eξ∼K(ξ) [ ξ>∇2zl(z) ∣∣ zi ξ ]
(19)
= l(zi) + 1
2 Eξ∼K(ξ)
[ ξ>H(i)ξ ] (20)
= l(zi) + 1
2 Eξ∼K(ξ) ∑ j ∑ k ξjξk ∂2l(z) ∂z(j)∂z(k) ∣∣∣∣ zi  (21) = l(zi) + 1
2 ∑ j Eξ [ ξ2j ] ∂2l(z) ∂z(j)∂z(j) ∣∣∣∣ zi + 1 2 ∑ j ∑ k 6=j Eξ [ξjξk] ∂2l(z) ∂z(j)∂z(k) ∣∣∣∣ zi
(22)
= l(zi) + η2
2 ∑ j ∂2l(z) ∂z(j)∂z(j) ∣∣∣∣ zi
(23)
= l(zi) + η2
2 tr(H(i)) (24) In that, l(zi) is the loss without noise and H(i) = ∇2zl(z) ∣∣ zi
the Hessian of l at zi. With z(j) we denote the elements of the column vector z.
B VANILLA CONDITIONAL MLE OBJECTIVE IS UNIFORMLY CONSISTENT
The objective function corresponding to a conditional M-projection. l(θ) = −Ep(x,y) [ log f̂θ(y|x) ] (25)
The sample equivalent:
l̂n(θ) = − 1
n n∑ i=1 log f̂θ(yi|xi) , (xi, yi) i.i.d∼ P (X,Y ) (26)
Corollary 1 Let Θ be a compact set and and f̂θ : Rl × Rm × Θ → (0,∞) continuous in θ for all (x, y) ∈ Rl × Rm such that Ep(x,y) [ supθ∈Θ log f̂θ(y|x) ] <∞. Then, as n→∞, we have
sup θ∈Θ ∣∣∣l̂n(θ)− l(θ)∣∣∣ a.s.−−→ 0 (27) Proof. The corollary follows directly from the uniform law of large numbers.

heading: C CONSISTENCY PROOFS
text: Lemma 1 Suppose for some > 0 there exists a constant B( )p such that∫ | log f̂θ(z)|1+ p(z)dz ≤ B( )p <∞ ∀θ ∈ Θ (28)
and there exists an n0 such that for all n > n0 there exists a constant B ( ) q̂ such that∫
| log f̂θ(z)|1+ q̂(hn)n (z)dz ≤ B ( ) q̂ <∞ ∀θ ∈ Θ (29)
almost surely. Then, the inequality
sup θ∈Θ ∣∣∣l(h)n (θ)− l(θ)∣∣∣ ≤ C (∫ |q̂(h)n (z)− p(z)|dz) 1+ (30) where C is a constant holds with probability 1 for all n > n0.
Proof of Lemma 1 Using Hoelder’s inequality and the nonnegativity of p and q̂(h)n , we obtain∣∣∣l(h)n (θ)− l(θ)∣∣∣ = ∣∣∣∣∫ log f̂θ(z)(q̂(h)n (z)− p(z))dz∣∣∣∣ ≤ ∫ | log f̂θ(z)| |q̂(h)n (z)− p(z)|dz
= ∫ | log f̂θ(z)| |q̂(h)n (z)− p(z)| 1 1+ |q̂(h)n (z)− p(z)| 1+ dz
≤ (∫ | log f̂θ(z)|1+ |q̂(h)n (z)− p(z)|dz ) 1 1+
(∫ |q̂(h)n (z)− p(z)|dz ) 1+
≤ (∫
| log f̂θ(z)|1+ q̂(h)n (z)dz
+ ∫ | log f̂θ(z)|1+ p(z)dz ) 1 1+ (∫ |q̂(h)n (z)− p(z)|dz ) 1+
Employing the regularity conditions (28) and (29) and writing C( ) = B( )p + B ( ) q̂ , it follows that ∃n0 such that ∀n > n0
sup θ∈Θ ∣∣∣l(h)n (θ)− l(θ)∣∣∣ ≤ (B( )p +B( )q̂ )(∫ |q̂(h)n (z)− p(z)|dz) 1+ = C( ) (∫ |q̂(h)n (z)− p(z)|dz ) 1+
with probability 1.
Lemma 1 states regularity conditions ensuring that the expectations in l(h)n (θ) and l(θ) are wellbehaved in the limit. In particular, (28) and (29) imply uniform and absolute integrability of the loglikelihoods under the respective probability measures induced by p and q̂(h)n . Since we are interested in the asymptotic behavior, it is sufficient for (29) to hold for n large enough with probability 1.
Inequality (30) shows that we can make |l(h)n (θ) − l(θ)| small by reducing the l1-distance between the true density p and the kernel density estimate q̂(h)n . There exists already a vast body of literature, discussing how to properly choose the kernel K and the bandwidth sequence (hn) so that∫ |q̂(hn)n (z) − p(z)|dz → 0. We employ the results in Devroye (1983) for our purposes, leading us to Proposition 1.
Proof of Proposition 1. Let A denote the event that ∃n0∀n > n0 inequality (30) holds for some constant C( ). From our regularity assumptions it follows that P(Ac) = 0. Given that A holds, we just have to show that ∫ |q̂(h)n (z)− p(z)|dz
a.s.−−→ 0. Then, the upper bound in (30) tends to zero and we can conclude our proposition.
For any δ > 0 let Bn denote the event∫ |q̂(h)n (z)− p(z)|dz ≤ δ (31)
wherein q̂(h)n (z) is a kernel density estimate obtained based on n samples from p(z). Under the conditions in (11) we can apply Theorem 1 of Devroye (1983), obtaining an upper bound on the probability that (31) does not hold, i.e. ∃u,m0 such that P(Bcn) ≤ e−un for all n > m0. Since we need both A and Bn for n → ∞ to hold, we consider the intersection of the events (A ∩Bn). Using a union bound argument it follows that ∃k0 such that ∀n > k0 : P((A ∩Bn)c) ≤ P(Ac) + P(Bcn) = 0 + e−un = e−un. Note that we can simply choose k0 = max{n0,m0} for this to hold. Hence, ∑∞ n=k0+1 P ((A ∩ Bn)c) < ∑∞ n=1 e
−un = 1eu−1 < ∞ and by the Borel-Cantelli lemma we can conclude that
lim n→∞ sup θ∈Θ ∣∣∣l(h)n (θ)− l(θ)∣∣∣ = 0 (32) holds with probability 1.
Proof of Theorem 1. The inequality in (10) implies that for any n > n0,
lim r→∞ sup θ∈Θ ∣∣∣l̂(h)n,r(θ)− l(θ)∣∣∣ ≤ lim r→∞ sup θ∈Θ ∣∣∣l̂(h)n,r(θ)− l(h)n (θ)∣∣∣+ sup θ∈Θ ∣∣∣l(h)n (θ)− l(θ)∣∣∣ (33) Let n > n0 be fixed but arbitrary and denote
Jn,r = sup θ∈Θ ∣∣∣l̂(h)n,r(θ)− l(h)n (θ)∣∣∣ r ∈ N, n > n0 (34) It is important to note that Jn,r is a random variable that depends on the samples Z(n) = (Z1, ..., Zn) as well as on the randomness inherent in Algorithm 1. We define I(r) = (I1, ...Ir) as the indices sampled uniformly from {1, ..., n} and Ξ(r) = (ξ1, ...xir) as the sequence of perturbation vectors sampled from K. Let P (Z(n)), P (I(r)) and P (Ξ(r)) be probability measures of the respective random sequences.
If we fix Z(n) to be equal to an arbitrary sequence Z(n), then q̂(h)n is fixed and we can treat Jn,r as the regular difference between a sample estimate and expectation under q̂(h)n . By the regularity condition (29), the compactness of Θ and the continuity of fθ in θ, we can invoke the uniform law of large numbers to show that
lim r→∞ Jn,r = lim r→∞ sup θ∈Θ ∣∣∣l̂(h)n,r(θ)− l(h)n (θ)∣∣∣ = 0 (35) with probability 1.
Now we want to show that (35) also holds with probability 1 for random training samples Z(n). First, we write Jn,r as a deterministic function of random variables:
Jn,r = J(Z (n), I(r),Ξ(r)) (36)
This allows us to restate the result in (35) as follows: PI(r),Ξ(r) ( ∀δ > 0 ∃r0 ∀r > r0 : J(Z(n) = Z(n), I(r),Ξ(r)) < δ ) = ∫ ∫ 1 ( ∀δ > 0 ∃r0 ∀r > r0 : J(Z(n) = Z(n), I(r),Ξ(r)) < δ ) dP (Ξ(r))dP (I(r))
= 1
(37)
In that 1(A) denotes an indicator function which returns 1 if A is true and 0 else. Next we consider the probability that the convergence in (35) holds for random Z(n):
PZ(n),I(r),Ξ(r) ( ∀δ > 0 ∃r0 ∀r > r0 : J(Z(n), I(r),Ξ(r)) < δ ) = ∫ ∫ ∫ 1 ( ∀δ > 0 ∃r0 ∀r > r0 : J(Z(n), I(r),Ξ(r)) < δ ) dP (Ξ(r))dP (I(r))dP (Z(n))
= ∫ dP (Z(n)) (∫ ∫ 1 ( ∀δ > 0 ∃r0 ∀r > r0 : J(Z(n), I(r),Ξ(r)) < δ ) dP (Ξ(r))dP (I(r)) ) ︸ ︷︷ ︸
=1
= 1
Note that we can dP (Z(n)) move outside of the inner integrals, since Z(n) is independent from I(r) and Ξ(r). Hence, we can conclude that (35) also holds, which we denote as eventA, with probability 1 for random training data.
From Proposition 1 we know, that
lim n→∞ sup θ∈Θ ∣∣∣l(h)n (θ)− l(θ)∣∣∣ = 0 (38) with probability 1. We denote the event that (38) holds as B. Since P (Ac) = P (Bc) = 0, we can use a union bound argument to show that P (A ∩B) = 1. From (35) and (33) it follows that for any n > n0,
lim r→∞ sup θ∈Θ ∣∣∣l̂(h)n,r(θ)− l(θ)∣∣∣ ≤ sup θ∈Θ ∣∣∣l(h)n (θ)− l(θ)∣∣∣ (39) with probability 1. Finally, we combine this result with (38), obtaining that
lim n→∞ lim r→∞ sup θ∈Θ ∣∣∣l̂(h)n,r(θ)− l(θ)∣∣∣ = 0 (40) almost surely, which concludes the proof.
Proof of Theorem 2. The proof follows the argument used in Theorem 1 of White (1989). In the following, we assume that (13) holds. From Theorem 1 we know that this is the case with probability 1. Respectively, we only consider realizations of our training data Z(n) and noise samples I(r), Ξ(r) for which the convergence in (13) holds (see proof of Theorem 1 for details on this notation).
For such realization, let (θ̂(h)n,r) be minimizers of l̂ (h) n,r. Also let (ni)i and for any i, (ri,j)j be increasing sequences of positive integers. Define vi,j := θ̂ (h) ni,ri,j and µi,j(θ) := l̂ (h) ni,ri,j (θ). Due to the compactness of Θ and the Bolzano-Weierstrass property thereof, there exists a limit point θ0 ∈ Θ and increasing subsequences (ik)k, (jk)k so that vik,jk → θ0 as k →∞. From the triangle inequality, it follows that for any > 0 there exists k0 so that ∀k > k0
|µik,jk(vik,jk)− l(θ0)| ≤ |µik,jk(vik,jk)− l(vik,jk)|+ |l(vik,jk)− l(θ0)| < 2 (41)
given the convergence established in Theorem 1 and the continuity of l in θ. Next, the result above is extended to
l(θ0)− l(θ) = [l(θ0)−µik,jk(vik,jk)] + [µik,jk(vik,jk)−µik,jk(θ)] + [µik,jk(θ)− l(θ)] ≤ 3 (42)
which again holds for k large enough. This due to (41), µik,jk(vik,jk)−µik,jk(θ) ≤ 0 since vik,jk is the minimizer of µik,jk , and µik,jk(θ)− l(θ) < by Theorem 1. Because can be made arbitrarily small, l(θ0) ≤ l(θ) as k → ∞. Because θ ∈ Θ is arbitrary, θ0 must be in Θ∗. In turn, since (ni)i, (ri,j)j and (ik)k, (jk)k were chosen arbitrarily, every limit point of a sequence (vik,jk)k must be in Θ∗.
In the final step, we proof the theorem by contradiction. Suppose that (14) does not hold. In this case, there must exist an > 0 and sequences (ni)i, (ri,j)j and (ik)k, (jk)k such that ||(vik,jk)k−θ̄||2 > for all k and θ̄ ∈ Θ∗. However, by the previous argument the limit point of the any sequence (vik,jk)k must be in Θ ∗. That is a contradiction to ||(vik,jk)k − θ̄||2 > ∀ k, θ̄ ∈ Θ∗. Since
the random sequences Z(n), I(r), Ξ(r) where chosen from a set with probability mass of 1, we can conclude our proposition that
lim n→∞ lim r→∞
d(θ̂(h)n,r,Θ ∗) = 0
almost surely.
Discussion of Theorem 2. Note that, similar to θ∗, θ̂(h)n,r does not have to be unique. In case there are multiple minimizers of l̂(h)n,r, we can chose one of them arbitrarily and the proof of the theorem still holds. Theorem 2 considers global optimizers over a set of parameters Θ, which may not be attainable in practical settings. However, the application of the theorem to the context of local optimization is straightforward when Θ is chosen as a compact neighborhood of a local minimum θ∗ of l (Holmstrom & Koistinen, 1992b). If we set Θ∗ = {θ∗} and restrict minimization over l̂(h)n,r to the local region, then θ̂(h)n,r converges to Θ∗ as n, r →∞ in the sense of Theorem 2.

heading: D CONDITIONAL DENSITY ESTIMATION MODELS

heading: D.1 MIXTURE DENSITY NETWORK
text: Mixture Density Networks (MDNs) combine conventional neural networks with a mixture density model for the purpose of estimating conditional distributions p(y|x) (Bishop, 1994). In particular, the parameters of the unconditional mixture distribution p(y) are outputted by the neural network, which takes the conditional variable x as input.
For our purpose, we employ a Gaussian Mixture Model (GMM) with diagonal covariance matrices as density model. The conditional density estimate p̂(y|x) follows as weighted sum of K Gaussians
p̂(y|x) = K∑ k=1 wk(x; θ)N ( y|µk(x; θ), σ2k(x; θ) ) (43)
wherein wk(x; θ) denote the weight, µk(x; θ) the mean and σ2k(x; θ) the variance of the k-th Gaussian component. All the GMM parameters are governed by the neural network with parameters θ and input x. The mixing weights wk(x; θ) must resemble a categorical distribution, i.e. it must hold that∑K k=1 wk(x; θ) = 1 and wk(x; θ) ≥ 0 ∀k. To satisfy the conditions, the softmax linearity is used for the output neurons corresponding to wk(x; θ). Similarly, the standard deviations σk(x) must be positive, which is ensured by a sofplus non-linearity. Since the component means µk(x; θ) are not subject to such restrictions, we use a linear output layer without non-linearity for the respective output neurons.
For the experiments in 5.2 and 5.1, we set K = 10 and use a neural network with two hidden layers of size 32.

heading: D.1.1 KERNEL MIXTURE NETWORK
text: While MDNs resemble a purely parametric conditional density model, a closely related approach, the Kernel Mixture Network (KMN), combines both non-parametric and parametric elements (Ambrogioni et al., 2017). Similar to MDNs, a mixture density model of p̂(y) is combined with a neural network which takes the conditional variable x as an input. However, the neural network only controls the weights of the mixture components while the component centers and scales are fixed w.r.t. to x. For each of the kernel centers, M different scale/bandwidth parameters σm are chosen. As for MDNs, we employ Gaussians as mixture components, wherein the scale parameter directly coincides with the standard deviation.
Let K be the number of kernel centers µk and M the number of different kernel scales σm. The KMN conditional density estimate reads as follows:
p̂(y|x) = K∑ k=1 M∑ m=1 wk,m(x; θ)N (y|µk, σ2m) (44)
As previously, the weights wk,m correspond to a softmax function. The M scale parameters σm are learned jointly with the neural network parameters θ. The centers µk are initially chosen by k-means clustering on the {yi}ni=1 in the training data set. Overall, the KMN model is more restrictive than MDN as the locations and scales of the mixture components are fixed during inference and cannot be controlled by the neural network. However, due to the reduced flexibility of KMNs, they are less prone to over-fit than MDNs.
For the experiments in 5.2 and 5.1, we set K = 50 and M = 2. The respective neural network has two hidden layers of size 32.

heading: D.2 NORMALIZING FLOW NETWORK
text: The Normalizing Flow Network (NFN) is similar to the MDN and KMN in that a neural network takes the conditional variable x as its input and outputs parameters for the distribution over y. For the NFN, the distribution is given by a Normalizing Flow (Rezende & Mohamed, 2015). It works by transforming a simple base distribution and an accordingly distributed random variable Z0 through a series of invertible, parametrized mappings f = fN ◦ · · · ◦ f1 into a successively more complex distribution p(f(Z0)). The PDF of samples zN ∼ p(f(Z0)) can be evaluted using the change-ofvariable formula:
log p(zN ) = log p(z0)− N∑ n=1 log ∣∣∣det ∂fn ∂zn−1 ∣∣∣ (45) The Normalizing Flows from Rezende & Mohamed (2015) were introduced in the context of posterior estimation in variational inference. They are optimized for fast sampling while the likelihood evaluation for externally provided data is comparatively slow. To make them useful for CDE, we invert the direction of the flows, defining a mapping from the transformed distribution p(ZN ) to the base distribution p(Z0) by setting f̂−1i (zi) = fi(zi).
We experimented with three types of flows: planar flows, radial flows as parametrized by Trippe & Turner (2018) and affine flows f−1(z) = exp(a)z+b. We have found that one affine flow combined with multiple radial flows performs favourably in most settings.
For the experiments in 5.2 and 5.1, we used a standard Gaussian as the base distribution that is transformed through one affine flow and ten radial flows. The respective neural network has two hidden layers of size 32.

heading: E SIMULATED DENSITIES AND DATASETS

heading: E.1 SKEWNORMAL
text: The data generating process (x, y) ∼ p(x, y) resembles a bivariate joint-distribution, wherein x ∈ R follows a normal distribution and y ∈ R a conditional skew-normal distribution (Anděl et al., 1984). The parameters (ξ, ω, α) of the skew normal distribution are functionally dependent on x. Specifically, the functional dependencies are the following:
x ∼ N ( · ∣∣∣∣µ = 0, σ = 12 ) (46)
ξ(x) = a ∗ x+ b a, b ∈ R (47) ω(x) = c ∗ x2 + d c, d ∈ R (48)
α(x) = αlow + 1
1 + e−x ∗ (αhigh − αlow) (49) y ∼ SkewNormal ( ξ(x), ω(x), α(x) ) (50)
Accordingly, the conditional probability density p(y|x) corresponds to the skew normal density function:
p(y|x) = 2 ω(x) N ( y − ξ(x) ω(x) ) Φ ( α(x) y − ξ(x) ω(x) ) (51)
In that, N (·) denotes the density, and Φ(·) the cumulative distribution function of the standard normal distribution. The shape parameter α(x) controls the skewness and kurtosis of the distribution.
We set αlow = −4 and αhigh = 0, giving p(y|x) a negative skewness that decreases as x increases. This distribution will allow us to evaluate the performance of the density estimators in presence of skewness, a phenomenon that we often observe in financial market variables. Figure 4a illustrates the conditional skew normal distribution.

heading: E.2 GAUSSIAN MIXTURE
text: The joint distribution p(x, y) follows a Gaussian Mixture Model in R4 with 5 Gaussian components, i.e. K = 5. We assume that x ∈ R2 and y ∈ R2 can be factorized, i.e.
p(x, y) = K∑ i=1 wk N (y|µy,k,Σy,k)N (x|µx,k,Σx,k) (52)
When x and y can be factorized as in (52), the conditional density p(y|x) can be derived in closed form:
p(y|x) = K∑ i=1 Wk(x) N (y|µy,k,Σy,k) (53)
wherein the mixture weights are a function of x:
Wk(x) = wk N (x|µx,k,Σx,k)∑K j=1 wk N (x|µx,j ,Σx,j)
(54)
For details and derivations we refer the interested reader to Guang Sung (2004) and Gilardi et al. (2002). The weightswk are sampled from a uniform distributionU(0, 1) and then normalized to sum to one. The component means are sampled from a spherical Gaussian with zero mean and standard deviation of σ = 1.5. The covariance matrices Σy,k) and Σy,k) are sampled from a Gaussian with mean 1 and standard deviation 0.5, and then projected onto the cone of positive definite matrices.
Since we can hardly visualize a 4-dimensional GMM, Figure 4b depicts a 2-dimensional equivalent, generated with the procedure explained above.

heading: E.3 EURO STOXX 50 DATA
text: The Euro Stoxx 50 data comprises 3169 trading days, dated from January 2003 until June 2015. The goal is to predict the conditional probability density of 1-day log-returns, conditioned on 14 explanatory variables. These conditional variables comprise classical return factors from finance as well as option implied moments. For details, we refer to Rothfuss et al. (2019). Overall, the target variable is one-dimensional, i.e. y ∈ Y ⊆ R, whereas the conditional variable x constitutes a 14-dimensional vector, i.e. x ∈ X ⊆ R14.

heading: E.4 NYC TAXI DATA
text: We follow the setup in Dutordoir et al. (2018). The dataset contains records of taxi trips in the Manhattan area operated in January 2016. The objective is to predict spatial distributions of the dropoff location, based on the pick-up location, the day of the week, and the time of day. In that, the two temporal features are represented as sine and cosine with natural periods. Accordingly, the target variable y is 2-dimensional (longitude and latitude of dropoff-location) whereas the conditional variable is 6-dimensional. From the ca. 1 million trips, we randomly sample 10,000 trips to serve as training data.

heading: E.5 UCI
text: Boston Housing Concerns the value of houses in the suburban area of Boston. Conditional variables are mostly socio-economic as well as geographical factors. For more details see https://archive.ics.uci.edu/ml/machine-learning-databases/housing/
Concrete The task is to predict the compressive strength of concrete given variables describing the conrete composition. For more details see https://archive.ics.uci.edu/ml/machine-learningdatabases/concrete/compressive/
Energy Concerns the energy efficiency of homes. The task is to predict the cooling load based on features describing the build of the respective house. For more details see https://archive.ics.uci.edu/ml/datasets/energy+efficiency

Review summaries:
The paper considers the problem of parametric conditional density estimation, i.e. given a set of points {(x_n, y_n)} drawn from a distribution $p(x,y)$, the task is to estimate the conditional distribution p(x|y). The paper considers parametric estimation where in given a parametrized family of distributions f_{theta} we wish to minimize the likelihood of seeing the given data over theta. The parametric family in a lot of applications consists of highly expressive families like neural networks, which leads to the issue of overfitting in small data regimes. This has been tackled via regularization over the parameter space which might be hard to interpret as the associated inductive bias is not well understood and depends on the parametric family under consideration. On the other hand the paper proposes to add explicit noise in the examples used during training, i.e. irrespective of the optimization procedure (which could be mini-bath sgd) the paper proposes to draw examples from the data set, explicitly add noise onto the examples and create a proxy objective over the augmented data set. 

The paper presents a regularization technique for conditional density estimation. The method is simple: adding noise to the data points, and training on the noisy data points. The paper also further gives an interpretation of the method, as a form of smoothing the curvature of the density function. It further proves the consistency of the method.