{'fraction_non_alphanumeric': 0.04068381688172927, 'fraction_numerical': 0.010136708583326484, 'mean_word_length': 4.2932134570765665, 'pattern_counts': {'":': 0, '<': 0, '<?xml version=': 0, '>': 0, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 4, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Gaussian Process Conditional Density Estimation

Conditional Density Estimation (CDE) models deal with estimating conditional distributions. The conditions imposed on the distribution are the inputs of the model. CDE is a challenging task as there is a fundamental trade-off between model complexity, representational capacity and overfitting. In this work, we propose to extend the model's input with latent variables and use Gaussian processes (GP) to map this augmented input onto samples from the conditional distribution. Our Bayesian approach allows for the modeling of small datasets, but we also provide the machinery for it to be applied to big data using stochastic variational inference. Our approach can be used to model densities even in sparse data regions, and allows for sharing learned structure between conditions. We illustrate the effectiveness and wide-reaching applicability of our model on a variety of real-world problems, such as spatio-temporal density estimation of taxi drop-offs, non-Gaussian noise modeling, and few-shot learning on omniglot images.


Introduction
Conditional Density Estimation (CDE) is the very general task of inferring the probability distribution p(f (x) | x), where f (x) is a random variable for each x. Regression can be considered a CDE problem, although the emphasis is on modeling the mapping rather than the conditional density. The conditional density is commonly Gaussian with parameters that depend on x. This simple model for data may be inappropriate if the conditional density is multi-modal or has non-linear associations.
Throughout this paper we consider an input x to be the condition, and the output y to be a sample from the conditional density imposed by x. For example, in the case of estimating the density of taxi drop-offs, the input or condition x could be the pick-up location and the output y would be the corresponding drop-off. In this context, we are more interested in learning the complete density over drop-offs rather than only a single point estimate, as we would expect the taxi drop-off to be multi-modal because passengers need to go to different places (e.g., airport/city center/suburbs). We would also expect the drop-off location to depend on the starting point and time of day: therefore, we are interested in conditional densities. In the experiment section we will return to this example.
In this work, we present a Gaussian process (GP) based model for estimating conditional densities, abbreviated as GP-CDE. While a vanilla GP used directly is unlikely to be a good model for conditional density estimation as the marginals are Gaussian, we extend the inputs to the model with latent variables to allow for modeling richer, non-Gaussian densities when marginalizing the latent variable. Fig. 1 shows a high-level overview of the model. The added latent variables are denoted by w. The latent variable w and condition x are used as input of the GP. A recognition/encoder network is used to amortize the learning of the variational posterior of the latent variables. The matrices A and P act as probabilistic linear transforms on the input and the output of the GP, respectively.
The GP-CDE model is closely related to both supervised and unsupervised, non-Bayesian and Bayesian models. We first consider the relationship to parametric models, in particular Variational A x x y NN GP w Decoder Encoder f P y Figure 1: Diagram of the GP-Conditional Density Estimator. The GP-CDE consists of an encoder (blue) and a decoder (orange). The observed variables x and y are respectively the condition on the distribution and a sample from the conditional distribution. The encoder consists of a neural network, and uses the variables x and y as inputs to produce the parameters for the posterior of the latent variable w. The decoder part is built out of a GP and two linear transformation matrices, A and P.
A is applied to the condition x to reduce the dimensionality of x before it is combined with the latent w and fed into the GP. The matrix P can be used to correlate the outputs of the latent function f .
Autoencoders (VAEs) [17,23] and their conditional counterparts (CVAE) [18,24]. (C)VAEs use a deterministic, parametrized neural network as decoder, whereas we adopt a Bayesian non-parametric GP. Using a GP for this non-linear decoder mapping offers two advantages. First, it allows us to specify our prior beliefs about the mapping, resulting in a model that can gracefully accommodate sparse or missing data. Note that even large datasets can be sparse, e.g., the omniglot images or taxi drop-off locations at the fringes of a city. As the CVAE has no prior on the decoder mapping, the model can overfit the training data and the latent variable posterior becomes over-concentrated, leading to poor test log-likelihoods. Second, the GP allows tractable uncertainty propagation through the decoder mapping (for certain kernels). This allows us to calculate the variational objective deterministically, and we can further exploit natural gradients for fast inference. Neural network decoders do not admit such structure, and are typically optimized with general-purpose tools.
A second perspective can be seen through the connection to Gaussian process models. By dropping the latent variables w we recover standard multiple-output GP regression with sparse variational inference. If we drop the known inputs x and use only the latent variables, we obtain the Bayesian GP-LVM [26,20]. Bayesian GP-LVMs are typically used for modeling complex distributions and non-linear mappings from a lower-dimensional latent variable into a high-dimensional space. By combining the GP-LVM framework with known inputs we create a model that outputs conditional samples in this high-dimensional space.
Our primary contribution is to show that our GP-CDE can be applied to a wide variety of settings, without the necessity for fine-tuning or regularization. We show that GP-CDE outperforms GP regression on regression benchmarks; we study the importance of accurate density estimation in high-dimensional spaces; and we deal with learning the correlations between conditions in a large spatio-temporal dataset. We achieve this through three specific contributions. (i) We extend the model of Wang and Neal [29] with linear transformations for the inputs and outputs. This allows us to deal with high-dimensional conditions and enables a priori correlations in the output. (ii) We apply natural gradients to address the difficulty of mini-batched optimization in [26]. (iii) We derive a free-form optimal posterior distribution over the latent variables. This provides a tighter bound and reduces the number of variational parameters to optimize.

Background: Gaussian Processes and Latent Variable Models
Gaussian Processes A Gaussian process (GP) is a Bayesian non-parametric model for functions. Bayesian models have two significant advantages that we exploit in this work: we can specify prior beliefs leading to greater data efficiency, and we can obtain uncertainty estimates for predictions. A GP is defined as a set of random variables {f (x 1 ), f (x 2 ), . . .}, any finite subset of which follow a multivariate Gaussian distribution [22]. When a stochastic function f : R D → R follows a GP it is fully specified by its mean m(·) and covariance function k(·, ·), and we write f ∼ GP(m(·), k(·, ·)).
The most common use of a GP is the regression task of inferring an unknown function f , given a set of N observations y = [y 1 , . . . , y N ] and corresponding inputs x 1 , . . . , x n . The likelihood p(y i | f ) generally is taken to depend on f (x i ) only, and the Gaussian likelihood N (y n | f (x n ), σ 2 ) is widely used as it results in analytical closed-form inference.

Conditional Deep Latent Variable Models Conditional Deep Latent
Variable Models (C-DLVMs) consist of two components: a prior distribution p(w n ) over the latent variables 2 which is assumed to factorize over the data, and a generator or decoder function g θ (x n , w n ) : R Dx+Dw → R Dy . The VAE and CVAE are examples where the generator function is a deep (convolutional) neural network with weights θ [17,24]. The outputs of the generator function are the parameters of a likelihood, commonly the Gaussian for continuous data or the Bernoulli for binary data. The joint distribution for a single data point is p θ (y n , w n | x n ) = p(w n )p θ (y n | x n , w n ). We assume the data to be i.i.d., then the the marginal likelihood in the Gaussian case is given by where X = {x} N n=1 , and likewise for Y and W. As g θ (x n , w n ) is a complicated non-linear function of its inputs, this integral cannot be calculated in closed-form. Kingma and Welling [17] and Rezende et al. [23] addressed this problem by using variational inference. Variational inference posits an approximate posterior distribution q φ (W), and finds the closest q φ to the true posterior, measured by KL divergence, i.e. arg min q φ KL [q φ (W) p(W | X, Y)]. It can be shown that this optimization objective is equal to the Evidence Lower Bound (ELBO) A mean-field distribution is typically used for the latent variables W with a multivariate Gaussian form for q φ (w n ) = N (w n | µ wn , Σ wn ). Rather than representing the Gaussian parameters µ wn and Σ wn for each data point directly, Kingma and Welling [17] and Rezende et al. [23] instead amortize these parameters into a set of global parameters φ, where φ parameterizes an auxiliary function h φ : (x n , y n ) → (µ wn , Σ wn ), referred to as the encoder/recognition network.

Conditional Density Estimation with Gaussian Processes
This section details our model and the inference scheme. Our contributions are threefold:(i) we derive an optimal free-form variational distribution q(W) (Section 3.2); (ii) we ease the burden of jointly optimizing q(f (·)) and q(W) using natural gradients (Section 3.2.1) for the variational parameters of q(f (·)); (iii) we extend the model to allow for the modeling of high-dimensional inputs and impose correlation on the outputs using linear transformations (Section 3.3).

Model
The key idea of our model is to substitute the neural network decoder in the C-DLVM framework with a GP, see Fig. 1. Treating the decoder in a Bayesian manner leads to several advantages. In particular, in the small-data regime, a probabilistic decoder will be advantageous to leverage prior assumptions and avoid over-fitting.
As we want to apply our model on both high-dimensional correlated outputs (e.g., images) and highdimensional inputs (e.g., one-hot encodings of omniglot labels), we introduce two matrices A and P. They are used for probabilistic linear transformations of the inputs and the outputs, respectively. The likelihood of the GP-CDE model is then given by where [·, ·] denotes concatenation. We assume the GP f (·) consists of L independent GPs f (·) for each output dimension . The latent variables are a priori independent for each data point and have a standard-normal prior distribution. We discuss priors for A and P in section 3.3.

Inference
In this section, we present our inference scheme, initially in the case without A and P to lighten the notation. We will return to these matrices in section section 3.3. We calculate an ELBO on the marginal likelihood, similarly to (1). Assuming a factorized posterior q(f (·), W) = q(f (·)) n q(w n ), where W = {w n } N n=1 , between the GP and the latent variables we get the ELBO Since the ELBO is a sum over the data we can calculate unbiased estimates of the bound using mini-batches. We follow Hensman et al. [15] and choose independent sparse GPs over the output dimensions q(f (·)) = p(f (·) | u ) q(u ) du , where = 1, . . . , L and u ∈ R M are inducing outputs corresponding to the inducing inputs z m , m = 1, . . . , M , so that , the integral can be calculated in closed-form. The result is a new sparse GP for each of the output dimensions q(f (·)) = GP(µ (·), σ (·, ·)) with closed-form mean and variance. See Appendix C for a detailed derivation.
Using the results of Matthews et al. [21], the KL-term over the multi-dimensional latent function f (·) simplifies to KL [q(u ) p(u )], which is closed-form, since q(u ) and p(u ) are both Gaussian.
The inner expectation over the variational posterior q(f (·)) in (2) can be calculated in closed-form (see Appendix C) as the likelihood is Gaussian. We define this analytically tractable quantity as Using this definition, and using the sparse variational posterior for q(f (·)) as described above, we write the bound in (2) as We consider two options for q(w n ): (i) we can either make a further Gaussian assumption, and have a variational posterior of the form q(w n ) = N (w n | µ wn , Σ wn ), or (ii) we can find the analytically optimal value of the bound for a free-form q(w n ).
(i) Gaussian q(w n ) First, a Gaussian q(w n ) implies that the KL over the latent variables is closedform, as both the prior and posterior are Gaussian. Therefore, we are left with the calculation of the first term in (4) E q(wn) L wn . We follow the approach of C-DLVMs, explained in Section 2, and use Monte Carlo sampling to estimate the expectation. To enable differentiability, we use the re-parameterization trick and write w n = µ wn + L wn ξ n with p(ξ n ) = N (0, I), independent for each data point, and L wn L wn = Σ wn . Note that this does not change the distribution of q(w n ), but now the expectation is over a parameterless distribution. We can then take a differentiable unbiased estimate of the bound by sampling from ξ n .
In practice, rather than represent the Gaussian parameters µ wn and L wn for each data-point directly, we instead amortize these parameters into a set of global parameters φ, where φ parameterizes an auxiliary function h φ , (or 'recognition network') of the data: (µ wn , L wn ) = h φ (x n , y n ). This is identical to the decoder component of C-DLVMs.
An alternative approach would be to use the kernel expectation results of Girard et al. [13] to evaluate E q(wn) L wn . Using these results we can evaluate the bound in closed-form, rather than an approximate using Monte Carlo. However, the computations involved in calculating the kernel expectations can be prohibitive, as it requires evaluating a N M 2 D y sized tensor. Furthermore, closed-form solutions for the kernel expectations only exist for RBF and polynomial kernels, which makes this approach less favorable in practice.
(ii) Analytically optimal q(w n ) So far, we assumed that the variational distribution q(w n ) is Gaussian. When q(·) is non-Gaussian, it is possible to integrate over w n with quadrature as we detail in the following. We first bound the conditional p(Y | X, W) and use the same sparse variational posterior for the GP as before, to obtain As shown in [15] and explained above we can calculate L wn analytically. By expressing the marginal likelihood as log p( where we exploited the monotonicity of the logarithm. We can compute this integral with quadrature when w n is low-dimensional (the dimensionality of x n does not matter here). Assuming we have sufficient quadrature points, this gives the analytically optimal bound for q(w n ). The analytical optimal approach does not resort to the Gaussian approximation for q(w n ), so it is a tighter bound. See Appendix D for a proof that this bound is necessarily tighter than the bound in the Gaussian case.

Natural Gradient
Optimizing q(u ) together with q(W) can be challenging due to problems of local optima and the strong coupling between the inducing outputs u and the latent variables W. One option is to analytically optimize the bound with respect to the variational parameters of u , but this prohibits the use of mini-batches and reduces the applicability to large-scale problems. Recall that the variational parameters of u are the mean and the covariance of the approximate posterior distribution q(u ) = N (m , S ) over the inducing outputs. We can use the natural gradient [3] to update the variational parameters m and S .
This approach has the attractive property of recovering exactly the analytically optimal solution for q(u ) in the full-batch case if the natural gradient step size is taken to be 1. While natural gradients have been used before in GP models [15], they have not been used in combination with uncertain inputs. Due to the quadratic form of the log-likelihood as a function of the kernel inputs X, W and Z, we can calculate the expectation w.r.t. q(W), which will still be quadratic in the inducing outputs u . Therefore, the expression is still conjugate, and the natural gradient step of size 1 recovers the analytic solution.
In practice, the natural gradient is used for the Gaussian variational parameters m and S , and ordinary gradients are used for the inducing inputs Z, the recognition network (if applicable) and other hyperparameters of the model (the kernel and likelihood parameters). The variational parameters of q(A) and the parameters of P are also updated using the ordinary gradient.

Probabilistic Linear Transformations
Input For high-dimensional inputs it may not be appropriate to define a GP directly in the augmented space [x n , w n ] ∈ R Dx+Dw . This might be the case if the input data is a one-hot encoding of many classes. We can extend our model with a linear projection to a lower-dimensional space before concatenating with the latent variables, [Ax n , w n ]. We denote this projection matrix by A, as shown in Fig. 1.
We use an isotropic Gaussian prior for elements of A and a Gaussian variational posterior that factorizes between A and the other variables in the model: term and include the mean and variance for Ax as the input of the GP. A similar approach was used in the regression case by Titsias and Lázaro-Gredilla [27].
Output We can move beyond the assumption of a priori independent outputs to a correlated model by using a linear transformation of the outputs of the GP. This model is equivalent to a 'multi-output' GP model with a linear model of covariance between tasks. In the multi-output GP framework [2], the D y outputs are stacked to a single vector of length N D y , and a single GP is used jointly with a structured covariance. In the simplest case, the covariance can be structured as R ⊗ K, where R can be any positive semi-definite matrix of size D y × D y , and K is an N × N matrix. By transforming the outputs with the matrix P we recover exactly this model with R = P P. Apart from the simplicity of implementation, another advantage is that we can handle degenerate cases (i.e., where the number of outputs is less than D y ) without having to deal with issues of ill-conditioning. It would be possible to use a Gaussian prior for P while retaining conjugacy, but in our experiments we use a non-probabilistic P and optimize it using MAP.

Related work
The GP-CDE model is closely related to both supervised an unsupervised Gaussian process based models. If we drop the latent variables W our approach recovers standard multiple-output Gaussian process regression with sparse variational inference. If we drop the known inputs X and use only the latent variables, we obtain a Bayesian GP-LVM [26]. Bayesian GP-LVMs are typically used for modeling complex distributions and non-linear mappings from a lower-dimensional variable into a high-dimensional space. By combining the GP-LVM framework with known inputs we create a model that outputs conditional samples in this high-dimensional space. Differently from Titsias and Lawrence 26, during inference we do not marginalize out the inducing variables u but rather treat them as variational parameters of the model. This scales our model to arbitrarily large datasets through the use of Stochastic Variational Inference (SVI) [15,16]. While the mini-batch extension to the Bayesian GP-LVM was suggested in [15], its notable absence from the literature may be due to the difficulty in the joint optimization of W and u . We found that the natural gradients were essential to alleviate this problem. A comparison demonstrating this is presented in the experiments.
Wang and Neal [29] proposed the Gaussian Process Latent Variable model (GP-LV), which is a special case of our model. The inference they employ is based on Metroplis sampling schemes and does not scale to large datasets or high dimensions. In this work, we extend their model using linear projection matrices on both input and output, and we present an alternative method of inference that scales to large datasets. Damianou and Lawrence [11] also propose a special case of our model, though they use it for missing data imputation rather than to induce non-Gaussian densities. They also use sparse variational inference, but they analytically optimize q(u ) so cannot use mini-batches. Depeweg et al. [12] propose a similar model, but use a Bayesian neural network instead of a GP.
The use of a recognition model, as in VAEs, was first proposed by Lawrence and Quiñonero-Candela [20] in the context of a GP-LVM, though it was motivated as a constraint on the latent variable locations rather than an amortization of the optimization cost. Recognitions models were later used by Bui and Turner [7] and by Dai et al. [9] for deep GPs.
A GP model with latent variables and correlated multiple outputs was recently proposed in Dai et al. [10]. In this model, the latent variables determine the correlations between outputs via a Kronecker-structured covariance, whereas we have a fixed between-output covariance. That is, in our model the covariance of the stacked outputs is (PP ) ⊗ (K X K W ), whereas in Dai et al. [10] the covariance is K W ⊗ K X . These models are complementary and perform different functions. [6] proposed a model that is also similar to ours, but with categorical variables in the latent space. Other approaches to non-parametric density estimation include modeling the log density directly with a GP [1], and using an infinite generalization of the exponential family [25] which was recently extended to the conditional case [4].

Experiments
Spatio-temporal density estimation We apply our model to a New York City taxi dataset to perform conditional spatial density estimation. The dataset holds records of more than 1.4 million taxi trips, which we filter to include trips that start and end within the Manhattan area. Our 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. The two temporal features are encoded as sine and cosine with the natural periods, giving 6-dimensional inputs in total 3 . Trippe and Turner [28] use a similar setup but predict a distribution over the taxi pick-up locations given the fare of the ride and the percent tip.  [26]. We evaluate the models using negative log predictive probability (NLPP) of the test set. The test sets are constructed by sequentially adding points that have greatest minimum distance from the testing set. In this way we cover as much of the input space as possible. We use a test set of 1000 points, and vary the number of training points to establish the utility of models in both sparse and dense data regimes. We use 1K, 5K and 1M randomly selected training points to evaluate the models in both sparse and dense data regimes.  We also evaluate several MDNs models with differing number of mixing components (MDN-k, where k is the number of components), using fully connected neural networks with 3 layers. The MDN model perform poorly except in the large data regime, where the model with the largest number of components is the best performing. The MDN with a large number of components can put mass at localized locations, which for this data is likely to be appropriate as the taxis are confined to streets.
We test three GP-based models: our GP-CDE model with 2-dimensional latent variables, and two special cases: one without conditioning (GP-LVM) and one without latent variables (GP). The GP-LVM [26] is our model without the conditioning, and does not perform well on this task as it has not access to the inputs and models all conditions identically. The GP model has no latent variables and independent Gaussian marginals, and so cannot model this data well as the drop-off location is quite strongly non-Gaussian. We added predictive probabilities for all models in Appendix F to illustrate these findings.
The GP-CDE performs best on this dataset for the small data regimes. For the large data case the MDN model is superior. We attribute this to the high density of data when 1 million training points are used. We used a 2D latent space and Gaussian q(W) for the latent variables, with a recognition network amortizing the inference. We use the RBF kernel and use Monte Carlo sampling to evaluate the bound, as described in Section 3.2. For training we use the Adam optimizer with a exponentially decaying learning rate starting at 0.01 for the hyperparameters, the inducing inputs and on the recognition network parameters. Natural gradient steps of size 0.05 are used for the GP's variational parameters. Fig. 2 shows the density of our GP-CDE model for two different conditions. Similar figures for the other methods are in Appendix F. Few-shot learning We demonstrate the GP-CDE model for the challenging task of few-shot conditional density estimation on the omniglot dataset. Our task is to obtain a density over the pixels given the class label. We use the training/test split from [19], using all the examples in the training classes and four samples from each of the test classes. The inputs are one-hot encoded (1623 classes) and the outputs are the pixel intensities, which we resize to 28 × 28.
We apply a linear transformation on both the input and output (see Section 3.3). We use a 1623 × 30 matrix A with independent standardnormal priors on each component to project the labels onto a 30-dimensional space. To prevent the model from overfitting, it is important to treat the A transformation in a Bayesian way and marginalize it.
To correlate the outputs a priori we use a linear transformation P of the GP outputs, which is equivalent to considering multiple outputs jointly in a single GP with a Kronecker-structured covariance. See Section 3.3. We use 400 GP outputs so P has shape 400 × 784. To initialize P we use a model of local correlation by using the Matérn-5 /2 kernel with unit lengthscale on the pixel coordinates and taking the first 400 eigenvectors scaled by the square-root eigenvalues. We then optimize the matrix P as a hyperparameter of the model. Learning the P matrix is a form of transfer learning: we update our prior in light of the training classes to make better inference about the few-shot test classes.
We obtain a log-likelihood of 7.2 × 10 −2 nat/pixel, averaging over the all the test images (659 classes with 16 images per class). We train for 25K iterations with the same training procedure as in the previous experiment. Samples from the posterior on a selection of test classes are shown in Fig. 3. For a larger selection, see Fig. 8 in Appendix F.
Heteroscedastic noise modeling We use 10 UCI regression datasets to compare two variants of our CDE model with a standard sparse GP and a CVAE. Since we model a 1D target we consider w n to be uni-dimensional, allowing us to use the quadrature method (Section 3.2) to obtain the bound for an analytically optimal q(w n ). We compare also to an amortized Gaussian approximation to q(w n ), where we use a three-layer fully connected neural network with tanh activations for the recognition model. In all three models we use a RBF kernel and 100 inducing points, optimizing for 20K iterations using Adam optimizer for the hyperparameters and a natural gradient optimizer with step size 0.1 for the Gaussian variational parameters. The quadrature model use Gauss-Hermite quadrature with 100 points. For the CVAE we use, given the modest size of the UCI datasets, a relatively small encoder and decoder network architecture together with dropout. See Appendix B for details. Fig. 3 shows the test log-likelihoods using 20-fold cross validation with 10% test splits. We normalize the inputs to have zero mean and unit variance. We see that the quadrature CDE model outperforms the standard GP and CVAE on many of the datasets. The optimal GP-CDE model performs better than the GP-CDE with Gaussian q(w) on all datasets. This can be attributed to three reasons: we impose fewer restrictions on the variational posterior, there is no amortization gap (i.e. the recognition network might not find the optimal parameters [8]), and problems of local optima are likely to be less severe as there are fewer variational parameters to optimize. Density estimation of image data In this experiment we compare the test log-likelihood of the GP-CDE and the CVAE [24] on the MNIST dataset for the image generation task. We train the models with N = 2, 4, 8, . . . , 512 images per class, to test their utility in different data regimes. The model's input is a one-hot encoding of the image label which we concatenate with a 2-dimensional latent variable. We use all 10, 000 test images to calculate the average test log-likelihood, which we estimate using Monte Carlo.
For the CVAE's encoder and decoder network architecture we follow Wu et al. [30] and regularize the network using dropout. Appendix B contains more details on the CVAE's setup. The GP-CDE has the same setup as in the few-shot learning experiment, except that we set the shape of the output mixing matrix P to 50 × 784. We reduce the size of P, compared to the omniglot experiment, as the MNIST digits are relatively easier to model. Since we are considering small datasets in this experiment the role of the mixing matrix becomes more important: it enables the encoding of prior knowledge about the structure in images.
Wu et al. [30] point out that when evaluating test-densities for generative models, the assumed noise variance σ 2 plays an important role, so for both models we compare two different cases: one with the likelihood variance parameter fixed and one where it is optimized. Table 2 shows that in low-data regimes the highly parametrized CVAE severely overfits to the data and underestimates the variance. The GP-CDE operates much more gracefully in these regimes: it estimates the variance correctly, even for N = 2 (where N is the number of training points), and the gap between train/test log-likelihood is considerably smaller. Necessity of natural gradients Natural gradients are a vital component of our approach. We demonstrate this with the simplest possible example modeling a dataset of 100 '1' digits, using an unconditional model with no projection matrices, no mini-batches and no recognition model (i.e. exactly the GP-LVM in Titsias and Lawrence [26]). We compare our natural gradient approach with step size of 0.1 against using the Adam optimizer (learning rate 0.001) directly for the variational parameters. We compare also to the analytic solution in Titsias and Lawrence [26], which is possible as we are not using mini-batches. We find that the analytic model and our natural gradient method obtain test log-likelihoods (using all the '1's in the standard testing set) of 1.02, but the ordinary gradient approach attains a test log-likelihood of only −0.13. See Fig. 9 in Appendix F for samples from the latent space, and Fig. 10 for the training curves. We see that the ordinary gradient model cannot find a good solution, even in a large number of iterations, but the natural gradient model performs similarly to the analytic case.

Conclusion
We presented a model for conditional density estimation with Gaussian processes. Our approach extends prior work in three significant ways. We perform Bayesian linear transformations on both input and output spaces to allow for the modeling of high-dimensional inputs and strongly-coupled outputs. Our model is able to operate in low and high data regimes. Compared with other approaches we have shown that our model does not over-concentrate its density, even with very few data. For inference, we derived an optimal posterior for the latent variable inputs and we demonstrated the usefulness of natural gradients for mini-batched training of GPs with uncertain inputs. These improvements provide us with a more accurate variational approximation, and allow us to scale to larger datasets than were previous possible. We applied the model in different settings across a wide range of dataset sizes and input/output domains, demonstrating its general utility.

A Notation
prior covariance function of the d th GP

B Network Architectures
In all neural network setups we apply dropout to the output of the hidden layers. The optimal dropout rate is found using grid-search over {0.2, 0.1, 0.01, 0.0}. We use the Adam optimizer for optimization and perform grid-search over {0.01, 0.001, 0.0001} to determine the optimal learning rate. The biases are initialized to zero and the weights using the Xavier distribution [14].
Heteroscedastic noise modeling on UCI datasets Given the modest size of the UCI datasets we choose a relatively small encoder and decoder architecture. The encoder has layers of the following size: D x + D y , 50, 100, 50 and 2 × D w . The decoder layers have size D x + D w , 10, 50, 50, 10 and D y . For this experiment we choose a unidimensional latent variable, D w = 1. The targets of the UCI datasets are also 1D. We use the tanh activation function for all hidden layers and the linear activation function for the final layer.
Density estimation on MNIST In this experiment we follow Wu et al. [30] for the network architecture. The decoder has 5 fully connected layers of size: D x + D w , 64, 256, 256, 256, 1024 and D y . The tanh activation function is applied to the outputs of the hidden layers and the sigmoid function on the final one. The encoder's fully connected layers have size: D x + D y , 256, 64, 2 × D w . We use the tanh activation function for the encoder's hidden layers and the linear activation function for the final layer. The inputs are one-hot encoding of the labels, D x = 10 and D y = 784.

C Inference details
In section 3.2 we follow Hensman et al. [15] and use sparse GPs to approximate the full GP. We do this by introducing M inducing or pseudo-inputs z m ∈ R Dx+Dw and collect them in the matrix Z = {x m } M m=1 . The inducing outputs u = f (Z) are the function f (·) evaluated at the inducing inputs. Similarly, we collect the inducing outputs over all dimensions in the matrix U = {u } L =1 . We assume a Gaussian prior for the inducing outputs p(u ) = N (0, K) and define a posterior distribution of the form q(u ) = N (m , S ). The sparse GP framework assumes u and f to be jointly Gaussian. We can now write f (·) conditioned on the inducing outputs u as where [k Z (·)] m = k(·, z m ) and [K ZZ ] mm = k(z m , z m ). Marginalizing with respect to u ∼ N (m , S ) leads to the variational posterior q(f (·)) As both, the variational posterior over the GP and the likelihood, are Gaussian, L wn (defined in (3)) can be calculated analytically. We start by using the fact that the likelihood factorizes over the output dimension E q(f (·)) log p(y n | f (·), x n , w n ) = E q(f (·)) log p(y n,l | f (·), x n , w n ) . We define a single term of the previous sum as L wn , which can be computed as L wn = E q(f (·)) − 1 2 log(2πσ 2 ) − 1 2σ 2 y 2 n, + f 2 n, − 2 y n, f n, = − 1 2 log(2πσ 2 ) − 1 2σ 2 y 2 n, + (σ 2 ([x n , w n ]) − µ 2 ([x n , w n ])) − 2 y n, µ ([x n , w n ]) .

D Proof of optimality for latent variable posterior
In Section 3.2 we argue that using the optimal free-form distribution q(w n ) leads to a tighter lower bound. The bound for the optimal and Gaussian form q(w n ) are, respectively      [26]. Middle: our natural gradient approach, with a step size of 0.1. Bottom: using the Adam optimizer to optimize the variational parameters m and S. We see that the ordinary gradient approach is prone to getting 'stuck' in poor local optima, due to the difficulty of optimization. ELBO optimal nat grads ordinary grads Figure 10: The training objective for the dataset of 100 '1's. The three plots show the same three curves with different ranges on the y-axis to highlight the similarities and differences. We can see that natural gradients provide a striking improvement. See Fig. 9 for the samples at the end of training