Year: 2021

Abstract:
Deep Ensemble is a flexible and effective alternative to Bayesian neural networks for uncertainty estimation in deep learning. However, Deep Ensemble is broadly criticized for lacking a proper Bayesian justification. Some attempts try to fix this issue, while they are typically coupled with the regression likelihood or rely on restrictive assumptions. In this work, we propose to define a Gaussian process (GP) approximate posterior with Deep Ensemble, based on which we perform variational inference directly in function space. We further develop a function-space posterior regularization mechanism to properly incorporate prior knowledge. We provide strategies to make the training feasible, and demonstrate the algorithmic benefits of variational inference in the GP family. As a result, our method consumes only marginally added training cost than the standard Deep Ensemble. Empirically, our approach achieves better uncertainty estimation than Deep Ensemble and its variants across diverse scenarios.

heading: 1 INTRODUCTION
text: Bayesian treatment of deep neural networks (DNNs) is promised to enjoy principled Bayesian uncertainty while unleashing the capacity of DNNs, with Bayesian neural networks (BNNs) as popular examples (MacKay, 1992; Hinton & Van Camp, 1993; Neal, 1995; Graves, 2011). Nevertheless, despite the surge of advance in BNNs (Louizos & Welling, 2016; Zhang et al., 2018), many of existing BNNs still face challenges in accurate and scalable inference (Sun et al., 2019), and exhibit limitations in uncertainty estimation and out-of-distribution robustness (Ovadia et al., 2019).
Alternatively, Deep Ensemble (DE) (Lakshminarayanan et al., 2017) trains multiple independent, randomly-initialized DNNs for ensemble, presenting higher flexibility and effectiveness than BNNs. However, it is hard to interpret DE as a Bayesian approach which seeks for the Bayesian posterior of a certain model, and there is no guarantee that the uncertainty estimates given by DE are reliable. To chase a Bayesian justification for DE, recent works like RMS (Lu & Van Roy, 2017; Osband et al., 2018; Pearce et al., 2020) and NTKGP (He et al., 2020) refine DE to be a “sample-then-optimize” approach (Matthews et al., 2017). Nonetheless, these works rely on the regression likelihood and often make strong assumptions like linearised and/or infinite-width models.
To build a Bayesian refinement of DE without reliance on restrictive assumptions, we propose to use DE to define a Gaussian process (GP) approximate posterior, which is dubbed as DE-GP for short. The imposition of a GP form on the approximate posterior is empowered by the equivalence between GPs and BNNs (Khan et al., 2019; Neal, 1996; Lee et al., 2018; Garriga-Alonso et al., 2018; Matthews et al., 2018; Novak et al., 2018; He et al., 2020), as well as the associated computational benefits when adopting a GP prior without loss of generality.
Concretely, we evaluate the empirical mean and covariance of the ensemble members and use them to specify DE-GP. This is inspired by the formula of neural network Gaussian process (NN-GP) (Novak et al., 2018). Nevertheless, we clarify that we define a parameteric and learnable GP kernel with finitely many and adaptive basis functions, in contrast to the NN-GP kernels built with infinitely many yet fixed basis functions. The flexibility of DE-GP is assured by the expressiveness of DNN ensemble members – by tuning them it is enabled to evolve over a rich variety of kernels.
Given these setups, we perform function-space variational inference (VI) for learning. We maximize the functional evidence lower bound (fELBO) (Sun et al., 2019; Rudner et al., 2021) to push DE-GP towards the true posterior over functions. We further present a posterior regularization scheme in function space to conveniently incorporate structural priors, and present a concrete example on how
to control function complexity. By VI, we are enabled to handle classification problems directly and exactly, without casting them into regression ones (He et al., 2020). In practice, DE-GP can be easily implemented on top of a standard DE, while only introducing minimally added cost.
The function-space Bayesian inference enables DE-GP to ameliorate the pathologies induced by the over-parameterization nature of DNNs (Sun et al., 2019; Wang et al., 2019). Thus DE-GP can offer calibrated uncertainty estimates, especially when adopting large networks (see Fig. 1). We empirically demonstrate that DE-GP outperforms strong baselines on various regression datasets, and presents superior uncertainty estimates and out-of-distribution robustness without compromising accuracy in challenging classification tasks. DE-GP also shows promise in solving contextual bandit problems, where the uncertainty plays a vital role in guiding exploration.

heading: 2 RELATED WORK
text: Bayesian neural networks (BNNs). Bayesian treatment of DNNs is an emerging topic in deep learning yet with a long history (Mackay, 1992; Hinton & Van Camp, 1993; Neal, 1995; Graves, 2011). BNNs can be learned by virtue of variational inference (Blundell et al., 2015; HernándezLobato & Adams, 2015; Louizos & Welling, 2016; Zhang et al., 2018; Khan et al., 2018; Deng et al., 2020), Laplace approximation (Mackay, 1992; Ritter et al., 2018), Markov chain Monte Carlo (Welling & Teh, 2011; Chen et al., 2014; Zhang et al., 2019), particle-optimization based variational inference (Liu & Wang, 2016), Monte Carlo dropout (Gal & Ghahramani, 2016), and other methods (Maddox et al., 2019; Izmailov et al., 2020). To avoid the difficulties of posterior inference in weight space, some recent works advocate performing Bayesian reasoning in function space (Sun et al., 2019; Rudner et al., 2021; Wang et al., 2019). In function space, BNNs of infinite or even finitely width equal to Gaussian processes (GPs) (Neal, 1996; Lee et al., 2018; Novak et al., 2018; He et al., 2020; Khan et al., 2019), which provides valuable insights for this work.
Deep Ensemble (DE) (Lakshminarayanan et al., 2017) is a qualified alternative to BNNs for uncertainty estimation (Ovadia et al., 2019), yet lacks a proper Bayesian interpretation. Wilson & Izmailov (2020) conceptually interpreted DE as a method that approximates the Bayesian posterior predictive distribution, but it is hard to judge whether the approximation is reliable or not in practice. Recently, a suite of works based on the notion of “sample-then-optimize” (Matthews et al., 2017) have been developed to make DE Bayesian. For example, RMS (Lu & Van Roy, 2017; Osband et al., 2018; Pearce et al., 2020) regularizes the ensemble members towards randomised priors to obtain posterior samples after training, while it typically assumes linear data likelihood which is impractical for deep models and classification tasks. He et al. (2020) proposed to add a randomised function to each ensemble member to realise a function-space Bayesian interpretation, but the method is asymptotically exact in the infinite width limit and is limited to regressions. By contrast, DE-GP works without restrictive assumptions. In parallel, D’Angelo & Fortuin (2021) proposed to add a repulsive term to DE to achieve Bayesian inference, which yet relies on unscalable gradient estimators.

heading: 3 BACKGROUND

heading: 3.1 BAYESIAN NEURAL NETWORKS
text: Consider learning on a dataset D = (X,Y) = {(xi,yi)}ni=1, with xi ∈ X as data and yi as C-dimensional targets. Let g(·,w) : X → RC denote the function represented by a DNN with weights w. Bayesian neural networks (BNNs) treat w as random variables and chase the posterior p(w|D) by imposing a prior p(w). Among the family, variational BNNs (VBNNs) are particularly appealing as they cast posterior inference as an optimization problem (Graves, 2011; Blundell et al., 2015; Louizos & Welling, 2016), where a variational distribution q(w) is introduced to approximate the true posterior by maximizing the evidence lower bound (ELBO):
Eq(w)[log p(D|w)]−DKL[q(w)‖p(w)]. (1)
BNNs marginalize over the posteriors to obtain the posterior predictive for new data x∗:
p(y|x∗,D) = Ep(w|D)p(y|x∗,w) ≈ Eq(w)p(y|x∗,w) ≈ S∑ s=1 p(y|x∗,ws), (2)
where ws ∼ q(w), s = 1, ..., S. Such a procedure is expected to propagate the embedded model uncertainty into the prediction. However in practice, most of the existing BNN approaches face obstacles in precise posterior inference due to non-trivial and convoluted posterior dependencies (Louizos & Welling, 2016; Zhang et al., 2018; Shi et al., 2018; Sun et al., 2019), and deliver unsatisfactory uncertainty estimation and out-of-distribution (OOD) robustness (Ovadia et al., 2019).

heading: 3.2 DEEP ENSEMBLE
text: As a workaround of BNNs, Deep Ensemble (DE) (Lakshminarayanan et al., 2017) deploys a set of M DNNs {g(·,wi)}Mi=1 for ensemble. {wi}Mi=1 are independently trained to interpret the data from different angles, under maximum likelihood estimation (MLE) principle in standard DE, or maximum a posteriori (MAP) principle in regularized DE (rDE):
max w1,...,wM
LDE = 1
M M∑ i=1 log p(D|wi), max w1,...,wM LrDE = 1 M M∑ i=1 [log p(D|wi) + log p(wi)]. (3)
Due to the randomness in network initialization and stochastic gradient descent (SGD), DE can effectively explore the non-convex, multimodal loss landscape of DNNs (Fort et al., 2019; Wilson & Izmailov, 2020), thus has shown promise in various uncertainty estimation tasks (Lakshminarayanan et al., 2017; Ovadia et al., 2019). However, due to the lack of a Bayesian justification, there is no guarantee that the uncertainty estimates given by DE are reliable.

heading: 3.3 RELATION OF THE BAYESIAN POSTERIORS OF DNNS TO GAUSSIAN PROCESSES
text: It is early shown that shallow BNNs converge to Gaussian processes (GPs) in the infinite width limit (Neal, 1996). Recently, there is ongoing effort to extend the result to multiple layers (Lee et al., 2018; Matthews et al., 2018), convolutional architectures (Garriga-Alonso et al., 2018; Novak et al., 2018), and beyond (Yang, 2019). Briefly, an infinitely wide BNN amount to the NN-GP:
GP(0, κ(x,x′)),with κ(x,x′) := Ew∼p(w)[g(x,w)g(x′,w)>], (4) where p(w) is a layerwise isotropic Gaussian and the DNN specifying g has infinite width. Notably, κ is a compositional, matrix-valued kernel, with values in the space of C × C matrices, and can be analytically estimated in some cases (Novak et al., 2018; 2019). Due to the i.i.d. weights, the kernel evaluation κ(x∗,x∗) of some data x∗ is a scaled identity matrix.
Based on the canonical results of GP regression, we know the posteriors of infinitely wide BNNs are still GPs (in regression tasks), corresponding to the Bayesian training of the last readout layer (Arora et al., 2019; Lee et al., 2019). Despite a golden standard for small regression problems, analytical posterior inference for NN-GP faces challenges when handling contemporary DNN architectures (Novak et al., 2018), large data (Shi et al., 2019) and classification problems.
Interestingly, finitely wide BNNs with Gaussian posteriors on weights can also be connected to GPs (Khan et al., 2019; Rudner et al., 2021). We then arrive at the hypothesis that the Bayesian posteriors of DNNs are closely related to GPs. So, we specify a GP approximate posterior with DE.

heading: 4 METHODOLOGY

heading: 4.1 DEEP ENSEMBLE AS A GAUSSIAN PROCESS APPROXIMATE POSTERIOR
text: Basically, we utilize DE’s ensemble members (referred to as basis functions hereinafter) to define a GP approximate posterior as follows:
q(f ;w1, ...,wM ) := GP(f |mq(x), kq(x,x′)),
mq(x) := 1 M ∑M i=1 gi(x), (5)
kq(x,x ′) := 1
M ∑M i=1 (gi(x)−mq(x)) (gi(x ′)−mq(x′))> + λIC ,
where gi refers to g(·,wi). Namely, the empirical mean and central covariance of finitely many, adaptive basis functions are leveraged to specify DE-GP. kq(x,x′) is a linear, matrix-valued kernel with a small scaled identity matrix λIC appended to avoid singularity. The variations in kq(x,x′) are confined to having up to M − 1 rank, echoing the recent investigations showing that low-rank approximate posteriors for DNNs conjoin effectiveness and efficiency (Maddox et al., 2019; Izmailov et al., 2020; Dusenberry et al., 2020). Note that the popular Nyström method (Williams & Seeger, 2001) also uses a low-rank matrix to approximate the original kernel matrix.
Akin to the kernels in (Wilson et al., 2016), the DE-GP kernels are highly flexible, and may automatically discover the underlying structures of high-dimensional data without manual participation.

heading: 4.2 VARIATIONAL INFERENCE IN FUNCTION SPACE
text: We do function-space variational inference to push DE-GP towards the true posterior over functions.
Prior We have high freedom to determine which prior to use attributed to the variational inference paradigm. Without loss of generality, we use the MC estimate of NN-GP (MC NN-GP) (Novak et al., 2018) as the prior due to its accessibility. Concretely, supposing a finitely wide DNN composed of a feature projector h(·,w) : X → RĈ and a linear readout layer with weight variance σ2w and bias variance σ2b , the MC NN-GP is defined as p(f) = GP(f |0, k(x,x′)), where
k(x,x′) = (σ2wk̂(x,x ′) + σ2b )IC ,with k̂(x,x ′) = 1
SĈ S∑ s=1 h(x,ws) >h(x′,ws). (6)
IC refers to the indentity matrix of size C × C and ws are i.i.d. samples from the Gaussian prior on weights p(w). Similar setups can be found in some related works like (Wang et al., 2019).
There may be a subtle difference between the MC NN-GP and the NN-GP, so what we are actually doing is not approximating the exact NN-GP posterior with DE-GP. Also of note that the NN-GP prior and the DE-GP posterior can be defined with various architectures.
fELBO Following Sun et al. (2019), we maximize the functional ELBO (fELBO) to optimize DE-GP:
max q(f)
Eq(f)[log p(D|f)]−DKL[q(f)‖p(f)]. (7)
Notably, there is a KL divergence between two GPs, which, on its own, is challenging to cope with. Fortunately, as proved by Sun et al. (2019), we can take the KL divergence between the marginal distributions of function evaluations as a substitute for it, giving rise to a more tractable objective:
max q(f) L = ∑ (xi,yi)∈D Eq(f)[log p(yi|f(xi))]−DKL[q(f X̃)‖p(f X̃)], (8)
where X̃ denotes a measurement set including all training inputs X, and f X̃ is the concatenation of the vectorized outputs of f for X̃, i.e., f X̃ ∈ R|X̃|C .

heading: 4.3 POSTERIOR REGULARIZATION IN FUNCTION SPACE
text: Though Bayesian learning in function space enables a direct imposition of prior knowledge on function properties like smoothness and periodicity (Sun et al., 2019), sometimes the prior knowledge cannot be trivially designed as a prior distribution over functions. For example, the functions are
periodical only over a limited input range, etc. In these cases, posterior regularization (Ganchev et al., 2010; Zhu et al., 2014) provides a principled workaround to impose structural constraints.
Briefly, we apply posterior regularization to functional variational inference by solving: max q(f) L s.t.: q(f) ∈ Q. (9)
Q = {q(f)|Eq(f)Ω(f) ≤ 0} is a valid set defined in terms of a functional Ω which delivers some statistic of interest of a function.1 For tractable optimization, we slack the constraint as a penalty:
max q(f) L′ = L − βmax{Eq(f)Ω(f), 0}, (10)
where β is a trade-off coefficient.
We next present an example on how to impose prior knowledge on function (hypothesis) complexity on the DE-GP approximate posterior using this paradigm.2 For the sake of brevity, we assume a binary classification scenario in the following discussion where y ∈ {−1, 1} and f, gi : X → R. We use 0-1 loss `(f(x), y) = 1y 6=sign(f(x)) to measure the classification error on one datum. We assume an underlying distribution µ = µ(x, y) supported on X × {−1, 1} for generating the training data D, based on which we can define the true risk of a function (hypothesis) f : R(f) := E(x,y)∼µ`(f(x), y). We set Eq(f)Ω(f) := Eq(f)R(f) in the seek of a posterior over functions that can generalize well.
By definition, a hypothesis sample f ∼ q(f) = GP(mq(x), kq(x,x′)) can be decomposed as f(x) = 1M ∑M i=1 gi(x) + ζ(x) with ζ(x) ∼ GP(0, kq(x,x′)). If sign(f(x)) 6= y, it is impossible that sign(g1(x)) = y, ..., sign(gM (x)) = y, and sign(ζ(x)) = y all hold. In other words,
`(f(x), y) ≤ M∑ i=1 [`(gi(x), y)] + `(ζ(x), y). (11)
We can further re-parameterize ζ(x) as ζ(x) = 1√ M ∑M i=1 i(gi(x) − mq(x)) + √ λ 0, i ∼ N (0, 1), i = 0, ...,M , which is essentially a real-valued random function symmetric around 0. Thus, for any (x, y) ∼ µ, we have Eq(f)`(ζ(x), y) = E 0,..., M `(ζ(x), y) = 1/2. As a result,
Eq(f)R(f) ≤ Eq(f)E(x,y)∼µ M∑ i=1 [`(gi(x), y)] + E(x,y)∼µ[1/2] = M∑ i=1 [R(gi)] + 1/2. (12)
Namely, the expected generalization error of the approximately posteriori functions can be bounded from above by those of the DNN basis functions.3 Recalling the theoretical and empirical results showing that DNNs’ generalization error R(gi) can be decreased by controlling model capacity in terms of norm-based regularization minwi ||wi||22 (Neyshabur et al., 2015; 2017; Bartlett et al., 2017; Jiang et al., 2019), we opt to maximize the following refined fELBO for learning DE-GP:
max w1,...,wM
LDE-GP = ∑
(xi,yi)∈D
Eq(f)[log p(yi|f(xi))]−DKL[q(f X̃)‖p(f X̃)]− β ∑ i ||wi||22. (13)
0 25 50 75 100 125 150 175 200 Epoch
5000
10000
15000
20000
25000
L 2 n
or m
o f w
ei gh ts DE-GP ( = 0.1) DE-GP ( = 0) DE rDE
Figure 2: Comparison on the L2 norm of weights. We traine models on CIFAR-10 with ResNet-20. DE-GP (β = 0) finds solutions with high complexity and poor test accuracy (see Table 1), yet DE-GP (β = 0.1) settles this.
Intuitively, the weight space of an over-parameterized DNN contains many redundant interpretations for the data, thus pushing the weights towards zero would not hurt data fitting but boost generalization. In practice, when processing with shallow architectures, β can be set as 0 to remove the constraint on model capacity; while when using deep architectures (e.g., ResNet (He et al., 2016)), we can set β according to commonly used weight decay coefficient for better generalization. Comparisons in Fig. 2 and Table 1 support this.

heading: 4.4 THE ALGORITHM
text: We outline the training procedure of DE-GP in Algorithm 1, and elaborate some details below. 1Here we assume one-dimensional outputs for Ω for notation compactness. 2The motivation is that, as observed, the fELBO cannot cause proper regularization effects on weightswi for deep architectures, thus the learned DE-GP may suffer from high-complexity basis functions. 3A similar conclusion can be drawn in multi-class classification scenario (see Appendix A).
Algorithm 1: Deep Ensemble as a Gaussian process posterior (DE-GP). 1: Input: D: dataset; {gi}Mi=1: M DNNs; k: prior kernel; α, β: trade-off coefficients; ν: distribution for
sampling extra measurement points; U : number of MC samples for estimating the expected log-likelihood 2: while not converged do 3: Ds = (Xs,Ys) ⊂ D, Xν ∼ ν, X̃s = {Xs,Xν} 4: gi = gi(X̃s) ∈ R|X̃s|C , i = 1, ...,M 5: mX̃sq = 1M ∑ i gi, k X̃s,X̃s q = 1 M ∑M i=1(gi −m X̃s q )(gi −mX̃sq )> + λI|X̃s|C 6: kX̃s,X̃s = k(X̃s, X̃s) ∈ R|Xs|C×|Xs|C
7: L1 = 1U ∑U i=1 ∑ (x,y)∈Ds log p(y|fi(x)), fi ∼ N (m X̃s q ,k X̃s,X̃s q )
8: L2 = DKL[N (mX̃sq ,kX̃s,X̃sq )‖N (0,kX̃s,X̃s)], L3 = ∑ i ||wi|| 2 2 9: wi = wi + η∇wi(L1 − αL2 − βL3), i = 1, ...,M

heading: 4.4.1 MINI-BATCH TRAINING
text: In deep learning scenarios, DE-GP should proceed by mini-batch training. At each step, we manufacture a stochastic measurement set with a mini-batch Ds = (Xs,Ys) from the training data D and random samples Xν from a continuous distribution ν (e.g., a uniform distribution) supported on X . Then, we adapt the objective defined in Eq. (13) to the following form:
LDE-GP = ∑
(xi,yi)∈Ds
Eq(f)[log p(yi|f(xi))]− αDKL[q(f X̃s)‖p(f X̃s)]− β ∑ i ||wi||22, (14)
where X̃s indicates the union of Xs and Xν . When α = 1 and β = 0, LDE-GP is a lower bound of log p(Ds) according to (Sun et al., 2019). Yet, given that we do not perform intensive hyper-parameter tuning for the NN-GP prior, it is reasonable to make α a tunable hyper-parameter to better trade off between the data fitting and the priori inductive bias.

heading: 4.4.2 AN EXACT AND EFFICIENT ESTIMATION OF THE MARGINAL KL DIVERGENCE
text: The marginal distributions of function evaluations are multivariate Gaussian by the definition of GP:
q(f X̃s) = N (f X̃s |mX̃sq ,kX̃s,X̃sq ), p(f X̃s) = N (f X̃s |0,kX̃s,X̃s), (15)
with the kernel matrices kX̃s,X̃sq ,k X̃s,X̃s ∈ R|X̃s|C×|X̃s|C as the joints of pair-wise outcomes.
Therefore, the marginal KL divergence and its gradients can be estimated exactly without resorting to some approximations (Rudner et al., 2021). What’s more, as discussed in Section 3.3, there is a simple structure in the prior kernel matrices, so we can write them in the form of Kronecker product:
kX̃s,X̃s ≈ (σ2wk̂X̃s,X̃s + σ2b )⊗ IC , (16)
where k̂X̃s,X̃s ∈ R|X̃s|×|X̃s| corresponds to the evaluation of kernel k̂. Hence we can exploit the property of Kronecker product to inverse kX̃s,X̃s in O(|X̃s|3) complexity.
Besides, as kX̃s,X̃sq is low-rank, we can leverage the Woodbury matrix identity (Woodbury, 1950) and matrix determinant lemma (Harville, 1998) to efficiently compute the inverse and determinant of kX̃s,X̃sq in O(|X̃s|CM2) time complexity given that usually M |X̃s|C (e.g., 10 256C).

heading: 4.5 DISCUSSION
text: Diversity. The diversity among the ensemble members in function space is explicitly encouraged by the KL divergence between variational and the prior in fELBO. Nonetheless, the expected loglikelihood in fELBO enforces each ensemble member to yield the same, correct outcomes for the training data. Thereby, the diversity mainly exists in the regions far away from the training data (see Fig. 1 and Fig. 9 in Appendix B.2). Yet, the diversity in DE does not have a clear theoretical support.
Efficiency. Compared to the overhead introduced by DNNs, the effort of estimating the KL divergence between Gaussians is negligible. The added cost of DE-GP primarily arises from the extra measurement points and the evaluation of the prior kernels. In practice, we use a small batch size
for the extra measurement points. We build the prior NN-GP kernels with cheap architectures and perform MC estimation in parallel. Eventually, DE-GP is only marginally slower than DE in training.
Weight sharing. DE-GP does not care about how the basis functions are parameterized, so we can perform weight sharing among the basis functions, for example, using a shared feature extractor and M independent MLP classifiers to construct M basis functions (Deng et al., 2021). With shared weights, DE-GP is still likely to be reliable because our learning principle induces diversity in function space regardless of the weights. Experiments in Section 5.3 validate this.
Limitations. Despite being Bayesian in principle, DE-GP is likely to be less flexible than DE and other variants because all the involved basis functions need to be updated simultaneously.

heading: 5 EXPERIMENTS
text: We perform extensive evaluation to prove that DE-GP yields better uncertainty estimates than the baselines, while preserving non-degraded predictive performance. The baselines include DE, rDE, NN-GP, RMS, etc. In all experiments, we estimate the prior kernel with 10 MC samples and set the sampling distribution for extra measurement points ν as the uniform distribution over the data region.

heading: 5.1 ILLUSTRATIVE REGRESSION
text: We build two regression problems with data from y = x3/4 + , ∼ N (0, 0.1) and y = x sin 5x+ , ∼ N (0, 0.2), respectively. We consider two architectures: MLPs with 2 hidden layers of 64 units and MLPs with 3 hidden layers of 512 units, where ReLU activation is used. For NN-GP, we analytically estimate the GP kernel and perform GP regression without training MLPs. For DE-GP, DE, rDE, and RMS, we train 50 MLPs. For MC dropout, we train a MLP with the commonly used 0.2 dropout rate. We use the default values for α and β, i.e., α = 1 and β = 0.
Fig. 1 and Appendix B.2 show the results on the two problems. We also provide the comparison on training speed in Fig. 1. As shown, DE-GP delivers calibrated uncertainty estimates across settings, on par with the non-parametric Bayesian baseline NN-GP. Yet, the baselines like DE, rDE, and RMS suffer from degeneracy issue as the dimension of weights increases. Though NN-GP outperforms other methods, the involved analytical GP regression may have scalability and effectiveness issues when facing modern architectures (Novak et al., 2018), while DE-GP does not suffer from them.

heading: 5.2 UCI REGRESSION
text: We then assess DE-GP on 5 UCI real-valued regression problems. The used architecture is a MLP with 2 hidden layers of 256 units and ReLU activation. 10 networks are trained for DE, DE-GP and other variants. For DE-GP, we set β = 0 and tune α according to validation sets.
We perform cross validation with 5 splits. Fig. 3 shows the results. DE-GP surpasses or approaches the baselines across scenarios in aspects of both test negative log-likelihood (NLL) and test root mean square error (RMSE). DE-GP even beats NN-GP, which is probably attributed to that the variational family specified by DE enjoys the beneficial inductive bias of practically sized SGD-trained DNNs, and DE-GP can flexibly trade off between the likelihood and the prior by tuning α.

heading: 5.3 IMAGE CLASSIFICATION ON FASHION-MNIST AND CIFAR-10
text: In the classification experiments, we augment the data log-likelihood (i.e., the first term in Eq. (14)) with a trainable temperature to tackle oversmoothing and avoid underconfidence. 4
Fashion-MNIST. We use a widened LeNet5 architecture with batch normalizations (BNs) (Ioffe & Szegedy, 2015) for the Fashion-MNIST dataset (Xiao et al., 2017). Considering the inefficiency of
4We can make the temperature a Bayesian variable, but it is unnecessary as our model is already Bayesian.
NN-GP, we mainly compare DE-GP to DE, rDE, and RMS. We set α as well as the regularization coefficients for rDE and RMS all as 0.1 according to validation accuracy. For DE-GP, we use β = 0 given the limited capacity of the architecture. The in-distribution performance is averaged over 3 runs. Fig. 4-(Left) and Fig. 4-(Middle) display how ensemble size impacts the test results. Surprisingly, the test error of DE-GP is even lower than the baselines.
Besides, to compare the quality of uncertainty estimates, we use the trained models to make prediction and quantify epistemic uncertainty for both the in-distribution test set and the out-of-distribution (OOD) MNIST test set. All predictions on OOD data are regarded as wrong. The epistemic uncertainty is estimated by the mutual information between the prediction and the variable function:
I(f, y|x,D) ≈ H
( 1
S S∑ s=1 p(y|fs(x)) ) − 1 S S∑ s=1 H (p(y|fs(x))) , s = 1, ..., S, (17)
where H indicates Shannon entropy, with fs = g(·,ws) for DE, rDE, and RMS, and fs ∼ q(f ;w1, ...,wM ) for DE-GP. We normalize the uncertainty estimates into [0, 1]. For each threshold τ ∈ [0, 1], we plot the average test error for data with ≤ τ uncertainty in Fig. 4-(Right). It is prominent that under various uncertainty thresholds, DE-GP makes fewer mistakes than the baselines, implying that DE-GP succeeds to assign relatively higher uncertainty for the OOD data.
CIFAR-10. Next, we apply DE-GP to the real-world image classification task CIFAR-10 (Krizhevsky et al., 2009). We consider the popular ResNet architectures (He et al., 2016) including ResNet-20 and ResNet-56. The ensemble size is fixed as 10. We split the data as training set, validation set, and test set of size 45000, 5000, 10000, respectively. We set β = 0.1, equivalent with the regularization coefficient on weight in rDE. We set α = 0.1 according to an ablation study in Appendix B.5. We use a lite ResNet-20 architecture without BNs and residual connections to set up the prior NN-GP kernel for both the ResNet-20 and ResNet-56 based variational posteriors.
We present the in-distribution test accuracy in Table 1 and the error versus uncertainty plots on the combination of CIFAR-10 and SVHN test sets in Fig. 5. It is noteworthy that DE-GP is on par with the practically-used, competing rDE in aspect of test accuracy. DE-GP (β = 0) shows unsatisfactory test accuracy, verifying the necessity of performing posterior regularization to penalize function complexity when using ultra-deep networks. The error versus uncertainty plots are similar to those for Fashion-MNIST, substantiating the universality of DE-GP.
We further test the trained methods on CIFAR-10 corruptions (Hendrycks & Dietterich, 2018), a challenging OOD generalization/robustness benchmark for deep models. As shown in Fig. 6 and Appendix B.3, DE-GP reveals smaller Expected Calibration Error (ECE) (Guo et al., 2017) and lower NLL at various levels of skew, reflecting its ability to make conservative predictions under corruptions and good OOD robustness.
More results for the deeper ResNet-110 architecture and the more challenging CIFAR-100 benchmark are provided in Appendix B.3 and Appendix B.4.
Weight Sharing. We build a ResNet-20 with 10 classification heads and a shared feature extraction module to evaluate the methods under weight sharing. We set a larger value for α for DE-GP to induce higher magnitudes of functional diversity. The test accuracy (over 6 trials) and error versus uncertainty plots on CIFAR-10 are illustrated in Fig. 7. We exclude RMS from the comparison as it assumes i.i.d. basis functions which may be incompatible with weight sharing. DE-GP benefits from the diversity in function space, hence performs better than DE and rDE, which purely hinge on the diversity in weight space.
5.4 CONTEXTUAL BANDIT
Finally, we apply DE-GP to contextual bandit, an important decision-making task where the uncertainty helps to guide exploration. Following (Osband et al., 2016), we use DE-GP to achieve efficient exploration inspired by Thompson sampling. We reuse most of the settings for UCI regression. We leverage the GenRL library to build two contextual bandit problems Covertype and Mushroom (Riquelme et al., 2018). The cumulative reward is depicted in Fig. 8. As desired, DE-GP offers better uncertainty estimates and hence beats the baselines by clear margins. The potential of DE-GP in more reinforcement learning and Bayesian optimization scenarios deserves future investigation.

heading: 6 CONCLUSION
text: In this work, we attempt to build a Bayesian refinement of Deep Ensemble by drawing inspiration from the connection between BNNs and GPs. We propose to leverage the ensemble members to specify an adaptive GP approximate posterior, and perform variational inference directly in function space. We inherit some theoretical results from existing works on functional variational inference, and further develop a posterior regularization scheme to conveniently induce prior knowledge on function properties. We demonstrate how to penalize the function complexity based on this scheme and empirically show the necessity of doing so. The whole algorithm can be implemented easily and efficiently. Extensive experiments validate the effectiveness of our method. We hope this work may shed light on the development of better Bayesian deep learning approaches.
ETHICS STATEMENT
This work proposes a Bayesian refinement of Deep Ensemble. Its potential positive impacts in the society are evident: its ability to enable better uncertainty estimation while maintaining predictive performance is crucial in industry, e.g., automatic driving, disease analysis, and financial applications. In this scenarios, the uncertainty estimates could be used to reject uncertain predictions, and raise the requirement of inviting humans into the decision process. As a fundamental research in machine learning, the negative consequences are not obvious. Though in theory any technique can be misused, it is not likely to happen at the current stage.
REPRODUCIBILITY STATEMENT
It is easy to reproduce the experiments given the code in the supplementary material.

heading: A PENALIZE FUNCTION COMPLEXITY IN MULTI-CLASS CLASSIFICATION
text: In the multi-class classification scenario where y ∈ {1, 2, ..., C} and f, gi : X → RC , we use the loss `(f(x), y) = 1f(x)[y]<maxy′ 6=y f(x)[y′] to measure prediction error where f(x)[j] denotes j-th coordinate of f(x). The distinct difference between this scenario and the binary classification scenario is that in this setting, ζ(x) is a vector-valued function:
ζ(x) = 1√ M M∑ i=1 i(gi(x)−mq(x)) + √ λ 0, (18)
where 0 ∼ N (0, IC) and i ∼ N (0, 1), i = 1, ...,M . We then make a mild assumption to simplify the analysis. Assumption 1. For any (x, y) ∈ µ, the elements on the diagonal of kq(x,x) have the same value.
This assumption implies that for any j, j′ ∈ {1, ..., C},
1
M M∑ i=1 (gi(x)[j]−mq(x)[j])2 + λ = 1 M M∑ i=1 (gi(x)[j ′]−mq(x)[j′])2 + λ. (19)
I.e., M∑ i=1 (gi(x)[j]−mq(x)[j])2 = M∑ i=1 (gi(x)[j ′]−mq(x)[j′])2. (20)
Therefore, ζ(x) possesses the same variance across its output coordinates and becomes a random guess classifier. Based on this, we have Eq(f)`(ζ(x), y) = E 0,..., M `(ζ(x), y) = (C − 1)/C. We can then derive a similar conclusion to that in Section 4.3.
The validity of Assumption 1. When the data dimension |X | is high and the number of training data n is finitely large, with zero probability the sampled data (x, y) ∼ µ resides in the training set. Therefore, only the KL divergence term of the fELBO explicitly affects the predictive uncertainty at x. Because the NN-GP prior possesses a diagonal structure, it is hence reasonable to make the above assumption.

heading: B MORE OF EXPERIMENTS
text: We provide more experimental details and results in this section. In all experiments, the number of MC samples for estimating the expected log-likelihood (i.e., U in Line 7 of Algorithm 1) is set as 256. Unless otherwise stated, we set the regularization constant λ as 0.05 times of the average eigenvalue of the central covariance matrices, and set the weight and bias variance for defining the NN-GP prior kernel at each layer as 2/fan_in and 0.01, where fan_in is the number of input features, as suggested by He initialization (He et al., 2015).
B.1 DETAILED SETTINGS
Illustrative regression. For the problem on y = x3/4 + , ∼ N (0, 0.1), we randomly sample 20 data points from [−1.5,−0.6] ∪ [0.6, 1.5]. For the problem on y = x sin 5x+ , ∼ N (0, 0.2), we use 6 data points {−0.8,−0.1, 0.02, 0.2, 0.6, 0.8} following (Pearce et al., 2020). For optimizing the DNN basis functions, we use an Adam (Kingma & Ba, 2015) optimizer with 0.003 learning rate. The optimization takes 300 iterations. We set hidden size as 64 for the architecture of the prior NN-GP kernel for speedup. The regularization constant λ is set as 1e− 4 times of the average eigenvalue of the central covariance matrices.
UCI regression. We pre-process the UCI data by standard normalization. We set the variance for data noise and the weight variance for the prior kernel following (Pearce et al., 2020). The batch size for stochastic training is 256. We use an Adam optimizer to optimize for 1000 epochs. The learning rate is initialized as 0.01 and decays by 0.99 every 5 epochs.
Fashion-MNIST classification. The used architecture is Conv(32, 3, 1)-BN-ReLU-MaxPool(2)Conv(64, 3, 0)-BN-ReLU-MaxPool(2)-Linear(256)-ReLU-Linear(10), where Conv(x, y, z) represents
a 2D convolution with x output channels, kernel size y, and padding z. The batch size for training data is 64. The batch size for extra measurement points is 0. We use an SGD optimizer to optimize for 24 epochs. The learning rate is initialized as 0.1 and follows a cosine decay schedule. We use an Adam optimizer with 1e− 3 learning rate to optimize the temperature. We use 1000 MC samples to estimate the posterior predictive and the epistemic uncertainty, because the involved computation is only the cheap softmax transformation on the sampled function values.
CIFAR-10 classification. We perform data augmentation including random horizontal flip and random crop. The batch size for training data is 128. The batch size for extra measurement points is 0. We use a SGD optimizer with 0.9 momentum to optimize for 200 epochs. The learning rate is initialized as 0.1 and decays by 0.1 at 100-th and 150-th epochs. We use an Adam optimizer with 1e− 3 learning rate to optimize the temperature. We use 1000 MC samples to estimate the posterior predictive and the epistemic uncertainty. Suggested by (Ovadia et al., 2019; He et al., 2020), we train models on CIFAR-10, and test them on the combination of CIFAR-10 and SVHN test sets. This is a standard benchmark for evaluating the uncertainty on OOD data.
Contextual bandit. We use MLPs with 2 hidden layers of 256 units. The batch size for training data is 512. The batch size for extra measurement points is 0. We update the model (i.e., the agent) for 100 epochs with an Adam optimizer every 50 rounds. We set α = 1 and β = 0 without tuning. DE, rDE, and RMS all randomly choose an ensemble member at per iteration, but our method randomly draws a sample from the GP for decision. This is actually emulating Thompson Sampling and advocated by Bootstrapped DQN (Osband et al., 2016). “Random” baseline corresponds to the Uniform algorithm.
B.2 MORE RESULTS ON ILLUSTRATIVE REGRESSION
We depict the results on y = x sin 5x+ , ∼ N (0, 0.2) problem in Fig. 9, which are consistently with those in Fig. 1.
B.3 MORE RESULTS ON CIFAR-10 CLASSIFICATION
We plot the negative log-likelihood and test accuracy on CIFAR-10 corruptions for models trained with ResNet-20 and ResNet-56 in Fig. 10 and Fig. 11. As shown, DE-GP outperforms the baselines in aspect of negative log-likelihood, but yields similar test accuracy to the baselines. Recapping the results in Fig. 6, DE-GP indeed has improved OOD robustness, but may still face problems in OOD generalization.
We then conduct experiments with the deeper ResNet-110 architecture. Due to resource constraint, we use 5 ensemble members. The other settings are roughly the same as those for ResNet-56. The results are offered in Fig. 12, which validate the effectiveness of DE-GP for large networks.
B.4 RESULTS ON CIFAR-100
We further perform experiments on the more challenging CIFAR-100 benchmark. We present the in-distribution test accuracy of DE-GP as well as the baselines in Table 2. We can see that DE-GP (β = 0.1) is still on par with rDE. We depict the error versus uncertainty plots on the combination of CIFAR-100 and SVHN test sets in Fig. 13. It is shown that the uncertainty estimates yielded by DE-GP for OOD data are more calibrated than the baselines. We further test the trained methods on CIFAR-100 corruptions (Hendrycks & Dietterich, 2018), and present the comparisons in aspects of test accuracy and NLL in Fig. 14. It is evident that DE-GP reveals lower NLL than the baselines at various levels of skew.
B.5 ABLATION STUDY ON α
We have conducted an ablation study on α (using ResNet-20 on CIFAR-10). The results are presented in Table 3. We can see that DE-GP is not sensitive to the value of α. We in practice set α = 0.01 in the CIFAR experiments. We did not use a smaller α as it may result in colder posteriors and in turn worse uncertainty estimates.
B.6 ABLATION STUDY ON THE ARCHITECTURE OF PRIOR KERNEL
We perform an ablation study on the architecture for defining the prior MC NN-GP kernel, with the results listed in Table 4. Surprisingly, using the cheap ResNet-20 architecture results in DE-GP with
better test accuracy. We deduce this is because a deeper prior architecture induces more complex, black-box correlation for the function, which may lead to over-regularization.

Review summaries:
The paper presents a novel perspective on deep ensembles (DE) that aims at providing a Bayesian justification to DE. The proposed approach builds a Gaussian process (GP) using the ensemble members and performs variational inference in the functional space. Moreover, the paper introduces a regularization method that works directly in the function space. The proposed algorithm DE-GP shows better results compared to DE an other well known Bayesian deep learning methods. 

The motivation for using fully Bayesian methods over ensemble methods has been a contentious topic in recent years - while ensemble methods are prized for their relative simplicity, they lack the theoretical framework that grounds fully Bayesian approaches. In this work, the authors propose an interpretation of Deep Ensemble models (DE) under a variational Bayesian framework. In particular, the authors demonstrate how VI can be carried out in the function space, and leverage posterior regularisation on functions to incorporate prior knowledge into the model architecture. The benefits of the proposed approach are verified via an extensive evaluation covering a variety of different problem settings, whereby it appears that the DE-GP consistently yields predictions having superior uncertainty calibration, and without compromising on predictive accuracy.

The authors propose using Deep Ensemble (DE) as basis functions to train a Gaussian Process (GP) with the primary motivation of making DEs more Bayesian. The Evidence Lower BOund (ELBO) being maximized is now optimized via variational inference (VI) over the functional space of bases (i.e. fELBO from [3] Sun et al.) defining the mean and covariance of a GP similar to NN-GP [4]. The paper evaluates the proposed method on standard UCI regression and small-scale MNIST / CIFAR10 image classification tasks in terms of accuracy, likelihood, and uncertainty estimation.

This paper presents a GP model with a deep kernel defined in terms of a (finite) deep ensemble (DE). A variational approximation and a regularization scheme were introduced to optimize the GP. The proposed methods were then demonstrated on several benchmark regression/classification datasets.   Here, the claimed contributions are (1) a DE-parameterized GP posterior that provides a Bayesian justification for DE; and (b) a regularization scheme that improves the generalization of the trained GP.