# Correlated Parameters to Accurately Measure Uncertainty in Deep Neural Networks

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.

## Authors

• 4 publications
• 4 publications
• ### Variational Inference to Measure Model Uncertainty in Deep Neural Networks

We present a novel approach for training deep neural networks in a Bayes...
02/26/2019 ∙ by Konstantin Posch, et al. ∙ 0

• ### Stochastic Bayesian Neural Networks

Bayesian neural networks perform variational inference over weights but ...
06/27/2020 ∙ by Abhinav Sagar, et al. ∙ 1

• ### Scalable Bayesian neural networks by layer-wise input augmentation

We introduce implicit Bayesian neural networks, a simple and scalable ap...
10/26/2020 ∙ by Trung Trinh, et al. ∙ 0

• ### A multilayer exponential random graph modelling approach for weighted networks

This paper introduces a new modelling approach to analyse weighted netwo...
11/16/2018 ∙ by Alberto Caimo, et al. ∙ 0

• ### Structured and Efficient Variational Deep Learning with Matrix Gaussian Posteriors

We introduce a variational Bayesian neural network where the parameters ...
03/15/2016 ∙ by Christos Louizos, et al. ∙ 0

• ### Bayesian Convolutional Neural Networks with Bernoulli Approximate Variational Inference

Convolutional neural networks (CNNs) work well on large datasets. But la...
06/06/2015 ∙ by Yarin Gal, et al. ∙ 0

• ### Local Critic Training for Model-Parallel Learning of Deep Neural Networks

This paper proposes a novel approach to train deep neural networks in a ...
05/03/2018 ∙ by Hojung Lee, et al. ∙ 0

##### This week in AI

Get the week's most popular data science and artificial intelligence research sent straight to your inbox every Saturday.

## 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 -divergence minimization. According to Amari [27] the -divergence between two densities and is given by . Thus, the -divergence converges for to the KL divergence [28] which is typically used in variational Bayes. It has been shown [29, 30] that an optimal choice of is task specific and that non-standard settings, i.e. settings with can lead to better prediction results and uncertainty estimates.

However, in this paper we do not propose an optimal choice of . Rather, for the classical case 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 denote the random vector covering all parameters (weights and biases) of a given neural net . Further, let denote the density used to define a priori knowledge regarding . According to the Bayesian theorem the posterior distribution of is given by the density

where denotes a set of training examples and holds the corresponding class labels. Note that the probability is given by the product in accordance with the classical assumptions on stochastic independence and modeling in deep learning for classification. The integral in the above representation of is commonly intractable due to its high dimension . Variational inference aims at approximating the posterior with another parametric distribution, the so-called variational distribution . To this end the so-called variational parameters

are optimized by minimizing the Kullback Leibler divergence (KL divergence)

 DKL(qϕ(w)||p(w|y,X))=Eqϕ(w)(lnqϕ(w)p(w|y,X))

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 cannot be minimized directly. Anyway, the minimization of is equivalent to the minimization of the so-called negative log evidence lower bound

 LVI =−Eqϕ(w)[lnp(y|w,X)]+DKL(qϕ(w)||p(w)) =−β∑i=1{Eqϕ(w)[lnf(xi;w)yi]}+DKL(qϕ(w)||p(w)).

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

. Inspired by stochastic gradient descent it is not unusual to approximate the expected values in

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 -th iteration the function to minimize is given by

 LkVI=−1mm∑i=1{lnf(~xi;wk)~yi]}+1βDKL(qϕ(w)||p(w))

where denotes a sample from , denotes the mini-batch size, and 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

(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

for a given example conditioned on the observed data . The distribution can be approximated via Monte Carlo integration

 p(y∗|x∗,y,X) =∫p(y∗|w,x∗)p(w|y,X) dw ≈∫p(y∗|w,x∗)qϕ(w) dw ≈1NN∑i=1f(x∗;wi)y∗

where denote samples from . Therefore, the class of an object 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 can directly be translated in uncertainty about the random network output

and thus the posterior probability of an object

belonging to class . Therefore, at first multiple samples are drawn from the variational distribution . Then the corresponding network outputs are determined. Finally, the empirical and quantiles of these outputs provide an estimate of the credible interval for the probability .

## 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 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 denote the random weights of layer . Further, let denote the corresponding random bias terms. The integers and denote the number of weights and the number of biases in layer , 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

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

 Σ=⎛⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜⎝a1c1c1a2c2c2a3⋱⋱⋱cn−1cn−1an⎞⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟⎠∈Rn×n

a symmetric tridiagonal matrix with positive diagonal entries. If

 c2i<14aiai+11cos2(πn+1)   i=1,...,n−1 (1)

then is positive definite. Consider now a matrix defined as follows:

 Σ :=LLT =⎛⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜⎝a1c1a2c2a3⋱⋱cn−1an⎞⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟⎠⎛⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜⎝a1c1a2c2a3⋱⋱cn−1an⎞⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟⎠ =⎛⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜⎝a21c1a1c1a1c21+a22c2a2c2a2c22+a23c3a3⋮cn−2an−2c2n−2+a2n−1cn−1an−1cn−1an−1c2n−1+a2n⎞⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟⎠

If satisfies condition and has positive diagonal entries, the matrix defines its Cholesky decomposition and further 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

 c2i−1+a2i :=τ2m2i   i=1,...,n (2) c0 :=0 (3)

where . 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

 ρ=ciai√τ2m2i√τ2m2i+1=ciaiτ2|mi||mi+1| (4) ⇔ ciai=ρτ2|mi||mi+1| (5)

for . By rearranging Equations and one obtains a recursive formula for the elements of the matrix :

 (5)⇔ai =ρτ2|mi||mi+1|ci (6)
 (2) ⇔c2i−1+ρ2τ4m2im2i+1c2i=τ2m2i (7) ⇔c2i=ρ2τ4m2im2i+1τ2m2i−c2i−1 (8)

Note that Equation for instance is satisfied for

 ci=ρτ2mimi+1√τ2m2i−c2i−1. (9)

By defining the this way one does not end up by the Cholseky decomposition which assumes the diagonal elements of to be positive and thus by also assumes the to be positive. Taking the absolute values of the according to Equation would result in the Cholesky decomposition, but is not necessary for our purposes and therefore not done. Thus, the elements of are recursively defined as

 c0 :=0 (10) ci :=ρτ2mimi+1√τ2m2i−c2i−1   i=1,...,n−1 (11) ai :=ρτ2|mi||mi+1|ci   i=1,...,n−1 (12) an :=√τ2m2n−c2n−1. (13)

Note that the matrix defined by the Equations satisfies condition iff

 (ρτ2|mi||mi+1|)2<14τ2m2iτ2m2i+11cos2(πn+1) (14) ⇔ ρ2<141cos2(πn+1) (15) ⇔ ρ∈⎛⎜ ⎜ ⎜ ⎜⎝−121√cos2(πn+1),121√cos2(πn+1)⎞⎟ ⎟ ⎟ ⎟⎠≈for large n(−12,12). (16)

Thus, provided that condition holds there exists a unique Cholesky decomposition of which again guarantees that the according to are well defined and, further, also that according to is well defined.

Using the considerations above the variational distribution can finally be defined as follows. In a first step lower bidiagonal matrices are specified for :

 Lj:=⎛⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜⎝aj1cj1aj2cj2aj3⋱⋱cj,Kj−1ajKj⎞⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟⎠ (17)
 cj0 :=0 (18) cji :=ρjτ2jmjimj,i+1√τ2jm2ji−c2j,i−1   i=1,...,Kj−1 (19) aji :=ρjτ2j|mji||mj,i+1|cji   i=1,...,Kj−1 (20) ajKj :=√τ2jm2jKj−c2j,Kj−1 (21)

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

 Qϕj(wj)=N(wj;mj,Σj)

with expected value and a tridiagonal covariance matrix . According to the considerations above, the variances of the normal distribution are given by () and the covariances are given by (). This again implies that the correlations are all the same and given by the parameter . Since the parameter 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

 τj =ln(1+exp(δj))>0. (22)

Moreover, the parameter should lie in the interval to ensure that the matrix is positive definite. In deep learning dimensions are commonly high, such that the approximation in Equation can be considered as valid. The following reparameterization ensures that the desired property holds:

 ρj =11+exp(−γj)−12∈(−12,12) (23)

In addition, the diagonal entries of have to be non-zero to ensure positive definiteness, which again implies that each component of 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, , , and can be summarized as the variational parameters corresponding to the weights of the -th network layer.

One can easily sample from a random vector belonging to this distribution using samples from a standard normal distribution since it can be written as

 Wj =mj+LjXj   with   Xj∼N(0Kj,IKj). (24)

Note that Equation can also be written as:

 Wj1 =mj1+aj1Xj1 (25) ⋮ Wji =mji+cj,i−1Xj,i−1+ajiXji (26) ⋮ WjKj =mjKj+cj,Kj−1Xj,Kj−1+ajKjXjKj (27)

The layer-wise variational distributions of the bias terms denoted by 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

 qϕ(w)=d∏j=1qϕj(wj)qϕbj(bj)

where , , denotes the density of , denotes the density of , and 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:

 p(w)=d∏j=1p(wj)p(bj)

where denotes the density of and denotes the density of .

### 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]):

 DKL(qϕ(w)||p(w)) =d∑j=1[DKL(qϕj(wj)||p(wj))] +d∑j=1[DKL(qϕbj(bj)||p(bj))]

Thus, computing the overall divergence can be reduced to compute for fixed , since the remaining divergences compute completely analogously (only the indices differ). According to Hershey and Olsen [36] the KL divergence between two -dimensional normal distributions, given by and , computes as

 DKL(H||G)=12 [ln|Σg||Σh|+tr(Σ−1gΣh)−p (28) +(μh−μg)TΣ−1g(μh−μg)].

Thus, the determinant of the covariance matrix is required for the computation of the Kullback Leibler divergence . Using basic properties of determinants computes as follows for fixed :

 (29)

 DKL(qϕj(wj)||p(wj))=12[ln|ζ2jIKj||Σj|+tr((ζ2jIþKj)−1Σj) −Kj+(mj−μj)T(ζ2jIKj)−1(mj−μj)] =12⎡⎣−ln⎛⎝Kj∏i=1a2ji⎞⎠+τ2jζ2j||mj||22+1ζ2j||mj−μj||22+c⎤⎦ =12⎡⎣−Kj∑i=1(lna2ji)+τ2jζ2j||mj||22+1ζ2j||mj−μj||22+c⎤⎦ (30)

where 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 according to our novel approach the partial derivatives of the approximation of the negative log evidence lower bound described in Section II with respect to the variational parameters , are required. In particular, the partial derivatives of the loss function 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, depends on the network itself, with parameters sampled from the variational distribution

. With the help of the multivariate chain rule the required partial derivatives of

can be computed based on the classical derivatives used in non-Bayesian deep learning:

 ∂L∂mj=(∂wj∂mj)T∂L∂wj ⇒   ∂L∂mji=∑l∂L∂wjl∂wjl∂mji (31) ⇒   ∂Lδj=∑l∂L∂wjl∂wjl∂δj (32) ⇒   ∂Lγj=∑l∂L∂wjl∂wjl∂γj (33)

Equations only deal with the derivatives of 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 from the variational distribution the layer-wise derivatives () are computed as in non-Bayesian deep learning. Thus, the problem of finding closed form expressions for the required derivatives of reduces to the problem of finding these expressions for the ’s and the ’s. Taking account of the Equations the needed derivatives of the weights can be expressed in terms of the corresponding derivatives of the and the

 ∂wji∂mjk =⎧⎪ ⎪⎨⎪ ⎪⎩∂cj,i−1∂mjkxj,i−1+∂aji∂mjkxjik≠i1+∂cj,i−1∂mjixj,i−1+∂aji∂mjixjik=i (34) wji∂δj =∂cj,i−1∂δjxj,i−1+∂aji∂δjxji (35) wji∂γj =∂cj,i−1∂γjxj,i−1+∂aji∂γjxji (36)

where the index lies in the set , while for a given the indices and lie in the set . The derivatives of the () with respect to the variational parameters are given by

 ∂cji∂mjk =⎧⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪⎨⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪⎩(ρjτ2jmjimj,i+1)u−32cj,i−1∂cj,i−1∂mjkki+1 (37) ∂cji∂δj =wδvmji[2√u−τju−12(τjm2ji−cj,i−1∂cj,i−1∂τj)]τju (38) ∂cji∂γj =wγτ2jmjimj,i+1[√u−ρju−12(−cj,i−1∂cj,i−1∂ρj)]u (39)

where , , , and , and obviously each derivative of equals zero. Moreover, the derivatives of the () with respect to the variational parameters are given by:

 ∂aji∂mjk =⎧⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪⎨⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪⎩(ρjτ2j|mji||mj,i+1|)(−1)c−2ji∂cji∂mjkki+1 (40) ∂aji∂δj =wδρjτj|mji||mj,i+1|[2cji−τj∂cji∂τj]c2ji (41) ∂aji∂γj =wγτ2j|mji||mj,i+1|[cji−ρj∂cji∂ρj]c2ji (42)

In addition, the derivatives of with respect to the variational parameters are given by

 ∂ajKj∂mjk =⎧⎪ ⎪ ⎪⎨⎪ ⎪ ⎪⎩(−1)y−12cj,Kj−1∂cj,Kj−1∂mjkk

where .

Finally, the partial derivatives of the KL divergence (abbreviated with ) with respect to the variational parameters are given by:

 ∂∂mjkDKL =−Kj∑i=1(1aji∂aji∂mjk)+τ2jζ2jmjk+1ζ2j(mjk−μjk) (46) ∂∂δjDKL =<