Noise Regularization for Conditional Density Estimation

TL;DR

This paper introduces a noise regularization technique for neural network-based conditional density estimation, improving model generalization and outperforming existing methods across multiple datasets.

Contribution

A novel, model-agnostic noise regularization method for CDE that enhances generalization and is proven to be asymptotically consistent.

Findings

Noise regularization outperforms other regularization methods.

Method is effective across seven datasets and three CDE models.

Neural CDE becomes preferable over non- and semi-parametric approaches.

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.

Tables

Table 1. 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	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 2. 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.

	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

Equations

$$\text{arg}\underset{\theta \in \Theta }{\max }\sum _{i=1}^n\log {\hat{f}}_{\theta }(x_i)=\text{arg}\underset{\theta \in \Theta }{\min }{\mathcal{D}}_{KL}(p_{\mathcal{D}}||{\hat{f}}_{\theta })\text{vspace}-4pt$$

$${\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]$$

$$\text{vspace}-4pt{\theta }^{\ast }=\text{arg}\underset{\theta }{\min }-\sum _{i=1}^n\log {\hat{f}}_{\theta }(y_i|x_i)$$

$${\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\text{vspace}-2pt$$

$${\mathbb{E}}_{K(\xi )}\left[l(z_i+\xi )\right]=$$

$$\approx$$

$$\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}$$

$$l(\theta )=-{\mathbb{E}}_{p(z)}\left[\log {\hat{f}}_{\theta }(z)\right]$$

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

$${\hat{l}}_{n,r}^{(h)}(\theta )=-\frac{1}{r}\sum _{j=1}^r\log {\hat{f}}_{\theta }({\tilde{z}}_j),{\tilde{z}}_j\sim {\hat{q}}_n^{(h)}(z)\text{vspace}-2pt$$

$$\begin{array}{rl}\underset{\theta \in \Theta }{\sup }\left|{\hat{l}}_{n,r}^{(h)}(\theta )-l(\theta )\right|\le & \underset{\theta \in \Theta }{\sup }\left|{\hat{l}}_{n,r}^{(h)}(\theta )-l_n^{(h)}(\theta )\right|\\ & +\underset{\theta \in \Theta }{\sup }\left|l_n^{(h)}(\theta )-l(\theta )\right|\end{array}\text{vspace}-2pt$$

$$\underset{n\to \infty }{\lim }{h}_n=0,\underset{n\to \infty }{\lim }n(h_n{)}^d=\infty$$

$$\underset{n\to \infty }{\lim }\underset{\theta \in \Theta }{\sup }\left|l_n^{(h)}(\theta )-l(\theta )\right|=0$$

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

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

$$l_{\mathcal{D}}(\mathcal{D})=\sum _{i=1}^{n}l(z_i)$$

$${\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$$

$$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})$$

$$\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]$$

$${\mathbb{E}}_{\xi \sim K(\xi )}\left[l(z_i+\xi )\right]$$

$$=l(z_i)+\frac{1}{2}{\mathbb{E}}_{\xi \sim K(\xi )}\left[{\xi }^{\top }{\mathbf{H}}^{(i)}\xi \right]$$

$$\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]$$

$$\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}}$$

$$\displaystyle=l(z_{i})+\frac{\eta^{2}}{2}\sum_{j}\frac{\partial^{2}l(z)}{\partial z^{(j)}\partial z^{(j)}}\bigg{|}_{z_{i}}$$

$$=l(z_i)+\frac{{\eta }^2}{2}\text{tr}({\mathbf{H}}^{(i)})$$

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

$${\hat{l}}_n(\theta )=-\frac{1}{n}\sum _{i=1}^n\log {\hat{f}}_{\theta }(y_i|x_i),(x_i,y_i)\stackrel{i.i.d}{\sim }P(X,Y)$$

$$\underset{\theta \in \Theta }{\sup }\left|{\hat{l}}_n(\theta )-l(\theta )\right|\stackrel{a.s.}{\to }0$$

$$\int |\log {\hat{f}}_{\theta }(z){|}^{1+ϵ}p(z)dz\le B_p^{(ϵ)}<\infty \forall \theta \in \Theta$$

$$\int |\log {\hat{f}}_{\theta }(z){|}^{1+ϵ}{\hat{q}}_n^{(h_n)}(z)dz\le B_{\hat{q}}^{(ϵ)}<\infty \forall \theta \in \Theta$$

$$\underset{\theta \in \Theta }{\sup }\left|l_n^{(h)}(\theta )-l(\theta )\right|\le C_{ϵ}{\left(\int |{\hat{q}}_n^{(h)}(z)-p(z)|dz\right)}^{\frac{ϵ}{1+ϵ}}$$

$$\left|l_n^{(h)}(\theta )-l(\theta )\right|$$

$$\le \int |\log {\hat{f}}_{\theta }(z)||{\hat{q}}_n^{(h)}(z)-p(z)|dz$$

$$=\int |\log {\hat{f}}_{\theta }(z)||{\hat{q}}_n^{(h)}(z)-p(z){|}^{\frac{1}{1+ϵ}}|{\hat{q}}_n^{(h)}(z)-p(z){|}^{\frac{ϵ}{1+ϵ}}dz$$

Full text

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}$:

[TABLE]

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)$:

[TABLE]

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

[TABLE]

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.

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

[TABLE]

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

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

[TABLE]

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

[TABLE]

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

[TABLE]

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

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:

[TABLE]

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:

[TABLE]

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:

[TABLE]

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

[TABLE]

Then,

[TABLE]

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,

[TABLE]

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

[TABLE]

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.

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.

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.

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}$:

[TABLE]

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.

[TABLE]

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

[TABLE]

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):

[TABLE]

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

[TABLE]

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.

[TABLE]

The sample equivalent:

[TABLE]

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

[TABLE]

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

[TABLE]

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

[TABLE]

almost surely. Then, the inequality

[TABLE]

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

[TABLE]

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

[TABLE]

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

[TABLE]

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

[TABLE]

holds with probability 1. $\Box$

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

[TABLE]

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

[TABLE]

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

[TABLE]

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:

[TABLE]

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

[TABLE]

In that $\mathbf{1}(A)$ denotes an indicator function which returns $1$ if $A$ is true and [math] else. Next we consider the probability that the convergence in (33) holds for random $\mathbf{Z}^{(n)}$:

[TABLE]

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

[TABLE]

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}$,

[TABLE]

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

[TABLE]

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

[TABLE]

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

[TABLE]

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

[TABLE]

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

[TABLE]

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:

[TABLE]

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:

[TABLE]

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:

[TABLE]

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

[TABLE]

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.

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.

[TABLE]

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

[TABLE]

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

[TABLE]

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