# Noise Regularization for Conditional Density Estimation

Jonas Rothfuss  Fabio Ferreira  Simon Boehm  Simon Walther  Maxim Ulrich
Tamim Asfour  Andreas Krause

###### Abstract

Capturing statistical relationships beyond the conditional mean is crucial in
many applications. To this end, conditional density estimation (CDE) aims to
learn the full conditional probability density from data. Though expressive,
neural network based CDE models can suffer from severe over-fitting when
trained with the maximum likelihood objective. Their particular structure
renders classical regularization in the parameter space ineffective. To
address this challenge, we propose a model-agnostic noise regularization
method for CDE that adds carefully controlled random perturbations to the data
during training. We prove that the proposed approach corresponds to a
smoothness regularization and establish its asymptotic consistency. Our
extensive experiments show that noise regularization consistently outperforms
other regularization methods across a range of neural CDE models. Furthermore,
we demonstrate the effectiveness of noise regularized neural CDE over
classical non- and semi-parametric methods, even when training data is scarce.

Machine Learning, ICML

##  1 Introduction

While regression analysis aims to describe the conditional mean
$\mathbb{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
observations 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 flexible 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 a Gaussian noise model, CDE uses expressive
distribution families to model deviations from the conditional mean. Hence,
the over-fitting problem tends to be even more severe in CDE than in
regression. Standard regularization of the neural network weights such as
weight decay (Pratt & Hanson, 1989) has been shown 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 hard to analyze.

The lack of an effective regularization scheme renders neural network based
CDE impractical in most scenarios where data is scarce. As a result, classical
non- and semi-parametric CDE tends to be the primary method in application
areas such as econometrics (Zambom & Dias, 2013).

To address this issue, we propose and analyze the use of noise regularization,
an approach well-studied in the context of regression and classification, for
the purpose of CDE. By adding small, carefully controlled, 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 a penalty on large
second derivatives in the training point locations which results in an
inductive bias towards smoother density estimates. Moreover, under mild
regularity conditions, we show that the proposed regularization scheme is
_consistent_ , converging to the unbiased maximum likelihood estimator. This
does not only support the soundness of the proposed method but also provides
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 consistently across various data sets. In a
comprehensive benchmark study, we demonstrate that, with noise regularization,
neural network based CDE is able to significantly improve upon state-of-the
art non-parametric estimators, even when only 400 training observations are
available. This is a relevant, and perhaps surprising, finding since non-
parametric CDE is considered one of the primary approaches for settings with
scarce and noisy data (Zambom & Dias, 2013). By using non-parametric
regularization for training parametric high-capacity models, we are able to
combine inductive biases from both worlds, making neural networks the
preferable CDE method, even for low-dimensional and small-scale tasks.

##  2 Background

##### Density Estimation.

Let $X$ be a random variable with probability density function (PDF) $p(x)$
defined over the domain $\mathcal{X}\subseteq\mathbb{R}^{d_{x}}$. Given a
collection $\mathcal{D}=\\{x_{1},...,x_{n}\\}$ of observations sampled from
$p(x)$, the goal is to find a good estimate $\hat{f}(x)$ of the true density
function $p$. In parametric estimation, the PDF $\hat{f}$ is assumed to belong
to a parametric family
$\mathcal{F}=\\{\hat{f}_{\theta}(\cdot)|\theta\in\Theta\\}$ where the density
function is described by a finite dimensional parameter $\theta\in\Theta$. The
standard method for estimating $\theta$ is maximum likelihood estimation
(MLE), wherein $\theta$ is chosen so that the likelihood of the data
$\mathcal{D}$ is maximized. This is equivalent to minimizing the Kullback-
Leibler divergence between the empirical data distribution
$p_{\mathcal{D}}(x)=\frac{1}{n}\sum_{i=1}^{n}\delta(||x-x_{i}||)$ (i.e.,
mixture of point masses in the observations $x_{i}$) and the parametric
distribution $\hat{f}_{\theta}$:

$\text{arg}\max_{\theta\in\Theta}\sum_{i=1}^{n}\log\hat{f}_{\theta}(x_{i})=\text{arg}\min_{\theta\in\Theta}\mathcal{D}_{KL}(p_{\mathcal{D}}||\hat{f}_{\theta})\vspace{-4pt}$
(1)

From a geometric perspective, (1) can be viewed as an orthogonal projection of
$p_{\mathcal{D}}(x)$ onto $\mathcal{F}$ w.r.t. the 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
$\hat{q}(x)=\frac{1}{nh^{d}}\sum_{i=1}^{n}K\left(\frac{x-x_{i}}{h}\right).$
One popular choice of $K(\cdot)$ is a Gaussian
$K(z)=(2\pi)^{-\frac{d}{2}}\exp{(-\frac{1}{2}z^{2}})$. Beyond the appropriate
choice of $K(\cdot)$, 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
$\mathcal{X}\subseteq\mathbb{R}^{d_{x}}$ and
$\mathcal{Y}\subseteq\mathbb{R}^{d_{y}}$ 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 $\mathcal{D}=\\{(x_{n},y_{n})\\}_{n=1}^{N}$ drawn from the joint
distribution $(x_{n},y_{n})\sim p(x,y)$, the aim of conditional density
estimation (CDE) is to find an estimate $\hat{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)$:

$\mathbb{E}_{p(x)}\left[\mathcal{D}_{KL}(p(y|x)||\hat{f}(y|x))\right]=\mathbb{E}_{p(x,y)}\left[\frac{\log
p(y|x)}{\log\hat{f}(y|x)}\right]\qquad\qquad$ (2)

Corresponding to (1), we refer to the minimization of (2) w.r.t. $\theta$ as
_conditional M-projection_. Given a dataset $\mathcal{D}$ drawn i.i.d. from
$p(x,y)$, the conditional MLE following from (2) can be stated as

$\vspace{-4pt}\theta^{*}=\text{arg}\min_{\theta}-\sum_{i=1}^{n}\log\hat{f}_{\theta}(y_{i}|x_{i})$
(3)

##  3 Related work

The first part of this section discusses 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 studies nonparametric kernel density
estimators (KDEs) (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)
(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 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 approximate 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 (Sietsma & Dow, 1991; Burges
& Schölkopf, 1996; Goodfellow et al., 2015; Yuan et al., 2017). Our approach
transforms the observations by adding small random perturbations. While this
form of regularization has been studied in the context of regression and
classification (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 $\theta$ is ill suited.
As can be observed in Figure 1, simply minimizing the negative log-likelihood
of the data leads to severe over-fitting. 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) 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.

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

###  4.1 The algorithm

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
$\tilde{x}_{i}=x_{i}+\xi_{x}$ and $\tilde{y}_{i}=y_{i}+\xi_{y}$ where the
perturbation vectors are sampled from noise distributions $K_{x}(\xi_{x})$ and
$K_{y}(\xi_{y})$. Further, we choose the noise to be zero centered as well as
i.i.d among the data dimensions, with standard deviation $h$:

$\mathbb{E}_{\xi\sim
K(\xi)}\left[\xi\right]=0~{}~{}~{}\text{and}~{}~{}~{}\mathbb{E}_{\xi\sim
K(\xi)}\left[\xi\xi^{\top}\right]=h^{2}I\vspace{-2pt}$ (4)

This can be seen as a form of 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 ${\mathcal{D}}$, and
adding a random perturbation vector to the respective data point. This
procedure is formalized in Algorithm 1.

For notational brevity, we set
${\mathcal{Z}}:={\mathcal{X}}\times{\mathcal{Y}}$,
$z:=(x^{\top},y^{\top})^{\top}$ and denote
$\hat{f}_{\theta}(z):=\hat{f}_{\theta}(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

0: ${\mathcal{D}}=\\{z_{1},...,z_{n}\\}$, noise intensity $h$

0: number of perturbed samples $r$,

1: for j = 1 to r do

2: Select $i\in\\{1,...,n\\}$ with equal prob.

3: Draw perturbation $\xi\sim K$

4: Set $\tilde{z}_{j}=z_{i}+h\xi$

5: end for

6: return
$\operatorname*{arg\,min}_{\theta\in\Theta}-\sum_{j=1}^{r}\log\hat{f}_{\theta}(\tilde{z}_{j})$

Algorithm 2 (Conditional) MLE with Noise Regularization - Mini-Batch Gradient
Descent

0: ${\mathcal{D}}=\\{z_{1},...,z_{n}\\}$, noise intensity $h$

0: learning rate $\alpha$, mini-batch size $m$

1: Initialize $\theta$

2: while $\theta$ not converged do

3: Sample minibatch $\\{z_{1},...,z_{m}\\}\subset{\mathcal{D}}$

4: for j = 1 to m do

5: Draw perturbation $\xi\sim K$

6: Set $\tilde{z}_{j}=z_{j}+h\xi$

7: end for

8:
$\theta\leftarrow\theta+\alpha\nabla_{\theta}\sum_{j=1}^{m}\log\hat{f}_{\theta}(\tilde{z}_{j})$

9: end while

10: return optimized parameter $\theta$

When considering highly flexible parametric families such as Mixture Density
Networks (MDNs) (Bishop, 1994), the maximum likelihood solution in line 6 of
Algorithm 1 is no longer tractable. In such case, one typically resorts to
stochastic optimization techniques such as mini-batch gradient descent and
variations thereof. The generic procedure in Algorithm 1 can be transformed
into a simple extension 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.2 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(\mathcal{D})$ be a loss function
over a set of data points $\mathcal{D}=\\{z_{1},...,z_{n}\\}$, which can be
partitioned into a sum of losses $l(\mathcal{D})=\sum_{i=1}^{n}l(z_{i})$,
corresponding to each data point $z_{i}$: The expected loss $l(z_{i}+\xi)$,
resulting from adding random perturbations, can be approximated by a second
order Taylor expansion around $z_{i}$. Using the assumption about $\xi$ in
(4), the expected loss an be written as

$\displaystyle\mathbb{E}_{K(\xi)}\left[l(z_{i}+\xi)\right]=$ $\displaystyle
l(z_{i})+\frac{1}{2}\mathbb{E}_{K(\xi)}\left[\xi^{\top}\mathbf{H}^{(i)}\xi\right]+\mathcal{O}(\xi^{3})$
$\displaystyle\approx$ $\displaystyle
l(z_{i})+\frac{h^{2}}{2}\text{tr}(\mathbf{H}^{(i)})$

where $l({\bm{z}}_{i})$ is the loss without noise and
$\mathbf{H}^{(i)}=\frac{\partial^{2}l}{\partial z^{2}}(z)\big{\rvert}_{z_{i}}$
the Hessian of $l$ w.r.t $z$, evaluated at $z_{i}$. Assuming that the noise
$\xi$ is small in its magnitude, $\mathcal{O}(\xi^{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
$\hat{f}_{\theta}(y|x)$, the loss function coincides with the negative
conditional log-likelihood
$l(y_{i},x_{i})=-\log\hat{f}_{\theta}(y_{i}|x_{i})$. Let the standard
deviation of the additive data noise $\xi_{x}$, $\xi_{y}$ be $h_{x}$ and
$h_{y}$ respectively. Maximum likelihood estimation (MLE) with data noise is
equivalent to minimizing the loss

$\displaystyle\begin{split}l(\mathcal{D})\approx&-\sum_{i=1}^{n}\log\hat{f}_{\theta}(y_{i}|x_{i})-\frac{h^{2}_{y}}{2}\sum_{i=1}^{n}\sum_{j=1}^{d_{y}}\frac{\partial^{2}\log\hat{f}_{\theta}(y|x)}{\partial
y^{(j)}\partial y^{(j)}}\bigg{|}_{\begin{subarray}{c}x_{i}\\\
y_{i}\end{subarray}}\\\
&-\frac{h^{2}_{x}}{2}\sum_{i=1}^{n}\sum_{j=1}^{d_{x}}\frac{\partial^{2}\log\hat{f}_{\theta}(y|x)}{\partial
x^{(j)}\partial x^{(j)}}\bigg{|}_{\begin{subarray}{c}x_{i}\\\
y_{i}\end{subarray}}\end{split}$

Hereby, the first term corresponds to the standard MLE objective, while the
other two terms constitute a form of smoothness regularization. The second
term of $l(\mathcal{D})$ penalizes concavity of the conditional log density
estimate $\log\hat{f}_{\theta}(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,
thus smoothing 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_{x}^{2}$ and $h_{y}^{2}$) 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.3 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 $\hat{f}_{\theta}(z):\mathbb{R}^{d_{z}}\times\Theta\rightarrow(0,\infty)$
a continuous function of $z$ and $\theta$. Moreover, we assume that the
parameter space $\Theta$ 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(\theta)=-\mathbb{E}_{p(z)}\left[\log\hat{f}_{\theta}(z)\right]$ (5)

As we typically just have a finite number of samples from $p(z)$, the
respective empirical estimate
$\hat{l}_{n}(\theta)=-\frac{1}{n}\sum_{i=1}^{n}\log\hat{f}_{\theta}(z_{i}),~{}z_{i}\stackrel{{\scriptstyle
i.i.d}}{{\sim}}p(z)$ is used as training objective. Note that we now define
the loss as function of $\theta$, and, for fixed $\theta$, treat
$l_{n}(\theta)$ 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_{\theta\in\Theta}|\hat{l}_{n}(\theta)-l(\theta)|\xrightarrow{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 $\\{\tilde{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
$\hat{q}_{n}^{(h)}(z)=\frac{1}{n}\sum_{i=1}^{n}\frac{1}{h^{d_{z}}}K\left(\frac{z-z_{i}}{h}\right).$
Hence, MLE with variable noise can be understood as

  1. 1.

forming a kernel density estimate $\hat{q}_{n}^{(h)}$ of the data,

  2. 2.

followed by a (conditional) M-projection of $\hat{q}_{n}^{(h)}$ onto the
parametric family.

Hereby, step 2 aims to find the $\theta^{*}$ that minimizes the following
objective:

$l_{n}^{(h)}(\theta)=-\mathbb{E}_{\hat{q}_{n}^{(h)}(z)}\left[\log\hat{f}_{\theta}(z)\right]\vspace{-2pt}$
(6)

Since (6) is generally intractable, $r$ samples are drawn from the kernel
density estimate, forming the following Monte Carlo approximation of (6) which
corresponds to the loss in line 6 Algorithm 1:

$\hat{l}_{n,r}^{(h)}(\theta)=-\frac{1}{r}\sum_{j=1}^{r}\log\hat{f}_{\theta}(\tilde{z}_{j})~{},\quad\tilde{z}_{j}\sim\hat{q}_{n}^{(h)}(z)\vspace{-2pt}$
(7)

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_{\theta\in\Theta}\left|\hat{l}_{n,r}^{(h)}(\theta)-l(\theta)\right|\xrightarrow{a.s.}0$
as $n,r\rightarrow\infty$. 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:

$\displaystyle\begin{split}\sup_{\theta\in\Theta}\left|\hat{l}_{n,r}^{(h)}(\theta)-l(\theta)\right|\leq&\sup_{\theta\in\Theta}\left|\hat{l}_{n,r}^{(h)}(\theta)-l_{n}^{(h)}(\theta)\right|\\\
&+\sup_{\theta\in\Theta}\left|l_{n}^{(h)}(\theta)-l(\theta)\right|\end{split}\vspace{-2pt}$
(8)

Note that $\hat{l}_{n,r}^{(h)}(\theta)$ is based on samples from the KDE,
which are obtained by adding random noise vectors $\xi\sim K(\cdot)$ 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_{\theta\in\Theta}|\hat{l}_{n,r}^{(h)}(\theta)-l_{n}^{(h)}(\theta)|$
arbitrary small by the uniform law of large numbers. In order to make
$\sup_{\theta\in\Theta}|l_{n}^{(h)}(\theta)-l(\theta)|$ small in the limit
$n\rightarrow\infty$, 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}(\theta)$ and $l(\theta)$ 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 (26) and (27) in Appendix C are satisfied,
and that

$\lim_{n\rightarrow\infty}~{}h_{n}=0,\qquad\lim_{n\rightarrow\infty}~{}n(h_{n})^{d}=\infty$
(9)

Then,

$\lim_{n\rightarrow\infty}\sup_{\theta\in\Theta}\left|l_{n}^{(h)}(\theta)-l(\theta)\right|=0$
(10)

almost surely.

In (9) 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.4 for
discussion). The result in (10) 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}^{(h)}(\theta)$ is
intractable and, thus, replaced by its sample estimate $\hat{l}_{n,r}^{(h)}$.
Since we can draw an arbitrary amount of samples from $\hat{q}_{n}^{(h)}$, we
can approximate $l_{n}^{(h)}(\theta)$ with arbitrary precision. Given a fixed
data set ${\mathcal{D}}$ of size $n>n_{0}$, this means that
$\lim_{r\rightarrow\infty}\sup_{\theta\in\Theta}\left|\hat{l}_{n,r}^{(h)}(\theta)-l_{n}^{(h)}(\theta)\right|=0$
almost surely, by (27) and the uniform law of large numbers. Since our
original goal was to also show consistency for $n\rightarrow\infty$, this
result is combined with Proposition 1, obtaining the following consistency
theorem.

###### Theorem 1

Suppose the regularity conditions (26) and (27) are satisfied, $h_{n}$
fulfills (9) and $\Theta$ is compact. Then,

$\lim_{n\rightarrow\infty}\underset{r\rightarrow\infty}{\overline{\lim}}\sup_{\theta\in\Theta}\left|\hat{l}_{n,r}^{(h)}(\theta)-l(\theta)\right|=0$
(11)

almost surely.

In that, $\overline{\lim}$ used to denote the limit superior (“lim sup”) of a
sequence. Training a (conditional) density model with noise regularization
means minimizing $\hat{l}_{n,r}^{(h)}(\theta)$ w.r.t. $\theta$. As result of
this optimization, one obtains a parameter vector $\hat{\theta}_{n,r}^{(h)}$,
which we hope is close to the minimizing parameter $\bar{\theta}$ of the ideal
objective function $l(\theta)$. In the following, we establish consistency
results, similar to Theorem 1, in the parameter space. For that, we first have
to formalize the concept of closeness and optimality in the parameter space.
Since a minimizing parameter $\bar{\theta}$ of $l(\theta)$ may not be unique,
we define $\Theta^{*}=\\{\theta^{*}~{}|~{}l(\theta^{*})\leq
l(\theta)~{}\forall\theta\in\Theta\\}$ as the set of global minimizers of
$l(\theta)$, and
$d(\theta,\Theta^{*})=\min_{\theta^{*}\in\Theta^{*}}\\{||\theta-\theta^{*}||_{2}\\}$
as the distance of an arbitrary parameter $\theta$ to $\Theta^{*}$. Based on
these definitions, it can be shown that Algorithm 1 is consistent in a sense
that the minimizer of $\hat{\theta}_{n,r}^{(h)}$ converges almost surely to
the set of optimal parameters $\Theta^{*}$.

###### Theorem 2

Suppose the regularity conditions (26) and (27) are satisfied, $h_{n}$
fulfills (9) and $\Theta$ is compact. For $r>0$ and $n>n_{0}$, let
$\hat{\theta}_{n,r}^{(h)}\in\Theta$ be a global minimizer of the empirical
objective $\hat{l}_{n,r}^{(h)}$. Then

$\lim_{n\rightarrow\infty}\underset{r\rightarrow\infty}{\overline{\lim}}d(\hat{\theta}_{n,r}^{(h)},\Theta^{*})=0$
(12)

almost surely.

Note that Theorem 2 considers global optimizers, but equivalently holds for
compact neighborhoods of a local minimum $\theta^{*}$ (see discussion in
Appendix C).

###  4.4 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_{n}^{(h)}(\theta)-l(\theta)|$ is minimized, which is practically not
feasible since $l(\theta)$ is intractable. Inequality (28) gives as an upper
bound on this quantity, suggesting to minimize $l_{1}$ distance between the
kernel density estimate $q_{n}^{(h)}$ 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\hat{\sigma}n^{-\frac{1}{4+d}}$ (Silverman,
1986). In that, $\hat{\sigma}$ denotes the estimated standard deviation of the
data, $n$ the number of data points and $d$ the dimensionality of
$\mathcal{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 (9) 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^{-\frac{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^{-\frac{1}{1+d}}$ may
be appropriate.

##  5 Experiments

We now perform 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 perturbations, i.e., $K(\xi)=\mathcal{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
($d_{x}=2,d_{y}=2$) and a Skew-Normal distribution whose parameters are
functionally dependent on $x$ ($d_{x}=1,d_{y}=1$). In terms of real-world
data, we use the following three data sources. EuroStoxx: Daily returns of the
Euro Stoxx 50 index conditioned on various stock return factors. NYC Taxi:
Drop-off locations of Manhattan taxi trips conditioned on the pickup location,
weekday and time. UCI datasets: Boston Housing , Concrete and Energy datasets
from the UCI machine learning repository (Dua & Graff, 2017). 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 on our supplementary
website111https://sites.google.com/view/noisereg/.

###  5.1 Noise intensity schedules

We complement the discussion in Section 4.4 with an empirical investigation of
different schedules of $h_{n}$. In particular, we compare a) the rule of thumb
$h_{n}\propto n^{-\frac{1}{4+d}}$ b) a square root decay schedule
$h_{n}\propto n^{-\frac{1}{1+d}}$ c) a constant bandwidth
$h_{n}=const.\in(0,\infty)$ 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.

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.

First, we observe that bandwidth rates that conform with the decay conditions
seem to converge in performance to the non-regularized maximum likelihood
estimator (red) as $n$ becomes large. This validates the theoretical result of
Theorem 1. Second, a fixed bandwidth across $n$ (green), violating (9),
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)222Note that an $l_{2}$
regularizer and weight decay are not equivalent since we use the adaptive
learning rate technique Adam. See Loshchilov & Hutter (2019) for details., 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.

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.  |  | Euro Stoxx | NYC Taxi | Boston | Concrete | Energy
---|---|---|---|---|---|---
MDN | noise (ours) | 3.94$\pm$0.03 | 5.25$\pm$0.04 | -2.49$\pm$0.11 | -2.92$\pm$0.08 | -1.04$\pm$0.09
weight decay | 3.78$\pm$0.06 | 5.07$\pm$0.04 | -3.29$\pm$0.32 | -3.33$\pm$0.14 | -1.21$\pm$0.10
l1 reg. | 3.19$\pm$0.19 | 5.00$\pm$0.05 | -4.01$\pm$0.36 | -3.87$\pm$0.29 | -1.44$\pm$0.22
l2 reg. | 3.16$\pm$0.21 | 4.99$\pm$0.04 | -4.64$\pm$0.52 | -3.84$\pm$0.26 | -1.55$\pm$0.26
Bayes | 3.26$\pm$0.43 | 5.08$\pm$0.03 | -3.46$\pm$0.47 | -3.19$\pm$0.21 | -1.25$\pm$0.23
KMN | noise (ours) | 3.92$\pm$0.01 | 5.39$\pm$0.02 | -2.52$\pm$0.08 | -3.09$\pm$0.06 | -1.62$\pm$0.06
weight decay | 3.85$\pm$0.03 | 5.31$\pm$0.02 | -2.69$\pm$0.15 | -3.15$\pm$0.06 | -1.79$\pm$0.12
l1 reg. | 3.76$\pm$0.04 | 5.39$\pm$0.02 | -2.75$\pm$0.13 | -3.25$\pm$0.07 | -1.82$\pm$0.10
l2 reg. | 3.71$\pm$0.05 | 5.37$\pm$0.02 | -2.66$\pm$0.13 | -3.18$\pm$0.07 | -1.79$\pm$0.13
Bayes | 3.33$\pm$0.02 | 4.47$\pm$0.02 | -3.40$\pm$0.11 | -4.08$\pm$0.05 | -3.65$\pm$0.07
NFN | noise (ours) | 3.90$\pm$0.01 | 5.20$\pm$0.03 | -2.48$\pm$0.11 | -3.03$\pm$0.13 | -1.21$\pm$0.08
weight decay | 3.82$\pm$0.06 | 5.19$\pm$0.03 | -3.12$\pm$0.39 | -3.12$\pm$0.14 | -1.22$\pm$0.16
l1 reg. | 3.50$\pm$0.10 | 5.12$\pm$0.05 | -12.6$\pm$12.8 | -3.91$\pm$0.52 | -1.29$\pm$0.16
l2 reg. | 3.50$\pm$0.09 | 5.13$\pm$0.05 | -14.2$\pm$9.60 | -3.99$\pm$0.66 | -1.34$\pm$0.19
Bayes | 3.34$\pm$0.33 | 5.10$\pm$0.03 | -5.99$\pm$2.45 | -3.55$\pm$0.46 | -1.11$\pm$0.22
Table 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).  | Euro Stoxx | NCY Taxi | Boston | Conrete | Energy
---|---|---|---|---|---
num. train obs. | 2536 | 8000 | 405 | 824 | 615
MDN (ours) | 4.00$\pm$0.03 | 5.41$\pm$0.02 | -2.39$\pm$0.02 | -2.89$\pm$0.03 | -1.04$\pm$0.05
KMN (ours) | 3.98$\pm$0.03 | 5.42$\pm$0.02 | -2.44$\pm$0.02 | -3.06$\pm$0.03 | -1.59$\pm$0.09
NFN (ours) | 4.00$\pm$0.03 | 5.12$\pm$0.03 | -2.40$\pm$0.04 | -2.93$\pm$0.02 | -1.23$\pm$0.06
LSCDE | 3.44$\pm$0.10 | 4.85$\pm$0.02 | -2.78$\pm$0.00 | -3.63$\pm$0.00 | -2.16$\pm$0.02
CKDE R.O.T. | 3.36$\pm$0.01 | 4.87$\pm$0.02 | -3.12$\pm$0.03 | -3.78$\pm$0.02 | -2.90$\pm$0.01
CKDE CV-ML | 3.87$\pm$0.01 | 5.27$\pm$0.06 | -2.76$\pm$0.26 | -3.35$\pm$0.13 | -1.14$\pm$0.02
NKDE R.O.T | 3.16$\pm$0.02 | 4.34$\pm$0.04 | -3.52$\pm$0.05 | -4.08$\pm$0.02 | -3.35$\pm$0.03
NKDE CV-ML | 3.41$\pm$0.02 | 4.93$\pm$0.08 | -3.34$\pm$0.13 | -3.93$\pm$0.05 | -2.21$\pm$0.12
Table 2: Comparison of conditional density estimators across 5 data sets. The
neural network based models (MDN, KMN, NFN) are trained with noise
regularization. Reported is the test log-likelihood and its standard deviation
across seeds.

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 yields similar performance
across the 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
methods mostly fail to do so.

Next, we consider real world data. We use 5-fold cross-validation on the
training set 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 respective data set. 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 question how such estimators
compete against popular non-parametric methods in small data regimes. In
particular, we compare to Conditional Kernel Density Estimation (CKDE) (Li &
Racine, 2007), $\epsilon$-Neighborhood Kernel Density Estimation (NKDE), and
Least-Squares Conditional Density Estimation (LSCDE) (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, the three neural network based methods outperform all
of the non- and semi-parametric methods. In summary, neural network based CDE
models, trained with noise regularization, work surprisingly well, even when
training data is scarce and fairly low-dimensional such as in case of the
Boston Housing data set $(d_{x}=13,d_{y}=1)$. We hypothesize this may be due
to the fact, that, when using non-parametric regularization for training
parametric high-capacity models, we combine favorable inductive biases of both
KDE and neural networks.

##  6 Conclusion

The paper proposes a regularization technique for neural CDE that adds random
perturbations to the data during training. It can be seamlessly integrated
into standard stochastic optimization and is model-agnostic. We show that the
proposed noise regularization inherits an inductive bias from non-parametric
KDE, effectively smoothing the estimated conditional density. Our experiments
demonstrate that it consistently outperforms other regularization methods
across models and datasets. Moreover, empirical results suggest that, when
trained with noise regularization, neural network based models are able to
significantly improve upon state-of-the art non-parametric estimators. This
makes neural CDE the preferable method, even in settings where training data
is low-dimensional and scarce.

#### Acknowledgments

This project has received funding from the European Research Council (ERC)
under the European Union’s Horizon 2020 research and innovation program grant
agreement No 815943. Simon Walter has been supported by the Konrad-Adenauer-
Stiftung. Finally, we thank Sara van der Geer for her advise regarding the
consistency proofs.

## References

  * Ambrogioni et al. (2017) Ambrogioni, L., Güçlü, U., van Gerven, M. A. J., and Maris, E.  The Kernel Mixture Network: A Nonparametric Method for Conditional Density Estimation of Continuous Random Variables.  2017\.  URL http://arxiv.org/abs/1705.07111.
  * Anděl et al. (1984) Anděl, J., Netuka, I., and Zvára, K.  On Threshold Autoregressive Processes.  _Kybernetica_ , 20(2):89–106, 1984.  URL https://dml.cz/bitstream/handle/10338.dmlcz/124493/Kybernetika_20-1984-2_1.pdf.
  * Bishop (1994) Bishop, C. M.  Mixture Density Networks.  1994\.
  * Bishop (1995) Bishop, C. M.  Training with Noise is Equivalent to Tikhonov Regularization.  _Neural Computation_ , 7(1):108–116, 1995.  ISSN 0899-7667.  doi: 10.1162/neco.1995.7.1.108.  URL http://www.mitpressjournals.org/doi/10.1162/neco.1995.7.1.108.
  * Blei et al. (2017) Blei, D. M., Kucukelbir, A., and McAuliffe, J. D.  Variational Inference: A Review for Statisticians.  _Journal of the American Statistical Association_ , 112(518):859–877, 2017.  doi: 10.1080/01621459.2017.1285773.  URL https://doi.org/10.1080/01621459.2017.1285773.
  * Burges & Schölkopf (1996) Burges, C. J. C. and Schölkopf, B.  Improving the accuracy and speed of support vector machines.  In _NIPS_ , pp. 375–381. MIT Press, 1996.  URL https://dl.acm.org/citation.cfm?id=2999034.
  * Cao et al. (1994) Cao, R., Cuevas, A., and González Manteiga, W.  A comparative study of several smoothing methods in density estimation.  _Computational Statistics & Data Analysis_, 17(2):153–176, 2 1994.  ISSN 0167-9473.  doi: 10.1016/0167-9473(92)00066-Z.  URL https://www.sciencedirect.com/science/article/pii/016794739200066Z?via%3Dihub.
  * Chen et al. (2014) Chen, T., Fox, E. B., and Guestrin, C.  Stochastic Gradient Hamiltonian Monte Carlo.  In _ICML_ , 2014.  URL https://arxiv.org/pdf/1402.4102.pdf.
  * Devroye (1983) Devroye, L.  The equivalence of weak, strong and complete convergence in L1 for kernel density estimates.  _Annals of Statistics_ , 11(3):896–904, 1983\.  URL https://pdfs.semanticscholar.org/71d9/b1c7a54cb48ab12bc3c8dcad626dc93d867b.pdf.
  * Devroye & Luc (1987) Devroye, L. and Luc.  _A course in density estimation_.  Birkhäuser, 1987.  ISBN 0817633650.  URL https://dl.acm.org/citation.cfm?id=27672.
  * Dinh et al. (2017) Dinh, L., Sohl-Dickstein, J., and Bengio, S.  Density estimation using Real NVP.  In _Proceedings of the International Conference on Learning Representations_ , 5 2017.  URL http://arxiv.org/abs/1605.08803.
  * Dua & Graff (2017) Dua, D. and Graff, C.  {UCI} Machine Learning Repository, 2017.  URL http://archive.ics.uci.edu/ml.
  * Dutordoir et al. (2018) Dutordoir, V., Salimbeni, H., Deisenroth, M. P., and Hensman, J.  Gaussian Process Conditional Density Estimation.  In _NeurIPS_ , 2018.
  * Gal & Uk (2016) Gal, Y. and Uk, Z. A.  Dropout as a Bayesian Approximation: Representing Model Uncertainty in Deep Learning Zoubin Ghahramani.  In _ICML_ , 2016.  URL http://yarin.co.
  * Gilardi et al. (2002) Gilardi, N., Bengio, S., and Kanevski, M.  Conditional Gaussian Mixture Models for Environmental Risk Mapping.  In _NNSP_ , 2002.  URL http://www.idiap.ch.
  * Goodfellow et al. (2015) Goodfellow, I. J., Shlens, J., and Szegedy, C.  Explaining and Harnessing Adversarial Examples.  In _ICLR_ , 2015.  URL http://arxiv.org/abs/1412.6572.
  * Guang Sung (2004) Guang Sung, H.  _Gaussian Mixture Regression and Classification_.  PhD thesis, 2004.  URL http://www.stat.rice.edu/~hgsung/thesis.pdf.
  * Hall et al. (1992) Hall, P., Marron, J., and Park, B. U.  Smoothed cross-validation.  _Probability Theory_ , 92:1–20, 1992.  URL https://link.springer.com/content/pdf/10.1007%2FBF01205233.pdf.
  * Hinton & Van Camp (1993) Hinton, G. E. and Van Camp, D.  Keeping Neural Networks Simple by Minimizing the Description Length of the Weights.  In _COLT_ , 1993.  URL http://www.cs.toronto.edu/~fritz/absps/colt93.pdf.
  * Holmstrom & Koistinen (1992a) Holmstrom, L. and Koistinen, P.  Using additive noise in back-propagation training.  _IEEE Transactions on Neural Networks_ , 3(1):24–38, 1992a.  ISSN 10459227.  doi: 10.1109/72.105415.  URL http://ieeexplore.ieee.org/document/105415/.
  * Holmstrom & Koistinen (1992b) Holmstrom, L. and Koistinen, P.  Using Additive Noise in Back Propagation Training.  _IEEE Transactions on Neural Networks_ , 3(1):24–38, 1992b.  ISSN 19410093.  doi: 10.1109/72.105415.  URL http://ieeexplore.ieee.org/document/105415/.
  * Hyndman et al. (1996) Hyndman, R. J., Bashtannyk, D. M., and Grunwald, G. K.  Estimating and Visualizing Conditional Densities.  _Journal of Computational and Graphical Statistics_ , 5(4):315, 12 1996.  ISSN 10618600.  doi: 10.2307/1390887.  URL https://www.jstor.org/stable/1390887?origin=crossref.
  * Kingma & Dhariwal (2018) Kingma, D. P. and Dhariwal, P.  Glow: Generative Flow with Invertible 1x1 Convolutions.  Technical report, 7 2018.  URL http://arxiv.org/abs/1807.03039.
  * Krogh & Hertz (1991) Krogh, A. and Hertz, J. A.  A Simple Weight Decay Can Improve Generalization.  In _NIPS_ , 1991.
  * Krogh & Hertz (1992) Krogh, A. and Hertz, J. A.  A Simple Weight Decay Can Improve Generalization.  Technical report, 1992.
  * Li & Racine (2007) Li, Q. and Racine, J. S.  _Nonparametric econometrics : theory and practice_.  Princeton University Press, 2007.
  * Loshchilov & Hutter (2019) Loshchilov, I. and Hutter, F.  Decoupled Weight Decay Regularization.  In _ICLR_ , 2019.  URL https://github.com/loshchil/AdamW-and-SGDW.
  * Maaten et al. (2013) Maaten, L., Chen, M., Tyree, S., and Weinberger, K.  Learning with marginalized corrupted features.  In _International Conference on Machine Learning_ , pp. 410–418, 2013.
  * Mackay (1992) Mackay, D. J. C.  A Practical Bayesian Framework for Backprop Networks.  _Neural Computation_ , 1992.  URL http://citeseerx.ist.psu.edu/viewdoc/download?doi=10.1.1.29.274&rep=rep1&type=pdf.
  * Mirza & Osindero (2014) Mirza, M. and Osindero, S.  Conditional Generative Adversarial Nets.  Technical report, 11 2014.  URL http://arxiv.org/abs/1411.1784.
  * Murphy (2012) Murphy, K. P.  _Machine Learning: A Probabilistic Perspective_.  2012\.
  * Murray & Edwards (1993) Murray, A. F. and Edwards, P. J.  Synaptic Weight Noise During MLP Learning Enhances Fault-Tolerance, Generalisation and Learning Trajectory.  In _NIPS_ , 1993.  URL https://pdfs.semanticscholar.org/b0fc/40f4a4e9db0a67bf644cd1d509044fd3c6c8.pdf.
  * Natarajan et al. (2013) Natarajan, N., Dhillon, I. S., Ravikumar, P. K., and Tewari, A.  Learning with noisy labels.  In _Advances in neural information processing systems_ , pp. 1196–1204, 2013.
  * Neal (2012) Neal, R. M.  _Bayesian learning for neural networks_ , volume 118.  Springer Science & Business Media, 2012.
  * Ng (2004) Ng, A. Y.  Feature selection, L 1 vs. L 2 regularization, and rotational invariance.  In _ICML_ , 2004.  URL http://ai.stanford.edu/~ang/papers/icml04-l1l2.pdf.
  * Nielsen (2018) Nielsen, F.  What is an information projection?  _Notices of the AMS_ , 2018.  doi: 10.1090/noti1647.  URL http://dx.doi.org/10.1090/noti1647.
  * Nowlan & Hinton (1992) Nowlan, S. J. and Hinton, G. E.  Simplifying Neural Networks by Soft Weight Sharing.  _Neural Computation_ , 1992.  URL http://www.cs.toronto.edu/~hinton/absps/sunspots.pdf.
  * Parzen (1962) Parzen, E.  On Estimation of a Probability Density Function and Mode.  _The Annals of Mathematical Statistics_ , 33(3):1065–1076, 9 1962.  ISSN 0003-4851.  doi: 10.1214/aoms/1177704472.  URL http://projecteuclid.org/euclid.aoms/1177704472.
  * Pratt & Hanson (1989) Pratt, L. Y. and Hanson, S. J.  Comparing Biases for Minimal Network Construction with Back-Propagation.  In _NIPS_ , 1989.
  * Rezende & Mohamed (2015) Rezende, D. J. and Mohamed, S.  Variational Inference with Normalizing Flows.  In _Proceedings of the 32nd International Conference on Machine Learning_ , 5 2015.  URL http://arxiv.org/abs/1505.05770.
  * Rosenblatt (1956) Rosenblatt, M.  Remarks on Some Nonparametric Estimates of a Density Function.  _The Annals of Mathematical Statistics_ , 27(3):832–837, 9 1956.  ISSN 0003-4851.  doi: 10.1214/aoms/1177728190.  URL http://projecteuclid.org/euclid.aoms/1177728190.
  * Rothfuss et al. (2019) Rothfuss, J., Ferreira, F., Walther, S., and Ulrich, M.  Conditional Density Estimation with Neural Networks: Best Practices and Benchmarks.  Technical report, 2019.  URL http://arxiv.org/abs/1903.00954.
  * Rudemo (1982) Rudemo, M.  Empirical Choice of Histograms and Kernel Density Estimators, 1982.  URL https://www.jstor.org/stable/4615859.
  * Scott & Wand (1991) Scott, D. W. and Wand, M. P.  Feasibility of Multivariate Density Estimates.  _Biometrika_ , 78(1):197, 3 1991.  ISSN 00063444.  doi: 10.2307/2336910.  URL https://www.jstor.org/stable/2336910?origin=crossref.
  * Sheather & Jones (1991) Sheather, S. J. and Jones, M. C.  A Reliable Data-Based Bandwidth Selection Method for Kernel Density Estimation.  _Journal of the Royal Statistical Society_ , 53:683–690, 1991.  doi: 10.2307/2345597.  URL https://www.jstor.org/stable/2345597.
  * Sietsma & Dow (1991) Sietsma, J. and Dow, R. J.  Creating artificial neural networks that generalize.  _Neural Networks_ , 4(1):67–79, 1 1991.  ISSN 0893-6080.  doi: 10.1016/0893-6080(91)90033-2.  URL https://www.sciencedirect.com/science/article/pii/0893608091900332.
  * Silverman (1982) Silverman, B.  On the estimation of a probability density function by the maximum penalized likelihood method.  _The Annals of Statistics_ , 10(3):795–810, 1982\.  URL https://projecteuclid.org/download/pdf_1/euclid.aos/1176345872.
  * Silverman (1986) Silverman, B.  Density estimation for statistics and data analysis.  _Monographs on Statistics and Applied Probability_ , 1986.
  * Sohn et al. (2015) Sohn, K., Lee, H., and Yan, X.  Learning Structured Output Representation using Deep Conditional Generative Models.  In _Advances in Neural Information Processing Systems_ , pp. 3483–3491, 2015.
  * Srivastava et al. (2014) Srivastava, N., Hinton, G., Krizhevsky, A., Sutskever, I., and Salakhutdinov, R.  Dropout: A Simple Way to Prevent Neural Networks from Overfitting.  _Journal of Machine Learning Research_ , 15:1929–1958, 2014\.  URL http://jmlr.org/papers/v15/srivastava14a.html.
  * Sugiyama & Takeuchi (2010) Sugiyama, M. and Takeuchi, I.  Conditional density estimation via Least-Squares Density Ratio Estimation.  In _Proceedings of the Thirteenth International Conference on Artificial Intelligence and Statistics_ , volume 9, pp. 781–788, 2010.  URL http://machinelearning.wustl.edu/mlpapers/paper_files/AISTATS2010_SugiyamaTSKHO10.pdf.
  * Tresp (2001) Tresp, V.  Mixtures of Gaussian Processes.  In _NIPS_ , 2001.  URL https://papers.nips.cc/paper/1900-mixtures-of-gaussian-processes.pdf.
  * Trippe & Turner (2018) Trippe, B. L. and Turner, R. E.  Conditional Density Estimation with Bayesian Normalising Flows.  Technical report, 2 2018.  URL http://arxiv.org/abs/1802.04908.
  * Wan et al. (2013) Wan, L., Zeiler, M., Zhang, S., Cun, Y. L., and Fergus, R.  Regularization of Neural Networks using DropConnect.  In _ICML_ , pp. 1058–1066, 2 2013.  URL http://proceedings.mlr.press/v28/wan13.html.
  * Webb (1994) Webb, A.  Functional approximation by feed-forward networks: a least-squares approach to generalization.  _IEEE Transactions on Neural Networks_ , 5(3):363–371, 5 1994.  ISSN 10459227.  doi: 10.1109/72.286908.  URL http://ieeexplore.ieee.org/document/286908/.
  * Welling & Teh (2011) Welling, M. and Teh, Y. W.  Bayesian Learning via Stochastic Gradient Langevin Dynamics.  In _ICML_ , 2011.  URL https://www.ics.uci.edu/~welling/publications/papers/stoclangevin_v6.pdf.
  * White (1989) White, H.  Learning in Artificial Neural Networks: A Statistical Perspective.  _Neural Computation_ , 1(4):425–464, 12 1989\.  ISSN 0899-7667.  doi: 10.1162/neco.1989.1.4.425.  URL http://www.mitpressjournals.org/doi/10.1162/neco.1989.1.4.425.
  * Yuan et al. (2017) Yuan, X., He, P., Zhu, Q., and Li, X.  Adversarial Examples: Attacks and Defenses for Deep Learning.  Technical report, 12 2017.  URL http://arxiv.org/abs/1712.07107.
  * Zambom & Dias (2013) Zambom, A. Z. and Dias, R.  A Review of Kernel Density Estimation with Applications to Econometrics.  _International Econometric Review (IER)_ , 5(1):20–42, April 2013.

##  Appendix A Derivation Smoothness Regularization

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

$l_{\mathcal{D}}(\mathcal{D})=\sum_{i=1}^{n}l(z_{i})$ (13)

Also, let each $z_{i}$ be perturbed by a random noise vector $\xi\sim K(\xi)$
with zero mean and i.i.d. elements, i.e.

$\mathbb{E}_{\xi\sim
K(\xi)}\left[\xi\right]=0~{}~{}~{}\text{and}~{}~{}~{}\mathbb{E}_{\xi\sim
K(\xi)}\left[\xi_{n}\xi_{j}^{\top}\right]=h^{2}I$ (14)

The resulting loss $l(z_{i}+\xi)$ can be approximated by a second order Taylor
expansion around $z_{i}$

$l(z_{i}+\xi)=l(z_{i})+\xi^{\top}\nabla_{z}l(z)\big{\rvert}_{z_{i}}+\frac{1}{2}\xi^{\top}\nabla_{z}^{2}l(z)\big{\rvert}_{z_{i}}\xi+\mathcal{O}(\xi^{3})$
(15)

Assuming that the noise $\xi$ is small in its magnitude,
$\mathcal{O}(\xi^{3})$ may be neglected. The expected loss under $K(\xi)$
follows directly from (15):

$\mathbb{E}_{\xi\sim
K(\xi)}\left[l(z_{i}+\xi)\right]=l(z_{i})+\mathbb{E}_{\xi\sim
K(\xi)}\left[\xi^{\top}\nabla_{x}l(z)\big{\rvert}_{z_{i}}\right]+\frac{1}{2}\mathbb{E}_{\xi\sim
K(\xi)}\left[\xi^{\top}\nabla_{x}^{2}l(z)\big{\rvert}_{z_{i}}\xi\right]$ (16)

Using the assumption about $\xi$ in (14) we can simplify (16) as follows:

$\displaystyle\mathbb{E}_{\xi\sim K(\xi)}\left[l(z_{i}+\xi)\right]$
$\displaystyle=l(z_{i})+\mathbb{E}_{\xi\sim
K(\xi)}\left[\xi\right]^{\top}\nabla_{z}l(z)\big{|}_{z_{i}}+\frac{1}{2}\mathbb{E}_{\xi\sim
K(\xi)}\left[\xi^{\top}\nabla_{z}^{2}l(z)\big{\rvert}_{z_{i}}\xi\right]$ (17)
$\displaystyle=l(z_{i})+\frac{1}{2}\mathbb{E}_{\xi\sim
K(\xi)}\left[\xi^{\top}\mathbf{H}^{(i)}\xi\right]$ (18)
$\displaystyle=l(z_{i})+\frac{1}{2}\mathbb{E}_{\xi\sim
K(\xi)}\left[\sum_{j}\sum_{k}\xi_{j}\xi_{k}\frac{\partial^{2}l(z)}{\partial
z^{(j)}\partial z^{(k)}}\bigg{|}_{z_{i}}\right]$ (19)
$\displaystyle=l(z_{i})+\frac{1}{2}\sum_{j}\mathbb{E}_{\xi}\left[\xi_{j}^{2}\right]\frac{\partial^{2}l(z)}{\partial
z^{(j)}\partial z^{(j)}}\bigg{|}_{z_{i}}+\frac{1}{2}\sum_{j}\sum_{k\neq
j}\mathbb{E}_{\xi}\left[\xi_{j}\xi_{k}\right]\frac{\partial^{2}l(z)}{\partial
z^{(j)}\partial z^{(k)}}\bigg{|}_{z_{i}}$ (20)
$\displaystyle=l(z_{i})+\frac{\eta^{2}}{2}\sum_{j}\frac{\partial^{2}l(z)}{\partial
z^{(j)}\partial z^{(j)}}\bigg{|}_{z_{i}}$ (21)
$\displaystyle=l(z_{i})+\frac{\eta^{2}}{2}\text{tr}(\mathbf{H}^{(i)})$ (22)

In that, $l(z_{i})$ is the loss without noise and
$\mathbf{H}^{(i)}=\nabla_{z}^{2}l(z)\big{\rvert}_{z_{i}}$ the Hessian of $l$
at $z_{i}$. With $z^{(j)}$ we denote the elements of the column vector $z$.

##  Appendix B Vanilla conditional MLE objective is uniformly consistent

The objective function corresponding to a conditional M-projection.

$l(\theta)=-\mathbb{E}_{p(x,y)}\left[\log\hat{f}_{\theta}(y|x)\right]$ (23)

The sample equivalent:

$\hat{l}_{n}(\theta)=-\frac{1}{n}\sum_{i=1}^{n}\log\hat{f}_{\theta}(y_{i}|x_{i})~{},\quad(x_{i},y_{i})\stackrel{{\scriptstyle
i.i.d}}{{\sim}}P(X,Y)$ (24)

###### Corollary 1

Let $\Theta$ be a compact set and and
$\hat{f}_{\theta}:\mathbb{R}^{l}\times\mathbb{R}^{m}\times\Theta\rightarrow(0,\infty)$
continuous in $\theta$ for all $(x,y)\in\mathbb{R}^{l}\times\mathbb{R}^{m}$
such that
$\mathbb{E}_{p(x,y)}\left[\sup_{\theta\in\Theta}\log\hat{f}_{\theta}(y|x)\right]<\infty$.
Then, as $n\rightarrow\infty$, we have

$\sup_{\theta\in\Theta}\left|\hat{l}_{n}(\theta)-l(\theta)\right|\xrightarrow{a.s.}0$
(25)

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

##  Appendix C Consistency Proofs

###### Lemma 1

Suppose for some $\epsilon>0$ there exists a constant $B^{(\epsilon)}_{p}$
such that

$\int|\log\hat{f}_{\theta}(z)|^{1+\epsilon}p(z)dz\leq
B^{(\epsilon)}_{p}<\infty\quad\forall\theta\in\Theta$ (26)

and there exists an $n_{0}$ such that for all $n>n_{0}$ there exists a
constant $B^{(\epsilon)}_{\hat{q}}$ such that

$\int|\log\hat{f}_{\theta}(z)|^{1+\epsilon}\hat{q}_{n}^{(h_{n})}(z)dz\leq
B^{(\epsilon)}_{\hat{q}}<\infty\quad\forall\theta\in\Theta$ (27)

almost surely. Then, the inequality

$\sup_{\theta\in\Theta}\left|l_{n}^{(h)}(\theta)-l(\theta)\right|\leq
C_{\epsilon}\left(\int|\hat{q}_{n}^{(h)}(z)-p(z)|dz\right)^{\frac{\epsilon}{1+\epsilon}}$
(28)

where $C_{\epsilon}$ 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
$\hat{q}^{(h)}_{n}$, we obtain

$\displaystyle\left|l_{n}^{(h)}(\theta)-l(\theta)\right|$
$\displaystyle=\left|\int\log\hat{f}_{\theta}(z)(\hat{q}_{n}^{(h)}(z)-p(z))dz\right|$
$\displaystyle\leq\int|\log\hat{f}_{\theta}(z)|~{}|\hat{q}_{n}^{(h)}(z)-p(z)|dz$
$\displaystyle=\int|\log\hat{f}_{\theta}(z)|~{}|\hat{q}_{n}^{(h)}(z)-p(z)|^{\frac{1}{1+\epsilon}}~{}|\hat{q}_{n}^{(h)}(z)-p(z)|^{\frac{\epsilon}{1+\epsilon}}dz$
$\displaystyle\leq\left(\int|\log\hat{f}_{\theta}(z)|^{1+\epsilon}~{}|\hat{q}_{n}^{(h)}(z)-p(z)|dz\right)^{\frac{1}{1+\epsilon}}$
$\displaystyle~{}\quad\left(\int|\hat{q}_{n}^{(h)}(z)-p(z)|dz\right)^{\frac{\epsilon}{1+\epsilon}}$
$\displaystyle\leq\bigg{(}\int|\log\hat{f}_{\theta}(z)|^{1+\epsilon}\hat{q}_{n}^{(h)}(z)dz$
$\displaystyle~{}\quad+\int|\log\hat{f}_{\theta}(z)|^{1+\epsilon}p(z)dz\bigg{)}^{\frac{1}{1+\epsilon}}\left(\int|\hat{q}_{n}^{(h)}(z)-p(z)|dz\right)^{\frac{\epsilon}{1+\epsilon}}$

Employing the regularity conditions (26) and (27) and writing
$C^{(\epsilon)}=B^{(\epsilon)}_{p}+B^{(\epsilon)}_{\hat{q}}$, it follows that
$\exists n_{0}$ such that $\forall n>n_{0}$

$\displaystyle\sup_{\theta\in\Theta}\left|l_{n}^{(h)}(\theta)-l(\theta)\right|$
$\displaystyle\leq\left(B^{(\epsilon)}_{p}+B^{(\epsilon)}_{\hat{q}}\right)\left(\int|\hat{q}_{n}^{(h)}(z)-p(z)|dz\right)^{\frac{\epsilon}{1+\epsilon}}$
$\displaystyle=C^{(\epsilon)}\left(\int|\hat{q}_{n}^{(h)}(z)-p(z)|dz\right)^{\frac{\epsilon}{1+\epsilon}}$

with probability 1. $\Box$

Lemma 1 states regularity conditions ensuring that the expectations in
$l^{(h)}_{n}(\theta)$ and $l(\theta)$ are well-behaved in the limit. In
particular, (26) and (27) imply uniform and absolute integrability of the log-
likelihoods under the respective probability measures induced by $p$ and
$\hat{q}_{n}^{(h)}$. Since we are interested in the asymptotic behavior, it is
sufficient for (27) to hold for $n$ large enough with probability 1.

Inequality (28) shows that we can make $|l_{n}^{(h)}(\theta)-l(\theta)|$ small
by reducing the $l_{1}$-distance between the true density $p$ and the kernel
density estimate $\hat{q}_{n}^{(h)}$. There exists already a vast body of
literature, discussing how to properly choose the kernel $K$ and the bandwidth
sequence $(h_{n})$ so that $\int|\hat{q}_{n}^{(h_{n})}(z)-p(z)|dz\rightarrow
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 $\exists n_{0}\forall
n>n_{0}$ inequality (28) holds for some constant $C^{(\epsilon)}$. From our
regularity assumptions it follows that $\mathbb{P}(A^{c})=0$. Given that $A$
holds, we just have to show that
$\int|\hat{q}_{n}^{(h)}(z)-p(z)|dz\xrightarrow{a.s.}0$. Then, the upper bound
in (28) tends to zero and we can conclude our proposition.

For any $\delta>0$ let $B_{n}$ denote the event

$\int|\hat{q}_{n}^{(h)}(z)-p(z)|dz\leq\delta$ (29)

wherein $\hat{q}_{n}^{(h)}(z)$ is a kernel density estimate obtained based on
$n$ samples from $p(z)$. Under the conditions in (9) we can apply Theorem 1 of
Devroye (1983), obtaining an upper bound on the probability that (29) does not
hold, i.e. $\exists u,m_{0}$ such that $\mathbb{P}(B_{n}^{c})\leq e^{-un}$ for
all $n>m_{0}$.

Since we need both $A$ and $B_{n}$ for $n\rightarrow\infty$ to hold, we
consider the intersection of the events $(A\cap B_{n})$. Using a union bound
argument it follows that $\exists k_{0}$ such that $\forall
n>k_{0}:\mathbb{P}((A\cap
B_{n})^{c})\leq\mathbb{P}(A^{c})+\mathbb{P}(B_{n}^{c})=0+e^{-un}=e^{-un}$.
Note that we can simply choose $k_{0}=\max\\{n_{0},m_{0}\\}$ for this to hold.
Hence, $\sum_{n={k_{0}+1}}^{\infty}P((A\cap
B_{n})^{c})<\sum_{n=1}^{\infty}e^{-un}=\frac{1}{e^{u}-1}<\infty$ and by the
Borel-Cantelli lemma we can conclude that

$\lim_{n\rightarrow\infty}\sup_{\theta\in\Theta}\left|l_{n}^{(h)}(\theta)-l(\theta)\right|=0$
(30)

holds with probability 1. $\Box$

Proof of Theorem 1. The inequality in (8) implies that for any $n>n_{0}$,

$\underset{r\rightarrow\infty}{\overline{\lim}}\sup_{\theta\in\Theta}\left|\hat{l}_{n,r}^{(h)}(\theta)-l(\theta)\right|\leq\underset{r\rightarrow\infty}{\overline{\lim}}\sup_{\theta\in\Theta}\left|\hat{l}_{n,r}^{(h)}(\theta)-l_{n}^{(h)}(\theta)\right|+\sup_{\theta\in\Theta}\left|l_{n}^{(h)}(\theta)-l(\theta)\right|$
(31)

Let $n>n_{0}$ be fixed but arbitrary and denote

$J_{n,r}=\sup_{\theta\in\Theta}\left|\hat{l}_{n,r}^{(h)}(\theta)-l_{n}^{(h)}(\theta)\right|\quad
r\in{\mathbb{N}},n>n_{0}$ (32)

It is important to note that $J_{n,r}$ is a random variable that depends on
the samples $\mathbf{Z}^{(n)}=(Z_{1},...,Z_{n})$ as well as on the randomness
inherent in Algorithm 1. We define $\mathbf{I}^{(r)}=(I_{1},...I_{r})$ as the
indices sampled uniformly from $\\{1,...,n\\}$ and
$\Xi^{(r)}=(\xi_{1},...xi_{r})$ as the sequence of perturbation vectors
sampled from $K$. Let $P(\mathbf{Z}^{(n)})$, $P(\mathbf{I}^{(r)})$ and
$P(\Xi^{(r)})$ be probability measures of the respective random sequences.

If we fix $\mathbf{Z}^{(n)}$ to be equal to an arbitrary sequence $Z^{(n)}$,
then $\hat{q}^{(h)}_{n}$ is fixed and we can treat $J_{n,r}$ as the regular
difference between a sample estimate and expectation under
$\hat{q}^{(h)}_{n}$. By the regularity condition $(\ref{eq:reg_cond_2})$, the
compactness of $\Theta$ and the continuity of $f_{\theta}$ in $\theta$, we can
invoke the uniform law of large numbers to show that

$\lim_{r\rightarrow\infty}J_{n,r}=\lim_{r\rightarrow\infty}\sup_{\theta\in\Theta}\left|\hat{l}_{n,r}^{(h)}(\theta)-l_{n}^{(h)}(\theta)\right|=0$
(33)

with probability 1.

Now we want to show that (33) also holds with probability 1 for random
training samples $\mathbf{Z}^{(n)}$. First, we write $J_{n,r}$ as a
deterministic function of random variables:

$J_{n,r}=J(\mathbf{Z}^{(n)},\mathbf{I}^{(r)},\Xi^{(r)})$ (34)

This allows us to restate the result in (33) as follows:

$\displaystyle\begin{split}&{\mathbb{P}}_{\mathbf{I}^{(r)},\Xi^{(r)}}\left(\forall\delta>0~{}\exists
r_{0}~{}\forall
r>r_{0}:J(\mathbf{Z}^{(n)}=Z^{(n)},\mathbf{I}^{(r)},\Xi^{(r)})<\delta\right)\\\
=&\int\int\mathbf{1}\left(\forall\delta>0~{}\exists r_{0}~{}\forall
r>r_{0}:J(\mathbf{Z}^{(n)}=Z^{(n)},\mathbf{I}^{(r)},\Xi^{(r)})<\delta\right)dP(\Xi^{(r)})dP(\mathbf{I}^{(r)})\\\
=&~{}1\end{split}$ (35)

In that $\mathbf{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
(33) holds for random $\mathbf{Z}^{(n)}$:

$\displaystyle\begin{split}&{\mathbb{P}}_{\mathbf{Z}^{(n)},\mathbf{I}^{(r)},\Xi^{(r)}}\left(\forall\delta>0~{}\exists
r_{0}~{}\forall
r>r_{0}:J(\mathbf{Z}^{(n)},\mathbf{I}^{(r)},\Xi^{(r)})<\delta\right)\\\
=&\int\int\int\mathbf{1}\left(\forall\delta>0~{}\exists r_{0}~{}\forall
r>r_{0}:J(\mathbf{Z}^{(n)},\mathbf{I}^{(r)},\Xi^{(r)})<\delta\right)dP(\Xi^{(r)})dP(\mathbf{I}^{(r)})dP(\mathbf{Z}^{(n)})\\\
=&\int
dP(\mathbf{Z}^{(n)})\underbrace{\left(\int\int\mathbf{1}\left(\forall\delta>0~{}\exists
r_{0}~{}\forall
r>r_{0}:J(\mathbf{Z}^{(n)},\mathbf{I}^{(r)},\Xi^{(r)})<\delta\right)dP(\Xi^{(r)})dP(\mathbf{I}^{(r)})\right)}_{=1}\\\
=&~{}1\end{split}$

Note that we can $dP(\mathbf{Z}^{(n)})$ move outside of the inner integrals,
since $\mathbf{Z}^{(n)}$ is independent from $\mathbf{I}^{(r)}$ and
$\Xi^{(r)}$. Hence, we can conclude that (33) also holds, which we denote as
event $A$, with probability 1 for random training data.

From Proposition 1 we know, that

$\lim_{n\rightarrow\infty}\sup_{\theta\in\Theta}\left|l_{n}^{(h)}(\theta)-l(\theta)\right|=0$
(36)

with probability 1. We denote the event that (36) holds as $B$. Since
$P(A^{c})=P(B^{c})=0$, we can use a union bound argument to show that $P(A\cap
B)=1$. From (33) and (31) it follows that for any $n>n_{0}$,

$\underset{r\rightarrow\infty}{\overline{\lim}}\sup_{\theta\in\Theta}\left|\hat{l}_{n,r}^{(h)}(\theta)-l(\theta)\right|\leq\sup_{\theta\in\Theta}\left|l_{n}^{(h)}(\theta)-l(\theta)\right|$
(37)

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

$\lim_{n\rightarrow\infty}\underset{r\rightarrow\infty}{\overline{\lim}}\sup_{\theta\in\Theta}\left|\hat{l}_{n,r}^{(h)}(\theta)-l(\theta)\right|=0$
(38)

almost surely, which concludes the proof. $\Box$

Proof of Theorem 2. The proof follows the argument used in Theorem 1 of White
(1989). In the following, we assume that (11) holds. From Theorem 1 we know
that this is the case with probability 1. Respectively, we only consider
realizations of our training data $\mathbf{Z}^{(n)}$ and noise samples
$\mathbf{I}^{(r)}$, $\Xi^{(r)}$ for which the convergence in (11) holds (see
proof of Theorem 1 for details on this notation).

For such realization, let $(\hat{\theta}_{n,r}^{(h)})$ be minimizers of
$\hat{l}_{n,r}^{(h)}$. Also let $(n_{i})_{i}$ and for any $i$, $(r_{i,j})_{j}$
be increasing sequences of positive integers. Define
$v_{i,j}:=\hat{\theta}_{n_{i},r_{i,j}}^{(h)}$ and
$\mu_{i,j}(\theta):=\hat{l}^{(h)}_{n_{i},r_{i,j}}(\theta)$. Due to the
compactness of $\Theta$ and the Bolzano-Weierstrass property thereof, there
exists a limit point $\theta^{0}\in\Theta$ and increasing subsequences
$(i_{k})_{k},(j_{k})_{k}$ so that $v_{i_{k},j_{k}}\rightarrow\theta^{0}$ as
$k\rightarrow\infty$.

From the triangle inequality, it follows that for any $\epsilon>0$ there
exists $k_{0}$ so that $\forall k>k_{0}$

$|\mu_{i_{k},j_{k}}(v_{i_{k},j_{k}})-l(\theta^{0})|\leq|\mu_{i_{k},j_{k}}(v_{i_{k},j_{k}})-l(v_{i_{k},j_{k}})|+|l(v_{i_{k},j_{k}})-l(\theta^{0})|<2\epsilon$
(39)

given the convergence established in Theorem 1 and the continuity of $l$ in
$\theta$. Next, the result above is extended to

$l(\theta^{0})-l(\theta)=[l(\theta^{0})-\mu_{i_{k},j_{k}}(v_{i_{k},j_{k}})]+[\mu_{i_{k},j_{k}}(v_{i_{k},j_{k}})-\mu_{i_{k},j_{k}}(\theta)]+[\mu_{i_{k},j_{k}}(\theta)-l(\theta)]\leq
3\epsilon$ (40)

which again holds for $k$ large enough. This due to (39),
$\mu_{i_{k},j_{k}}(v_{i_{k},j_{k}})-\mu_{i_{k},j_{k}}(\theta)\leq 0$ since
$v_{i_{k},j_{k}}$ is the minimizer of $\mu_{i_{k},j_{k}}$, and
$\mu_{i_{k},j_{k}}(\theta)-l(\theta)<\epsilon$ by Theorem 1. Because
$\epsilon$ can be made arbitrarily small, $l(\theta^{0})\leq l(\theta)$ as
$k\rightarrow\infty$. Because $\theta\in\Theta$ is arbitrary, $\theta^{0}$
must be in $\Theta^{*}$. 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 $\Theta^{*}$.

In the final step, we proof the theorem by contradiction. Suppose that (12)
does not hold. In this case, there must exist an $\epsilon>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}-\bar{\theta}||_{2}>\epsilon$ for all $k$ and
$\bar{\theta}\in\Theta*$. However, by the previous argument the limit point of
the any sequence $(v_{i_{k},j_{k}})_{k}$ must be in $\Theta^{*}$. That is a
contradiction to
$||(v_{i_{k},j_{k}})_{k}-\bar{\theta}||_{2}>\epsilon~{}\forall~{}k,\bar{\theta}\in\Theta*$.
Since the random sequences $\mathbf{Z}^{(n)}$, $\mathbf{I}^{(r)}$, $\Xi^{(r)}$
where chosen from a set with probability mass of 1, we can conclude our
proposition that

$\lim_{n\rightarrow\infty}\underset{r\rightarrow\infty}{\overline{\lim}}d(\hat{\theta}_{n,r}^{(h)},\Theta^{*})=0$

almost surely. $\Box$

Discussion of Theorem 2. Note that, similar to $\theta^{*}$,
$\hat{\theta}_{n,r}^{(h)}$ does not have to be unique. In case there are
multiple minimizers of $\hat{l}_{n,r}^{(h)}$, 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 $\Theta$, which may not be
attainable in practical settings. However, the application of the theorem to
the context of local optimization is straightforward when $\Theta$ is chosen
as a compact neighborhood of a local minimum $\theta^{*}$ of $l$ (Holmstrom &
Koistinen, 1992b). If we set $\Theta^{*}=\\{\theta^{*}\\}$ and restrict
minimization over $\hat{l}_{n,r}^{(h)}$ to the local region, then
$\hat{\theta}_{n,r}^{(h)}$ converges to $\Theta^{*}$ as $n,r\rightarrow\infty$
in the sense of Theorem 2.

##  Appendix D Conditional density estimation models

###  D.1 Mixture Density Network

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
$\hat{p}(y|x)$ follows as weighted sum of $K$ Gaussians

$\hat{p}(y|x)=\sum_{k=1}^{K}w_{k}(x;\theta)\mathcal{N}\left(y|\mu_{k}(x;\theta),\sigma_{k}^{2}(x;\theta)\right)$
(41)

wherein $w_{k}(x;\theta)$ denote the weight, $\mu_{k}(x;\theta)$ the mean and
$\sigma_{k}^{2}(x;\theta)$ the variance of the k-th Gaussian component. All
the GMM parameters are governed by the neural network with parameters $\theta$
and input $x$.

The mixing weights ${w_{k}(x;\theta)}$ must resemble a categorical
distribution, i.e. it must hold that $\sum_{k=1}^{K}w_{k}(x;\theta)=1$ and
$w_{k}(x;\theta)\geq 0~{}\forall k$. To satisfy the conditions, the softmax
linearity is used for the output neurons corresponding to $w_{k}(x;\theta)$.
Similarly, the standard deviations $\sigma_{k}(x)$ must be positive, which is
ensured by a sofplus non-linearity. Since the component means
$\mu_{k}(x;\theta)$ 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 $\hat{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 $\sigma_{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 $\mu_{k}$ and $M$ the number of
different kernel scales $\sigma_{m}$. The KMN conditional density estimate
reads as follows:

$\hat{p}(y|x)=\sum_{k=1}^{K}\sum_{m=1}^{M}w_{k,m}(x;\theta)\mathcal{N}(y|\mu_{k},\sigma_{m}^{2})$
(42)

As previously, the weights $w_{k,m}$ correspond to a softmax function. The $M$
scale parameters $\sigma_{m}$ are learned jointly with the neural network
parameters $\theta$. The centers $\mu_{k}$ are initially chosen by k-means
clustering on the $\\{y_{i}\\}_{i=1}^{n}$ 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}\circ\dots\circ f_{1}$ into a successively more complex distribution
$p(f(Z_{0}))$. The PDF of samples ${\bm{z}}_{N}\sim p(f(Z_{0}))$ can be
evaluted using the change-of-variable formula:

$\log p({\bm{z}}_{N})=\log
p({\bm{z}}_{0})-\sum_{n=1}^{N}\log\Big{|}\det\frac{\partial
f_{n}}{\partial{\bm{z}}_{n-1}}\Big{|}$ (43)

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
$\hat{f}^{-1}_{i}({\bm{z}}_{i})=f_{i}({\bm{z}}_{i})$.

We experimented with three types of flows: planar flows, radial flows as
parametrized by Trippe & Turner (2018) and affine flows
$f^{-1}({\bm{z}})=\exp(a){\bm{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.

##  Appendix E Simulated densities and datasets

###  E.1 SkewNormal

The data generating process $(x,y)\sim p(x,y)$ resembles a bivariate joint-
distribution, wherein $x\in\mathbb{R}$ follows a normal distribution and
$y\in\mathbb{R}$ a conditional skew-normal distribution (Anděl et al., 1984).
The parameters $(\xi,\omega,\alpha)$ of the skew normal distribution are
functionally dependent on $x$. Specifically, the functional dependencies are
the following:

$\displaystyle x$
$\displaystyle\sim\mathcal{N}\left(~{}\cdot~{}\bigg{|}\mu=0,\sigma=\frac{1}{2}\right)$
(44) $\displaystyle\xi(x)$ $\displaystyle=a*x+b\qquad a,b\in\mathbb{R}$ (45)
$\displaystyle\omega(x)$ $\displaystyle=c*x^{2}+d\qquad c,d\in\mathbb{R}$ (46)
$\displaystyle\alpha(x)$
$\displaystyle=\alpha_{low}+\frac{1}{1+e^{-x}}*(\alpha_{high}-\alpha_{low})$
(47) $\displaystyle y$ $\displaystyle\sim
SkewNormal\big{(}\xi(x),\omega(x),\alpha(x)\big{)}$ (48)

Accordingly, the conditional probability density $p(y|x)$ corresponds to the
skew normal density function:

$p(y|x)=\frac{2}{\omega(x)}\mathcal{N}\left(\frac{y-\xi(x)}{\omega(x)}\right)\Phi\left(\alpha(x)\frac{y-\xi(x)}{\omega(x)}\right)$
(49)

In that, $\mathcal{N}(\cdot)$ denotes the density, and $\Phi(\cdot)$ the
cumulative distribution function of the standard normal distribution. The
shape parameter $\alpha(x)$ controls the skewness and kurtosis of the
distribution. We set $\alpha_{low}=-4$ and $\alpha_{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 4(a) illustrates the conditional skew normal distribution.

(a) SkewNormal

(b) GaussianMixture

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

###  E.2 Gaussian Mixture

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

$p(x,y)=\sum_{i=1}^{K}w_{k}~{}\mathcal{N}(y|\mu_{y,k},\Sigma_{y,k})\mathcal{N}(x|\mu_{x,k},\Sigma_{x,k})$
(50)

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

$\displaystyle
p(y|x)=\sum_{i=1}^{K}W_{k}(x)~{}\mathcal{N}(y|\mu_{y,k},\Sigma_{y,k})$ (51)

wherein the mixture weights are a function of $x$:

$W_{k}(x)=\frac{w_{k}~{}\mathcal{N}(x|\mu_{x,k},\Sigma_{x,k})}{\sum_{j=1}^{K}w_{k}~{}\mathcal{N}(x|\mu_{x,j},\Sigma_{x,j})}$
(52)

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 $\sigma=1.5$. The covariance matrices $\Sigma_{y,k})$ and
$\Sigma_{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 4(b) depicts a
2-dimensional equivalent, generated with the procedure explained above.

###  E.3 Euro Stoxx 50 data $(d_{x}=14,d_{y}=1)$

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\in\mathcal{Y}\subseteq\mathbb{R}$,
whereas the conditional variable $x$ constitutes a 14-dimensional vector, i.e.
$x\in\mathcal{X}\subseteq\mathbb{R}^{14}$.

###  E.4 NYC Taxi data $(d_{x}=6,d_{y}=2)$

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 drop-off 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.

###  E.5 UCI

##### Boston Housing $(d_{x}=13,d_{y}=1)$

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 $(d_{x}=8,d_{y}=1)$

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-learning-
databases/concrete/compressive/

##### Energy $(d_{x}=9,d_{y}=1)$

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