Correlated Parameters to Accurately Measure Uncertainty in Deep Neural Networks

TL;DR

This paper introduces a Bayesian deep learning method using correlated variational parameters to improve uncertainty estimation and robustness against overfitting, demonstrated on MNIST and CIFAR-10 datasets.

Contribution

It proposes a novel variational inference approach with correlated parameters for Bayesian neural networks, reducing complexity while enhancing uncertainty quantification.

Findings

Effective uncertainty estimation demonstrated on benchmark datasets

Robustness to overfitting shown in experiments

Fewer parameters needed compared to traditional Bayesian methods

Abstract

In this article a novel approach for training deep neural networks using Bayesian techniques is presented. The Bayesian methodology allows for an easy evaluation of model uncertainty and additionally is robust to overfitting. These are commonly the two main problems classical, i.e. non-Bayesian, architectures have to struggle with. The proposed approach applies variational inference in order to approximate the intractable posterior distribution. In particular, the variational distribution is defined as product of multiple multivariate normal distributions with tridiagonal covariance matrices. Each single normal distribution belongs either to the weights, or to the biases corresponding to one network layer. The layer-wise a posteriori variances are defined based on the corresponding expectation values and further the correlations are assumed to be identical. Therefore, only a few additional parameters need to be optimized compared to non-Bayesian settings. The novel approach is successfully evaluated on basis of the popular benchmark datasets MNIST and CIFAR-10.

Figures

Tables

Table 1. TABLE I : Test errors of the trained models. The predictions of the Bayesian models are based on 200 200 200 samples from the corresponding variational distributions per test image.

Model	Test error
Frequentist (100 neurons)	0.84%
Bayesian (100 neurons)	0.7%
Frequentist (250 neurons)	0.7%
Bayesian (250 neurons)	0.61%

Table 2. TABLE II : Parameters ρ j subscript 𝜌 𝑗 \rho_{j} and τ j subscript 𝜏 𝑗 \tau_{j} corresponding to the variational distributions of the Bayesian models with 100 100 100 and 250 250 250 neurons in the first fully connected layer ( j = 1 , … , 4 𝑗 1 … 4 j=1,...,4 ).

Layer	Neurons	${\tau }_j$	${\tau }_{bj}$	${\rho }_j$	${\rho }_{bj}$
Convolutional 1	100	0.04	0.05	-0.44	-0.04
Convolutional 2	100	0.52	0.05	-0.25	-0.01
Fully connected 1	100	0.86	0.05	-0.15	0.01
Fully connected 2	100	0.05	0.04	-0.21	0.01
Convolutional 1	200	0.03	0.05	-0.44	0.03
Convolutional 2	200	0.35	0.05	-0.21	-0.01
Fully connected 1	200	2.02	0.06	-0.15	-0.01
Fully connected 2	200	0.06	0.05	-0.18	0.01

Table 3. TABLE III : Overview of quite certain and uncertain prediction results.

Model	Prediction result	Quite certain	Uncertain
100 neurons	correct	9668	262
100 neurons	wrong	14	56
250 neurons	correct	9587	352
250 neurons	wrong	8	53

Table 4. TABLE IV : Test errors of the trained models. The predictions of the Bayesian model are based on 200 200 200 samples from the corresponding variational distribution per test image.

Model	Test error
Frequentist	20.9%
Bayesian	20.32%

Table 5. TABLE V : Parameters ρ j subscript 𝜌 𝑗 \rho_{j} and τ j subscript 𝜏 𝑗 \tau_{j} of the variational distribution ( j = 1 , … , 4 𝑗 1 … 4 j=1,...,4 ).

Layer	${\tau }_j$	${\tau }_{bj}$	${\rho }_j$	${\rho }_{bj}$
Convolutional 1	0.04	0.05	-0.47	0.01
Convolutional 2	0.46	0.05	-0.49	-0.01
Convolutional 3	0.37	0.05	-0.41	0.01
Fully connected	0.14	0.05	-0.29	0.01

Table 6. TABLE VI : Overview of quite certain and uncertain prediction results for α = 95 % 𝛼 percent 95 \alpha=95\% and α = 99 % 𝛼 percent 99 \alpha=99\% .

Prediction result	$\alpha$	Quite certain	Uncertain
correct	95%	5894	2074
wrong	95%	403	1629
correct	99%	5281	2687
wrong	99%	250	1782

Equations

$$p(\mathbf{w}|\mathbf{y},\mathbf{X})=\frac{p(\mathbf{y}|\mathbf{w},\mathbf{X})p(\mathbf{w})}{\int p(\mathbf{y}|\mathbf{w},\mathbf{X})p(\mathbf{w})d\mathbf{w}}$$

$$D_{KL}(q_{\mathbf{\varphi }}(\mathbf{w})||p(\mathbf{w}|\mathbf{y},\mathbf{X}))={\mathbb{E}}_{q_{\mathbf{\varphi }}(\mathbf{w})}\left(\ln \frac{q_{\mathbf{\varphi }}(\mathbf{w})}{p(\mathbf{w}|\mathbf{y},\mathbf{X})}\right)$$

$$L_{VI}$$

$$=-\sum _{i=1}^{\beta }\left\{{\mathbb{E}}_{q_{\mathbf{\varphi }}(\mathbf{w})}\left[\ln \mathbf{f}({\mathbf{x}}_i;\mathbf{w}{)}_{y_i}\right]\right\}+D_{KL}(q_{\mathbf{\varphi }}(\mathbf{w})||p(\mathbf{w})).$$

$$L_{VI}^k=-\frac{1}{m}\sum _{i=1}^m\left\{\ln \mathbf{f}({\tilde{\mathbf{x}}}_i;{\mathbf{w}}_k{)}_{{\tilde{y}}_i}\right]\}+\frac{1}{\beta }{D}_{KL}(q_{\mathbf{\varphi }}(\mathbf{w})||p(\mathbf{w}))$$

$$p(y^{\ast }|{\mathbf{x}}^{\ast },\mathbf{y},\mathbf{X})$$

$$\approx \int p(y^{\ast }|\mathbf{w},{\mathbf{x}}^{\ast })q_{\mathbf{\varphi }}(\mathbf{w})d\mathbf{w}$$

$$\approx \frac{1}{N}\sum _{i=1}^N\mathbf{f}({\mathbf{x}}^{\ast };{\mathbf{w}}_i{)}_{y^{\ast }}$$

$$\mathbf{\Sigma }=\left(\begin{array}{ccccc}{a}_1& c_1& & & \\ c_1& a_2& c_2& & \\ & c_2& a_3& \ddots & \\ & & \ddots & \ddots & c_{n-1}\\ & & & c_{n-1}& a_n\end{array}\right)\in {\mathbb{R}}^{n\times n}$$

$$c_i^2<\frac{1}{4}{a}_i{a}_{i+1}\frac{1}{{\cos }^2\left(\frac{\pi }{n+1}\right)}i=1,...,n-1$$

$$\mathbf{\Sigma }$$

$$=\left(\begin{array}{ccccc}{a}_1& & & & \\ c_1& a_2& & & \\ & c_2& a_3& & \\ & & \ddots & \ddots & \\ & & & c_{n-1}& a_n\end{array}\right)\left(\begin{array}{ccccc}{a}_1& c_1& & & \\ & a_2& c_2& & \\ & & a_3& \ddots & \\ & & & \ddots & c_{n-1}\\ & & & & a_n\end{array}\right)$$

$$=\left(\begin{array}{ccccc}{a}_1^2& c_1{a}_1& & & \\ c_1{a}_1& c_1^2+a_2^2& c_2{a}_2& & \\ & c_2{a}_2& c_2^2+a_3^2& c_3{a}_3& \\ & & ⋮& & \\ & & c_{n-2}{a}_{n-2}& c_{n-2}^2+a_{n-1}^2& c_{n-1}{a}_{n-1}\\ & & & c_{n-1}{a}_{n-1}& c_{n-1}^2+a_n^2\end{array}\right)$$

$$c_{i-1}^2+a_i^2$$

$$c_0$$

$$\rho =\frac{c_i{a}_i}{\sqrt{{\tau }^2{m}_i^2}\sqrt{{\tau }^2{m}_{i+1}^2}}=\frac{c_i{a}_i}{{\tau }^2|m_i||m_{i+1}|}$$

$$\iff$$

$$(5)\iff a_i$$

$$(2)$$

$$\iff c_i^2=\frac{{\rho }^2{\tau }^4{m}_i^2{m}_{i+1}^2}{{\tau }^2{m}_i^2-c_{i-1}^2}$$

$$c_i=\frac{\rho {\tau }^2{m}_i{m}_{i+1}}{\sqrt{{\tau }^2{m}_i^2-c_{i-1}^2}}.$$

$$c_0$$

$$c_i$$

$$a_i$$

$$a_n$$

$$(\rho {\tau }^2|m_i||m_{i+1}|{)}^2<\frac{1}{4}{\tau }^2{m}_i^2{\tau }^2{m}_{i+1}^2\frac{1}{{\cos }^2\left(\frac{\pi }{n+1}\right)}$$

$$\iff$$

$${\mathbf{L}}_j:=\left(\begin{array}{ccccc}{a}_{j1}& & & & \\ c_{j1}& a_{j2}& & & \\ & c_{j2}& a_{j3}& & \\ & & \ddots & \ddots & \\ & & & c_{j,K_j-1}& a_{j{K}_j}\end{array}\right)$$

$$c_{j0}$$

$$c_{ji}$$

$$a_{ji}$$

$$a_{j{K}_j}$$

$$Q_{{\mathbf{\varphi }}_j}({\mathbf{w}}_j)=\mathcal{N}({\mathbf{w}}_j;{\mathbf{m}}_j,{\mathbf{\Sigma }}_j)$$

$${\tau }_j$$

$${\rho }_j$$

$${\mathbf{W}}_j$$

$$W_{j1}$$

Full text

I Introduction

Nowadays, due to excellent results obtained in many fields of applied machine learning including computer vision and natural language processing [1] the popularity of deep learning is increasing rapidly. One of the reasons can surely be found in the fact that Krizhevsky et al. [2] outperformed the competitors in the ImageNet Large Scale Visual Recognition Challenge 2012 by proposing a convolutional neural network (CNN) named AlexNet. While AlexNet includes eight layers, more recent architectures for image classification go even deeper [3, 4]. It is well known, that a feed-forward network with merely one hidden layer can approximate a broad class of functions abitrarily well. A mathematical more profound formulation of this statement, the so-called universal approximation theorem, was proven by Hornik et al. [5]. However, Liang and Srikant [6] could show that deep nets require exponential less parameters than shallow ones in order to achive a given degree of approximation. Possible applications of deep nets for computer vision include medical imaging, psychology, the automotive industry, finance and life sciences [7, 8, 9, 10, 11, 12].

Despite the large and ever-increasing number of real world use cases deep learning comes along with two restrictions which still limit its areas of application. The first restriction is that deep networks require a large amount of training data, otherwise they are prone to overfiting. The reason for this is the huge amount of parameters neural nets hold. Although deep nets require exponential less parameters than shallow ones, the remaining number is nevertheless very high. Thus, in many potential fields of application, where such an amount cannot be provided, deep learning is of limited use, or often even cannot be used. To counteract this problem commonly diverse regularization techniques are applied. Besides classical approaches, such as the penalization of the L2 norm or the L1 norm, stochastic regularization methods gain increasing attention. For instance, dropout [13] and dropconnect [14] count to these stochastic techniques. The first one randomly sets the activation of non-output neurons to zero during network training and the second one randomly sets network weights to zero. While dropout classically is interpreted as an efficient way of performing model averaging with neural networks, Gal and Ghahramani [15] as well as Kingma et al. [16] recently showed that it can also be considered as an application of Bayesian statistics. The second restriction deep networks struggle with is that prediction uncertainty cannot be measured. Especially, in the medical field or for self-driving vehicles it is essential that the prediction uncertainty can be determined [17]. In these areas of application a model which predicts on average quite well is not good enough. One has to know if the model is certain in its predictions or not, such that in the case of high uncertainty a human can decide instead of the machine. Please be aware of the fact that the probabilites obtained when running a deep net for a classification task should not be interpreted as the confidence of the model. As a matter of fact a neural net can guess randomly while returning a high class probability [18].

A possible strategy to overcome the restrictions classical deep learning has to deal with is applying Bayesian statistics. In so-called Bayesian deep learning the network parameters are treated like random variables and are not considered to be fixed deterministic values. In particular, an a priori distribution is assigned to them and updating the prior knowledge after observing traning data results in the so-called posterior distribution. The uncertainty in the network parameters can be directly translated in uncertainty about predictions. Further, Bayesian methods are robust to overfitting because of the built-in regularization due to the prior. Buntine and S. Weigend [19] were one of the first who presented approximate Bayesian methods for neural nets. Two years later Hinton and van Camp [20] already proposed the first variational methods. Variational methods try to approximate the true posterior distribution with another parameteric distribution, the so-called variational distribution. The approximation takes place due to an optimization of the parameters of the variational distribution. They followed the idea that there should be much less information in the weights than in the output vectors of the training cases in order to allow for a good generalization of neural networks. Denker and Lecun [21] as well as J.C. MacKay [22] used Laplace approximation in order to investigate the posterior distributions of neural nets. Neal [23] proposed and investigated hybrid Monte Carlo training for neural networks as a less limited alternative to the Laplace approximation. However, the approaches mentioned up to now are often not scalable for modern applications which go along with highly parameterized networks. Graves [24] was the first to show how variational inference can be applied to modern deep neural networks due to Monte Carlo integration. He used a Gaussian distribution with a diagonal covariance matrix as variational distribution. Blundell et al. [25] extended and improved the work of Graves [24] and also used a diagonal Gaussian to approximate the posterior. As already mentioned, Gal and Ghahramani [15] as well as Kingma et al. [16] showed that using the regularization technique dropout can also be considered as variational inference.

The indepenece assumptions going along with variational inference via diagonal Gaussians (complete independence of network parameters), or also going along with variational inference according to dropout (independence of neurons) are restrictive. Permitting an exchange of information between different parts of neural network architectures should lead to more accurate uncertainty estimates. Louizos and Welling [26] used a distribution over random matrices in order to define the variational distribution. Thus, they could reduce the amount of variance-related parameters that have to be estimated and further allow for an information sharing between the network weights. Note that in the diagonal Gaussian approach assigning one variance term to each random weight and one variance term to each random bias doubles the amount of parameters to optimize in comparison to frequentist deep learning. Consequently, network training becomes more complicated and computationally expensive.

It should be mentioned that variational Bayes is just a specific case of local $\alpha$-divergence minimization. According to Amari [27] the $\alpha$-divergence between two densities $p(\mathbf{w})$ and $q(\mathbf{w})$ is given by $D_{\alpha}(p(\mathbf{w})||q(\mathbf{w}))=\frac{1}{\alpha(1-\alpha)}\left(1-\int p(\mathbf{w})^{\alpha}q(\mathbf{w})^{(1-\alpha)}~{}d\mathbf{w}\right)$. Thus, the $\alpha$-divergence converges for $\alpha\rightarrow 0$ to the KL divergence [28] which is typically used in variational Bayes. It has been shown [29, 30] that an optimal choice of $\alpha$ is task specific and that non-standard settings, i.e. settings with $\alpha\neq 0$ can lead to better prediction results and uncertainty estimates.

However, in this paper we do not propose an optimal choice of $\alpha$. Rather, for the classical case $\alpha=0$ we will propose a good and easy to interpret variational distribution. For this task recent work from Posch et al. [31] is extended. They used a product of Gaussian distributions with specific diagonal covariance matrices in order to define the variational distribution. In particular, the a posteriori uncertainty of the network parameters is represented per network layer and depending on the estimated parameter expectation values. Therefore, only few additional parameters have to be optimized compared to classical deep learning and the parameter uncertainty itself can easily be analyzed per network layer. We extened this distribution by allowing network parameters to be correlated with each other. In particular, the diagonal covariance matrices are replaced with tridiagonal ones. Each tridiagonal matrix is defined in such a way that the correlations between neighboured parameters are identical. This way of treating network layers as units in terms of dependence allows for an easy analysis of the dependence between network parameters. Moreover, again only few additional parameters compared to classical deep learning need to be optimized, which guarantees that the difficulty of the network optimization does not increase significantly. Note that our extension allows for an exchange of information between different parts of the network and therefore should lead to more reliable uncertainty estimates. We have evaluated our approach on basis of the popular benchmark datasets MNIST [32] and CIFAR-10 [33]. The promising results can be found in Section IV.

II Background

In this section based on previous work [20, 24, 18, 34, 31] we briefly discuss how variational inference can be applied in deep learning. Since we are particularly interested in image classification the focus will be on networks designed for classification tasks. Note that the methodology also holds for regression models after some slight modifications.

Let $\mathbf{W}$ denote the random vector covering all parameters (weights and biases) of a given neural net $\mathbf{f}$. Further, let $p(\mathbf{w})$ denote the density used to define a priori knowledge regarding $\mathbf{W}$. According to the Bayesian theorem the posterior distribution of $\mathbf{W}$ is given by the density

[TABLE]

where $\mathbf{X}=\{\mathbf{x}_{1},...,\mathbf{x}_{\beta}\}$ denotes a set of training examples and $\mathbf{y}=(y_{1},...,y_{\beta})^{T}$ holds the corresponding class labels. Note that the probability $p(\mathbf{y}|\mathbf{w},\mathbf{X})$ is given by the product $\prod_{i=1}^{\beta}\mathbf{f}(\mathbf{x}_{i};\mathbf{w})_{y_{i}}$ in accordance with the classical assumptions on stochastic independence and modeling in deep learning for classification. The integral in the above representation of $p(\mathbf{w}|\mathbf{y},\mathbf{X})$ is commonly intractable due to its high dimension $\beta$. Variational inference aims at approximating the posterior with another parametric distribution, the so-called variational distribution $Q_{\boldsymbol{\phi}}(\mathbf{w})$. To this end the so-called variational parameters $\boldsymbol{\phi}$ are optimized by minimizing the Kullback Leibler divergence (KL divergence)

[TABLE]

between the variational distribution and the posterior. Although the KL divergence is no formal distance measure (does not satisfy some of the requested axioms) it is a common choice to measure the similarity of probability distributions.

Since the posterior distribution is unknown the divergence $D_{KL}(q_{\boldsymbol{\phi}}(\mathbf{w})||p(\mathbf{w}|\mathbf{y},\mathbf{X}))$ cannot be minimized directly. Anyway, the minimization of $D_{KL}(q_{\boldsymbol{\phi}}(\mathbf{w})||p(\mathbf{w}|\mathbf{y},\mathbf{X}))$ is equivalent to the minimization of the so-called negative log evidence lower bound

[TABLE]

Thus, the optimization problem reduces to the minimization of the sum of the expected negative log likelihood and the KL divergence between the variational distribution and the prior with respect to $\boldsymbol{\phi}$. Inspired by stochastic gradient descent it is not unusual to approximate the expected values in $L_{VI}$ via Monte Carlo integration with one sample during network training. Note that the re-sampling in each training iteration guarantees that a sufficient amount of samples is drawn overall. Commonly, mini-batch gradient descent is used for optimization in deep learning. To take account of the resulting reduction of the number of training examples used in each iteration of the optimization the objective function has to be rescaled. Thus, in the $k$-th iteration the function to minimize is given by

[TABLE]

where $\mathbf{w}_{k}$ denotes a sample from $q_{\boldsymbol{\phi}}(\mathbf{w})$, $m$ denotes the mini-batch size, and $\tilde{\mathbf{x}}_{1},...,\tilde{\mathbf{x}}_{m},\tilde{y}_{1},...,\tilde{y}_{m}$ denote the mini-batch itself.

Summing up, optimizing a Bayesian neural net is quite similar to the optimization of a classical one. While in frequentist deep learning it is common to penalize the Euclidean norm of the network parameters in terms of regularization in Bayesian deep learning deviations of the variational distribution from the prior are penalized. In principal, the same loss function $L$ (cross entropy loss) is minimized, but with the crucial difference that network parameters have to be sampled since they are random.

In Bayesian deep learning predictions are based on the posterior predictive distribution, i.e. the distribution of a class label $y^{*}$ for a given example $\mathbf{x}^{*}$ conditioned on the observed data $\mathbf{y},\mathbf{X}$. The distribution can be approximated via Monte Carlo integration

[TABLE]

where $\mathbf{w}_{1},...,\mathbf{w}_{N}$ denote samples from $q_{\boldsymbol{\phi}}(\mathbf{w})$. Therefore, the class of an object $\mathbf{x}^{*}$ is predicted by computing multiple network outputs with parameters sampled from the variational distribution. Averaging the output vectors results in an estimate of the posterior predictive distribution, such that the a posteriori most probable class finally serves as prediction.

The a posteriori uncertainty in the network parameters $\mathbf{W}$ can directly be translated in uncertainty about the random network output $\mathbf{f}(\mathbf{x}^{*};\mathbf{W})$ and thus the posterior probability of an object $\mathbf{x}^{*}$ belonging to class $y^{*}$. Therefore, at first multiple samples $\mathbf{w}_{1},...,\mathbf{w}_{N}$ are drawn from the variational distribution $Q_{\boldsymbol{\phi}}(\mathbf{w})$. Then the corresponding network outputs $\mathbf{f}(\mathbf{x}^{*};\mathbf{w}_{1})_{y^{*}},...,\mathbf{f}(\mathbf{x}^{*};\mathbf{w}_{n})_{y^{*}}$ are determined. Finally, the empirical $\frac{\alpha}{2}$ and $(1-\frac{\alpha}{2})$ quantiles of these outputs provide an estimate of the $(1-\alpha)$ credible interval for the probability $p(y^{*}|\mathbf{x}^{*},\mathbf{y},\mathbf{X})$.

III Methodology

In this section we describe our novel approach. At first we give a formal definition of the variational distribution we use to approximate the posterior and, additionally, we propose a normal prior. Moreover, we report the derivatives of the approximation $L_{VI}^{k}$ of the negative log evidence lower bound (see Section II) with respect to the variational parameters, i.e. the learnable parameters. Finally, we present the pseudocode which is the basis of our implementation of the proposed method.

III-A Variational Distribution and Prior

Let $\mathbf{W}_{j}=(W_{j1},...,W_{jK_{j}})^{T}$ denote the random weights of layer $j$ $(j=1,...,d)$. Further, let $\mathbf{B}_{j}=(B_{j1},...,B_{jk_{j}})^{T}$ denote the corresponding random bias terms. The integers $K_{j}$ and $k_{j}$ denote the number of weights and the number of biases in layer $j$, respectively.

As already mentioned in Section I, we define the variational distribution as a product of multivariate normal distributions with tridiagonal covariance matrices. Applying variational inference to Bayesian deep learning presupposes that samples from the variational distribution can be drawn during network training as well as at the stage in which new samples are used for prediction. Especially at the training phase, it is essential that the random sampling can be reduced to the sampling from a (multivariate) standard normal distribution and appropriate affine-linear transformation of the drawn samples based on the learnable parameters. A direct sampling from the non-trivial normal distributions would mask the variational parameters and thus make it impossible to optimize them by gradient descent. Provided that a covariance matrix is positive definite (in general covariance matrices are only positive semidefinte) there exists a unique Cholesky decomposition of the matrix which can be used for this task. Note that for each real-valued symmetric positive-definite square matrix a unique decomposition (Cholesky decomposition) of the form $\mathbf{\Sigma}=\mathbf{L}\mathbf{L}^{T}$ exists, where L is a lower triangular matrix with real and positive diagonal entries. Thus, we are interested in symmetric tridiagonal matrices which always stay positive definite no matter how the corresponding learnable parameters are adjusted during network training. The first thing required is a criterion for the positive definitness of for our purposes approriate tridiagonal matrices. Andelić and Fonseca [35] gave the following sufficient condition for positive definiteness of tridiagonal matrices: Let

[TABLE]

a symmetric tridiagonal matrix with positive diagonal entries. If

[TABLE]

then $\mathbf{\Sigma}$ is positive definite. Consider now a matrix $\mathbf{\Sigma}\in\mathbb{R}^{n\times n}$ defined as follows:

[TABLE]

If $\mathbf{\Sigma}$ satisfies condition $(1)$ and has positive diagonal entries, the matrix $\mathbf{L}$ defines its Cholesky decomposition and further $\mathbf{\Sigma}$ is a valid covariance matrix since every real, symmetric, and positive semidefinite square matrix defines a valid covariance matrix. As in the work of Posch et al. [31] we define the variances as multiples of the corresponding expectation values, denoted by $m_{1},...,m_{n}\in\mathbb{R}\setminus\{0\}$

[TABLE]

where $\tau\in\mathbb{R}^{+}$. Defining the variances proportional to the expectation values allows for a useful specification of them. This specification requires, besides the expectation values, only one additional variational parameter. Moreover, we want the correlations to be identical, which leads to the following covariances $c_{i}a_{i}$

[TABLE]

for $i=1,...,n-1$. By rearranging Equations $(2)$ and $(5)$ one obtains a recursive formula for the elements of the matrix $L$:

[TABLE]

[TABLE]

Note that Equation $(8)$ for instance is satisfied for

[TABLE]

By defining the $c_{i}$ this way one does not end up by the Cholseky decomposition which assumes the diagonal elements $a_{i}$ of $L$ to be positive and thus by $(6)$ also assumes the $c_{i}$ to be positive. Taking the absolute values of the $c_{i}$ according to Equation $(9)$ would result in the Cholesky decomposition, but is not necessary for our purposes and therefore not done. Thus, the elements of $\mathbf{L}$ are recursively defined as

[TABLE]

Note that the matrix $\mathbf{\Sigma}$ defined by the Equations $(10-13)$ satisfies condition $(1)$ iff

[TABLE]

Thus, provided that condition $(16)$ holds there exists a unique Cholesky decomposition of $\mathbf{\Sigma}$ which again guarantees that the $c_{i}$ according to $(11)$ are well defined and, further, also that $a_{n}$ according to $(13)$ is well defined.

Using the considerations above the variational distribution can finally be defined as follows. In a first step lower bidiagonal matrices $\mathbf{L}_{j}$ are specified for $j=1,...,d$:

[TABLE]

[TABLE]

The variational distribution of the weights of the $j$-th layer is then defined as a multivariate normal distribution

[TABLE]

with expected value $\mathbf{m}_{j}$ and a tridiagonal covariance matrix $\mathbf{\Sigma}_{j}=\mathbf{L}_{j}\mathbf{L}_{j}^{T}$. According to the considerations above, the variances of the normal distribution are given by $\tau_{j}^{2}m_{ji}^{2}$ ($i=1,...,K_{j}$) and the covariances are given by $\rho_{j}\tau_{j}^{2}|m_{ji}||m_{j,i+1}|$ ($i=1,...,K_{j}-1$). This again implies that the correlations are all the same and given by the parameter $\rho_{j}$. Since the parameter $\tau_{j}$ regulates the variances of the distribution it should not take negative values during optimization. To guarantee this it is reparameterized with help of the softplus function

[TABLE]

Moreover, the parameter $\rho_{j}$ should lie in the interval $(-\frac{1}{2},\frac{1}{2})$ to ensure that the matrix $\mathbf{\Sigma}_{j}$ is positive definite. In deep learning dimensions are commonly high, such that the approximation in Equation $(16)$ can be considered as valid. The following reparameterization ensures that the desired property holds:

[TABLE]

In addition, the diagonal entries of $\mathbf{\Sigma}_{j}$ have to be non-zero to ensure positive definiteness, which again implies that each component of $\mathbf{m}_{j}$ has to be non-zero. We decide to set values which are not significantly different from zero to small random numbers in our implementation instead of introducing another reparameterization. Finally, $\mathbf{m}_{j}\in\mathbb{R}^{K_{j}}\setminus\{\mathbf{0}\}$, $\delta_{j}þ\in þ\mathbb{R}$, and $\gamma_{j}\in\mathbb{R}$ can be summarized as the variational parameters $\boldsymbol{\phi}_{j}$ corresponding to the weights of the $j$-th network layer.

One can easily sample from a random vector $\mathbf{W}_{j}$ belonging to this distribution using samples from a standard normal distribution $\mathcal{N}(0,1)$ since it can be written as

[TABLE]

Note that Equation $(24)$ can also be written as:

[TABLE]

The layer-wise variational distributions of the bias terms denoted by $q_{\boldsymbol{\phi}_{bj}}(\mathbf{b}_{j})$ are defined completely analogous to those of the weights. Assuming independence of the layers as well as independence between weights and biases, the overall variational distribution is given by

[TABLE]

where $\boldsymbol{\phi}_{j}=\{\mathbf{m}_{j}$, $\delta_{j},\gamma_{j}\},\boldsymbol{\phi}_{bj}=\{\mathbf{m}_{bj},\delta_{bj},\gamma_{bj}\}$, $q_{\boldsymbol{\phi}_{j}}(\mathbf{w}_{j})$ denotes the density of $\mathcal{N}(\mathbf{m}_{j},\mathbf{\Sigma}_{j})$, $q_{\boldsymbol{\phi}_{bj}}(\mathbf{b}_{j})$ denotes the density of $\mathcal{N}(\mathbf{m}_{bj},\mathbf{\Sigma}_{bj})$, and $\mathbf{w}$ is a vector including all weights and all biases.

We define the a priori distribution completely analogous to Posch et al. [31]. In particular, its density is given by:

[TABLE]

where $p(\boldsymbol{w}_{j})$ denotes the density of $N(\boldsymbol{\mu}_{j},\zeta_{j}^{2}\mathbf{I}_{K_{j}})$ and $p(\boldsymbol{b}_{j})$ denotes the density of $N(\boldsymbol{\mu}_{bj},\zeta_{bj}^{2}\mathbf{I}_{k_{j}})$.

III-B Kullback Leibler Divergence

The fact that the variational distribution as well as the prior factorize simplifies the computation of the Kullback Leibler divergence. Indeed, the overall divergence is given by the sum of the layer-wise divergences (for further details refer to Posch et al. [31]):

[TABLE]

Thus, computing the overall divergence can be reduced to compute $D_{KL}(q_{\boldsymbol{\phi}_{j}}(\mathbf{w}_{j})||p(\mathbf{w}_{j}))$ for fixed $j\in\{1,...,K_{j}\}$, since the remaining divergences compute completely analogously (only the indices differ). According to Hershey and Olsen [36] the KL divergence between two $p$-dimensional normal distributions, given by $H(\mathbf{x})=\mathcal{N}(\mathbf{x};\boldsymbol{\mu}_{h},\mathbf{\Sigma}_{h})$ and $G(\mathbf{x})=\mathcal{N}(\mathbf{x};\boldsymbol{\mu}_{g},\mathbf{\Sigma}_{g})$, computes as

[TABLE]

Thus, the determinant of the covariance matrix $\mathbf{\Sigma}_{j}$ is required for the computation of the Kullback Leibler divergence $D_{KL}(q_{\boldsymbol{\phi}_{j}}(\mathbf{w}_{j})||p(\mathbf{w}_{j}))$. Using basic properties of determinants $|\mathbf{\Sigma}_{j}|$ computes as follows for fixed $j$:

[TABLE]

Using $(28)$ and $(29)$ $D_{KL}(q_{\boldsymbol{\phi}_{j}}(\mathbf{w}_{j})||p(\mathbf{w}_{j}))$ then reads

[TABLE]

where $c$ always denotes an additive constant.

III-C Derivatives

Commonly neural networks are optimized via mini-batch gradient descent. Thus, in order to train a neural net $\mathbf{f}(\cdot,\mathbf{w})$ according to our novel approach the partial derivatives of the approximation $L_{VI}^{k}$ of the negative log evidence lower bound described in Section II with respect to the variational parameters $\boldsymbol{\phi}_{j}=\{\mathbf{m}_{j}$, $\delta_{j},\gamma_{j}\},\boldsymbol{\phi}_{bj}=\{\mathbf{m}_{bj},\delta_{bj},\gamma_{bj}\}$ are required. In particular, the partial derivatives of the loss function $L$ typically used in deep learning and the partial derivatives of the Kullback Leibler divergence between prior and variational distribution have to be computed. Note, that the loss function equals the negative log likelihood of the data and is given by the cross-entropy loss in the case of classification and by the Euclidean loss in the case of regression. Thus, $L$ depends on the network $\mathbf{f}$ itself, with parameters sampled from the variational distribution $q_{\boldsymbol{\phi}}(\mathbf{w})$. With the help of the multivariate chain rule the required partial derivatives of $L$ can be computed based on the classical derivatives used in non-Bayesian deep learning:

[TABLE]

Equations $(31-33)$ only deal with the derivatives of $L$ with respect to the variational parameters belonging to the network weights. Completely analogous equations hold for the bias terms. In the sequel we focus on the derivatives of the weights, since the derivatives for the biases are obviously of the same form. Note that for a given sample $\mathbf{w}$ from the variational distribution the layer-wise derivatives $\frac{\partial L}{\partial\mathbf{w}_{j}},\frac{\partial L}{\partial\mathbf{b}_{j}}$ ($j=1,...,d$) are computed as in non-Bayesian deep learning. Thus, the problem of finding closed form expressions for the required derivatives of $L$ reduces to the problem of finding these expressions for the $\mathbf{w}_{j}$’s and the $\mathbf{b}_{j}$’s. Taking account of the Equations $(25-27)$ the needed derivatives of the weights can be expressed in terms of the corresponding derivatives of the $c_{ji}$ and the $a_{ji}$

[TABLE]

where the index $j$ lies in the set $\{1,...,d\}$, while for a given $j$ the indices $i$ and $k$ lie in the set $\{1,...,K_{j}\}$. The derivatives of the $c_{ji}$ ($i=1,...,K_{j}-1$) with respect to the variational parameters are given by

[TABLE]

where $u:=\tau_{j}^{2}m_{ji}^{2}-c_{j,i-1}^{2}$, $v:=\rho_{j}\tau_{j}^{2}m_{j,i+1}$, $w_{\delta}:=\operatorname{exp}(\delta_{j})/[1+\operatorname{exp}(\delta_{j})]$, and $w_{\gamma}:=\operatorname{exp}(-\gamma_{j})/[1+\operatorname{exp}(-\gamma_{j})]^{2}$, and obviously each derivative of $c_{j0}$ equals zero. Moreover, the derivatives of the $a_{ji}$ ($i=1,...,K_{j}-1$) with respect to the variational parameters are given by:

[TABLE]

In addition, the derivatives of $a_{jK_{j}}$ with respect to the variational parameters are given by

[TABLE]

where $y:=\tau_{j}^{2}m_{jK_{j}}^{2}-c_{j,K_{j}-1}^{2}$.

Finally, the partial derivatives of the KL divergence $D_{KL}(q_{\boldsymbol{\phi}}(\mathbf{w})||p(\mathbf{w}))$ (abbreviated with $D_{KL}$) with respect to the variational parameters are given by:

[TABLE]

For reasons of brevity the calculation process corresponding to the derivatives presented is not reported.

III-D Implementation and Pseudocode

We have implemented the proposed apporach by modifying and extending the popular open-source deep learning framework Caffe, see [37]. In particular, our implementation includes a Bayesian version of the classical inner product layer as well as the convolutional layer. This allows for the training of (deep) multilayer perceptrons (MLPs) and, moreover, convolutional neural networks (CNNs) according to our novel approach. Up to now we have not parallelized our code such that it can run on GPU. This is left for future research and will enable the training of state-of-the-art CNNs for image classification in a reasonable amount of time. The pseudocode which was starting point of our implementation is presented below. The code shows how a classical, i.e. frequentist, inner product, or convolutional layer can be extended in order to fit with the methodology presented.

Besides parameters, for which obviously initial values are required, the presented pseudo code also asks for the two additional parameters $\nu$ and $\kappa$. Empirically, we discovered that in practical applications it can be helpful to use another penalization strength of the KL divergence than the fixed one proposed in Section II. In particular, the derivatives of the KL divergence do not suffer from the vanishing gradients problem in contrast to the derivatives of the loss function. For this reason, reducing the penalization strength of the divergence by a fixed factor, which is equivalent to reducing the learning rate of the divergence, might be a good decision. Moreover, again, empirically we discovered that the derivatives with respect to the $\gamma_{j}$ are quite small which slows down the learning process. To overcome this problem the learning rate multiplier $\kappa$ can be used.

In the forward pass at first variational parameters which hold “critical” values are assigned similar but less “critical” ones. This way of proceeding should guarantee that no numerical issues occur during the training process. Thus, very small expectation values, variances (smaller than $0.01$ times abs(corresponding expectation value) in absolute values), and correlations (smaller than $0.01$ in absolute values) are replaced, as well as correlations which are nearly given by $0.5$, or $-0.5$. Recall that the correlations have to stay in the interval $(-\frac{1}{2},\frac{1}{2})$ in order for the covariance matrices to be positive definite. After the checkings for numerical stability the weights and the biases of the network are sampled from the variational distribution. The sampled parameters are then used in place of the deterministic ones in non-Bayesian deep learning in order to perform the classical forward pass.

In the backward pass at first the derivatives of the loss function with respect to the sampled weights and biases are computed. This is done completely analogous as in classical deep learning, but the sampled parameters are used in place of the deterministic ones. In the next steps the derivatives of the sampled network parameters with respect to the variational parameters are computed, and the derivatives of the KL divergence with respect to the variational parameters are calculated. Finally, an appropriate merging of all the computed derivatives results in the derivatives of the approximation $L_{VI}^{k}$ of the negative log evidence lower bound. These derivatives are used to update the variational parameters according to some learning schedule.

IV Experiments

This section investigates how well the proposed approach performs on real world datasets. On the one hand the prediction accuracy and on the other hand the quality of the uncertainty information provided are of particular interest. The popular benchmark datasets MNIST [32] and CIFAR-10 [33] form basis of the evaluations.

IV-A MNIST

In this section, the performance of the proposed method is evaluated based on the MNIST dataset. This dataset contains $70,000$ grayscale images of handwritten digits ($0-9$). The complete dataset is partitioned in a training dataset which counts $60,000$ images and a test set of $10,0000$ examples. Each image is of size $28\times 28$.

The architecture chosen for the performance evaluation is the popular LeNet proposed by Lecun et al. [38]. This CNN mainly consists of two convolutional layers and, further, two fully connected layers. The special version of LeNet used in this paper has the following additional specifications: The ReLU activation function is assigned to the first fully connected layer, while the other layers simply use the identity function as activation function. Moreover, the first convolutional layer includes $20$ and the second one includes $50$ $5\times 5$ kernels. Max-pooling with kernel size $2\times 2$ and stride of $2$ is applied after both convolutional layers. The number of neurons of the second fully connected layer is determined by the number of possible classes and thus given by $10$. In contrast to the second fully connected layer the number of neurons in the first fully connected layer can be freely specified. In our experiments we set this number once to $100$ and once to $250$.

To get an idea how well the proposed Bayesian approach predicts compared to the frequentist one, we train LeNet according to both approaches. Mini-batch gradient descent with a batch size of $64$ is chosen for optimization. Further, a decreasing learning rate policy is selected. In particular, the learning rate in the $k$-th iteration is specified as $0.01\cdot(1+0.0001\cdot i)^{−0.75}$. Moreover, momentum is applied with a value of $0.9$. The total amount of iterations is set to $100,000$. In terms of regularization in the frequentist approach a penalization of the Euclidean norm with a penalization strength of $0.0005$ takes place. Further, dropout is applied after the first fully connected layer with a dropping rate given by $0.5$. These regularization techniques are not applied in the Bayesian approach which is naturally regularized due to the sampling from the variational distribution and the penalization of the KL divergence during network training. However, in Bayesian deep learning the parameters of the prior have to be specified. We set the expectation values of all network parameters to [math] and the variances to $1$. There is no a priori knowledge available such that the prior can merely be used to guarantee that network parameters do not diverge. Moreover, for each network layer we set the parameters $\nu$ and $\kappa$ described in Section III-D to $1/(60,000\cdot 100)$ and $50$, respectively. Thus, the penalization strength of the KL divergence is reduced by a factor of $100$ due to the reasons described in Section III-D.

For the model with $100$ neurons in the first fully connected layer both the Bayesian and the frequentist training process are visualized in the Figures 2 and 2. One can see that the loss decreases in a quite similar way in both processes. This is an interesting result since in the work of Posch et al. [31] the Bayesian loss fluctuates heavily compared to the classical one. Thus, the learning of correlations between network parameters (which is our extension of their work) enables a smoother network training. Note that the test error plotted in Figure 2 is just a rough approximation of the true one, which is based on one sample from the variational distribution per test image. Usually, the corresponding true errors are significantly lower. However, the rough approximation suffices to monitor network training. A more accurate and thus also computationally more expensive approximation based on multiple samples can be made at the end of the training process.

The finally achieved test errors are given in Table I. Note that the predictions of the Bayesian models are based on $200$ samples from the corresponding variational distributions per test image. One can see that the Bayesian models perform a little bit better than their frequentist analogues. The reason for the superior results of the Bayesian models can be found in their natural robustness against overfitting. It can be hard to regularize frequentist models in such a way that they obtain a robustness comparable to Bayesian models.

In Table II the correlations $\rho_{j}$ and the variance determining parameters $\tau_{j}$ of the two Bayesian models can be found ($j=1,...,4$). One can see that there is a low a posteriori uncertainty about the bias terms and further that they are nearly uncorrelated. This is a plausible result since the number of bias terms is small compared to the number of weights. The a posteriori uncertainty of the weights differs significantly from layer to layer. While the first convolutional layer and the second fully connected layer go along with a low a posteriori uncertainty about the network weights, the other layers show a high uncertainty. For the network with $250$ neurons in the first fully connected layer the standard deviation of the posterior distribution of a weight is given by two times the corresponding expectation value. Note that layers with low uncertainty are exactly the ones which directly act on the network input and the network output. The correlations of the weights are all negative.

As already mentioned in Section II the uncertainties in the network parameters of Bayesian nets can directly be translated in uncertainty about the predictions. In particular, credible intervals for the probability that an image shows an object of a given class can be estimated by computing multiple neural network outputs with weights sampled from the variational distribution. Figures 3, 4, and 5 show boxplots of respectively $200$ random network outputs corresponding to three representative images from the test dataset. The network which was used to produce these figures is the one with $100$ neurons in the first fully connected layer. In Figure 3 there are no boxes at all since the network is very sure about its correct prediction. This is the case for most of the correctly classified images. The Figures 4 and 5 reflect the common behavior of the net in case of incorrect classification results. Either the network has difficulties to decide between two classes where one of the two is the true one, or it is completely uncertain what class to predict. To quantify the quality of the prediction uncertainty of the Bayesian models we estimate $95\%$ credible intervals of the component-wise network outputs for all test images. The estimates are based on $200$ random network outputs, respectively. A prediction result is considered as quite certain if the credible interval of the predicted class does not overlap with the intervals of the other classes. In the other case a prediction result is considered as uncertain. Table III summarizes the corresponding results. One can see that the models are most of the times quite certain about correct prediction results while they are usually uncertain in wrong prediction results. This shows that the uncertainty information obtained from the proposed Bayesian approach is very good. In Figure 6 all the test images can be found where the model with $250$ neurons in the first fully connected layer was quite certain about its wrong prediction.

IV-B CIFAR-10

In this section, the performance of the proposed method is evaluated based on the CIFAR-10 dataset. This dataset consists of $60,000$ RGB images in $10$ classes (airplane, automobile, bird, cat, deer, dog, frog, horse, ship, truck). The complete dataset is partitioned into a training dataset which counts $50,000$ images and a test set of $10,000$ examples. Each image is of size $32\times 32$. The architecture chosen for the performance evaluation is the one included under the name CIFAR10_full in the Caffe framework [37]. This CNN mainly consists of three convolutional layers followed by one fully connected layer. For further details take a look at the model definition available in the Caffe framework.

As in Section IV-A we train the network according to the frequentist approach and according to the proposed Bayesian approach. Mini-batch gradient descent with a batch size of $100$ is chosen for optimization. Further, a fixed learning rate policy with a learning rate of $0.001$ is selected. The total amount of training iterations is given by $100,000$ for the frequentist approach and by $40,000$ for the Bayesian approach. A penalization of the Euclidean norm with a penalization strength of $0.004$ for the convolutional layers and a penalization strength of $1$ for the fully connected layer is used for regularization of the frequentist net. In the Bayesian case the prior distribution is used for regularization. Therefore, we assign independent normal distributions with zero mean (see Section III-A) to the network parameters. The standard deviations of these distributions are specified as $1$ for the convolutional layers and specified as $0.05$ for the fully connected layer. Thus, the fully connected layer gets stronger regularized than the other ones, as in the frequentist case. Moreover, for each network layer we set the parameters $\nu$ and $\kappa$ described in Section III-D to $1/(50,000\cdot 10)$ and $50$. Thus, the penalization strength of the KL divergence is reduced by a factor of $10$ due to the reasons described in Section III-D.

The Figures 8 and 8 show the frequentist and the Bayesian training process, respectively.

The achieved test errors are given in Table IV. Note that the predictions of the Bayesian model are based on $200$ samples from the corresponding variational distribution per test image. The Bayesian model performs slightly better than the frequentist one due to the natural robustness of Bayesian models to overfitting.

In Table V the correlations $\rho_{j}$ and the variance determining parameters $\tau_{j}$ of the Bayesian model can be found ($j=1,...,4$). The parameters can be interpreted as in Section IV-A.

To quantify the quality of the prediction uncertainty of the Bayesian model we estimate the $\alpha=95\%$ credible intervals and further the $\alpha=99\%$ credible intervals of the component-wise network outputs for all test images. The estimates are based on $200$ random network outputs, respectively. The terms quite certain and uncertain are to be understood as in Section IV-A. Then Table III summarizes in how many of the correct classification results the model is certain about its prediction and further summarizes the same for the wrong classification results. The quality of the uncertainty information is good, since for the correct classifications the majority of the predictions is considered as quite certain in contrast to the incorrect classifications where the opposite is the case. The amount of uncertainty can be regulated by the parameter $\alpha$. Assigning a higher value to $\alpha$ results in higher uncertainty and thus also in a lower number of quite certainly deemed wrong predictions. Thus, an optimal choice of $\alpha$ is dependent on the application of interest. For instance, for self-driving cars a value of $\alpha$ near one would be appropriate.

V Conclusion

We presented the first Bayesian approach to deep learning that allows to learn correlations between network parameters while introducing only few additional parameters to be optimized. In particular, we approximated the intractable posterior of the network parameters with a product of Gaussian distributions with tridiagonal covariance matrices. These distributions are defined in such a way, that the variances are multiples of the expectation values and the correlations belonging to a given distribution are identical. The novel approach was evaluated on basis of the popular benchmark datasets MNIST and CIFAR-10. Superior prediction accuracies compared to well regularized frequentist models show that the proposed method is more robust to overfitting than classical ones. Further, we showed that accurate uncertainty information about network predictions can be computed. This possibility of measuring prediction uncertainty should have a significant impact for real world applications such as self-driving vehicles. Finally, network parameter uncertainties and dependencies can readily be interpreted per layer, due to the fact that only few additional parameters are required by our method.

Acknowledgments

The publication is one of the results of the project iDev40 (www.idev40.eu). The iDev40 project has received funding from the ECSEL Joint Undertaking (JU) under grant agreement No 783163. The JU receives support from the European Union’s Horizon 2020 research and innovation programme. It is co-funded by the consortium members, grants from Austria, Germany, Belgium, Italy, Spain and Romania.