Year: 2021

Abstract:
Ensembles of independently trained neural networks are a state-of-the-art approach to estimate predictive uncertainty in Deep Learning, and can be interpreted as an approximation of the posterior distribution via a mixture of delta functions. The training of ensembles relies on non-convexity of the loss landscape and random initialization of their individual members, making the resulting posterior approximation uncontrolled. This paper proposes a novel and principled method to tackle this limitation, minimizing an f -divergence between the true posterior and a kernel density estimator in a function space. We analyze this objective from a combinatorial point of view, and show that it is submodular with respect to mixture components for any f . Subsequently, we consider the problem of greedy ensemble construction, and from the marginal gain of the total objective, we derive a novel diversity term for ensemble methods. The performance of our approach is demonstrated on computer vision out-of-distribution detection benchmarks in a range of architectures trained on multiple datasets. The source code of our method is made publicly available.

heading: 1 INTRODUCTION
text: Estimation of predictive uncertainty is one of the most important challenges to solve in Deep Learning (DL). Applications in finance, medicine and self-driving cars are examples where reliable uncertainty estimation may help to avoid substantial financial losses, improve patient outcomes, or prevent fatal accidents (Gal, 2016). However, to date, despite rapid progress, there is a lack of principled methods that reliably estimate the predictive uncertainty of deep neural networks (DNNs).
Numerous attempts have been made to develop Bayesian techniques for uncertainty estimation in DL (Gal & Ghahramani, 2016; Lakshminarayanan et al., 2017; Wilson & Izmailov, 2020; Maddox et al., 2019; Izmailov et al., 2020; Van Amersfoort et al., 2020; Wenzel et al., 2020b). One of the most practical and empirically best-performing approaches is based on training a series of independent DNNs (Lakshminarayanan et al., 2017; Wilson & Izmailov, 2020; Ashukha et al., 2020; Lu et al., 2020; Wenzel et al., 2020b). The main method in this category, Deep Ensembles (DE) (Lakshminarayanan et al., 2017), is used as a reference approach in the context of this paper.
Recent studies, e.g. Wilson & Izmailov (2020) interpret ensembles as an approximation of predictive posterior. While this interpretation is correct from a Bayesian point of view, obtaining individual ensemble members via maximum a posteriori probability (MAP) estimation, as e.g. done in DE (Lakshminarayanan et al., 2017), may not lead to obtaining good coverage of the full support of the posterior distribution, and has arbitrary bad approximation guarantees. For example, the resulting approximation can be poor in the case when the true posterior distribution is unimodal. In this work, we argue that enforcing coverage of posterior (i.e. ensemble diversity) is non-trivial, and needs a specialized principled approach. We highlight this graphically in Figure 1.
Another important line of work in modern Bayesian DL (BDL) is a paradigm of performing Bayesian inference in the weight space. While distributions over weights induce distributions over functions (Wilson & Izmailov, 2020), it is rather unclear what the properties of such functions are, and whether Bayesian posteriors obtained in the weight space yield good quality approximations in the function space. For example, it is known that diverse weights do not necessary yield diverse functions, thus sampling from weight-based posteriors may yield poor quality uncertainty estimation, e.g. in detecting out-of-distribution (OOD) data (Garipov et al., 2018; Hafner et al., 2020).
Recent studies (Wang et al., 2019; D’Angelo & Fortuin, 2021) show the promise of particle variational inference (POVI) done in the function space, however, the performance of those methods is not state-of-the-art due to the use of BNN priors. In practice, a weight decay prior combined with batch normalization, data augmentation and dropout are often employed due to empirical improvement, and theoretical guarantees are traded for practical performance. Furthermore, particle-based function space VI requires training an ensemble of BNNs simultaneously, which is not only difficult to implement, but also requires extensive resources for parallelization.
Summary of the contributions. In this paper, we propose a novel and principled methodology for approximate function space posterior inference for DNNs. Contrary to the mainstream approach in BDL, which is based on defining a posterior distribution over the model parameters (Maddox et al., 2019), we take a functional view, which allows us to treat the problem of training ensembles as optimization over sets. Namely, selecting a set of functions to approximate the Bayesian posterior is fundamentally similar to problems such as facility location and set cover, which are frequently solved using submodular optimization (Krause et al., 2008). Specifically, our contributions are:
1. We show that fitting a kernel density estimator to a distribution using an f -divergence is a cardinality-fixed non-monotone submodular maximization problem.
2. Inspired by the Random Greedy algorithm for submodular maximization (Buchbinder et al., 2014), we design a new method for the function space Bayesian posterior approximation via greedy training of ensembles with a justified coverage-promoting diversity term.
3. We demonstrate the effectiveness and competitiveness of our approach compared to DE in the OOD detection task on MNIST, CIFAR, and SVHN, benchmarks on different architectures and ensemble sizes. Furthermore, we show that our method can use state-of-the-art training technqiues, compared to the existing Bayesian approaches, and yeilds state-of-theart performance.

heading: 2 PRELIMINARIES

heading: 2.1 PROBLEM STATEMENT
text: Consider an ensemble to be parameterized by a set of functions Z = {zm}Mm=1 ⊂ F , where F is a class of continuous functions, zm : Rd → Rc, with d the dimensionality of the input data, and c the dimensionality of the output. When training ensembles, we generally want to solve the following optimization problem:
min Z,|Z|=M R(Z)− ΩλM (Z), (1)
whereR(Z) = 1N ∑N i=1 ` ( 1 M ∑M m=1 zm(xi), yi ) is the empirical risk of the ensemble, ` : Y×Y → R+ is a loss function, D = {xi, yi}Ni=1 is a training dataset of size N , and ΩλM (Z) is some diversitypromoting term, with diversity regularization strength λM .
Empirical observations in the earlier works on ensembles (Lakshminarayanan et al., 2017; Wilson & Izmailov, 2020; Fort et al., 2019) have shown that one can simply ignore ΩλM (Z), and rely on
non-convexity of the loss landscape, optimizing the risks of individual ensemble members. From a variational inference (VI) perspective (Zhang et al., 2019), this can be seen as a mode-seeking method, and it has been shown experimentally that every ensemble member may discover different modes of the posterior distribution in the function space p(z|D). Approximation of p(z|D) is the ultimate aim of this paper, and we argue that randomization-based mode-seeking is insufficient to obtain a good quality approximation of p(z|D), as this procedure does not maximize the coverage of the support of the posterior. By contrast, we aim to find an ΩλM (Z) such that the posterior coverage is also maximized. Taking a VI perspective again (Zhang et al., 2019), ΩλM (Z) needs to enforce that minZ,|Z|=M R(Z)− ΩλM (Z) has rather mean-seeking behavior, while still discovering the high density modes.

heading: 2.2 A COMBINATORIAL VIEW OF ENSEMBLE CONSTRUCTION
text: Having now defined the main criteria for (1), we highlight that the problem of constructing an ensemble can be seen from a combinatorial point of view. We therefore treat ensemble construction as subset selection from some ground set of functions, and introduce the main notions of submodular analysis, a powerful tool that enables the analysis of the optimization of set functions (Bach, 2013; Fujishige, 2005).
Definition 1 (Submodularity). A set function g : 2V → R, for the power set of a base set V , is submodular if for all A ⊆ B ⊂ V and x ∈ V \B
g(A ∪ {x})− g(A) ≥ g(B ∪ {x})− g(B). (2)
Definition 2 (Supermodularity and modularity). A set function is called supermodular if its negative is submodular, and modular if it is both submodular and supermodular.
Consider now problem (1). Assuming that the loss function ` is convex, we can derive an upper-bound on the riskR(Z) using Jensen’s inequality, and obtain a method, which generalizes DE
min Z,|Z|=M
1
M M∑ m=1 R(zm)− ΩλM (Z). (3)
If M is fixed during optimization, 1M ∑M m=1R(zm) contributes a positive modular term to the overall objective. Adding a positive modular function to any set function does not change its submodularity or supermodularity, thus we focus on ΩλM (Z). A trivial approach would be to enforce pair-wise diversity by computing a norm of the pairwise differences between functions, i.e. setting ΩλM (Z) = λM ∑ i6=j ‖zi − zj‖2∗. However, this is a cardinality-fixed submodular minimization problem. It is known that it is strongly NP-hard, i.e. there exists no general polynomial time approximation algorithm for it (Svitkina & Fleischer, 2011). The poor quality approximation of this approach is highlighted in Figure 1. We therefore conclude that the choice of ΩλM (Z) has a direct impact on the approximabilty of the objective.
3 SUBMODULAR ANALYSIS OF f -DIVERGENCES
3.1 f -DIVERGENCES ARE SUPERMODULAR FUNCTIONS
Main result. We now consider the problem of approximating a Bayesian posterior via minimization of an f -divergence. Here, we specifically aim our optimization procedure to have both mode and mean-seeking behaviors, i.e. cover the posterior distribution as much as possible, ending up in its mode. We furthermore aim to obtain a polynomial time algorithm that yields a good quality approximation guarantee. In this paper, we leverage classic definitions of approximation algorithm and approximation guarantees.
Definition 3 (Approximation algorithm and guarantees (Williamson & Shmoys, 2011)). A γapproximation algorithm for an optimization problem is a polynomial-time algorithm that for all instances of the problem produces a solution whose value is within a factor γ of the optimal solution. γ in this case is called approximation guarantee.
Let us now formally introduce f -divergences. Definition 4 (f -divergence). Let f : R+ → R be a convex function such that f(1) = 0. The f -divergence between distributions Pz and Qz (Pz being absolutely continuous with respect to Qz), admitting densities p(z) and q(z) is defined as
Df (Pz||Qz) = ∫ f ( p(z)
q(z)
) q(z)dz. (4)
Consider some density p(z) over continuous functions. We define qM (z) = 1M ∑M m=1K(d(z, zj)), Kj(z) := K(z, zj) is a kernel centered at zj used to approximate the modes of p(z). Theorem 1. Any f -divergence
Df (p||qM ) = ∫ f ( p(z)
1 M ∑M j=1Kj(z)
) 1
M M∑ m=1 Km(z)dz (5)
between a distribution p(z) and a normalized mixture of M kernels with equal weights is supermodualar in a cardinality-fixed setting, assuming that ∀zmaxqM Df (p(z)||qM (z)) <∞.
Proof. The proof is shown in Appendix A.1.
Minimization of (5) is equivalent to a cardinality-constrained maximization of a non-monotone submodular function of Z = {z1, . . . , zM}. Approximation guarantees for problems of this form are given for non-negative submodular functions (Buchbinder et al., 2014).
One can convert (5) to a non-negative function by defining: F (Z) := −Df (p||qM ) + C, (6)
where C = maxZ D(p||qM ) is a pre-defined constant, which is important for understanding approximation guarantees, but does not need to be computed in practice. After the described transformation, which leads to (6), we obtain F (Z), a non-negative non-monotone submodular function.
Approximation guarantees for f -divergences. In the context of submodular functions, there exists an inapproximability result, obtained in (Gharan & Vondrák, 2011), which states that no general polynomial time algorithm with guarantees better than 0.491 exists to solve a submodular maximization problem with constrained cardinality. We, however, that in practice, it might still be possible to obtain a better approximation factor than 0.491 for some specific types of divergences, as these guarantees are defined for all instances of the optimization problem (Definition 3).
Let us consider the approximation guarantees for (6), denoting by q∗M the optimal solution, and by q̂M a solution found by some algorithm. For an approximation factor γ, we have −Df (p||q̂M ) + C ≥ γ(−Df (p||q∗M ) + C). (7) Simple algebra shows that this implies
Df (p||q̂M ) ≤ γ min qM∈F Df (p||qM ) + (1− γ) max qM∈F Df (p||qM ), (8)
where F is a family of distributions over the spaces of candidate functions. The derived result indicates that the upper bound on the approximate solution found by minimizing an f -divergence can be substantially dominated by (1− γ) maxqM∈F Df (p||qM ) if F is chosen poorly. The impact of the choice of function space, determined largely by the neural network architecture, is therefore unavoidable when designing algorithms for approximating the Bayesian posterior.
3.2 GREEDY MINIMIZATION OF f -DIVERGENCES
Random greedy algorithm. Although submodular optimization has natural parallel extensions and associated approximation guarantees, due to the simplicity of presentation, we focus in this paper on forward greedy selection and use Algorithm 1 for optimizing submodular functions. This algorithm has approximation guarantee of 1/e in general. The only required step for this greedy algorithm is a computation of the marginal gain on the objective function F (Z), i.e. ∆(zk|Z) = F (Z ∪ {zk}) − F (Z). We argue that due to randomized initialization in neural networks, one can use a naive greedy algorithm instead of the sampling procedure in lines 7-8, and we show the implementation of the resulting training procedure in Appendix B.1.
Algorithm 1 Random Greedy algorithm 1: Input: V – Ground set 2: Input: F – Arbitrary submodular function 3: Input: M – Cardinality of the solution 4: Z ← ∅ 5: for m = 1 to M do 6: R← arg maxT⊂V \Z:|T |=M ∑ z′∈T ∆(z
′|Z) 7: ui ← Uniform(R) 8: Z ← Z ∪ {ui} 9: end for
10: return Z
Marginal gain. At each step of a greedy algorithm, a marginal gain ∆(zk|Z) = F (Z ∪ {zk})− F (Z) of adding a new element zk to an existing mixture ∑k−1 j=1 Kj(z) is maximized. For f - divergences, we thus formulate the following proposition: Proposition 1. Consider C = maxDf (p||qM ), where Df (p||qM ) is an arbitrary f -divergence between some distribution p(z) and a mixture of kernels qM (z) = 1M ∑M j=1Kj(z), and Df (p||qM ) < ∞. Then, maximization of a marginal gain for set function
F (Z) = − ∫ f ( p(z)
1 M ∑M j=1Kj(z)
) 1
M M∑ m=1 Km(z)dz + C, (9)
at a step k of a greedy algorithm corresponds to
arg max zk ∆(zk|Z) = arg min zk Ez∼Kk(z)f
( p(z)
1 M ∑k j=1Kj(z)
) . (10)
Proof. The proof is shown Appendix A.2.
The case of reverse KL divergence. Having mean-seeking behavior is useful to fit the kernel density estimator using forward divergences. However, if one wants to optimize marginal gains, they need to have a mode-seeking behavior, which is achieved via optimizing reverse divergences (Zhang et al., 2019).
If we consider the generator for the reverse KL-divergence, f(x) = − log x, the minimization simplifies to
min zk Ez∼Kk(z) − log p(z) + log  1 M k∑ j=1 Kj(z)  . (11) In the sequel, we leverage the standard reverse KL-divergence generator, which aims to discover the modes of p(z), but also enforces the coverage of the approximation. Finally, we note that it is costly to compute the expectation in practice; we thus use the following point estimate of (11):
min zk − log p(zk) + log  1 M k−1∑ j=1 Kj(zk)  . (12)

heading: 4 GREEDY APPROXIMATION OF THE BAYESIAN POSTERIOR
text: Objective function. In this section, we consider parametric functions zθ : Rd → Rc, where θ denotes the parameters determining the function. We now use Bayes’ theorem and a mean-field approximation to express the posterior of an ensemble of such functions as
p(zθ|D) ∝ p(zθ1 , . . . , zθM |D) ∝ M∏ m=1 p(D|zθm)p(zθm), (13)
where θm are parameters, p(D|zθ) the likelihood and p(zθ) the prior; p(D|zθ) ∝∏n i=1 exp(−`(zθ(xi), yi)), and p(zθ) ∝ exp(−λ‖θ‖22).
To solve the problem defined by (12) in the context of Bayesian posterior approximation, we define the kernel density components via generalized exponential kernels Kj(zθ) ∝ exp(−λMd(zθ, zθj )2), where λM is proportional to the kernel width. Therefore, at a kth greedy step, we minimize
J(θk) = E(x,y)∼p(x,y)`(zθk(x), y) + λ‖θk‖22︸ ︷︷ ︸ Marginal gain on R(Z) + log
k−1∑ j=1 exp ( −λM M d(zθk , zθj ) 2 ) ︸ ︷︷ ︸
Marginal gain on ΩλM (Z)
, (14)
which is similar to the marginal gain on our originally defined high-level objective (1), except that the diversity term Ω has a different form based on our submodular f -divergence optimization.
Sampling-based approximation of the diversity term. When minimizing (14), one needs to be able to compute the diversity term Ω = log ∑k−1 j=1 Kj(−λMd(zθ, zθj )2), which is derived from a kernel density estimator we aim to fit to the true posterior. We note that this needs to be done in the function space, which makes this computation non-trivial.
We earlier defined Kj(z) to be an individual kernel in a mixture 1M ∑M j=1Kj(z). In order to be able to use the f -divergence, the individual components Kj(z) must be density functions centered at zj , which implies the need of a notion of similarity or the existence of a function norm, which are known to be NP-hard to compute for any neural network with depth greater than 3 (Rannen-Triki et al., 2019). We note that this intrinsic hardness result applies to all methods that define a meaningful posterior distribution in function space through a kernel density estimator. We thus use here a method from (Rannen-Triki et al., 2019) to approximate ‖z‖22, via i.i.d. samples xi ∼ P ∗, where P ∗ is a weighting distribution, which is required to ensure that Monte Carlo integration yields a reasonable approximation. This leads to the following sampling-based approximation of the diversity term:
log k−1∑ j=1 exp ( −λM M Ex∼p∗(x)‖zθk(x)− zθj (x)‖22 ) . (15)
A note on practical implementation. To this point, we defined all the main components of our method, except the weighting distribution in (15). The desirable behavior, which we expect an ensemble to exhibit, is that it must be uncertain on the OOD data and certain in the regions where the training data are available. This implies that p∗(x) must include OOD samples. One can use OOD data in training explicitly, however, we resort to a setting when OOD data are unknown, and use a simple heuristic, which fits a Gaussian to every data dimension with the variance ×5 larger than the variance of the data. We specify further details about the full algorithm in Appendix B.1.

heading: 5 RELATED WORK
text: Randomization-based ensembles. Generally, ensembles have been studied in Machine Learning over several decades for different classes of models (Hansen & Salamon, 1990; Breiman, 1996; Freund & Schapire, 1997; Lakshminarayanan et al., 2017). One can distinguish several main methods for diverse ensemble construction (Pearce et al., 2020): randomization of training initialization and hyperparameters (Hansen & Salamon, 1990; Wenzel et al., 2020b; Wilson & Izmailov, 2020; Zaidi et al., 2020), bagging (Breiman, 1996; 2001), boosting (Freund & Schapire, 1997), and explicit diversity training (Kuncheva & Whitaker, 2003; Ross et al., 2020; Yang et al., 2020; Brown et al., 2005; Kariyappa & Qureshi, 2019; Sinha et al., 2020; Melville & Mooney, 2005). Randomizationbased ensemble construction has shown good results in in-domain uncertainty (Ashukha et al., 2020) estimation, but also in the detection of OOD data (Lakshminarayanan et al., 2017).
Submodular ensemble pruning and greedy ensemble construction Submodularity in ensemble learning has previously been discussed in the context of ensemble pruning (Sha et al., 2014). The
goal of ensemble pruning is to trim a large ensemble of models so that the accuracy of the ensemble remains the same. We note that recently, inspired by the submodular pruning approach, a greedy algorithm has been applied to randomization based ensembles (Wenzel et al., 2020b; Zaidi et al., 2020). However, we also note that neither of these works approached the problem of ensemble construction from a Bayesian posterior approximation point of view. However, they provide extensive experimental evidence on the plausibility of greedy approach. Our work sheds light on why in Wenzel et al. (2020b); Zaidi et al. (2020) ensembles constructed greedily worked well in OOD detection tasks.
Diversity-promoting regularization for ensemble training. The ensemble literature contains a line of work focusing on explicitly promoting regularization in ensembles (Kuncheva & Whitaker, 2003; Ross et al., 2020; Yang et al., 2020; Brown et al., 2005; Kariyappa & Qureshi, 2019; Sinha et al., 2020; Melville & Mooney, 2005; Havasi et al., 2021). In terms of promoting diversity outside training data, the closest work to ours is (Ross et al., 2020; Rame & Cord, 2021), and in terms of the form of diversity regularization it is (Kariyappa & Qureshi, 2019). However, none of these works takes the perspective of approximating the Bayesian posterior. Another limitation of most of these approaches is that they use either out-of-distribution data in training, adversarial examples, or expensive generative models, thus making those methods difficult to scale to large datasets.
POVI A problem of learning diverse ensembles can be seen from a POVI perspective: in particular Stein Variational Gradient Descent (SVGD) (Wang & Liu, 2019). SVGD aims to learn a diverse set of functions Z, approximating arbitrary distributions via a set of particles using a reverse KL divergence. This approach is similar to ours, however, we tackle the problem of optimizing a general f -divergence in the function space. Furthermore, to our knowledge, the present work is the first that takes a submodular minimization perspective in BDL.
Greedy POVI Another line of work in the POVI family also considers greedy approximations (Futami et al., 2019; Jerfel et al., 2021). In both of these works, the authors perform particle based inference by also optimizing the kernel widths of the KDE, which is different to our method. Furthermore, these works do not tell what the diversity term for ensemble training should look like, and how to compute it in practice.
Function space POVI Conventionally, Bayesian inference in DL is thought of in the weight space (MacKay, 1992; Blundell et al., 2015; Gal & Ghahramani, 2016; Izmailov et al., 2020; Wilson & Izmailov, 2020; Pearce et al., 2020; Wenzel et al., 2020a). However, recent studies point out that despite the fact that simple priors over weights may imply complex posteriors in the function space, the connection between two of these is difficult to establish (Sun et al., 2019; Hafner et al., 2020).
Recent papers on function-space POVI (Wang et al., 2019; D’Angelo & Fortuin, 2021) point out that one can do Bayesian inference in the function space by optimizing an objective function with a repulsive term. Notably, they draw connection to the reverse KL-divergence minimization, and we thus consider our method to be closely connected to function space POVI.

heading: 6 EXPERIMENTS

heading: 6.1 SETUP
text: Datasets and models. We ran our main experiments on CIFAR10, CIFAR100 (Krizhevsky, 2009) and SVHN (Netzer et al., 2011) in-distribution datasets. Our OOD detection benchmark included CIFAR10, CIFAR100, DTD (Cimpoi et al., 2014), SVHN (Netzer et al., 2011), LSUN (Yu et al., 2015), TinyImageNet (Le & Yang, 2015), Places 365 (Zhou et al., 2017), Bernoulli noise images, Gaussian noise, random blobs image, and uniform noise images. The composition of the benchmark was inspired by the work of Hendrycks et al. (2018). We excluded the in-distribution datasets for each of the settings, resulting in a total of 10 OOD datasets for each in-distribution dataset. The full description of the benchmark is shown in Appendix B.2.
The experiments were conducted using ResNet164 (pre-activated version; denoted as PreResNet164) (He et al., 2016), VGG16 (with batch normalization (Ioffe & Szegedy, 2015); denoted as VGG16BN) (Simonyan & Zisserman, 2015), and WideResNet28x10 (Zagoruyko & Komodakis,
2016). All our models in the ensembles were trained for 100 epochs using PyTorch (Paszke et al., 2019), each ensemble on a single NVIDIA V100 GPU. In the case of CIFAR and SVHN experiments, we trained ensembles of size M = 11, and report the results across 5 different random seeds. For the CIFAR experiments, we selected λM ∈ {0.001, 0.005, 0.01, 0.05, 0.1, 0.5, 1, 3, 5, 7, 10}. In addition to the CIFAR and SVHN experiments, we used MNIST (LeCun et al., 1998) with ResNet8. The details of those experiments are shown in Appendix B.2.
Model selection and metrics. We used mutual information (MI) between the distribution of the predicted label ŷ for the point x̂ and the posterior distribution over functions p(f |D), to evaluate the epistemic uncertainty (Malinin & Gales, 2018; Depeweg et al., 2018) and reported the area under the ROC curve (AUC) and area under the precision-recall (PR) curve, i.e. average precision (AP) to quantify the OOD detection performance. Furthermore, we computed the false positive rate at 95% true positive rate (FPR95). Details on the computation of epistemic uncertainty can be found in Appendix B.3.

heading: 6.2 RESULTS
text: Illustrative examples. Figure 2 illustrates how our method performs on the two moons dataset. Here, we used a two-layer fully-connected network with ReLU activations (Fukushima, 1988). Having a high λM is important to obtain good uncertainty estimation. As expected, compared to our method, the DE method does not explicitly maximize the coverage of the posterior, and thus fails to be uncertain outside the training data.
Out-of-distribution detection. We present aggregated results for all the models and in-distribution datasets in Table 1. It is clear that on average (across OOD datasets), our method is substantially better than DE. This holds for all the architectures and in-distribution datasets. We show the expanded version of all the OOD detection results in Appendix C.3. An example of these results is shown in Figure 3 for all the models trained on CIFAR 100. Here, one can see that our method is at least similar to DE, and substantially better overall. Finally, some examples of OOD detection by both DE and our method are shown in Figure 4. Here, we computed optimal thresholds for each of the methods by optimizing the trade-off between the true positive and true negative rates.

heading: 7 DISCUSSION
text: In this paper, we have introduced a novel paradigm for Bayesian posterior approximation in Deep Learning using greedy ensemble construction via submodular optimization. We have proven a new general theoretical result, which shows that minimization of an f -divergence between some distribution and a kernel density estimator has approximation guarantees, and can be done greedily. We then derived a novel coverage promoting diversity term for ensemble construction. The results presented in this paper, as well as in Appendix C.4, demonstrate that our method outperforms the
state-of-the-art approach for ensemble construction, DE (Lakshminarayanan et al., 2017), on a range of benchmarks.
This study has some limitations, which outline several directions for the future work. Firstly, we did not compare our approach to a variety of existing methods for ensemble generation, e.g. snapshot ensembles (Huang et al., 2017), batch ensembles (Wen et al., 2020) or hyperparameter ensembles (Wenzel et al., 2020b). However, these methods are heuristic, and as discussed in the related work, our method can be used in conjunction with them to make them more principled. Furthermore, we note that the main contribution of this paper is novel theory.
The second limitation of this work is that it does not compare to f-POVI Wang et al. (2019); D’Angelo & Fortuin (2021). However, as noted earlier, those methods are in a different class, as they focus on training models without state-of-the-art training techniques such as batch normalization (Ioffe & Szegedy, 2015) and data augmentation.
The third limitation is that we did not fully explore different techniques of generating weighting distributions for the diversity term. However, we still evaluated the possibility of using real images for this purpose in an outlier exposure (Appendix C.4).
The final limitation is related to our method, as it lacks parallelization possibilities compared to DE. We note, however, that there is a class of parallel submodular optimization algorithms that enable parallelization and retain approximation guarantees (Ene & Nguyen, 2020).
To conclude, this paper provides a novel foundational framework for Bayesian Deep Learning. We hope for the wide adaption of the proposed method by practitioners across the fields. Our code is available at https://anonymous.4open.science/r/greedy_ensembles_ training-0CD6/.

heading: A PROOFS

heading: A.1 PROOF OF THEOREM 1
text: To prove this proposition, we make use of the following theorem:
Theorem A1 (Theorem 1.4 from (Mitrinovic & Vasic, 1970)). A function f : D → R is convex on D = [a, b] if and only if ∀x1 < x2 < x3 ∈ D∣∣∣∣∣x1 f(x1) 1x2 f(x2) 1x3 f(x3) 1
∣∣∣∣∣ ≥ 0. (16) Theorem 1. Any f -divergence
Df (p||qM ) = ∫ f ( p(z)
1 M ∑M j=1Kj(z)
) 1
M M∑ m=1 Km(z)dz (17)
between a distribution p(z) and a mixture of M kernels with equal weights is supermodualar in a cardinality-fixed setting, assuming that ∀zmaxqM Df (p(z)||qM (z)) <∞.
Proof. In order to prove that f -divergences are supermodular, we need to show that ∀α > 0, xf(αx ) is convex, because qM (z) is a positive modular function (Bach, 2019, Proposition 6.1), due to M being fixed.
From Theorem A1, ∀x1 < x2 < x3 in [a, b] f(x) is convex if an only if
∣∣∣∣∣x1 f(x1) 1x2 f(x2) 1x3 f(x3) 1 ∣∣∣∣∣ ≥ 0. (18)
We know that x1 ≤ x2 ≤ x3 and x1, α > 0. We divide each ith row by xi:∣∣∣∣∣∣ 1 1x1 f(x1) 1 x1 1 1x2 f(x2) 1 x2
1 1x3 f(x3) 1 x3 ∣∣∣∣∣∣ ≥ 0. (19) After the division, we denote new variables y1 = 1x3 , y2 = 1 x2
, and y3 = 1x2 . One can see that y1 < y2 < y3, because x1 < x2 < x3. We then get
∣∣∣∣∣∣ 1 y3f( 1 y3 ) y3 1 y2f( 1 y2
) y2 1 y1f(
1 y1 ) y1 ∣∣∣∣∣∣ ≥ 0. (20) Changing the first and the third row of the determinant will change the sign. Changing the third and the first columns will also change the sign. Therefore
∣∣∣∣∣∣ y1 y1f( 1 y1 ) 1 y2 y2f( 1 y2 ) 1
y3 y3f( 1 y3 ) 1 ∣∣∣∣∣∣ ≥ 0, (21) and thus we get that
∀x ∈ [a, b], xf ( 1
x
) is convex ⇐⇒ f(x) is convex. (22)
Consider now a mapping 1αy, which preserves convexity:
∣∣∣∣∣∣ αy1 αy1f( α y1 ) 1 αy2 αy2f( α y2 ) 1
αy3 αy3f( α y3 ) 1 ∣∣∣∣∣∣ ≥ 0. (23) Division of the first and the second column by α does not change the sign of the determinant, therefore
xf (α x ) ⇐⇒ f(x) convex, (24)
which concludes the proof.

heading: A.2 PROOF OF PROPOSITION 1
text: Proposition 1. Consider C = maxDf (p||qM ), where Df (p||qM ) is an arbitrary f -divergence between some distribution p(z) and a mixture of kernels qM (z) = 1M ∑M j=1Kj(z), and Df (p||qM ) < ∞. Then, maximization of a marginal gain for set function
F (Z) = − ∫ f ( p(z)
1 M ∑M j=1Kj(z)
) 1
M M∑ m=1 Km(z)dz + C, (25)
at a step k of a greedy algorithm corresponds to
arg max zk ∆(zk|Z) = arg min zk Ez∼Kk(z)f
( p(z)
1 M ∑k j=1Kj(z)
) . (26)
Proof. We aim to derive a marginal gain of adding an element defined by zk to 1M ∑k−1 j=1 Kj(z) 1.
Let us denote G(z) = f (
Mp(z)∑k−1 j=1 Kj(z)
) − f ( Mp(z)∑k j=1Kj(z) ) . Then
− ∫ f ( Mp(z)∑k j=1Kj(z) ) 1 M k∑ m=1 Km(z)dz + C + ∫ f ( Mp(z)∑k−1 j=1 Kj(z) ) 1 M k−1∑ m=1 Km(z)dz − C =
∫ G(z)
M k−1∑ m=1 Km(z)dz − ∫ f ( Mp(z)∑k j=1Kj(z) ) Kk(z) M dz. (27)
1Note: 1 M is a constant, which remains unchanged at all iterations of the greedy algorithm.
One can observe that the first term of (27) is upper-bounded by a constant, which is not dependent on zk: ∫
G(z)
M k−1∑ m=1 Km(z)dz ≤ ∫ f ( kp(z)∑k−1 j=1 Kj(z) ) 1 M k−1∑ j=1 Kj(z)dz = const, (28)
therefore, to maximize the marginal gain, one needs to maximize the second term. Consequently, we write the objective corresponding to a marginal gain as
∆(zk|Z \ zk) = − ∫ f ( p(z)
1 M ∑k j=1Kj(z)
) Kk(z)dz, (29)
maximization of which is equivalent to
arg min zk Ez∼Kk(z)f
( p(z)
1 M ∑k j=1Kj(z)
) , (30)
which concludes the proof.
B IMPLEMENTATION DETAILS

heading: B.1 PRACTICAL IMPLEMENTATION OF THE ALGORITHM
text: Weighting distribution for the diversity term We propose the following simple heuristic, defining p∗(x) as a normal distribution N (µD, α · ΣD) of dimensionality, corresponding to the training data. The covariance ΣD for this distribution is set to be diagonal, such that the variance for every dimension j is ΣD[j, j] = (α · σj)2, where α > 1 is a scaling parameter, and σ2j is a variance of the dimension j computed from samples of the training dataset D. Similarly, µD, the vector of expected values for every dimension, is also computed from the training data. Finally, the hyperparameter α = 5 was found to work well, and we thus report all the experimental results with it fixed. We note that a similar technique, but for in-distribution data generation has been used earlier in (Melville & Mooney, 2005).
The resulting algorithm The resulting, computationally tractable optimization algorithm for ensembles, which minimizes marginal gains (10), is shown in Algorithm B2. For simplicity, we omit the snapshot selection step, i.e. early stopping.
We note that contrary to the general Random Greedy method, shown in the main text, we can resort to a method with complexity of O(k). This is achieved through the fact that a uniform selection of the elements maximizing the marginal gain can be avoided, since at each greedy step, we initialize the new models randomly before maximizing the marginal gain. Another performance improvement can be gained by storing the evaluations zj(xi) ∀j = 1, . . . , k − 1 in memory before executing each kth step.
We report here also one important practical trick, which we found important during the training. Specifically, freezing the batch normalization layers (Ioffe & Szegedy, 2015) before computing the diversity term turned out to help the convergence substantially. We anticipate that the diversity term weighting distribution approximated as a simple multivariate Gaussian with diagonal covariance may be corrupting the batch norm statistics. We thus think that using other, more sophisticated techniques for generating the weighting distribution samples might provide better results.

heading: B.2 OOD DETECTION BENCHMARKS
text: MNIST The MNIST dataset benchmark included 6 datasets: Fashion MNIST Xiao et al. (2017), DTD Cimpoi et al. (2014), Omniglot Lake et al. (2015), Gaussian noise, Bernoulli noise, and uniform noise datasets (see Table B1). We re-scaled all the images to the range of [0, 1]. Subsequently, we applied the same mean and standard deviation as we have applied to the original images before feeding them to the network.
Algorithm B2 O(k) Random Greedy algorithm for training ensembles of neural networks. 1: Input: D = {(xi, yi)}ni=1 – Dataset 2: Input: M – Size of the ensemble 3: Input: N – Number of iterations 4: Input: α – Variance parameter for p∗(x) 5: Input: Nb – Mini-batch size 6: Z ← ∅ 7: D∗ ← {(x∗, y∗i )}ni=1 ∼ N (µD, α · ΣD) 8: while |Z| < M do 9: k ← |Z|; 10: Randomly initialize zθk ; 11: for i = 1 to N do 12: Di ← {(xb, yb)}Nbb=1 ∼ D; 13: D∗i ← {x∗b} Nb b=1 ∼ D∗; 14: L← E(x̂,ŷ)∼Di`(zθk(x̂), ŷ) + λ‖θk‖22 15: Ω← 0 16: if k > 1 then 17: for j = 1 to k − 1 do 18: dij ← Ex∗∼D∗i ‖zθk(x
∗)− zθj (x∗)‖22 19: end for 20: Ω← log ∑k−1 m=1 exp(− λM M dim) 21: end if 22: Update θ using∇θ (L+ Ω) 23: end for 24: Z ← Z ∪ {zθk} 25: end while 26: return Z
CIFAR and SVHN The CIFAR and SVHN OOD benchmark included 10 different datasets. Here, we also DTD, Bernoulli noise, Gaussian noise and uniform noise datasets, and added Places 365 (Zhou et al., 2017), Tiny ImageNet (Le & Yang, 2015; Deng et al., 2009), and LSUN (Yu et al., 2015) datasets to the benchmark. For CIFAR10 as in-domain data, we added CIFAR100 and SVHN (Netzer et al., 2011) to the benchmark. For CIFAR100 – CIFAR10 and SVHN. Finally, for SVHN, we added CIFAR10 and CIFAR100 as OOD datasets, making a total of 10 OOD datasets per 1 in-distribution dataset. The details about each of the datasets are shown in Table B1.
Before feeding the images to the network, we applied re-scaling similarly to the MNIST setting. We also used the same mean and standard deviation normalization procedure, as for the in-domain data.
Table B1: Description of the datasets used in all the exp. R indicates real images, S – synthetic.
Dataset Type # samples Comment Uniform
S 25, 000 N/A Gaussian 25, 000 Generated once, used in all experiments Blobs 25, 000 N/A Bernoulli 25, 000 N/A Omniglot 13, 181 Evaluation images
CIFAR10
R 10, 000 Test set (not used in training) CIFAR100 10, 000 Test set (not used in training) SVHN 73, 257 Test set (not used in training) Places 365 10, 000 First 10, 000 images from the test set (sorted alphabetically) TinyImageNet 10, 000 Original validation set images DTD 5, 640 Release 1.0.1 LSUN 10, 000 Test set Fashion MNIST 10, 000 Test set

heading: B.3 EPISTEMIC UNCERTAINTY COMPUTATION
text: We used the epistemic uncertainty, i.e. mutual information (MI) between the distribution of predicted label ŷ for the point x̂ and the posterior distribution over functions p(f |D), to evaluate the uncertainty (Malinin & Gales, 2018; Depeweg et al., 2018). As a distribution over weights induces a distribution over functions, we approximate the MI as:
I(ŷ; f |x̂,D) = H [ Ep(θ|x̂,D)p(ŷ|θ, x̂,D) ] − Ep(θ|x̂,D)H [p(ŷ|θ, x̂,D)] , (31)
where H[·] denotes the entropy. One can see that this metric can be efficiently computed from the predictions of an ensemble.

heading: C EXPERIMENTS

heading: C.1 EXPERIMENTAL DETAILS
text: Model selection Contrary to the commonly used practice, we did not use CIFAR10/100 and SVHN test set sets for model selection. Neither did we use any OOD data. Instead, we used validation set accuracy (10% of the training data; randomly chosen stratified split) to select the models when optimizing the marginal gain. The best snapshot was found using the validation data, was then selected for final testing. When selecting the models for evaluation on OOD data, we first evaluated ensembles on the in-distribution test set (Appendix C.2). Subsequently, we selected the highest λM that did not harm the test set (in-domain) performance (no overlap of confidence intervals defined as mean ± standard error). To provide additional information, we also analyzed adaptive calibration error (ACE) with 30 bins (Nixon et al., 2019).
Two moons dataset For the synthetic data experiments, we used scikit-learn (Pedregosa et al., 2011), and generated a two-moons dataset with 300 points in total, having the noise parameter fixed to 0.3. Here, we used a two-layer neural network with ReLU (Krizhevsky, 2009) activations and hidden layer size of 128.
CIFAR10/100 and SVHN The main training hyper-parameters were adapted from (Maddox et al., 2019) (see Table C2), but with additional modifications inspired by (Malinin & Gales, 2018; Smith & Topin, 2019), which helped to train the CIFAR models to state-of-the-art performance in only 100 epochs. As such, we first employed a warm-up of the learning rate (LR) from a value 10 times lower than the initial LR (LRinit in Table C2) for 5 epochs. Subsequently, after 50% of the training budget, we linearly annealed the LR to the value of LR× lrscale until 90% of the training budget is reached, after which we kept the value of LR constant.
All models were trained using stochastic gradient descent with momentum of 0.9 and a total batch size of 128. We employed standard training augmentations – horizontal flipping, reflective padding to 34× 34, and random crop to 34× 34 pixels.
Model LRinit Nesterov Weight Decay lrscale PreResNet164 0.1 Yes 0.0001 0.01 VGG16BN 0.05 No 0.0005 0.01 WideResNet28x10 0.1 No 0.0005 0.001 Table C2: Main hyperparameters of all the models used in the CIFAR and SVHN in-domain experiments.
MNIST In addition to the CIFAR10/100 experiments, we also trained our method on MNIST (LeCun et al., 1998) with PreResNet8 architecture. As OOD, we used FashionMNIST (Xiao et al., 2017) and Omniglot (Lake et al., 2015) datasets. We also tested other architectures, such as PreResNet20, but the models with higher depth than 8 already gave nearly perfect scores on MNIST.
Hyper-parameter-wise, we trained all the models for 20 epochs without warmup with the batch size of 256 using plain SGD with momentum. The weight decay was set to 1e−5. We used LR annealing similarly as for CIFAR experiments, but used lrscale = 0.0001. No data augmentations were used in any of the MNIST experiments. λM was searched in range {0.0001, 0.001, 0.01, 0.11, 7} for
0.0 2.5 5.0 7.5 10.0 M
93.0
93.5
94.0
94.5
95.0
95.5
96.0
96.5
97.0
Ac cu
ra cy
[% ]
PreResNet164 VGG16BN WideResNet28x10
(a)
0.0 2.5 5.0 7.5 10.0 M
0
1
2
3
4
5
6
AC E
[% ]
(b)
0.0 2.5 5.0 7.5 10.0 M
76
78
80
82
84
Ac cu
ra cy
[% ]
(c)
0.0 2.5 5.0 7.5 10.0 M
0.0
0.2
0.4
0.6
0.8
1.0
AC E
[% ]
(d)
Figure C1: Relationship between accuracy, ACE, and λM (M = 11). Subplots (a) and (b) show the results for CIFAR10. Subplots (c) and (d) show the results for CIFAR100.
M ∈ {3, 5, 9, 15}. This series of experiments was re-run 3 times, as the MNIST dataset is rather simple, and the test scores have low variance between the runs.

heading: C.2 CIFAR10/100 AND SVHN IN-DOMAIN PERFORMANCE
text: CIFAR10/100 in-distribution performance vs. diversity. Figure C1 provides an illustration of how the test set performance changes with λM on CIFAR data. One can see a general trend that when λM approaches M , the models lose the ability to make accurate predictions, which results in lower accuracy and poorer calibration. Interestingly, performance on the VGG model degrades much slower with λM compared to other architectures. Similar findings were also obtained for the SVHN dataset. Based on the test performance, we selected the models for further evaluation on OOD benchmark.
CIFAR and SVHN: best models’ performance Table C3 shows the results of all the trained models on the in-domain data. One can see that the results between Deep Ensembles (DE) (Lakshminarayanan et al., 2017) do not differ significantly. We trained all these models according to the earlier specified hyper-parameters and the learning rate schedule. Models selected in Table C3 are used to report the results in the main experiments.

heading: C.3 DETAILED CIFAR, SVHN RESULTS
text: Detalized versions of the results presented in the main text are shown in Table C4. The corresponding λM coefficients are the same as in Table C3.

heading: C.3.1 MNIST DETAILED RESULTS
text: In-domain performance Model selection scheme on MNIST was exactly the same as for the CIFAR10/100. We illustrate the relationship between λM , accuracy, and the Adaptive Calibration error (ACE) in Figure C2. One can see that λM remains the same even when M is increasing. Table C7 shows the detailed in-domain performance for the best λM , equal to 0.1.
0.0 0.5 1.0 1.5 2.0 M
98.00
98.25
98.50
98.75
99.00
99.25
99.50
99.75
100.00
Ac cu
ra cy
[% ]
M=3 M=5 M=9 M=15
(a)
0.0 0.5 1.0 1.5 2.0 M
0
5
10
15
20
25
30
NL L
×1 00
(b)
0.0 0.5 1.0 1.5 2.0 M
0.0
0.1
0.2
0.3
0.4
0.5
AC E
[% ]
(c)
Figure C2: Relationship between accuracy, negative loglikelihood, ACE, and λM for different M on MNIST (LeCun et al., 1998). Subplots (a) and (b) show the results for PreResNet8. Experiments were re-run 3 times with different seeds. We found that the best results were obtained with λM = 0.1.
Table C3: In-domain performance on the test sets of CIFAR10/100 and SVHN for all the models used in the experiments (M = 11). We report mean and standard error over 5 random seeds for each of the models. λM = 0 indicates Deep Ensembles (Lakshminarayanan et al., 2017). Standard errors are reported if they are more than 0.01 across runs. DE indicates Deep Ensembles.

heading: Architecture Dataset Method Accuracy (%) NLL ×100 ACE (%)
text: PreResNet164
C10 DE 95.70±0.02 13.28±0.06 0.18Ours (λM = 3) 95.66±0.02 13.18±0.09 0.16
C100 DE 79.97±0.04 73.44±0.17 0.08Ours (λM = 5) 79.93±0.04 74.16±0.91 0.09±0.01
SVHN DE 99.46±0.01 2.34±0.04 0.25Ours (λM = 1) 99.38±0.01 3.16±0.13 0.81±0.12
VGG16BN
C10 DE 94.55±0.02 17.59±0.05 0.23Ours (λM = 5) 94.48±0.06 17.84±0.15 0.23±0.01
C100 DE 76.32±0.09 91.07±0.35 0.11Ours (λM = 5) 76.78±0.07 89.10±0.33 0.10
SVHN DE 99.40±0.01 2.71±0.02 0.18±0.01Ours (λM = 1) 99.25±0.01 3.94±0.10 0.95±0.12
WideResNet28x10
C10 DE 96.56±0.02 10.76±0.06 0.16±0.01Ours (λM = 1) 96.54±0.01 10.99±0.03 0.18±0.01
C100 DE 83.08±0.09 62.05±0.18 0.09Ours (λM = 1) 83.02±0.06 62.20±0.13 0.08
SVHN DE 99.45±0.01 2.53±0.03 0.43±0.01Ours (λM = 1) 99.38 2.87±0.04 0.46±0.01
Out of distribution detection The final OOD benchmark results on MNIST are summarized in Table C8. One can see that on MNIST, our method yields a performance boost for any size of an ensemble, even for the very small ones (M = 3).

heading: C.4 ADDITIONAL RESULTS
text: Does the choice weighting distribution for the diversity term affect the results? We considered PreResNet164 trained on CIFAR100 with CIFAR10 as a weighting distribution dataset. This approach can be seen as a form of outlier exposure technique, earlier proposed by Hendrycks et al. (2018). The results on OOD benchmark (excluding CIFAR10) are shown in Table C9. One can observe that the relationship between changing to a dataset of real images and the OOD detection performance is non-trivial. While for some datasets outlier exposure did bring benefit, for some other datasets, it did not. Notably, the datasets on which the benefit of outlier exposure was the best were datasets of real images, and we think that in case unlabeled images with non-overlapping class distribution are available, using them can bring benefit.
Effect of model capacity While the main experiments in the paper were conducted using large models, we also investigated whether ensembles of smaller models can benefit from our method. We followed the same λM selection procedure as for the main experiments in the paper, and report results for PreResNet20 in Table C10 for an ensemble of size M = 11. We found that the optimal λM in this case is smaller compared to λM for PreResNet164, and our method here does not yield any substantial boost over Deep Ensembles on these data. We found that all the models with small capacity trained on CIFAR diverged with high λM
Robustness to distribution shift As an additional evaluation, we investigated whether our method performs on par with DE under the distribution shift, to make sure that introduction of additional regularization did not affect the robustness properties. We thus use a corrupted version of CIFAR10 test set released by (Hendrycks & Dietterich, 2019), and it has been recently shown that DE outperform many other methods on this benchmark (Ovadia et al., 2019). Here, we report the results for
Table C4: CIFAR10 results. We report mean and standard error over 5 random seeds for each of the models. Standard errors are reported if they are more than 0.01 across runs. DE indicates Deep Ensembles.
Architecture OOD dataset DE Ours AUC (↑) AP (↑) FPR95 (↓) AUC (↑) AP (↑) FPR95 (↓)
PreResNet164
bernoulli 0.98±0.01 0.97±0.01 0.04±0.01 1.00 1.00 0.00 blobs 0.96 0.98 0.12±0.01 0.96 0.98 0.12±0.01 cifar100 0.90 0.87 0.30 0.90 0.88 0.30 dtd 0.93 0.83 0.19±0.01 0.96 0.93 0.14 gaussian 0.93±0.01 0.94±0.01 0.16±0.01 0.96±0.02 0.97±0.02 0.11±0.03 lsun 0.93 0.89 0.20 0.95 0.94 0.18 places 0.92 0.89 0.21 0.94 0.93 0.19 svhn 0.94 0.99 0.16 0.95 0.99 0.14 tiny imagenet 0.91 0.88 0.28 0.92 0.89 0.26 uniform 0.98±0.01 0.97±0.01 0.04±0.01 1.00 1.00 0.00
VGG16BN
bernoulli 0.94±0.01 0.95±0.01 0.11±0.02 1.00 1.00 0.00 blobs 0.96 0.98 0.16 0.96 0.98 0.15 cifar100 0.89 0.86 0.34 0.89 0.86 0.34 dtd 0.90 0.77±0.01 0.25±0.01 0.96 0.92 0.18 gaussian 0.95 0.97 0.14±0.01 0.99 0.99 0.05±0.01 lsun 0.93 0.91 0.24 0.95 0.94 0.21 places 0.91 0.90 0.27±0.01 0.94 0.93 0.23 svhn 0.87 0.97 0.28±0.01 0.87 0.97 0.27±0.01 tiny imagenet 0.90 0.88 0.33 0.90 0.88 0.32 uniform 0.90±0.02 0.89±0.01 0.16±0.02 1.00 1.00 0.00
WideResNet28x10 bernoulli 1.00 1.00 0.00 1.00 1.00 0.00 blobs 0.96 0.97 0.11±0.01 0.97 0.98 0.10±0.01 cifar100 0.92 0.89 0.27 0.91 0.89 0.27 dtd 0.93 0.86±0.01 0.22±0.01 0.97 0.94 0.14±0.01 gaussian 0.96±0.01 0.97±0.01 0.09±0.01 1.00 1.00 0.00 lsun 0.93 0.91 0.20 0.95 0.95 0.17±0.01 places 0.93 0.91 0.21 0.95 0.94 0.18 svhn 0.96 0.99 0.12 0.95 0.99 0.13±0.01 tiny imagenet 0.92 0.90 0.27 0.93 0.91 0.25 uniform 1.00 1.00 0.00 1.00 1.00 0.00
VGG16BN and PreResNet164, as they yielded OOD performance gain in both SVHN and LSUN datasets. One can see from Figure C3 that our method performs on-par with DE, as no statistical significance in difference between methods can be concluded from this plot. This further supports our claims that the developed greedy ensemble training approach works the same or on par with DE.
Effect of ensemble size on CIFAR We ran our experiments using PreResNet164 with on CIFAR10/100, having M ∈ {3, 5, 7} and λM{0.1, 0.5, 0.8, 1, 1.5, 2, 3}. Both in-domain accuracy, and the OOD detection on LSUN and SVHN are shown in Table C11. The results in that table show that with small ensemble size, our method may yield better performance than DE on SVHN, abut even with an ensemble of size M = 3, it has a substantial boost over DE in detecting LSUN.
Qualitative results on CIFAR100. In Figure C4, we further illustrate the capabilities of uncertainty estimation of our method for the PreResNet164 trained on CIFAR100 (M = 11). Our method has a better true positive rate (as can also be seen from the histograms), and slightly better precision when the threshold for the recall is high. We note that the histograms of epistemic uncertainties still overlap significantly, however, with our method, the model does not have a high number of overconfident predictions coming from the OOD data anymore.
Table C5: CIFAR100 results. We report mean and standard error over 5 random seeds for each of the models. Standard errors are reported if they are more than 0.01 across runs. DE indicates Deep Ensembles.
Architecture OOD dataset DE Ours AUC (↑) AP (↑) FPR95 (↓) AUC (↑) AP (↑) FPR95 (↓)
PreResNet164
bernoulli 0.81±0.03 0.82±0.03 0.24±0.04 1.00 1.00 0.00 blobs 0.92±0.01 0.95 0.25±0.02 0.92±0.02 0.95±0.01 0.23±0.03 cifar10 0.80 0.76 0.57 0.78±0.01 0.75 0.67±0.02 dtd 0.76±0.01 0.61±0.01 0.64±0.01 0.80±0.01 0.75±0.01 0.78±0.05 gaussian 0.80±0.01 0.83±0.01 0.37±0.02 0.95±0.02 0.97±0.01 0.16±0.05 lsun 0.86 0.81 0.45 0.87 0.85 0.47±0.01 places 0.82 0.77 0.52 0.83 0.81 0.60±0.03 svhn 0.80±0.01 0.96 0.55±0.01 0.83±0.01 0.96 0.50±0.01 tiny imagenet 0.82 0.79 0.53 0.82 0.79 0.58±0.01 uniform 0.51±0.09 0.66±0.05 0.55±0.09 1.00 1.00 0.00
VGG16BN
bernoulli 0.87±0.02 0.88±0.02 0.24±0.03 1.00 1.00 0.00 blobs 0.95 0.97 0.16±0.01 0.97 0.98 0.12±0.02 cifar10 0.78 0.73 0.63 0.78 0.74 0.64±0.01 dtd 0.73 0.53 0.62±0.01 0.86 0.80±0.01 0.50±0.01 gaussian 0.87±0.02 0.90±0.01 0.35±0.04 0.95±0.01 0.97±0.01 0.15±0.03 lsun 0.85 0.82 0.47 0.90 0.89 0.40 places 0.82 0.78 0.55 0.86 0.85 0.51 svhn 0.76±0.01 0.95 0.67±0.02 0.76±0.01 0.95 0.72±0.04 tiny imagenet 0.81 0.78 0.56 0.83 0.79 0.55 uniform 0.86±0.02 0.89±0.02 0.28±0.04 1.00 1.00 0.00
WideResNet28x10 bernoulli 0.97±0.02 0.96±0.02 0.04±0.02 1.00 1.00 0.00 blobs 0.95 0.97 0.16±0.01 0.98±0.01 0.99 0.09±0.02 cifar10 0.80 0.74 0.53 0.81 0.76 0.54±0.01 dtd 0.84±0.01 0.75±0.01 0.54±0.01 0.88±0.01 0.84±0.01 0.49±0.02 gaussian 0.72±0.08 0.78±0.05 0.39±0.09 0.94±0.03 0.97±0.01 0.20±0.07 lsun 0.87 0.81±0.01 0.35 0.91 0.90±0.01 0.32±0.01 places 0.84 0.79±0.01 0.44±0.01 0.88 0.87±0.01 0.41±0.01 svhn 0.80 0.95 0.52±0.01 0.80±0.01 0.95 0.52±0.01 tiny imagenet 0.84 0.78 0.46 0.85 0.81 0.46 uniform 0.95±0.01 0.96±0.01 0.12±0.03 1.00 1.00 0.00
1 2 3 4 5 Corruption severity
20
40
60
80
Ac cu
ra cy
(% )
Method deepens greedy
(a)
1 2 3 4 5 Corruption severity
0
2
4
6
8
10
12
AC E
(% )
(b)
1 2 3 4 5 Corruption severity
20
40
60
80
Ac cu
ra cy
(% )
(c)
1 2 3 4 5 Corruption severity
0
2
4
6
8
10
12
AC E
(% )
(d)
Figure C3: CIFAR10 Robustness benchmark results (M = 11) (Hendrycks & Dietterich, 2019). Subplots (a) and (b) show the results for PreResNet164 trained with λM = 3. Subplots (c) and (d) show the results for VGG16BN trained with λM = 5. The results have been averaged over 5 seeds.
Table C6: SVHN results (averaged across 5 seeds)
Architecture OOD dataset DE Ours AUC (↑) AP (↑) FPR95 (↓) AUC (↑) AP (↑) FPR95 (↓)
PreResNet164
bernoulli 1.00 0.99 0.01 1.00 1.00 0.00 blobs 1.00 0.99 0.01 0.99 0.98 0.02 cifar10 0.99 0.97 0.02 0.99 0.97 0.02 cifar100 0.99 0.96 0.02 0.99 0.96 0.04 dtd 0.99 0.94 0.02 1.00 0.98 0.01 gaussian 1.00 0.99 0.01 1.00 1.00 0.00 lsun 0.99 0.96 0.02 1.00 0.98 0.01 places 0.99 0.96 0.02 1.00 0.98 0.01 tiny imagenet 0.99 0.96 0.02 1.00 0.97 0.02 uniform 1.00 0.99 0.01 1.00 1.00 0.00
VGG16BN
bernoulli 1.00 0.99 0.01 1.00 1.00 0.00 blobs 1.00 0.98 0.01 1.00 0.99 0.01 cifar10 0.99 0.95 0.02 0.99 0.96 0.03 cifar100 0.99 0.94 0.02 0.99 0.95 0.05 dtd 0.99 0.93 0.02 1.00 0.97 0.01 gaussian 1.00 0.99 0.01 1.00 1.00 0.00 lsun 0.99 0.95 0.02 1.00 0.99 0.01 places 0.99 0.96 0.02 1.00 0.98 0.01 tiny imagenet 0.99 0.95 0.02 0.99 0.97 0.03 uniform 1.00 0.99 0.01 1.00 1.00 0.00
WideResNet28x10 bernoulli 1.00 0.99±0.01 0.02±0.01 1.00 1.00 0.00 blobs 0.99 0.98 0.02 1.00 0.99 0.01 cifar10 0.99 0.96 0.02 1.00 0.97 0.02 cifar100 0.99 0.95 0.03 0.99 0.97 0.02 dtd 0.99 0.91±0.01 0.04 1.00 0.99 0.00 gaussian 1.00 0.98 0.02 1.00 1.00 0.00 lsun 0.99 0.95 0.03 1.00 0.99 0.00 places 0.99 0.95 0.03 1.00 0.99 0.00 tiny imagenet 0.99 0.96 0.02 1.00 0.98 0.01 uniform 0.99 0.98±0.01 0.03±0.01 1.00 1.00 0.00
Size Method Accuracy (%) NLL ×100 ECE (%)
3 DE 99.44±0.02 1.72±0.03 0.23±0.01 Ours 99.43±0.03 1.74±0.07 0.24±0.04
5 DE 99.43±0.01 1.65±0.04 0.22±0.01 Ours 99.46±0.01 1.72±0.02 0.32±0.01
7 DE 99.44±0.01 1.61±0.03 0.25±0.01 Ours 99.41±0.01 1.66±0.06 0.28
9 DE 99.44±0.01 1.65±0.04 0.28±0.01 Ours 99.47±0.01 1.62±0.03 0.31±0.01
15 DE 99.47±0.01 1.60±0.03 0.29±0.01 Ours 99.49 1.57±0.03 0.30±0.02
Table C7: Test set results (in-domain performance) of our method on PreResNet8 trained on MNIST with different ensemble sizes M . We report the means over 3 random seeds for each of the models. λM = 0.1 Was used for all the experiments when training our method. Standard errors are reported if they are more than 0.01 across runs. DE indicates Deep Ensembles.
Table C8: Out-of-distribution detection results of our method with PreResNet8 trained on MNIST with different ensemble sizes M . We report the means over 3 random seeds for each of the models. λM = 0.1 Was used for all the experiments when training our method. Standard errors are reported if they are more than 0.01 across runs. DE indicates Deep Ensembles.
Size OOD dataset DE Ours AUC (↑) AP (↑) FPR95 (↓) AUC (↑) AP (↑) FPR95 (↓)
3
bernoulli 0.09±0.05 0.52±0.01 0.98 1.00 1.00 0.00 dtd 0.41±0.16 0.47±0.13 0.97±0.01 1.00 0.99 0.00 fashion mnist 0.89±0.03 0.92±0.02 0.71±0.19 0.99±0.01 0.99 0.07±0.03 gaussian 0.98±0.01 0.99±0.01 0.06±0.02 0.99±0.01 1.00 0.03±0.02 omniglot 0.98 0.98 0.08 0.98 0.98 0.08 uniform 0.23±0.17 0.58±0.07 0.90±0.06 1.00 1.00 0.00
5
bernoulli 0.02 0.50 0.98 1.00 1.00 0.00 dtd 0.48±0.17 0.55±0.14 0.94±0.03 1.00 1.00 0.00 fashion mnist 0.87±0.05 0.92±0.03 0.67±0.24 0.99 0.99 0.04±0.01 gaussian 1.00 1.00 0.02±0.01 1.00 1.00 0.01±0.01 omniglot 0.98 0.99 0.06 0.98 0.99 0.07 uniform 0.30±0.22 0.62±0.10 0.77±0.18 1.00 1.00 0.00
7
bernoulli 0.33±0.24 0.66±0.13 0.77±0.18 1.00 1.00 0.00 dtd 0.75±0.12 0.76±0.11 0.67±0.25 1.00 1.00 0.00 fashion mnist 0.95±0.02 0.97±0.02 0.29±0.18 0.99 0.99 0.05±0.01 gaussian 1.00 1.00 0.02±0.01 1.00 1.00 0.02±0.01 omniglot 0.99 0.99 0.06 0.99 0.99 0.06±0.01 uniform 0.49±0.22 0.71±0.12 0.64±0.26 1.00 1.00 0.00
9
bernoulli 0.45±0.17 0.70±0.08 0.90±0.05 1.00 1.00 0.00 dtd 0.83±0.06 0.83±0.06 0.68±0.23 1.00 1.00 0.00 fashion mnist 0.97±0.01 0.97 0.16±0.04 1.00 1.00 0.01 gaussian 1.00 1.00 0.01±0.01 1.00 1.00 0.00 omniglot 0.99 0.99 0.06 0.99 0.99 0.06 uniform 0.72±0.15 0.82±0.09 0.50±0.21 1.00 1.00 0.00
15 bernoulli 0.18±0.06 0.55±0.02 0.99±0.01 1.00 1.00 0.00 dtd 0.82±0.04 0.80±0.04 0.84±0.09 1.00 1.00 0.00 fashion mnist 0.97 0.97 0.17±0.03 1.00 1.00 0.01 gaussian 1.00 1.00 0.01±0.01 1.00 1.00 0.01±0.01 omniglot 0.99 0.99 0.05 0.99 0.99 0.05 uniform 0.70±0.13 0.80±0.07 0.53±0.19 1.00 1.00 0.00
Table C9: Outlier exposure results for PreResNet164 trained on CIFAR100 with CIFAR10 as a weighting distribution for the diversity term. We found λM = 1 to be the best one for this experiment in the outlier exposure. The most optimal setting for our model was λM = 3 and λM = 5, respectively.
OOD dataset DE Ours Ours w. Outlier Exposure AUC (↑) AP (↑) FPR95 (↓) AUC (↑) AP (↑) FPR95 (↓) AUC (↑) AP (↑) FPR95 (↓)
bernoulli 0.81±0.03 0.82±0.03 0.24±0.04 1.00 1.00 0.00 1.00 1.00 0.00 blobs 0.92±0.01 0.95 0.25±0.02 0.92±0.02 0.95±0.01 0.23±0.03 0.96±0.01 0.98 0.14±0.01 dtd 0.76±0.01 0.61±0.01 0.64±0.01 0.80±0.01 0.75±0.01 0.78±0.05 0.81±0.01 0.71±0.01 0.57 gaussian 0.80±0.01 0.83±0.01 0.37±0.02 0.95±0.02 0.97±0.01 0.16±0.05 0.94±0.01 0.96±0.01 0.19±0.03 lsun 0.86 0.81 0.45 0.87 0.85 0.47±0.01 0.89 0.88±0.01 0.40±0.01 places 0.82 0.77 0.52 0.83 0.81 0.60±0.03 0.86 0.84±0.01 0.49±0.01 svhn 0.80±0.01 0.96 0.55±0.01 0.83±0.01 0.96 0.50±0.01 0.78±0.01 0.95 0.56±0.02 tiny imagenet 0.82 0.79 0.53 0.82 0.79 0.58±0.01 0.83 0.79 0.52 uniform 0.51±0.09 0.66±0.05 0.55±0.09 1.00 1.00 0.00 0.98±0.01 0.98±0.02 0.05±0.02
Table C10: Out of distribution detection for a small-capacity model – PreResNet20 (M = 11). Results were averaged over 5 random seeds. Standard errors are not reported if they less than 0.01 across runs. DE indicates Deep Ensembles.
Dataset OOD dataset DE Ours AUC (↑) AP (↑) FPR95 (↓) AUC (↑) AP (↑) FPR95 (↓)
C10
bernoulli 0.99±0.01 0.99±0.01 0.05±0.02 1.00 1.00 0.00 blobs 0.97 0.98 0.12 0.97 0.98 0.12±0.01 cifar100 0.88 0.86 0.37 0.89 0.86 0.36 dtd 0.93 0.86±0.01 0.22±0.01 0.94±0.01 0.89±0.01 0.20±0.01 gaussian 0.94±0.01 0.96±0.01 0.18±0.03 0.95±0.01 0.97±0.01 0.15±0.03 lsun 0.93 0.92 0.25 0.94 0.93 0.23±0.01 places 0.92 0.90 0.27 0.93 0.91 0.26 svhn 0.91±0.01 0.98 0.23±0.01 0.91 0.98 0.23±0.01 tiny imagenet 0.90 0.88 0.33 0.90 0.88 0.32 uniform 0.99 0.99 0.03±0.01 1.00 1.00 0.00
C100
bernoulli 0.91±0.05 0.93±0.04 0.19±0.08 1.00 1.00 0.00 blobs 0.95 0.97 0.17±0.01 0.97 0.98 0.12 cifar10 0.78 0.73 0.63 0.77 0.72 0.63 dtd 0.81±0.01 0.72±0.01 0.66±0.03 0.83±0.01 0.77±0.01 0.63±0.03 gaussian 0.98±0.01 0.99 0.07±0.01 0.99 0.99 0.04±0.01 lsun 0.88 0.87 0.43±0.01 0.89 0.88 0.41±0.01 places 0.84 0.82 0.55±0.01 0.85 0.83 0.52±0.01 svhn 0.85 0.97 0.48±0.01 0.84±0.01 0.97 0.49±0.02 tiny imagenet 0.82 0.78 0.56 0.82 0.79 0.54 uniform 0.90±0.03 0.92±0.02 0.23±0.07 1.00 1.00 0.00

heading: M Dataset Method Accuracy (%) NLL ×100 ACE (%) SVHN LSUN
text: 0 1 2 Epistemic unc.
0
500
1000
1500
2000
2500
3000
Co un
t
LSUN CIFAR100
(a)
0 1 2 Epistemic unc.
0
500
1000
1500
2000
2500
3000
Co un
t
LSUN CIFAR100
(b)
0.0 0.5 1.0 FPR
0.0
0.2
0.4
0.6
0.8
1.0
TP R
M = 5 M = 0
(c)
0.0 0.5 1.0 Recall
0.4
0.6
0.8
1.0
Pr ec
isi on
M = 5 M = 0
(d)
0 1 2 Epistemic unc.
0
500
1000
1500
2000
2500
3000
Co un
t
SVHN CIFAR100
(e)
0 1 2 Epistemic unc.
0
500
1000
1500
2000
2500
3000
Co un
t
SVHN CIFAR100
(f)
0.00 0.25 0.50 0.75 1.00 FPR
0.0
0.2
0.4
0.6
0.8
1.0
TP R
M = 5 M = 0
(g)
0.00 0.25 0.50 0.75 1.00 Recall
0.0
0.2
0.4
0.6
0.8
1.0
Pr ec
isi on
M = 5 M = 0
(h)
Figure C4: Uncertainty estimation quality on PreResNet164 (He et al., 2016) trained on CIFAR100 and evaluated on CIFAR100 test set vs LSUN and SVHN, respectively. Histograms (a) and (e) indicate Deep Ensembles (Lakshminarayanan et al., 2017). Histograms (b) and (f) show our method trained with λM = 5. Subplots (c) and (d) show the ROC and PR curves for the LSUN dataset, respectively. Subplots (g) and (h) show the ROC and PR curves for the SVHN dataset, respectively The curves were computed using average epistemic uncertainty per sample (5 seeds). Standard errors were < 0.01.

Review summaries:
The paper presents a combinatorial approach to controlling and improving the effects of diversity; a relevant and interesting feature in Bayesian networks. The approach is general in the sense that it involves minimization of an f-divergence, for any f. The approach is mainly compared to standard deep ensembles, i.e. without controlling diversity. Additionally, the paper provides a theoretical result with approximation guarantees, which might be useful outside this specific application.

This proposes studies ensemble construction from a combinatorial point of view. Analysis are conducted in function space, where the intuition are from and algorithm are developed. Specifically, it starts by formally defining the sub-modularity then proving f-divergence for a distribution and a mixture of kernels to be super-modular function. A greedy maximization of f-divergence is given and its implementation for approximating Bayesian posterior are presented. The experiments on synthetic data show the ability of covering the posterior of the proposed method, compared with the deep ensemble approach. It also conducts experiments on real-world data for OOD detection tasks. 

The paper proposes a novel and principled method to make Deep Ensemble Bayesian. It minimizes an f-divergence between the true posterior and a kernel density estimator fitting the functions associated with the ensemble candidates. The authors formulate the learning as a submodular problem and propose to use a greedy approach to solve it. The resultant learning objective consists of a novel diversity term. Some out-of-distribution detection results in computer vision demonstrate the efficacy of the proposed method.