Title: NOISE REGULARIZATION FOR CONDITIONAL DENSITY ESTIMATION

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 non-and semi-parametric approaches, even when training data is scarce.

1. INTRODUCTION:

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 {(x n , y n )} N n=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.

2. BACKGROUND:

Density Estimation. Let X be a random variable with probability density function (PDF) p(x) defined over the domain X ⊆ R dx . Given a collection D = {x 1 , ..., x n } 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 p D (x) = 1 n n i=1 δ(||x -x i ||) (i.e., mixture of point masses in the observations x i ) and the parametric distribution fθ :

EQUATION

From a geometric perspective, (1) can be viewed as an orthogonal projection of p D (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 x n (Rosenblatt, 1956; Parzen, 1962) . The resulting density estimate reads as q(x) =

1 nh d n i=1 K x-xi h . One popular choice of K(•) is a Gaussian K(z) = (2π) -d 2 exp (-1 2 z 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 ⊆ R dx and Y ⊆ R dy 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 = {(x n , y n )} N n=1 drawn from the joint distribution (x n , y n ) ∼ 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):

EQUATION

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θ (y i |x i ) (3)

3. RELATED WORK:

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) , l 1 /l 2 -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 & Schö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.

4. NOISE REGULARIZATION:

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; Kukač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 (x i , y i ) by random variables xi = x i + ξ x and ỹi = y i + ξ y where the perturbation vectors are sampled from noise distribu- tions K x (ξ x ) and K y (ξ 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(ξ) ξξ = h 2 I (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 = {z 1 , ..., z n }, noise intensity h Require: number of perturbed samples r, 1: for j = 1 to r do 2:

Select i ∈ {1, ..., n} with equal prob. Set zj = z j + 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).

4.1. VARIABLE NOISE AS SMOOTHNESS REGULARIZATION:

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 = {z 1 , ..., z n }, which can be partitioned into a sum of losses l(D) = n i=1 l(z i ), corresponding to each data point z i : The expected loss l(z i +ξ), resulting from adding random perturbations, can be approximated by a second order Taylor expansion around z i . Using the assumption about ξ in (4), the expected loss an be written as

E ξ∼K(ξ) [l(z i + ξ)] = l(z i ) + 1 2 E ξ∼K(ξ) ξ H (i) ξ + O(ξ 3 ) ≈ l(z i ) + h 2 2 tr(H (i) ) (5)

where l(z i ) is the loss without noise and H (i) = ∂ 2 l ∂z 2 (z) zi the Hessian of l w.r.t z, evaluated at z i . 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(y i , x i ) = -log fθ (y i |x i ). Let the standard deviation of the additive data noise ξ x , ξ y be h x and h y respectively. Maximum likelihood estimation (MLE) with data noise is equivalent to minimizing the loss

l(D) ≈ - n i=1 log fθ (y i |x i ) - h 2 y 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 y i , 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 (h 2

x 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.

4.2. CONSISTENCY OF NOISE REGULARIZATION:

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) : R dz × Θ → (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

EQUATION

As we typically just have a finite number of samples from p(z), the respective empirical estimate

ln (θ) = -1 n n i=1 log fθ (z i ), z i i.i.d

∼ p(z) is used as training objective. Note that we now define the loss as function of θ, and, for fixed θ, treat l n (θ) 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 θ∈Θ | ln (θ) -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 {z i }. 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 {z 1 , ..., z n } 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 h dz 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 (θ) = -E q(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:

EQUATION

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:

EQUATION

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 h n 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→∞ h n = 0, lim n→∞ n(h n ) 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 (h n ). 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

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 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, h n fulfills (11) and Θ is compact. For r > 0 and n > n 0 , let θ(h) n,r ∈ Θ be a global minimizer of the empirical objective l(h) n,r . Then lim

EQUATION

)

almost surely.

Note that Theorem 2 considers global optimizers, but equivalently holds for compact neighborhoods of a local minimum θ * (see discussion in Appendix C).

4.3. CHOOSING THE NOISE INTENSITY:

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 l 1 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 h n 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 h n may decay too slowly and a faster decay rate such as n -1 1+d may be appropriate.

5. EXPERIMENTS:

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) UCI datasets: Standard data sets from the UCI machine learning repository (Dua & Graff, 2017) . In particular, Boston Housing (d x = 13, d y = 1), Concrete (d x = 8, d y = 1), Energy (d x = 9, d y = 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

5.1. NOISE INTENSITY SCHEDULES:

We complement the discussion in 4.3 with an empirical investigation of different schedules of h n . In particular, we compare a) the rule of thumb h n ∝ n -1 4+d b) a square root decay schedule h n ∝ n -1 1+d c) a constant bandwidth h n = const. ∈ (0, ∞) and d) no noise regularization, i.e. h n = 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 h n 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.

5.2. REGULARIZATION COMPARISON:

We now investigate how the proposed noise regularization scheme compares to classical regularization techniques. In particular, we consider an l 1 and l 2 -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 1 : Comparison of various regularization methods for three neural network based CDE models across 5 data sets. We report the test log-likelihood and its respective standard deviation (higher log-likelihood values are better). 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.

5.3. CONDITIONAL DENSITY ESTIMATOR BENCHMARK STUDY:

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.

6. CONCLUSION:

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.

A DERIVATION SMOOTHNESS REGULARIZATION:

Let l(D) be a loss function over a set of data points D = {z 1 , ..., z N }, which can be partitioned into a sum of losses corresponding to each data point x n :

EQUATION

Also, let each z i 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 = h 2 I (16) The resulting loss l(z i + ξ) can be approximated by a second order Taylor expansion around

EQUATION

Assuming that the noise ξ is small in its magnitude, O(ξ 3 ) may be neglected. The expected loss under K(ξ) follows directly from ( 17):

EQUATION

Using the assumption about ξ in ( 16) we can simplify (18) as follows:

EQUATION

In that, l(z i ) is the loss without noise and H (i) = ∇ 2 z l(z) zi the Hessian of l at z i . 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.

EQUATION

The sample equivalent:

EQUATION

Corollary 1 Let Θ be a compact set and and fθ :

R l × R m × Θ → (0, ∞) continuous in θ for all (x, y) ∈ R l × R m such that E p(x,y) sup θ∈Θ log fθ (y|x) < ∞. Then, as n → ∞, we have sup θ∈Θ ln (θ) -l(θ) a.s.

--→ 0

Proof. The corollary follows directly from the uniform law of large numbers.

C CONSISTENCY PROOFS:

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 n 0 such that for all n > n 0 there exists a constant B ( )

EQUATION

almost surely. Then, the inequality

EQUATION

)

where C is a constant holds with probability 1 for all n > n 0 .

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 ∃n 0 such that ∀n > n 0 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 l 1 -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 (h n ) so that |q 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) : P Z (n) ,I (r) ,Ξ (r) ∀δ > 0 ∃r 0 ∀r > r 0 : J(Z (n) , I (r) , Ξ (r) ) < δ = 1 ∀δ > 0 ∃r 0 ∀r > r 0 : J(Z (n) , I (r) , Ξ (r) ) < δ dP (Ξ (r) )dP (I (r) )dP (Z (n) )

= dP (Z (n) ) 1 ∀δ > 0 ∃r 0 ∀r > r 0 : 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 event A, with probability 1 for random training data.

From Proposition 1 we know, that

EQUATION

with probability 1. We denote the event that (38) holds as B. Since P (A c ) = P (B c ) = 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 > n 0 , lim

EQUATION

with probability 1. Finally, we combine this result with (38), obtaining that

EQUATION

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 (n i ) i and for any i, (r i,j ) j be increasing sequences of positive integers. Define v i,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

(i k ) k , (j k ) k so that v i k ,j k → θ 0 as k → ∞.

From the triangle inequality, it follows that for any > 0 there exists k 0 so that ∀k > k 0

|µ i k ,j k (v i k ,j k ) -l(θ 0 )| ≤ |µ i k ,j k (v i k ,j k ) -l(v i k ,j k )| + |l(v i k ,j k ) -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 ) -µ i k ,j k (v i k ,j k )] + [µ i k ,j k (v i k ,j k ) -µ i k ,j k (θ)] + [µ i k ,j k (θ) -l(θ)] ≤ 3 (42)

which again holds for k large enough. This due to (41),

µ i k ,j k (v i k ,j k ) -µ i k ,j k (θ) ≤ 0 since v i k ,j k is

the minimizer of µ i k ,j k , and µ i k ,j k (θ) -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 (n i ) i , (r i,j ) j and (i k ) k , (j k ) k were chosen arbitrarily, every limit point of a sequence (v i k ,j k ) 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 (n i ) i , (r i,j ) j and

(i k ) k , (j k ) k such that ||(v i k ,j k ) k -θ|| 2 >

for all k and θ ∈ Θ * . However, by the previous argument the limit point of the any sequence r) where chosen from a set with probability mass of 1, we can conclude our proposition that lim 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 (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.

EQUATION

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

EQUATION

wherein w k (x; θ) denote the weight, µ k (x; θ) the mean and σ 2 k (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 w k (x; θ) must resemble a categorical distribution, i.e. it must hold that K k=1 w k (x; θ) = 1 and w k (x; θ) ≥ 0 ∀k. To satisfy the conditions, the softmax linearity is used for the output neurons corresponding to w k (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.

D.1.1 KERNEL MIXTURE NETWORK:

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:

EQUATION

As previously, the weights w k,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 {y i } n i=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.

D.2 NORMALIZING FLOW NETWORK:

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 Z 0 through a series of invertible, parametrized mappings f = f N • • • • • f 1 into a successively more complex distribution p(f (Z 0 )). The PDF of samples z N ∼ p(f (Z 0 )) can be evaluted using the change-ofvariable formula:

EQUATION

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(Z N ) to the base distribution p(Z 0 ) by setting f -1 i (z i ) = f i (z i ). 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.

E SIMULATED DENSITIES AND DATASETS E.1 SKEWNORMAL:

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 (Anděl et al., 1984) . The parameters (ξ, ω, α) of the skew normal distribution are functionally dependent on x. Specifically, the functional dependencies are the following: 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.

EQUATION

E.2 GAUSSIAN MIXTURE:

The joint distribution p(x, y) follows a Gaussian Mixture Model in R 4 with 5 Gaussian components, i.e. K = 5. We assume that x ∈ R 2 and y ∈ R 2 can be factorized, i.e. 

When x and y can be factorized as in (52), the conditional density p(y|x) can be derived in closed form:

EQUATION

wherein the mixture weights are a function of x:

EQUATION

For details and derivations we refer the interested reader to Guang Sung ( 2004) and Gilardi et al. (2002) . The weights w k are sampled from a uniform distribution U (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.

E.3 EURO STOXX 50 DATA:

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 ⊆ R 14 .

:

Note that an l2 regularizer and weight decay are not equivalent since we use the adaptive learning rate technique Adam. SeeLoshchilov & Hutter (2019) for details.

Figures and Tables:

Figure 1: Conditional MDN density estimate (red) and true conditional density (green) for different noise regularization intensities h z ∈ {0.0, 0.05, 0.2}. The MDN has been fitted with 3000 samples drawn from a conditional Gaussian.

zj = z i + 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 = {z 1 , ..., z n }, noise intensity h Require: learning rate α, mini-batch size m 1: Initialize θ 2: while θ not converged do 3: Sample minibatch {z 1 , ..., z m } ⊂ D

with arbitrary precision. Given a fixed data set D of size n > n 0 , this means that lim r→∞ 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, h n fulfills (11) and Θ is compact. Then, lim n→∞ lim r→∞ sup θ∈Θ l(h) n,r (θ) -l(θ) = 0 (13) almost surely.

Figure 2: Comparison of different noise intensity schedules h n and their implications on the performance of various CDE models across different training set sizes. NYC Taxi: Drop-off locations of Manhattan taxi trips conditioned on the pickup location, weekday and time (d x = 6, d y = 2).

Figure 3: Comparison of various regularization methods for three neural network based CDE models. The models are trained with simulated data sets of different sizes.

-p(z)|dz → 0. We employ the results in Devroye (1983) for our purposes, leading us to Proposition 1.

-x * (α high -α low ) (49) y ∼ SkewN ormal ξ(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.

Figure 4: Conditional density simulation models. Conditional probability densities corresponding to the different simulation models. The coloured graphs represent the probability densities p(y|x), conditioned on different values of x.

(y|µ y,k , Σ y,k )N (x|µ x,k , Σ x,k )

. Comparison of conditional density estimators across 5 data sets. Reported is the test loglikelihood and its respective standard deviation (higher log-likelihood values are better).

Table content: <table><tr><td/><td colspan="2">Euro Stoxx NCY Taxi</td><td>Boston</td><td>Conrete</td><td>Energy</td></tr><tr><td>num. train obs.</td><td>2536</td><td>8000</td><td>405</td><td>824</td><td>615</td></tr><tr><td>MDN</td><td colspan="5">4.00±0.03 5.41±0.02 -2.39±0.02 -2.89±0.03 -1.04±0.05</td></tr><tr><td>KMN</td><td colspan="5">3.98±0.03 5.42±0.02 -2.44±0.02 -3.06±0.03 -1.59±0.09</td></tr><tr><td>NFN</td><td colspan="5">4.00±0.03 5.12±0.03 -2.40±0.04 -2.93±0.02 -1.23±0.06</td></tr><tr><td>LSCDE</td><td colspan="5">3.44±0.10 4.85±0.02 -2.78±0.00 -3.63±0.00 -2.16±0.02</td></tr><tr><td>CKDE R.O.T.</td><td colspan="5">3.36±0.01 4.87±0.02 -3.12±0.03 -3.78±0.02 -2.90±0.01</td></tr><tr><td colspan="6">CKDE CV-ML 3.87±0.01 5.27±0.06 -2.76±0.26 -3.35±0.13 -1.14±0.02</td></tr><tr><td>NKDE R.O.T</td><td colspan="5">3.16±0.02 4.34±0.04 -3.52±0.05 -4.08±0.02 -3.35±0.03</td></tr><tr><td colspan="6">NKDE CV-ML 3.41±0.02 4.93±0.08 -3.34±0.13 -3.93±0.05 -2.21±0.12</td></tr></table>

). If we set Θ * = {θ * } and restrict minimization over l(h)

Table content: <table><tr><td/><td>n,r</td></tr><tr><td>to the local region, then</td><td>θ(h) n,r converges to Θ  *  as n, r → ∞ in the sense of Theorem 2.</td></tr><tr><td colspan="2">D CONDITIONAL DENSITY ESTIMATION MODELS</td></tr><tr><td colspan="2">D.1 MIXTURE DENSITY NETWORK</td></tr><tr><td colspan="2">Mixture Density Networks (MDNs) combine conventional neural networks with a mixture density</td></tr><tr><td colspan="2">model for the purpose of estimating conditional distributions p(y|x)</td></tr></table>

Peer review (iclr2020):
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. Pros: (i) The paper is well written. (ii) The method is simple to implement. (iii) The authors demonstrate a clear intuition of what the method does, i.e., as a form of smoothing. Cons: (i) The method itself is not novel. Adding noise to for regularization is a quite common technique used in many different applications. (ii) The experiments performed in the paper are all very small scaled. Not very convincing. (iii) Continuing the last point, it is unclear whether such technique can be scaled up. The method in another view is to replace the empirical distribution with a kernel density estimate. We know kernel methods don*t scale. It*s hence questionable whether the method can give any benefits in large scale dataset. (iv) If the use of noise regularization is to smooth the density function, how does it compare to methods that enforces Lipschitz continuity? e.g., gradient penalty, or construct a lipschtz constrained network such as in [1]. [1] Sorting out lipschitz function approximation. Anil et al.