Reliable training and estimation of variance networks

arXiv:1906.03260v2 [stat.ML] 4 Nov 2019

Nicki S. Detlefsen∗ †nsde@dtu.dk

Martin Jørgensen* †marjor@dtu.dk

Søren Hauberg †sohau@dtu.dk

Abstract
We propose and investigate new complementary methodologies for estimating
predictive variance networks in regression neural networks. We derive a locallyaware mini-batching scheme that results in sparse robust gradients, and we show
how to make unbiased weight updates to a variance network. Further, we formulatea heuristic for robustly fitting both the mean and variance networks post hoc. Finally,
we take inspiration from posterior Gaussian processes and propose a networkarchitecture with similar extrapolation properties to Gaussian processes. The
proposed methodologies are complementary, and improve upon baseline methods
individually. Experimentally, we investigate the impact of predictive uncertainty on
multiple datasets and tasks ranging from regression, active learning and generativemodeling. Experiments consistently show significant improvements in predictive
uncertainty estimation over state-of-the-art methods across tasks and datasets.

1

Introduction

The quality of mean predictions has dramatically increased in the last decade with the rediscovery ofneural networks [LeCun et al., 2015]. The predictive variance, however, has turned out to be a moreelusive target, with established solutions being subpar. The general finding is that neural networks
tend to make overconfident predictions [Guo et al., 2017] that can be harmful or offensive [Amodei
et al., 2016]. This may be explained by neural networks being general function estimators that does
not come with principled uncertainty estimates. Another explanation is that variance estimation is afundamentally different task than mean estimation, and that the tools for mean estimation perhaps do
not generalize. We focus on the latter hypothesis within regression.
To illustrate the main practical problems in variance estimation, we
consider a toy problem where data is generated as y = x · sin(x) +0.3·1 +0.3·x·2 , with 1 , 2 ∼ N (0, 1) and x is uniform on [0, 10](Fig. 1). As is common, we do maximum likelihood estimation of
N (µ(x), σ 2 (x)), where µ and σ 2 are neural nets. While µ provides
an almost perfect fit to the ground truth, σ 2 shows two problems: σ 2
is significantly underestimated and σ 2 does not increase outside thedata support to capture the poor mean predictions.
Figure 1: Max. likelihood fit
These findings are general (Sec. 4), and alleviating them is the mainof N (µ(x), σ 2 (x)) to data.
purpose of the present paper. We find that this can be achieved by a
combination of methods that 1) change the usual mini-batching to be location aware; 2) only optimize
variance conditioned on the mean; 3) for scarce data, we introduce a more robust likelihood function;and 4) enforce well-behaved interpolation and extrapolation of variances. Points 1 and 2 are achieved
through changes to the training algorithm, while 3 and 4 are changes to model specifications. We
empirically demonstrate that these new tools significantly improve on state-of-the-art across datasets
in tasks ranging from regression to active learning, and generative modeling.
∗†

Equal contribution
Section for Cognitive Systems, Technical University of Denmark

33rd Conference on Neural Information Processing Systems (NeurIPS 2019), Vancouver, Canada.

2

Related work

Gaussian processes (GPs) are well-known function approximators with built-in uncertainty estimators [Rasmussen and Williams, 2006]. GPs are robust in settings with a low amount of data, andcan model a rich class of functions with few hyperparameters. However, GPs are computationallyintractable for large amounts of data and limited by the expressiveness of a chosen kernel. Advances
like sparse and deep GPs [Snelson and Ghahramani, 2006, Damianou and Lawrence, 2013] partiallyalleviate this, but neural nets still tend to have more accurate mean predictions.
Uncertainty aware neural networks model the predictive mean and variance as two separate neuralnetworks, often as multi-layer perceptrons. This originates with the work of Nix and Weigend [1994]
and Bishop [1994]; today, the approach is commonly used for making variational approximations[Kingma and Welling, 2013, Rezende et al., 2014], and it is this general approach we investigate.
Bayesian neural networks (BNN) [MacKay, 1992] assume a prior distribution over the network
parameters, and approximate the posterior distribution. This gives direct access to the approximate
predictive uncertainty. In practice, placing an informative prior over the parameters is non-trivial.
Even with advances in stochastic variational inference [Kingma and Welling, 2013, Rezende et al.,2014, Hoffman et al., 2013] and expectation propagation [Hernández-Lobato and Adams, 2015], it isstill challenging to perform inference in BNNs.
Ensemble methods represent the current state-of-the-art. Monte Carlo (MC) Dropout [Gal andGhahramani, 2016] measure the uncertainty induced by Dropout layers [Hinton et al., 2012] arguingthat this is a good proxy for predictive uncertainty. Deep Ensembles [Lakshminarayanan et al.,2017] form an ensemble from multiple neural networks trained with different initializations. Bothapproaches obtain ensembles of correlated networks, and the extent to which this biases the predictive
uncertainty is unclear. Alternatives include estimating confidence intervals instead of variances[Pearce et al., 2018], and gradient-based Bayesian model averaging [Maddox et al., 2019].
Applications of uncertainty include reinforcement learning, active learning, and Bayesian optimization [Szepesvári, 2010, Huang et al., 2010, Frazier, 2018]. Here, uncertainty is the crucial elementthat allows for systematically making a trade-off between exploration and exploitation. It has also
been shown that uncertainty is required to learn the topology of data manifolds [Hauberg, 2018].
The main categories of uncertainty are epistemic and aleatoric uncertainty [Kiureghian andDitlevsen, 2009, Kendall and Gal, 2017]. Aleatoric uncertainty is induced by unknown or unmeasured features, and, hence, does not vanish in the limit of infinite data. Epistemic uncertainty
is often referred to as model uncertainty, as it is the uncertainty due to model limitations. It is thistype of uncertainty that Bayesian and ensemble methods generally estimate. We focus on the overall
predictive uncertainty, which reflects both epistemic and aleatoric uncertainty.

3

Methods

The opening remarks (Sec. 1) highlighted two common problems that appear when µ and σ 2 areneural networks. In this section we analyze these problems and propose solutions.
Preliminaries. We assume that datasets D = {xi , yi }N
i=1 contain i.i.d. observations yi ∈ R, xi ∈RD . The targets yi are assumed to be conditionally Gaussian, pθ (y|x) = N (y|µ(x), σ 2 (x)), where2 }. The maximum likelihood estimateµ and σ 2 are continuous functions parametrized by θ = {θµ , θσP12(MLE) of the variance of i.i.d. observations {yi }Nisi=1
i (yi − µ̂) , where µ̂ is the sample
N −1mean. This MLE does not exist based on a single observation, unless the mean µ is known, i.e. themean is not a free parameter. When yi is Gaussian, the residuals (yi − µ)2 are gamma distributed.
3.1

A local likelihood model analysis

By assuming that both µ and σ 2 are continuous functions, we are implicitly saying that σ 2 (x) iscorrelated with σ 2 (x + δ) for sufficiently small δ, and similar for µ. Consider the local likelihood
estimation problem [Loader, 1999, Tibshirani and Hastie, 1987] at a point xi ,NXlog p̃θ (yi |xi ) =wj (xi ) log pθ (yj |xj ),
(1)
j=1

2

where wj is a function that declines as kxj − xi k increases, implying that the local likelihood at xi
is dependent on the points nearest to xi . Notice p̃θ (yi |xi ) = pθ (yi |xi ) if wj (xi ) = 1i=j . Consider,with this w, a uniformly drawn subsample (i.e. a standard mini-batch) of the data {xk }M
k=1 and its
corresponding stochastic gradient of Eq. 1 with respect to θσ2 . If for a point, xi , no points near it
are in the subsample, then no other point will influence the gradient of σ 2 (xi ), which will point in
the direction of the MLE, that is highly uninformative as it does not exist unless µ(xi ) is known.
Local data scarcity, thus, implies that while we have sufficient data for fitting a mean, locally wehave insufficient data for fitting a variance. Essentially, if a point is isolated in a mini-batch, allinformation it carries goes to updating µ and none is present for σ 2 .
If we do not use mini-batches, we encounter that gradients wrt. θµ and θσ2 will both be scaled with1
2σ 2 (x) meaning that points with small variances effectively have higher learning rates [Nix and
Weigend, 1994]. This implies a bias towards low-noise regions of data.
3.2

Horvitz-Thompson adjusted stochastic gradients

We will now consider a solution to this problem within the local likelihood framework, which will
give us a reliable, but biased, stochastic gradient for the usual (nonlocal) log-likelihood. We will thenshow how this can be turned into an unbiased estimator.
If we are to add some local information, giving more reliable gradients, we should choose a w in Eq.1
that reflects this. Assume for simplicity that wj (xi ) = 1kxi −xj k<d for some d > 0. The gradientof log p̃θ (y|xi ) will then be informative, as more than one observation will contribute to the local
variance if d is chosen appropriately. Accordingly, we suggest a practical mini-batching algorithm
that samples a random point xj and we let the mini-batch consist of the k nearest neighbors of xj .3
In order to allow for more variability in a mini-batch, we suggest sampling m points uniformly,
and then sampling n points among the k nearest neighbors of each of the m initially sampled
points. Note that this is a more informative sample, as all observations in the sample are likely to
influence the same subset of parameters in θ, effectively increasing the degrees of freedom4 , hence
the quality of variance estimation. In other words, if the variance network is sufficiently expressive,
our Monte Carlo gradients under this sampling scheme are of smaller variation and more sparse. In
the supplementary material, we empirically show that this estimator yields significantly more sparse
gradients, which results in improved convergence. Pseudo-code of this sampling-scheme, can befound in the supplementary material.
While such a mini-batch would give rise to an informative stochastic gradient, it would not be an
unbiased stochastic gradient of the (nonlocal) log-likelihood. This can, however, be adjusted by using
the Horvitz-Thompson (HT) algorithm [Horvitz and Thompson, 1952], i.e. rescaling the log-likelihood
contribution of each sample xj by its inclusion probability πj . With this, an unbiased estimate of the
log-likelihood (up to an additive constant) becomes


N 
X
X 1  11(yi − µ(xi ))2(yj − µ(xj ))22− log(σ 2 (xi )) −
≈
log(σ
(x
))
−
(2)
−
j
2
2σ 2 (xi )πj2
2σ 2 (xj )
i=1
xj ∈O

where O denotes the mini-batch. With the nearest neighbor mini-batching, the inclusion probabilities
can be calculated as follows. The probability that observation j is in the sample is n/k if it is among
the k nearest neighbors of one of the initial m points, which are chosen with probability m/N , i.e.
N

πj =

mXn
1j∈Ok (i) ,N i=1 k(3)

where Ok (i) denotes the k nearest neighbors of xi .
Computational costs The proposed sampling scheme requires an upfront computational cost ofO(N 2 D) before any training can begin. We stress that this is pre-training computation and not3

By convention, we say that the nearest neighbor of a point is the point itself.
Degrees of freedom here refers to the parameters in a Gamma distribution – the distribution of varianceestimators under Gaussian likelihood. Degrees of freedom in general is a quite elusive quantity in regression
problems.
43

updated during training. The cost is therefore relative small, compared to training a neural network
for small to medium size datasets. Additionally, we note that the search algorithm does not have tobe precise, and we could therefore take advantage of fast approximate nearest neighbor algorithms[Fu and Cai, 2016].
3.3

Mean-variance split training

The most common training strategy is to first optimize θµ assuming a constant σ 2 , and then proceed
to optimize θ = {θµ , θσ2 } jointly, i.e. a warm-up of µ. As previously noted, the MLE of σ 2 does
not exist when only a single observation is available and µ is unknown. However, the MLE does
exist when µ is known, in which case it is σ̂ 2 (xi ) = (yi − µ(xi ))2 , assuming that the continuity ofσ 2 is not crucial. This observation suggests that the usual training strategy is substandard as σ 2 isnever optimized assuming µ is known. This is easily solved: we suggest to never updating µ and σ 2simultaneously, i.e. only optimize µ conditioned on σ 2 , and vice versa. This reads as sequentially
optimizing pθ (y|θµ ) and pθ (y|θσ2 ), as we under these conditional distributions we may think of µ
and σ 2 as known, respectively. We will refer to this as mean-variance split training (MV).
3.4

Estimating distributions of variance

When σ 2 (xi ) is influenced by few observations, underestimation is still likely due to the left skewnessof the gamma distribution of σ̂i2 = (yi − µ(xi ))2 . As always, when in a low data regime, it is sensible
to be Bayesian about it; hence instead of point estimating σ̂i2 we seek to find a distribution. Notethat we are not imposing a prior, we are training the parameters of a Bayesian model. We choose
the inverse-Gamma distribution, as this is the conjugate prior of σ 2 when data is Gaussian. Thismeans θσ2 = {θα , θβ } where α, β > 0 are the shape and scale parameters of the inverse-Gamma
respectively. So the log-likelihood is now calculated by integrating out σ 2
Zlog pθ (yi ) = log N (yi |µi , σi2 )dσi2 = log tµi ,αi ,βi (yi ),(4)where σi2 ∼ I NV-G AMMA(αi , βi ) and αi = α(xi ), βi = β(xi ) are modeled as neural networks.
Having an inverse-Gamma prior changes the predictive distribution to a located-scaled5 Student-t
distribution, parametrized with µ, α and β. Further, the t-distribution is often used as a replacementof the Gaussian when data is scarce and the true variance is unknown and yields a robust regression[Gelman et al., 2014, Lange et al., 1989]. We let α and β be neural networks that implicitly determine
the degrees of freedom and the scaling of the distribution. Recall the higher the degrees of freedom,
the better the Gaussian approximation of the t-distribution.
3.5

Extrapolation architecture

If we evaluate the local log-likelihood (Eq. 1) at a point x0 far away from all data points, then
the weights wi (x0 ) will all be near (or exactly) zero. Consequently, the local log-likelihood isapproximately 0 regardless of the observed value y(x0 ), which should be interpreted as a large
entropy of y(x0 ). Since we are working with Gaussian and t-distributed variables, we can recreate
this behavior by exploiting the fact that entropy is only an increasing function of the variance. We canre-enact this behavior by letting the variance tend towards an a priori determined value η if x0 tends
Daway from the training data. Let {ci }L
i=1 be points in R that represent the training data, akin to
inducing points in sparse GPs [Snelson and Ghahramani, 2006]. Then define δ(x0 ) = mini kci −x0 kandσ̂ 2 (x0 ) = 1 − ν(δ(x0 )) σ̂θ2 + ην(δ(x0 )),(5)where ν : [0, ∞) 7→ [0, 1] is a surjectively increasing function. Then the variance estimate will go toη as δ → ∞ at a rate determined by ν. In practice, we choose ν to be a scaled-and-translated sigmoid
function: ν(x) = sigmoid((x + a)/γ), where γ is a free parameter we optimize during training anda ≈ −6.9077γ to ensure that ν(0) ≈ 0. The inducing points ci are initialized with k-means andoptimized during training. This choice of architecture is similar to that attained by posterior Gaussian
processes when the associated covariance function is stationary. It is indeed the behavior of theseestablished models that we aim to mimic with Eq. 5.
5This means y ∼ F , where F = µ + σt(ν). The explicit density can be found in the supplementary material.

44

Experiments

4.1

Regression

To test our methodologies we conduct multiple experiments in various settings. We compare ourmethod to state-of-the-art methods for quantifying uncertainty: Bayesian neural network (BNN)[Hernández-Lobato and Adams, 2015], Monte Carlo Dropout (MC-Dropout) [Gal and Ghahramani,
2016] and Deep Ensembles (Ens-NN) [Lakshminarayanan et al., 2017]. Additionally we compare totwo baseline methods: standard mean-variance neural network (NN) [Nix and Weigend, 1994] andGPs (sparse GPs (SGP) when standard GPs are not applicable) [Rasmussen and Williams, 2006]. We
refer to our own method(s) as Combined, since we apply all the methodologies described in Sec. 3.
Implementation details and code can be found in the supplementary material. Strict comparisonsof the models should be carefully considered; having two seperate networks to model mean andvariance seperately (as NN, Ens-NN and Combined) means that all the predictive uncertainty, i.e. bothaleatoric and episteminc, is modeled by the variance networks alone. BNN and MC-Dropout have ahigher emphasis on modeling epistemic uncertainty, while GPs have the cleanest separation of noiseand model uncertainty estimation. Despite the methods quantifying different types of uncertainty,
their results can still be ranked by test set log-likelihood, which is a proper scoring function.
Toy regression. We first return to the toy problem of Sec. 1, where we consider 500 points fromy = x · sin(x) + 0.3 · 1 + 0.3 · x · 2 , with 1 , 2 ∼ N (0, 1). In this example, the variance isheteroscedastic, and models should estimate larger variance for larger values of x. The results6 canbe seen in Figs. 2 and 3. Our approach is the only one to satisfy all of the following: capture theheteroscedasticity, extrapolate high variance outside data region and not underestimating within.
7UXHVWG
*3
11
0&'URSRXW
(QV11
%11
&RPELQHGVWG[








Figure 2: From top left to bottom right: GP, NN,
BNN, MC-Dropout, Ens-NN, Combined.

 





x

 

Figure 3: Standard deviation estimatesas a function of x.

Variance calibration. To our knowledge, no benchmark for quantifying variance estimation exists.
We propose a simple dataset with known uncertainty information. More precisely, we consider
weather data from over 130 years.7 Each day the maximum temperature is measured, and theuncertainty is then given as the variance in temperature over the 130 years. The fitted models canbe seen in Fig. 4. Here we measure performance by calculating the mean error in uncertainty:
PN
22
Err = N1 i=1 |σtrue(xi ) − σest(xi )|. The numbers are reported above each fit. We observe that our
Combined model achieves the lowest error of all the models, closely followed by Ens-NN and GP.
Both NN, BNN and MC-Dropout all severely underestimate the uncertainty.
Ablation study. To determine the influence of each methodology from Sec. 3, we experimented
with four UCI regression datasets (Fig. 5). We split our contributions in four: the locality sampler(LS), the mean-variance split (MV), the inverse-gamma prior (IG) and the extrapolating architecture(EX). The combined model includes all four tricks. The results clearly shows that LS and IG
methodologies has the most impact on test set log likelihood, but none of the methodologies perform
worse than the baseline model. Combined they further improves the results, indicating that the
proposed methodologies are complementary.
67

The standard deviation plotted for Combined, is the root mean of the inverse-Gamma.
https://mrcc.illinois.edu/CLIMATE/Station/Daily/StnDyBTD2.jsp

5Err=7.6

(UU 

(UU 

(UU 

(UU 










60











40









20



















(a)x GP

 







(b) xNN

0 050



100 150 200 250 300 350 x(c) BNN









 x(d) MC-Dropoutyy

80

y

100



y

yy(UU 
















x(e) Ens-NN

 









x(f) Combined

Figure 4: Weather data with uncertainties. Dots are datapoints, green lines are the true uncertainty,
blue curves are mean predictions and the blue shaded areas are the estimated uncertainties.

NN
NN+LS
NN+MV
NN+EX
NN+IG
Combined

Concrete

Boston

NN
NN+LS
NN+MV
NN+EX
NN+IG
Combined

3.75 3.50 3.25 3.00 2.75 2.50 2.25 2.00

logp(x)Wine (red)

NN
NN+LS
NN+MV
NN+EX
NN+IG
Combined

4.0

3.5

3.0

logp(x)

2.5

2.0

NN
NN+LS
NN+MV
NN+EX
NN+IG
Combined

1.0

0.5

0.0

logp(x)0.5

1.0

1.5 1.5

Wine (white)

1.0

0.5

0.0

logp(x)0.5

1.0

Figure 5: The complementary methodologies from Sec. 3 evaluated on UCI benchmark datasets.

UCI benchmark. We now follow the experimental setup from Hernández-Lobato and Adams[2015], by evaluating models on a number of regression datasets from the UCI machine learning
database. Additional to the standard benchmark, we have added 4 datasets. Test set log-likelihood canbe seen in Table 1, and the corresponding RMSE scores can be found in the supplementary material.
Our Combined model performs best on 10 of the 13 datasets. For the small Boston and Yacht datasets,
the standard GP performs the best, which is in line with the experience that GPs perform well whendata is scarce. On these datasets our model is the best-performing neural network. On the Energy
and Protein datasets Ens-NN perform the best, closely followed by our Combined model. One clear
advantage of our model compared to Ens-NN is that we only need to train one model, whereas
Ens-NN need to train 5+ (see the supplementary material for training times for each model). The
worst performing model in all cases is the baseline NN model, which clearly indicates that the usual
tools for mean estimation does not carry over to variance estimation.
Active learning. The performance of active learning depends on predictive uncertainty [Settles,
2009], so we use this to demonstrate the improvements induced by our method. We use the same
network architectures and datasets as in the UCI benchmark. Each dataset is split into: 20% train, 60%
pool and 20% test. For each active learning iteration, we first train a model, evaluate the performanceon the test set and then estimate uncertainty for all datapoints in the pool. We then select the n points
with highest variance (corresponding to highest entropy [Houlsby et al., 2012]) and add these to the
N
D
GP
SGP
NN
BNN
MC-Dropout
Ens-NN
Combined
Boston506
13 −1.76 ± 0.3 −1.85 ± 0.25 −3.64 ± 0.09 −2.59 ± 0.11 −2.51 ± 0.31 −2.45 ± 0.25 −2.09 ± 0.09
Carbon10721 7
- 3.74 ± 0.53 −2.03 ± 0.14 −1.1 ± 1.76 −1.08 ± 0.05 −0.44 ± 7.28 4.35 ± 0.16
Concrete
1030 8 −2.13 ± 0.14 −2.29 ± 0.12 −4.23 ± 0.07 −3.31 ± 0.05 −3.11 ± 0.12 −3.06 ± 0.32 −1.78 ± 0.04
Energy768
8 −1.85 ± 0.34 −2.22 ± 0.15 −3.78 ± 0.04 −2.07 ± 0.08 −2.01 ± 0.11 −1.48 ± 0.31 −1.68 ± 0.13
Kin8nm8192 8- 2.01 ± 0.02 −0.08 ± 0.02 0.95 ± 0.08 0.95 ± 0.15 1.18 ± 0.03 2.49 ± 0.07
Naval11934 16
- 3.47 ± 0.21 3.71 ± 0.05 3.80 ± 0.09 5.55 ± 0.05 7.27 ± 0.13
Power plant 9568 4- −1.9 ± 0.03 −4.26 ± 0.14 −2.89 ± 0.01 −2.89 ± 0.14 −2.77 ± 0.04 −1.19 ± 0.03
Protein45730 9- −2.95 ± 0.09 −2.91 ± 0.00 −2.93 ± 0.14 −2.80 ± 0.02 −2.83 ± 0.05
Superconduct 21263 81- −4.07 ± 0.01 −4.92 ± 0.10 −3.06 ± 0.14 −2.91 ± 0.19 −3.01 ± 0.05 −2.43 ± 0.05
Wine (red) 1599 11 0.96 ± 0.18 −0.08 ± 0.01 −1.19 ± 0.11 −0.98 ± 0.01 −0.94 ± 0.01 −0.93 ± 0.09 1.21 ± 0.23
Wine (white) 4898 11- −0.14 ± 0.05 −1.29 ± 0.09 −1.41 ± 0.17 −1.26 ± 0.01 −0.99 ± 0.06 0.40 ± 0.42
Yacht
308
7 0.16 ± 1.22 −0.38 ± 0.32 −4.12 ± 0.17 −1.65 ± 0.05 −1.55 ± 0.12 −1.18 ± 0.21 −0.07 ± 0.05
Year515345 90- −5.21 ± 0.87 −3.97 ± 0.34 −3.78 ± 0.01 −3.42 ± 0.02 −3.01 ± 0.14

Table 1: Dataset characteristics and tests set log-likelihoods for the different methods. A - indicatesthe model was infeasible to train. Bold highlights the best results.

6

training set. We set n = 1% of the initial pool size. This is repeated 10 times, such that the last modelis trained on 30%. We repeat this on 10 random training-test splits to compute standard errors.
Fig. 6,show the evolution of average RMSE for each method during the data collection process for
the Boston, Superconduct and Wine (white) datasets (all remaining UCI datasets are visualized in the
supplementary material). In general, we observe two trends. For some datasets we observe that our
Combined model outperforms all other models, achieving significantly faster learning. This indicatesthat our model is better at predicting the uncertainty of the data in the pool set. On datasets where thesampling process does not increase performance, we are on par with other models.

Figure 6: Average test set RMSE and standard errors in active learning. The remaining datasets areshown in the supplementary material.
4.2

Generative models

To show a broader application of our approach, we also explore it in the context of generativemodeling. We focus on variational autoencoders (VAEs) [Kingma and Welling, 2013, Rezende et al.,2014] that are popular deep generative models. A VAE model the generative process:Zp(x) = pθ (x|z)p(z)dz, pθ (x|z) = N x|µθ (z), σθ2 (z) or pθ (x|z) = B x|µθ (z) , (6)where p(z) = N (0, Id ).
This is trained by introducing a variational approximationqφ (z|x) = N (z|µφ (x), σφ2 (x)) and then jointly training pθ and qφ . For our purposes, it
suffcient to note that a VAE estimates both a mean and a variance function. Thus using standard
training methods, the same problems arise as in the regression setting. Mattei and Frellsen [2018]
have recently shown that estimating a VAE is ill-posed unless the variance is bounded from below.
In the literature, we often find that
1. Variance networks are avoided by using a Bernoulli distribution, even if data is not binary.
2. Optimizing VAEs with a Gaussian posterior is considerably harder than the Bernoulli case. To
overcome this, the variance is often set to a constant e.g. σ 2 (z) = 1. The consequence is that the
log-likelihood reconstruction term in the ELBO collapses into an L2 reconstruction term.
3. Even though the generative process is given by Eq. 6, samples shown in the literature are often
reduced to x̃ = µ(z), z ∼ N (0, I). This is probably due to the wrong/meaningless variance term.
We aim to fix this by training the posterior variance σθ2 (z) with our Combined method. We do notchange the encoder variance σφ2 (x) and leave this to future study.
Artificial data. We first evaluate the benefits of more reliable variance networks in VAEs on
artificial data. We generate data inspired by the two moon dataset8 , which we map into four
dimensions. The mapping is thoroughly described in the supplementary material, and we emphasizethat we have deliberately used mappings that MLP’s struggle to learn, thus with a low capacity
network the only way to compensate is to learn a meaningful variance function.
In Fig. 7 we plot pairs of output dimensions using 5000 generated samples. For all pairwise
combinations we refer to the supplementary material. We observe that samples from our Comb-VAE
capture the data distribution in more detail than a standard VAE. For VAE the variance seems to be8

https://scikit-learn.org/stable/modules/generated/sklearn.datasets.make_moons.
html

7

0

2

0
1

3
2
1
0
1
2

0

3

2

2

1

1

0

0

1

1

0.6
0.5
0.4
0.3
0.2
0.1

35

3
2
1
0
1
2

30
3
25
2
20
15
1
10
5
0
1
Figure 8: Variance estimates in latent space for2
Groundtruth
00
2
2 02 0 2 standard
2 Comb-VAE
22
2
020 VAE02
0Comb-VAE
22 VAE
2 0 (top)
0 2 and2
our(bottom).
1
0 1 2 31 2 3
2 3
3 Figure 7:0 The
3Blue points are the encoded
ground truth and generated distributions.
training data.
2 Top: x1 vs. x2 . Bottom:
2 x2 vs x3 .
2
1
1
1
MNIST
FashionMNIST CIFAR10SVHN
00 2053.01 ± 1.60 1506.31
0
VAE± 2.71 1980.84 ± 3.32
3696.35 ± 2.94
ELBO
Comb-VAE
2152.31±
3.32
1621.29±7.23
2057.32±8.13
3701.41 ± 5.84
1
1
1
VAE
1914.77 ± 2.15 1481.38 ± 3.68 1809.43 ± 10.32 3606.28 ± 2.75
2 log p(x) Comb-VAE
2 2018.37 ± 4.35 1567.23
2± 4.82 1891.39 ± 20.21 3614.39 ± 7.91

2

3
2
1
0
1
2

3

1

2
1

1

0
1
2

Table 2: Generative modeling of 4 datasets. For each dataset we report training ELBO and test set
2
0 The2standard2 errors0are calculated
22
2 3 trained
02 models
022 with
2 random
02 initialization.
202 0 2log-likelihood.
over

0

0 1

0 1 2

1 2 3

underestimated, which is similar to the results from regression. The poor sample quality of a standard
VAE can partially be explained by the arbitrariness of decoder variance function σ 2 (z) away fromPD
data. In Fig. 8, we calculated the accumulated variance j=1 σj2 (z) over a grid of latent points.
We clearly see that for the standard VAE, the variance is low where we have data and arbitrary awayfrom data. However, our method produces low-variance region where the two half moons are anda high variance region away from data. We note that Arvanitidis et al. [2018] also dealt with the
problem of arbitrariness of the decoder variance. However their method relies on post-fitting of thevariance, whereas ours is fitted during training. Additionally, we note that [Takahashi et al., 2018]
also successfully modeled the posterior of a VAE as a Student t-distribution similar to our proposed
method, but without the extrapolation and different training procedure.
Image data. For our last set of experiments we fitted a standard VAE and our Comb-VAE tofour datasets: MNIST, FashionMNIST, CIFAR10, SVHN. We want to measure whether there is an
improvement to generative modeling by getting better variance estimation. The details about networkarchitecture and training can be found in the supplementary material. Training set ELBO and test
set log-likelihoods can be viewed in Table 2. We observe on all datasets that, on average tighter
bounds and higher log-likelihood are achieved, indicating that we better fit the data distribution. We
quantitatively observe (see Fig. 9) that variance has a more local structure for Comb-VAE and thatthe variance reflects the underlying latent structure.

5

2

0

2

20

02

3

3

0

3

2

2

3
2

3
2

3
2
1
0
1
2

3
2
1
0
1
2

3

3

2
1

3
2
1
0
1
2

1
0

1 0
0 1
1 2

1
0

1 0
0 1
1 2

Discussion & Conclusion

While variance networks are commonly used for modeling the predictive uncertainty in regression
and in generative modeling, there have been no systematic studies of how to fit these to data. We
have demonstrated that tools developed for fitting mean networks to data are subpar when applied to
8

2 2

2 3

Figure 9: Generated MNIST images on a grid in latent space using the standard variance network(left) and proposed variance network (right).
variance estimation. The key underlying issue appears to be that it is not feasible to estimate both amean and a variance at the same time, when data is scarce.
While it is beneficial to have separate estimates of both epistemic and aleatoric uncertainty, we havefocused on predictive uncertainty, which combine the two. This is a lesser but more feasible goal.
We have proposed a new mini-batching scheme that samples locally to ensure that variances are better
defined during model training. We have further argued that variance estimation is more meaningfulwhen conditioned on the mean, which implies a change to the usual training procedure of joint
mean-variance estimation. To cope with data scarcity we have proposed a more robust likelihood thatmodel a distribution over the variance. Finally, we have highlighted that variance networks need to
extrapolate differently from mean networks, which implies architectural differences between suchnetworks. We specifically propose a new architecture for variance networks that ensures similarvariance extrapolations to posterior Gaussian processes from stationary priors.
Our methodologies depend on algorithms that computes Euclidean distances. Since these often break
down in high dimensions, this indicates that our proposed methods may not be suitable for highdimensional data. Since we mostly rely on nearest neighbor computations, that empirical are known
to perform better in high dimensions, our methodologies may still work in this case. Interestingly, thevery definition of variance is dependent on Euclidean distance and this may indicate that varianceis inherently difficult to estimate for high dimensional data. This could possible be circumvented
through a learned metric.
Experimentally, we have demonstrated that proposed methods are complementary and provide
significant improvements over state-of-the-art. In particular, on benchmark data we have shownthat our method improves upon the test set log-likelihood without improving the RMSE, which
demonstrate that the uncertainty is a significant improvement over current methods. Another indicatorof improved uncertainty estimation is that our method speeds up active learning tasks compared
to state-of-the-art. Due to the similarities between active learning, Bayesian optimization, andreinforcement learning, we expect that our approach carries significant value to these fields as well.
Furthermore, we have demonstrated that variational autoencoders can be improved through better
generative variance estimation. Finally, we note that our approach is directly applicable alongsideensemble methods, which may further improve results.

9

Acknowledgements. This project has received funding from the European Research Council (ERC)
under the European Union’s Horizon 2020 research and innovation programme (grant agreement
no 757360). NSD, MJ and SH were supported in part by a research grant (15334) from VILLUM
FONDEN. We gratefully acknowledge the support of NVIDIA Corporation with the donation of
GPU hardware used for this research.

References
M. Abadi, A. Agarwal, P. Barham, E. Brevdo, Z. Chen, C. Citro, G. S. Corrado, A. Davis, J. Dean, M. Devin,
S. Ghemawat, I. Goodfellow, A. Harp, G. Irving, M. Isard, Y. Jia, R. Jozefowicz, L. Kaiser, M. Kudlur, J. Levenberg, D. Mané, R. Monga, S. Moore, D. Murray, C. Olah, M. Schuster, J. Shlens, B. Steiner, I. Sutskever,K. Talwar, P. Tucker, V. Vanhoucke, V. Vasudevan, F. Viégas, O. Vinyals, P. Warden, M. Wattenberg, M. Wicke,
Y. Yu, and X. Zheng. TensorFlow: Large-scale machine learning on heterogeneous systems, 2015. URL
https://www.tensorflow.org/. Software available from tensorflow.org.
D. Amodei, C. Olah, J. Steinhardt, P. Christiano, J. Schulman, and D. Mané. Concrete problems in ai safety.
arXiv preprint arXiv:1606.06565, 2016.
G. Arvanitidis, L. K. Hansen, and S. Hauberg. Latent space oddity: on the curvature of deep generative models.
In International Conference on Learning Representations, 2018.
C. M. Bishop. Mixture density networks. Technical report, Citeseer, 1994.
A. Damianou and N. D. Lawrence. Deep gaussian processes. Proceedings of the 16th International Conferenceon Artificial Intelligence and Statistics (AISTATS), 2013.
P. I. Frazier. A tutorial on bayesian optimization. arXiv preprint arXiv:1807.02811, 2018.
C. Fu and D. Cai. Efanna : An extremely fast approximate nearest neighbor search algorithm based on knn
graph. 09 2016.
Y. Gal and Z. Ghahramani. Dropout as a bayesian approximation: Representing model uncertainty in deeplearning. In international Conference on Machine Learning, pages 1050–1059, 2016.
A. Gelman, J. B. Carlin, H. S. Stern, D. B. Dunson, A. Vehtari, and D. B. Rubin. Bayesian Data Analysis. CRC
Press, 2014.
GPy. GPy: A gaussian process framework in python. http://github.com/SheffieldML/GPy, since 2012.
C. Guo, G. Pleiss, Y. Sun, and K. Q. Weinberger. On calibration of modern neural networks. In Proceedings ofthe 34th International Conference on Machine Learning-Volume 70, pages 1321–1330. JMLR. org, 2017.
S. Hauberg. Only bayes should learn a manifold (on the estimation of differential geometric structure from data).
arXiv preprint arXiv:1806.04994, 2018.
J. M. Hernández-Lobato and R. Adams. Probabilistic backpropagation for scalable learning of bayesian neuralnetworks. In International Conference on Machine Learning, pages 1861–1869, 2015.
G. Hinton, N. Srivastava, A. Krizhevsky, I. Sutskever, and R. R. Salakhutdinov. Improving neural networks by
preventing co-adaptation of feature detectors. arXiv preprint, arXiv, 07 2012.
M. D. Hoffman, D. M. Blei, C. Wang, and J. Paisley. Stochastic variational inference. The Journal of Machine
Learning Research, 14(1):1303–1347, 2013.
D. G. Horvitz and D. J. Thompson. A generalization of sampling without replacement from a finite universe.
Journal of the American Statistical Association, 47(260):663–685, 1952.
N. Houlsby, F. Huszar, Z. Ghahramani, and J. M. Hernández-lobato. Collaborative gaussian processes forpreference learning. In F. Pereira, C. J. C. Burges, L. Bottou, and K. Q. Weinberger, editors, Advances in
Neural Information Processing Systems 25, pages 2096–2104. Curran Associates, Inc., 2012.
S.-J. Huang, R. Jin, and Z.-H. Zhou. Active learning by querying informative and representative examples. In
Advances in neural information processing systems, pages 892–900, 2010.
A. Kendall and Y. Gal. What uncertainties do we need in bayesian deep learning for computer vision? In
Advances in neural information processing systems, pages 5574–5584, 2017.

10

D. P. Kingma and J. Ba. Adam: A method for stochastic optimization. In 3rd International Conference on
Learning Representations, ICLR 2015, San Diego, CA, USA, May 7-9, 2015, Conference Track Proceedings,
2015. URL http://arxiv.org/abs/1412.6980.
D. P. Kingma and M. Welling. Auto-encoding variational bayes. ICLR, 12 2013.
A. D. Kiureghian and O. Ditlevsen. Aleatory or epistemic? does it matter? Structural Safety, 31(2):105 – 112,
2009. Risk Acceptance and Risk Communication.
B. Lakshminarayanan, A. Pritzel, and C. Blundell. Simple and scalable predictive uncertainty estimation using
deep ensembles. In Advances in Neural Information Processing Systems, pages 6402–6413, 2017.
K. L. Lange, R. J. A. Little, and J. M. G. Taylor. Robust statistical modeling using the t distribution. Journal of
the American Statistical Association, 84(408):881–896, 1989.
Y. LeCun, Y. Bengio, and G. E. Hinton. Deep learning. Nature, 521(7553):436–444, 2015.
C. Loader. Local Regression and Likelihood. Springer, New York, 1999.
D. J. C. MacKay. A Practical Bayesian Framework for Backpropagation Networks. Neural Comput., 4(3):448–472, may 1992.
W. Maddox, T. Garipov, P. Izmailov, D. Vetrov, and A. G. Wilson. A Simple Baseline for Bayesian Uncertainty
in Deep Learning. CoRR, feb 2019.
P.-A. Mattei and J. Frellsen. Leveraging the exact likelihood of deep latent variable models. In Proceedings ofthe 32Nd International Conference on Neural Information Processing Systems, NIPS’18, pages 3859–3870,
USA, 2018. Curran Associates Inc.
D. Nix and A. Weigend. Estimating the mean and variance of the target probability distribution. In Proc. 1994
IEEE Int. Conf. Neural Networks, pages 55–60 vol.1. IEEE, 1994.
A. Paszke, S. Gross, S. Chintala, G. Chanan, E. Yang, Z. DeVito, Z. Lin, A. Desmaison, L. Antiga, and A. Lerer.
Automatic differentiation in pytorch. 2017.
T. Pearce, M. Zaki, A. Brintrup, and A. Neely. High-Quality Prediction Intervals for Deep Learning: A
Distribution-Free, Ensembled Approach. In Proceedings of the 35th International Conference on Machine
Learning, feb 2018.
C. E. Rasmussen and C. Williams. Gaussian Processes for Machine Learning. University Press Group Limited,2006.
D. J. Rezende, S. Mohamed, and D. Wierstra. Stochastic backpropagation and approximate inference in deep
generative models. In E. P. Xing and T. Jebara, editors, Proceedings of the 31st International Conference
on Machine Learning, volume 32 of Proceedings of Machine Learning Research, pages 1278–1286, Bejing,
China, 22–24 Jun 2014. PMLR.
B. Settles. Active learning literature survey. Technical report, University of Wisconsin-Madison Department ofComputer Sciences, 2009.
E. Snelson and Z. Ghahramani. Sparse gaussian processes using pseudo-inputs. In Advances in neural
information processing systems, pages 1257–1264, 2006.
C. K. Sønderby, T. Raiko, L. Maalø e, S. r. K. Sø nderby, and O. Winther. Ladder variational autoencoders. In
D. D. Lee, M. Sugiyama, U. V. Luxburg, I. Guyon, and R. Garnett, editors, Advances in Neural Information
Processing Systems 29, pages 3738–3746. Curran Associates, Inc., 2016. URL http://papers.nips.cc/
paper/6275-ladder-variational-autoencoders.pdf.
C. Szepesvári. Algorithms for reinforcement learning. Synthesis lectures on artificial intelligence and machinelearning, 4(1):1–103, 2010.
H. Takahashi, T. Iwata, Y. Yamanaka, M. Yamada, and S. Yagi. Student-t variational autoencoder for robust density estimation. In Proceedings of the Twenty-Seventh International Joint Conference on Artificial Intelligence,
IJCAI-18, pages 2696–2702. International Joint Conferences on Artificial Intelligence Organization, 7 2018.
R. Tibshirani and T. Hastie. Local likelihood estimation. Journal of the American Statistical Association, 82(398):559–567, 1987.

11

A

Further results

Gradient experiments In Fig. 10 we have plotted the sparsity index and variance of gradient for
both the mean (top row) and variance function (bottom row). We do this for both normal minibatching and our proposed locality sampler. Sparsity is measured as `00.001 (∇) = {j, ∇j ≤ 0.001}`0(∇)

and the sparsity index is then given by SI = 0.001. We observe for the mean function, that the|∇|sparsity index and variance is similar for the two methods, indicating that our locality sampler does
not improve on the fitting of the mean function, as expected. However for the variance function we
see a clear gap in sparsity and variance, indicating that our locality sampler gives more local andstable updates to variance networks.


PLQLEDWFK
ORFDOVDPSOHU

Var( (x))

6SDUVLW\LQGH[









,WHUDWLRQ











PLQLEDWFK
ORFDOVDPSOHUPLQLEDWFK
ORFDOVDPSOHU













,WHUDWLRQ

2(x))



Var(

6SDUVLW\LQGH[










PLQLEDWFK
ORFDOVDPSOHU








,WHUDWLRQ











,WHUDWLRQ



Figure 10: Left: Sparsity index for mean (top) and variance (bottom) network. Right: Variance of
gradient for mean (top) and variance (bottom) network. The variance network was disabled for the
first 2500 iterations, to warm up the mean function to secure convergence.
UCI benchmark (RMSE) In Table 3 the test set RMSE results for the UCI regression benchmark
can be seen. We clearly observe that all neural network based methods achieve nearly identical
RMSE for all datasets, indicating that the mean function is similarly trained for all the methods.
BostonCarbonConcrete
EnergyKin8nmNavalPower plantProtein
Superconduct
Wine (red)Wine (white)
Yacht
Year

N
506
10721
1030
768
8192
11934
9568
45730
21263
1599
4898
308
515345

D
13
7
8
8
8
16
4
9
81
11
11
7
90

GP
2.79 ± 0.52
6.03 ± 0.59
1.98 ± 0.76
0.88 ± 0.06
0.42 ± 0.21-

SGP
2.98 ± 0.55
1.01 ± 0.01
6.45 ± 0.64
2.12 ± 0.56
0.08 ± 0.0
4.65 ± 0.12
11.32 ± 0.38
0.65 ± 0.04
0.65 ± 0.03
0.72 ± 0.21-

NN
4.45 ± 1.41
0.41 ± 0.0
7.71 ± 1.32
1.67 ± 0.44
0.21 ± 0.01
0.01 ± 0.0
4.23 ± 0.33
4.38 ± 0.07
11.73 ± 0.46
0.66 ± 0.06
0.67 ± 0.04
1.63 ± 0.61
12.47 ± 0.96

BNN
3.45 ± 0.87
0.18 ± 0.1
5.78 ± 0.21
1.89 ± 0.04
0.18 ± 0.02
0.01 ± 0.0
4.12 ± 0.45
4.67 ± 0.94
11.07 ± 1.7
0.69 ± 0.41
0.68 ± 0.32
1.05 ± 0.11
9.01 ± 0.45

MC-Dropout
3.01 ± 0.99
0.29 ± 0.0
5.33 ± 0.65
1.69 ± 0.19
0.12 ± 0.01
0.01 ± 0.0
4.13 ± 0.13
4.19 ± 0.08
11.44 ± 0.39
0.64 ± 0.06
0.71 ± 0.04
1.11 ± 0.48
8.92 ± 0.23

Ens-NN3.33 ± 1.33
0.41 ± 0.0
5.65 ± 0.55
2.13 ± 0.46
0.01 ± 0.01
0.01 ± 0.0
4.11 ± 0.21
4.36 ± 0.07
11.63 ± 0.49
0.67 ± 0.06
0.78 ± 0.04
1.58 ± 0.58
8.88 ± 0.13

Combined
3.11 ± 0.35
0.35 ± 0.01
5.75 ± 0.41
1.70 ± 0.21
0.12 ± 0.01
0.01 ± 0.0
4.12 ± 0.13
4.52 ± 0.19
11.65 ± 0.65
0.68 ± 0.11
0.72 ± 0.09
1.27 ± 0.11
8.85 ± 0.05

Table 3: Dataset characteristics and RMSE for the different methods. A - indicates the models was
infeasible to train.

12

Timings of models In Table 4 we show the average computation time for each model. The
experiments was conducted with an Intel Xeon E5-2620v4 CPU and Nvidia GTX TITAN X GPU.
We note that our model suffers from long computations for very large datasets, mainly due to thecomputation of the neighborhood graph in the locality sampler. This could be reduced by using
fast approximative method for k-nearest-neighbor and by reducing data dimensionality e.g. PCA
dimensionality reduction.
BostonCarbonConcrete
EnergyKin8nmNavel
Power
Protein
Superconduct
Wine (red)Wine (white)
Yacht
Year

N
506
10721
1030
768
8192
11934
9568
45730
21263
1599
4898
308
515345

D
13
7
8
8
8
16
4
9
81
11
11
7
90

GP
8.37 +- 2.92
7.48 +- 1.2
10.48 +- 4.12
35.33 +- 19.09
0.93 +- 0.32-

SGP
91.86 +- 30.12
192.57 +- 72.28
173.73 +- 4.0
121.39 +- 52.2
1526.29 +- 20.28
9.79 +- 0.23
1267.82 +- 783.28
313.9 +- 2.9
262.69 +- 10.14
781.13 +- 8.2
22.74 +- 11.98-

NN
94.04 +- 2.24
90.05 +- 3.4
92.91 +- 4.17
92.91 +- 2.14
92.22 +- 4.65
91.2 +- 3.28
92.23 +- 3.28
138.49 +- 2.51
95.27 +- 3.43
93.51 +- 2.09
91.97 +- 2.41
92.82 +- 3.47
139.45 +- 6.15

BNN
81.32 +- 1.83
80.95 +- 2.04
81.02 +- 2.04
80.92 +- 2.14
80.65 +- 2.11
82.97 +- 1.84
81.29 +- 1.98
124.72 +- 1.86
90.71 +- 1.86
80.67 +- 2.04
80.94 +- 1.75
81.0 +- 1.92
643.96 +- 4.17

MC-Dropout
98.08 +- 2.0
98.62 +- 1.85
97.94 +- 1.84
97.86 +- 2.13
97.85 +- 2.87
98.64 +- 1.89
98.26 +- 1.87
140.73 +- 3.32
90.25 +- 1.67
97.78 +- 2.58
97.61 +- 2.73
98.27 +- 1.68
77.96 +- 2.09

Ens-NN
479.1 +- 12.49
439.45 +- 23.01
468.94 +- 10.81
475.04 +- 12.61
459.81 +- 27.66
482.74 +- 8.54
472.39 +- 9.73
707.69 +- 11.02
477.57 +- 11.91
416.02 +- 23.9
451.65 +- 19.04
422.9 +- 35.85
725.27 +- 20.08

Combined
93.39 +- 1.82
123.61 +- 2.84
97.65 +- 6.4
93.01 +- 1.29
123.15 +- 2.86
136.25 +- 3.29
118.26 +- 1.86
658.63 +- 11.75
235.72 +- 3.8
130.61 +- 7.91
99.44 +- 1.33
105.33 +- 2.32
2453.62 +- 18.7

Table 4: Timings(s) for the different models evaluated on the UCI benchmark

Ablation study (RMSE) In Fig. 11 we have plotted the RMSE for different combinations of our
methodologies. Since the RMSE is only influenced by how well µ(x) is fitted, the difference in log
likelihood that we observe between the models (see paper) must be explained by how well σ 2 (x) is
fitted.

Boston

NN
NN+LS
NN+MV
NN+LS+MV
Combined

Concrete

NN
NN+LS
NN+MV
NN+LS+MV
Combined

2.5

3.0

RMSE

3.5

4.0

3.0

NN
NN+LS
NN+MV
NN+LS+MV
Combined

3.5

4.0

4.5

RMSE

5.0

5.5

6.0

Wine (red)

NN
NN+LS
NN+MV
NN+LS+MV
Combined

Wine (white)

0.60 0.62 0.64 0.66 0.68 0.70 0.72 0.74 0.64 0.66 0.68 0.70 0.72 0.74 0.76

RMSE

RMSE

Figure 11: Our contributions evaluated on four different UCI benchmark datasets.

Active learning In Figs. 12 and 13 we respective shows the progress of RMSE and log likelihoodon all 13 dataset. We observe that for some of the datasets (Boston, Superconduct, Power) our
proposed Combine model achieves faster learning than other methods. On all other datasets we are
equally good as the best performing model.
Generative modeling toy data In Fig. 15 we show marginal distribution, pairwise pointplots andpairwise joint distribution for our artificially dataset used in the generative setting. Top row show
the ground true data, middle row shows reconstructions and samples from standard VAE model andbottom row show reconstructions and samples from our proposed Comb-VAE model. We observethat reconstruction are similar for the two models, but the quality of the generative samples are muchbetter for Comb-VAE.
Generative modeling of image data In Fig. 14 we show a meshgrid of samples from VAE andComb-VAE on the MNIST dataset. We clearly see how proper extrapolation of variance in CombVAE can be used to "mask" when we are inside our data region and when we are outside. For standard
VAE we observe a near constant variance added to the images.

B

On the parametrization of the t-distributed predictive distribution

We parametrize the Student-t distribution by letting the variance σ 2 have an inverse-Gamma distribution with shape and scale parameters α and β. We use that if σ 2 ∼ I NV-G AMMA(α, β) then
13

0.1

100.0

0.1

6.0

4.0

15.0

8.0
RMSE

6.0
5.0

RMSE

14.0

7.0

13.0

0 1 2 3 4 5 6 7 8 9
AL iter

40.0

0.1

20.0

0.1

Wine (red)0.1

0 1 2 3 4 5 6 7 8 9
AL iter

Wine (white)

10.0

0.0

0 1 2 3 4 5 6 7 8 9
AL iter

Yacht

1.5

14.0

40.0

3.0

1.2

12.0

35.0

1.0

10.0

30.0

8.0

45.0

6.0

20.0

1.5

0.2

4.0

15.0

1.0

0.0

2.0

0 1 2 3 4 5 6 7 8 9
AL iter

0 1 2 3 4 5 6 7 8 9
AL iter

0 1 2 3 4 5 6 7 8 9
AL iterYear
Ens-NN
Combined
MC-Dropout
NN
BNN

25.0

0.5

0 1 2 3 4 5 6 7 8 9
AL iter

6.0
4.0

0.0
0 1 2 3 4 5 6 7 8 9
AL iter

16.0

0.8

8.0

0.0

1.8

2.5

Power

0.1

3.5

2.0

12.0

4.0

0.1

Navel
0.1

0.1

60.0

0.0

0 1 2 3 4 5 6 7 8 9
AL iterSuperconduct

9.0

80.0

Kin8nmRMSE

Protein

4.0

0 1 2 3 4 5 6 7 8 9
AL iter

RMSE

0.4

0 1 2 3 4 5 6 7 8 9
AL iter

RMSE

120.0

0.4

3.0

3.0

0.1

RMSE

8.0

0.4

0.1

RMSE

10.0

0.4

Energy
140.0

RMSE

4.0

0.5
RMSE

5.0

RMSE

RMSE

6.0

Concrete
12.0

RMSE

Carbon0.5

RMSE

Boston10.0
0 1 2 3 4 5 6 7 8 9
AL iter

0 1 2 3 4 5 6 7 8 9
AL iter

Figure 12: Average test RMSE and standard errors in the active learning experiments for all datasets.

-6.0

-12.5

-60.0

Superconduct
-2.5

-3.5

-5.0

-4.5

-6.5

-5.0
-5.5

0.0

logp(x)

0.0

-7.5
-10.0

-2.0

-4.0

-3.0

-6.0

-4.0

0 1 2 3 4 5 6 7 8 9
AL iter

0 1 2 3 4 5 6 7 8 9
AL iterlogp(x)logp(x)

-8.0

-16.0

Ens-NN
Combined
MC-Dropout
NN
BNN

-6.0

-10.0

-8.0
-10.0
-12.0
-14.0

0 1 2 3 4 5 6 7 8 9
AL iter

0 1 2 3 4 5 6 7 8 9
AL iter

Year

-4.0

-14.0

-7.0

-2.0

-5.0

-7.0
0 1 2 3 4 5 6 7 8 9
AL iter

Yacht

-12.0

-6.0

-15.0

0.0
-2.0

-5.0

-12.5

0 1 2 3 4 5 6 7 8 9
AL iter

Wine (white)

-4.0

-6.0

-6.0
0 1 2 3 4 5 6 7 8 9
AL iter

-1.0

-3.0

0.0

-4.0

-50.0
0 1 2 3 4 5 6 7 8 9
AL iter

-2.0

-2.0

-35.0

Wine (red)

-3.0

-30.0

Power

-1.0

2.0

-40.0

-30.0
0 1 2 3 4 5 6 7 8 9
AL iter

-2.5

-4.0

-20.0
-25.0

-70.0
0 1 2 3 4 5 6 7 8 9
AL iterlogp(x)logp(x)4.0-20.0
logp(x)

-10.0

-50.0

-6.0

logp(x)-7.5

Protein0 1 2 3 4 5 6 7 8 9
AL iter

-10.0

-15.0

-30.0

-4.0

-7.0

-10.0

-20.0

-3.5

-5.5

6.0

-5.0

-40.0

Navel

Kin8nm
0.0-2.5

0 1 2 3 4 5 6 7 8 9
AL iter-5.0

Energy
-5.0

-15.0

-4.5

Concretelogp(x)

-5.0

0.0
-10.0

logp(x)-4.0

logp(x)logp(x)-3.0

Carbon0.0logp(x)2.5logp(x)Boston-2.0

0 1 2 3 4 5 6 7 8 9
AL iter

Figure 13: Average test log likelihood and standard errors in the active learning experiments for all
datasets.
1
σ 2 ∼ Γ(α, β). Then

Z ∞

βα
β
N (y|µ, σ 2 )(σ 2 )−(α+1) exp(− 2 )dσ 2
Γ(α)
σ
0
Z ∞
1
βα
1
1
√
=
(σ 2 )−(α+1)− 2 exp(− 2 (β + (y − µ)2 ))dσ 2σ2Γ(α) 2π 0Γ(α + 21 )βα√ =α+ 21 ,Γ(α) 2π12
β + 2 (y − µ)p(y|µ, α, β) =where we substituted the variable σ 2 with σ12 and used that the remaining was a Gamma integral.

C

Parameters of the locality sampler

In Algorithm 1 a pseudoimplementation of our proposed locality sampler can be seen. The two
important parameters in this algorithm are the primary sampling units (psu) and secondary sampling
units (ssu). In Fig. 16 and 17 we visually show the effect of these two parameters.
For all our experiments we set psu=3 and ssu=40 when we are training the mean function and forvariance function we set psu=1 and ssu=10.
14

(a) VAE(b) Comb-VAE

Figure 14: Meshgrid of latent space. For each subplot we sampled a mesh grid [−4, 4] × [−4, 4] of[50, 50] points, which we used to generate samples from.
Algorithm 1 Locality-sampler
Input N datapoints, a metric d on feature space RD , integers m, n, k.
1: For each datapoint calculate the k nearest neighbors under the metric d.
2: Sample m primary sampling units with uniform probability without replacement among all N
units.
3: For each of the primary sampling units sample n secondary sampling units among the primary
sampling units k nearest neighbors with uniform probability without replacement.
Output All secondary sampling units is a sample of at most m · n points. If a new sample is needed
repeat from Step 2.

D

Implementation details for regression experiments

All neural network based models were implemented in Pytorch [Paszke et al., 2017], except for theBaysian neural network which was implemented in Tensorflow [Abadi et al., 2015]. GP models
where implemented in GPy [GPy, since 2012]. Below we have stated details for all models:
GP Fitted using a ARD kernel and with default settings of GPy.
SGP Fitted using a ARD kernel and with default settings of GPy. Number of inducing points was
set to min(500, |Dtrain |).
NN Model use two networks, one for the predictive mean and one for the predictive variance. Model
was trained to optimize the log-likelihood of data.
BNN We use a standard factored Gaussian prior for the weights and use the Flipout approximation[?] for the variational approximation.
MC-Dropout Model use a single network, where we place dropout on all weights. The dropout
weight was set to 0.05. The model was trained to optimize the RMSE of data.
Ens-NN Model consist of an ensemble of 5 individual NN models, each modeled as two individual
networks. Each are trained to optimize the log-likelihood of data. Only difference between
ensemble models is initialization.
Combined Model use three networks, one for the mean function, one for the α parameter and one
for the β parameter. We set the number of inducing points to min(500, |Dtrain |). For the γin the scaled-and-translated sigmoid function ν(x) we initialize it to 1.5, and try to minimize
it during training. Model was trained to optimize the log-likelihood of data.
15

0

2
1
0
1
2
3

1

2
1
0

3

2

1
3
2
1
0
1
2
3
2
1
0
1
2
2

0

2

0

0

1

2

2

0

2

2

2

0

2

3

0

2
1
0
1
2

2
1
0
1
2
3

2

2

1

1

1

3

0

1

1

3
2
1
0
1
2

3
2
1
0
1
2

2

0

3
2
1
0
1
2

3

3

2

1

0

(a) Ground true

2

0

0

2

0

1

2

2

0

2

2

2

0

3

3
2
1
0
1
2

2

2

0

0

2

0

0

2
1
0
1
2

2

1

1

1

3

2

0

1

1

3
2
1
0
1
2

3
2
1
0
1
2

2

0

3
2
1
0
1
2
0

0

2

0

1

2

2

2

0

2

2

2

0

2

2

2

0

2

3

2
1
0
1
2

3

2

2

(c) VAE (S)

3

3

2

1

0

(b) VAE (R)

1

0

2

(d) Comb-VAE (R)

2

2

0

3

2

3
2
1
0
1
2
2

0

0

2

0

1

2

2

0

2

3

(e) Comb-VAE (S)Figure 15: Pairwise plots between all sets of variables for our artificially dataset in the generative
setting.
16

Figure 16: The effect of changing the size of the primary sampling unit. From top to bottom:
psu = [1, 2, 3]. Each column corresponds to a sample from the locality sampler.

Figure 17: The effect of changing the size of the secondary sampling unit. From top to bottom:
ssu = [10, 50, 100]. Each column corresponds to a sample from the locality sampler.

For each neural network based approach, we follow the experimental setup [Hernández-Lobato andAdams, 2015]. All individual networks was modeled a single hidden layer MLPs with 50 neurons forall other datasets than "Protein" and "Year" where we use 100 neurons. Except for the output of each
network, the activation function used is ReLU. For the output of the mean networks, no activation
function is applied. For the output of the variance network, the Softplus activation function is used tosecure positive variance. All neural network based models where trained using the Adam optimizer[Kingma and Ba, 2015] with a learning rate of 10−1 using a batch size of 256. All models weretrained for 10.000 iterations.
The code can found here: https://github.com/SkafteNicki/john.

E

Generative network architecture

Pixel values of the images were scaled to the interval [0,1]. Each pixel is assumed to be i.i.d. Gaussian
distributed. For the encoders and decoders we use multilayer perceptron networks, see table below.
The numbers corresponds to the size of the layer and the parenthesis states the used activation
function and whether or not batch normalization was used. D indicates the size of the images i.e.
17µencoder2
σencoderµdecoder2σdecoderLayer 1512 (BN + ReLU)
512 (BN + ReLU)
128 (BN + ReLU)
128 (BN + ReLU)

Layer 2256 (BN + ReLU)
256 (BN + ReLU)
256 (BN + ReLU)
256 (BN + ReLU)

Layer 3128 (BN + ReLU)
128 (BN + ReLU)
512 (BN + ReLU)
512 (BN + ReLU)

Layer 4d (BN + Linear)d (Softplus)
D (ReLU)
D (Softplus)

D = width × height × channels. For MNIST and FashionMNIST these are 28,28,1 and for
CIFAR10 and SVHN these are 32,32,1. d indicates the size of the latent space. For MNIST and
FashionMNIST we set d = 2 for visualization purpose and for CIFAR10 and SVHN we set d = 10
to be able to capture the higher complexity of these datasets.
To train the networks we used the Adam optimizer [Kingma and Ba, 2015] with learning rate 10−3
and a batch size of 512. We train for 20000 iterations without early stopping. Additionally, we useitwarm-up for the KL term [Sønderby et al., 2016], by scaling it with w = min 1, warmup where itis the current iteration number and warmup was set to half the number of iterations. This securesthat we converge to stable reconstructions before introducing too much structure in the latent space.

F

Artificial Data

The data considered in Section 4 are generated in the following way: first we sample points in R2 in
a two-moon type way. See details in Algorithm 2. We generate 500 points in this way to establish a’known’ latent space. We then map these to four dimensions (v1 , v2 , v3 , v4 ) byp
v1 (z1 , z2 ) = z1 − z2 +  · 0.03 + 0.05 · (3 + z1 ),
(7)
p
1
(8)
v2 (z1 , z2 ) = z12 + z2 +  · 0.03 + 0.03 · kzk2 ,2pv3 (z1 , z2 ) = z1 z2 − z1 +  · 0.03 + 0.05 · kzk2 ,
(9)
s
0.03
v4 (z1 , z2 ) = z1 + z2 +  · 0.03 +,(10)
0.2 + kzk2
where all  ∼ N (0, 1) and independent. A typical dataset from this procedure is shown in Figure 18.
Algorithm 2 Two moon sampler
1: Sample U ∼ Bernoulli(0.5).
2: if U = 1 then3:
Set c = (0.5, 0) and sample α1 ∼ unif[0, π].
4:Let z = c + (cos(α1 ), sin(α1 )), and sample α2 ∼ unif[0, 2π] and u ∼ unif[0, 1].
5:
Let z = z + u4 · (cos(α2 ), sin(α2 )).
6: else7:
Set c = (−0.5, 0) and sample α1 ∼ unif[π, 2π].
8:
Let z = c + (cos(α1 ), sin(α1 )), and sample α2 ∼ unif[0, 2π] and u ∼ unif[0, 1].
9:
Let z = z + u4 · (cos(α2 ), sin(α2 )).

return z

18

0

1

2

3

−2

−1

0

1

2

1

2

−1

0

2

0.5

3

−2

−1

1.0

V1

1

2

0.0

−1

0

1

V2

V4

−2

−1

0

−1.0

1

2

−2

−1

−0.5

0

V3

−1.5

−1.0

−0.5

0.0

0.5

1.0

1.5
−2

−1

0

1

2

−2

−1

0

1

2

Figure 18: An example of the two moon data, and its transformation into R4 . Shown as pairwise
scatterplots.

19