{'fraction_non_alphanumeric': 0.03367362486667683, 'fraction_numerical': 0.014109401188480877, 'mean_word_length': 4.776447808484422, 'pattern_counts': {'":': 0, '<': 0, '<?xml version=': 0, '>': 0, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 2, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Good Initializations of Variational Bayes for Deep Models

Stochastic variational inference is an established way to carry out approximate Bayesian inference for deep models. While there have been effective proposals for good initializations for loss minimization in deep learning, far less attention has been devoted to the issue of initialization of stochastic variational inference. We address this by proposing a novel layer-wise initialization strategy based on Bayesian linear models. The proposed method is extensively validated on regression and classification tasks, including Bayesian DeepNets and ConvNets, showing faster convergence compared to alternatives inspired by the literature on initializations for loss minimization.


INTRODUCTION
Deep Neural Networks (dnns) and Convolutional Neural Networks (cnns) have become the preferred choice to tackle various supervised learning problems, such as regression and classification, due to their ability to model complex problems and the mature development of regularization techniques to control overfitting (Le-Cun et al., 2015;Srivastava et al., 2014). There has been a recent surge of interest in the issues associated with their overconfidence in predictions, and proposals to mitigate these (Guo et al., 2017;Lakshminarayanan et al., 2017). Bayesian techniques offer a natural framework to deal with such issues, but they are characterized by computational intractability (Bishop, 2006;Ghahramani, 2015).
A popular way to recover tractability is to use variational inference (Jordan et al., 1999). In variational inference, an approximate posterior distribution is introduced and its parameters are adapted by optimizing a variational objective, which is a lower bound to the marginal likelihood. The variational objective can be written as the sum of an expectation of the log- likelihood under the approximate posterior and a regularization term which is the negative Kullback-Leibler (kl) divergence between the approximating distribution and the prior over the parameters. Stochastic Variational Inference (svi) offers a practical way to carry out stochastic optimization of the variational objective. In svi, stochasticity is introduced with a doubly stochastic approximation of the expectation term, which is unbiasedly approximated using Monte Carlo and by selecting a subset of the training points (minibatching) (Graves, 2011;Kingma and Welling, 2014).
While svi is an attractive and practical way to perform approximate inference for dnns, there are limitations. For example, the form of the approximating distribution can be too simple to accurately approximate complex posterior distributions (Ha et al., 2016;Ranganath et al., 2015;Rezende and Mohamed, 2015). Furthermore, svi increases the number of optimization parameters compared to optimizing model parameters through, e.g., loss minimization; for example, a fully factorized Gaussian posterior over model parameters doubles the number of parameters in the optimization compared to loss minimization. This has motivated research into other ways to perform approximate Bayesian inference for dnns by establishing connections between variational inference and dropout (Gal and Ghahramani, 2016a,b;. The development of a theory to fully understand the optimization landscape of dnns and cnns is still in its infancy (Dziugaite and Roy, 2017) and most works have focused on the practical aspects characterizing the optimization of their parameters (Duchi et al., 2011;Kingma and Ba, 2015;Srivastava et al., 2014). If this lack of theory is apparent for optimization of model parameters, this is even more so for the understanding of the optimization landscape of the objective in variational inference, where variational parameters enter in a nontrivial way in the objective (Graves, 2011;Rezende et al., 2014). However, the use of stochastic gradient optimization, which is guaranteed to asymptotically converge to a local optimum, motivates us to investigate ways to position the optimizer closer to a local optimum from the outset. Initialization in svi plays a huge role in the convergence of the learning process; the illustrative example in Figure 1 shows how a poor initialization can prevent svi to converge to good solutions in short amount of time even for simple problems.
In this work, we focus on this issue affecting svi for dnns and cnns. While there is an established literature on ways to initialize model parameters of dnns when minimizing its loss (Glorot and Bengio, 2010;Saxe et al., 2013;Mishkin and Matas, 2015), to the best of our knowledge, there is no study on the equivalent for svi for Bayesian dnns and cnns. Inspired by the literature on residual networks (He et al., 2016) and greedy initialization of dnns (Bengio et al., 2006;Mishkin and Matas, 2015), we propose Iterative Bayesian Linear Modeling (i-blm), which is an initialization strategy for svi grounded on Bayesian linear modeling. Iterating from the first layer, i-blm initializes the posteriors at layer l by learning Bayesian linear models which regress from the input, propagated up to layer l, to the labels.
We show how i-blm can be applied in a scalable way and without considerable overhead to regression and classification problems, and how it can be applied to initialize svi not only for dnns but also for cnns. Through a series of experiments, we demonstrate that i-blm leads to faster convergence compared to other initalizations inspired by the work on loss minimization for dnns. Furthermore, we show that i-blm makes it possible for svi with a Gaussian approximation applied to cnns to compete with Monte Carlo Dropout (mcd; Gal and Ghahramani (2016b)), which is currently the state-of-the-art method to perform approximate inference for cnns. In all, thanks to the proposed initialization, we make it possible to reconsider Gaussian svi for dnns and cnns as a valid competitor to mcd, as well as highlight the limitations of svi with a Gaussian posterior in applications involving cnns.
In summary, in this work we make the following contributions: (1) we propose a novel way to initialize svi for dnns based on Bayesian linear models; (2) we show how this can be done for regression and classification; (3) we show how to apply our initialization strategy to cnns; (4) we empirically demonstrate that our proposal allows us to achieve performance superior to other initializations of svi inspired by the literature on loss minimization; (5) for the first time, we achieve state-of-the-art performance with Gaussian svi for cnns.

RELATED WORKS
The problem of initialization of weights and biases in dnns for gradient-based loss minimization has been extensively tackled in the literature since early breakthroughs in the field (Rumelhart et al., 1986;Baldi and Hornik, 1989). LeCun (1998) is one of the seminal papers discussing practical tricks to achieve an efficient loss minimization through back-propagation.
More recently, Bengio et al. (2006) propose a greedy layer-wise unsupervised pre-training that proved to help optimization and generalization. A justification can be found in Erhan et al. (2010), where the authors show that pre-training can act as regularization; by initializing the parameters in a region corresponding to a better basin of attraction for the optimization procedure, the model can reach a better local minimum and increase its generalization capabilities. Glorot and Bengio (2010) propose a simple way to estimate the variance for random initialization of weights that makes it possible to avoid saturation both in forward and back-propagation steps. Another possible strategy can be found in Saxe et al. (2013), that investigate the dynamics of gradient descend optimization, and propose an initialization based on random orthogonal initial conditions. This algorithm takes a weight matrix filled with Gaussian noise, decomposes it to orthonormal basis using a singular value decomposition and replaces the weights with one of the components. Building on this work, Mishkin and Matas (2015) propose a data-driven weight initialization by scaling the orthonormal matrix of weights to make the variance of the output as close as possible to one.
Variational inference addresses the problem of intractable Bayesian inference by reinterpreting inference as an optimization problem. Its origins can be tracked back to early works in MacKay (1992); Hinton and van Camp (1993); Neal (1997). More recently, Graves (2011) proposes a practical way to carry out variational inference using stochastic optimization (Duchi et al., 2011;Zeiler, 2012;Sutskever et al., 2013;Kingma and Ba, 2015). Kingma and Welling (2014) propose a reparameterization trick that allows for the optimization of the variational lower bound through automatic differentiation. To decrease the variance of stochastic gradients, which impacts convergence speed, this work was extended using the so-called local reparameterization trick, where the sampling from the approximate posterior over model parameters is replaced by the sampling from the resulting distribution over the dnn units .
In the direction of finding richer posterior families for variational inference, we mention the works on Normalizing Flows, a series of transformations that, starting from simple distributions, can capture complex structures in the variational posterior (Rezende and Mohamed, 2015;Kingma et al., 2016;Louizos and Welling, 2017;Huang et al., 2018). Alternatives can be found in Stein variational inference (Liu and Wang, 2016), quasi-Monte Carlo variational inference (Buchholz et al., 2018) and variational boosting (Miller et al., 2017).
To the best of our knowledge, there is no study that either empirically or theoretically addresses the problem of initialization of parameters for variational inference, and this work aims to fill this gap.

PRELIMINARIES
In this section we introduce some background material on Bayesian dnns and svi.

Bayesian Deep Neural Networks
Bayesian dnns are statistical models whose parameters (weights and biases) are assigned a prior distribution and inferred using Bayesian inference techniques. Bayesian dnns inherit the modeling capacity of dnns while allowing for quantification of uncertainty in model parameters and predictions. Considering an input x ∈ R Din and a corresponding output y ∈ R Dout , the relation between inputs and outpus can be seen as a composition of nonlinear vector-valued functions f (l) for each hidden layer (l) Let W be a collection of all model parameters (weights and biases) W (l) at all hidden layers. Each neuron computes its output as where φ(·) denotes a so-called activation function which introduces a nonlinearity at each layer. Note that we absorbed the biases in w.
Given a prior over W, the objective of Bayesian inference is to find the posterior distribution over all model parameters W using the available input data X = {x 1 , . . . , x n } associated with labels Y = {y 1 , . . . , y n } Bayesian inference for dnns is analytically intractable and it is necessary to resort to approximations. One way to recover tractability is through the use of variational inference techniques as described next.

Stochastic Variational Inference
In variational inference, we introduce a family of distributions q θ (W), parameterized through θ, and attempt to find an element of this family which is as close to the posterior distribution of interest as possible (Jordan et al., 1999). This can be formulated as a minimization with respect to θ of the kl divergence (Kullback, 1959) between the elements of the family q θ (W) and the posterior: Simple manipulations allow us to rewrite this expression as the negative lower bound (nelbo) to the marginal likelihood of the model (see supplementary material) where the first term is the expected negative log- , and the second term acts as regularizer, penalizing distributions q θ that deviate too much from the prior. When the likelihood factorizes across data points, we can unbiasedly estimate the expectation term randomly selecting a mini-batch B of m out of n training points Each term in the sum can be further unbiasedly estimated using N MC Monte Carlo samples as where W i ∼ q θ (W). Following Kingma and Welling (2014), each sample W i is constructed using the reparameterization trick, which allows to obtain a deterministic dependence of the nelbo w.r.t. θ. Alternatively, it is possible to determine the distribution of the dnn units f i before activation from q θ (W). This trick, known as the local reparameterization trick, allows one to considerably reduce the variance of the stochastic gradient w.r.t. θ and achieve faster convergence as shown by . Figure 2: Visual representation of the proposed method for initialization. In (a) and (b), we learn two Bayesian linear models, whose outputs are used in (c) to infer the following layer.

PROPOSED METHOD
In this section, we introduce our proposed Iterative Bayesian Linear Model (i-blm) initialization for svi. We first introduce i-blm for regression with dnns, and we then show how this can be extended to classification and to cnns.

Initialization of DNNs for Regression
In order to initialize the weights of dnns, we proceed iteratively as follows. Given a fully-factorized Gaussian approximation q(W (l) ), starting from the first layer, we can set the parameters of q(W (0) ) = N (W (0) |µ (0) , Σ (0) ) by running Bayesian linear regression with inputs X and labels Y . After this, we initialize the approximate posterior over the weights at the second layer q(W (1) ) = N (µ (1) , Σ (1) ) by running Bayesian linear regression with inputs X = Φ(XW (0) ) and labels Y . Here, Φ(·) denotes the elementwise application of the activation function to the argument, whereasW (0) is a sample from q(W (0) ). We then proceed iteratively in the same way up to the last layer.
The intuition behind i-blm is as follows. If one layer is enough to capture the complexity of a regression task, we expect to be able to learn an effective mapping right after the initialization of the first layer. In this case, we also expect that the mapping at the next layers implements simple transformations, close to the identity. Learning a set of weights with these characteristics starting from a random initialization is extremely hard; this motivated the work in He et al. (2016) that proposed the residual network architecture. Our i-blm initialization takes this observation as an inspiration to perform initialization for svi for general deep models.
In practice, i-blm proceeds as follows. Before applying the nonlinearity through the activation function, each layer in a Bayesian dnn can be seen as multivariate Bayesian linear regression model. Denoting by h (l) the number of output neurons at layer (l), this is equivalent to h (l) univariate Bayesian linear models. Instead of using the entire training set to learn the linear models, each one of these is inferred based on a random mini-batch of data, whose inputs are prop-agated through the previous layers. The complexity of i-blm is linear in the batch size and cubic in the number of neurons to be initialized. Figure 2 gives an illustration of the proposed method for a simple architecture.

Fully factorized Gaussian posterior
In the case of a fully factorized approximate posterior over the weights q(W (l) ) = ij q ij (W (l) ), we are interested in obtaining the best approximate to the solution of the Bayesian linear regression, which in general yields a posterior that is not factorized. For simplicity of notation, let w be the parameters of interest in Bayesian linear regression for a given output y = Y ·i . We can formulate this problem by minimizing the kl divergence from q(w) = N (w|m, diag(s 2 )) to the actual posterior p(w|X, y). This results in the mean being equal to the mean of p(w|X, y) and the variances (s 2 i ) −1 = Σ −1 ii ; see the supplementary material for the full derivation.

Initialization for Classification
In this section we show how our proposal can be extended to k-class classification problems. We assume a one-hot encoding of the labels, so that Y is an n × k matrix of zeros and ones (one for each row of Y ). Recently, it has been shown that it is possible to obtain an accurate modeling of the posterior over classification functions by applying regression on a transformation of the labels (Milios et al., 2018). This is interesting because it allows us to apply Bayesian linear regression as before in order to initialize svi for dnns.
The transformation of the labels is based on the formalization of a simple intuition, which is the inversion of the softmax transformation. One-hot encoded labels are viewed as a set of parameters of a degenerate Dirichlet distribution. We resolve the degeneracy of the Dirichlet distribution by adding a small regularization, say α = 0.01, to the parameters. At this point, we leverage the fact that Dirichlet distributed random variables can be constructed as a ratio of Gamma ran- We can then approximate the Gamma random variables with log-Normals by moment matching. By doing so, we obtain a representation of the labels which allows us to use standard regression with a Gaussian likelihood, and which retrieves an approximate Dirichlet when mapping predictions back using the softmax transformation (8). As a result, the latent functions obtained represent probabilities of class labels.
The only small complication is that the transformation imposes a different noise level for labels that are 0 or 1, and this is because of the non-symmetric nature of the transformation. Nevertheless, it is a simple matter to extend Bayesian linear regression to handle heteroscedasticity; see the supplement for a derivation of heteroscedastic Bayesian linear regression and (Milios et al., 2018) for more insights on the transformation to apply regression on classification problems.

Initialization of CNNs
The same method can be also applied on cnns. Convolutional layers are commonly implemented as matrix multiplication (e.g. as a linear model) between a batched matrix of patches and a reshaped filter matrix (Jia, 2014). Rather than using the outputs of the previous layer as they are, for convolutional layers each Bayesian linear model learns the mapping from patches to output features.

RESULTS
In this section, we compare different initialization algorithms for svi to prove the effectiveness of i-blm.
We propose a number of competitors inspired from the literature developed for loss minimization in dnns and cnns. In the case of cnns, we also compare with Monte Carlo Dropout (mcd; Gal and Ghahramani (2016a)), which is the state-of-the-art for inference in Bayesian cnns.
As we discussed, in the case of a fully factorized Gaussian approximation we double the number of parameters compared to loss minimization, as we have a mean and a variance parameter for each weight. At layer (l), given the fully-factorized variational distribution q(W (l) ) = i,j N (w (l) i,j |µ i,j , σ 2 i,j ), we initialize µ i,j and σ 2 i,j with the following methods. Uninformative The optimization of the posterior starts from the prior distribution q(W ) = i,j N (w i,j |0, 1). Note that this yields an initial kl divergence in the nelbo equal to zero.
Random Heuristic An extension to commonly used heuristic with µ i,j = 0 and σ 2 i,j = 1 Din , with D in the number of input features at layer (l).
Xavier Normal Originally proposed by Glorot and Bengio (2010), it samples all weights independently from a Gaussian distribution with zero mean and σ 2 = 2 Din+Dout . This variance-based scaling avoids issues with vanishing or exploding gradients. Although this work considers only the case with linear activations, this initialization has been shown to work well in many applications. In this case, it is straightforward to extend it to the case of svi; indeed, instead of sampling, we directly set µ i,j = 0 and σ 2 i,j = 2 Din+Dout , knowing that the sampling is performed during the Monte Carlo estimate of the log-likelihood.
Orthogonal Starting from an analysis of learning dynamics of dnns with linear activations, Saxe et al. (2013) propose an initialization scheme with orthonormal weight matrices. The idea is to decompose a Gaussian random matrix onto an orthonormal basis, and use the resulting orthogonal matrix. We adapt this method for svi by initializing the mean matrix with the orthogonal matrix and σ 2 i,j = 1 Din . For our experiments, we use the implementation in pytorch (Paszke et al., 2017) provided by the Authors, which uses a QR-decomposition.
Layer-Sequential Unit-Variance (LSUV) Starting from orthogonal initialization, Mishkin and Matas (2015) propose to perform a layer sequential variance scaling of the weight matrix. By implementing a datadriven greedy initialization, it generalizes the results to any nonlinear activation function and even to any type of layers that can impact the variance of the activations. We implement Layer-Sequential Unit-Variance (lsuv) for the means, while the variances are set to σ 2 i,j = 1 Din .

Experiments
Throughout the experiments, we use adam optimizer (Kingma and Ba, 2015) with learning rate 10 −3 , batch size 64, and 16 Monte Carlo samples at training time and 128 at test time. All experiments are run on a workstation equipped with two 16c/32t Intel Xeon CPU and four NVIDIA Tesla P100, with a maximum time budget of 24 hours (never reached). To better understand the effectiveness of different initialization, all learning curves are plotted w.r.t. training iteration rather than wall-clock time. In all experiments, even the ones involving millions of parameters, the initialization costs only few seconds (a couple of minutes in the case of AlexNet), so we do not shift the curves as this cost is negligible in the whole training procedure. Regression with a shallow architecture In this experiment we compare initialization methods for a shallow dnn architecture on two datasets. The architecture used in these experiments has one single hidden layer with 100 hidden neurons and relu activations. We impose that the approximate posterior has fully factorized covariance. Figure 4 shows the learning curves on the powerplant (n = 9568, d = 4) and protein (n = 45730, d = 9) datasets, repeated over five different train/test splits. i-blm allows for a better initialization compared to the competitors, leading to a lower root mean square error (rmse) and lower mean negative log-likelihood (mnll) on the test for a given computational budget. We refer the reader to the supplementary material for a more detailed analysis of the results.
Regression with a deeper architecture Similar considerations hold when increasing the depth of the model, keeping the same experimental setup. Figure 5 shows the progression of the rmse and mnll error metrics when using svi to infer parameters of a dnn with five hidden layers and 100 hidden neurons per layer, and relu activations. Again, the proposed initialization allows svi to converge faster than when using other initializations.
Classification with a deep architecture Using the same deep dnn architecture as in the last experiment (five hidden layers with 100 neurons), we tested i-blm with classification problems on mnist (n = 70000, d = 784), eeg (n = 14980, d = 14), credit (n = 1000, d = 24) and spam (n = 4601, d = 57). Interestingly, with this architecture, some initialization strategies struggled to converge, e.g., uninformative on mnist and lsuv on eeg. The gains offered by i-blm achieves the most striking results on mnist. After less than 1000 training steps (less than an epoch), it can  already reach a test accuracy greater than 95%; other initalizations reach such performance much later during training. Even after 100 epochs, svi inference initialized with i-blm provides on avarage an increase up to 14% of accuracy at test time. Full results are reported in the supplementary material.
Experiments on CNNs For this experiment, we implemented a Bayesian version of the original LeNet-5 architecture proposed by  with two convolutional layers of 6 and 16 filters, respectively and relu activations applied after all convolutional layers and fully-connected layers. We tested our framework on mnist and on cifar10. The only initialization strategies that achieve convergence are orthogonal and lsuv, along with i-blm; the other methods did not converge, meaning that they push the posterior back Step I-BLM (this work) Orthogonal LSUV Figure 7: Progression of error rate and mnll over training iterations for different initialization strategies using LeNet-5 on mnist and cifar10.
to the prior. Figure 7 reports the progression of error rate and mnll. For both mnist and cifar10, i-blm places the parameters where the network can consistently deliver better performance both in terms of error rate and mnll throughout the entire learning procedure.
Comparison with Monte Carlo Dropout Monte Carlo Dropout (mcd; Gal and Ghahramani (2016b)) offers a simple and effective way to perform approximate Bayesian cnn inference, thanks to the connection that the Authors have established between dropout and variational inference. In this experiment, we aim to compare and discuss benefits and disadvantages of using a Gaussian posterior approximation with respect to the Bernoulli approximation that characterizes mcd. For a fair comparison, we implemented Step I-BLM (this work) Montecarlo Dropout

I-BLM Initialization
Test on mnist Test on not-mnist Figure 9: Entropy distribution while testing on mnist and not-mnist the same LeNet-5 architecture and the same learning procedure in Gal and Ghahramani (2016b) 1 . In particular, for mnist, the two convolutional layers have 20 and 50 filters, respectively. Dropout layers are placed after every convolutional and fully-connected layers with a dropout probability of 0.5. To replicate the results in Gal and Ghahramani (2016b), we used the same learning rate policy base-lr × (1 + ξ × iter) −p with ξ = 0.0001, p = 0.75, base-lr = 0.01 and weight decay of 0.0005. Figure 8 shows the learning curves. Monte Carlo Dropout achieves state-of-art error rate but the form assumed by mcd for the posterior is reflected on an higher mnll compared to svi with a Gaussian posterior. Provided with a nontrivial initialization, Gaussian svi can better fit the model and deliver a better quantification of uncertainty.
Analysis of out-of-sample uncertainty One of the advantages of Bayesian inference is the possibility to reason about uncertainty. With this experiment, we aim to demonstrate that svi with a Gaussian approximate posterior is competitive with mcd in capturing uncertainty in predictions. To show this, we focus on a cnn with the LeNet-5 architecture. We run mcd and svi with a Gaussian approximate posterior with the proposed initialization on mnist. At test time, we carry out predictions on both mnist and not-mnist; the latter is a dataset equivalent to mnist in input dimensions (1 × 28 × 28) and number of classes, but it represents letters rather than numbers 2 . This experimental setup is often used to check that the entropy of the predictions on not-mnist are actually higher than the entropy of the predictions on mnist. We report the entropy of the prediction on mnist and not-mnist in Figure 9. mcd and svi behave similarly on mnist, but on not-mnist the the histogram of the entropy indicates that svi yields a slightly higher uncertainty compared to mcd.
Deeper CNNs -AlexNet We report here error rate and mnll for svi with i-blm and mcd on AlexNet (Krizhevsky et al., 2012). The cnn is composed by a stack of five convolutional layers and three fully-connected layers for a total of more than 1M parameters (2M for svi). In this experiment, we have experienced the situation in which, due to the overparameterization of the model, the nelbo is completely dominated by the kl divergence. Therefore, the prior has a large influence on the optimization, so we decided to follow the approach in Graves (2011), allowing for a phase of optimization of the variances of the prior over the parameters. The results are reported in Figure 8. Once again, we show that svi with i-blm provides a lower negative log-likelihood with respect to Bernoulli approximation in mcd. In light of these results and considerations, we can consider Gaussian svi as a competitive method to carry out Bayesian inference on cnns.

CONCLUSIONS
This work fills an important gap in the literature of Bayesian deep learning, that is how to effectively initialize variational parameters in svi. We proposed a novel way to do so, i-blm, which is based on an iterative layer-wise initialization based on Bayesian linear models. Through a series of experiments, including regression and classification with dnns and cnns, we demonstrated the ability of our approach to consistently initialize the optimization in a way that makes convergence faster than alternatives inspired from the state-of-the-art in loss minimization for deep learning.
Thanks to i-blm, it was possible to carry out an effective comparison with mcd, which is currently the state-of-the-art method to carry out approximate inference for dnns and cnns. This suggests a number of directions to investigate to improve on svi and Bayesian cnns. We found that the choice of the prior plays an important role in the behavior of the optimization, so we are investigating ways to define sensible priors for these models. Furthermore, we are looking into ways to extend our initialization strategy to more complex posterior distributions beyond Gaussian and to other deep models, such as Deep Gaussian Processes.
A Full derivation of variational lower bound

B Bayesian linear regression
We express the likelihood and the prior on the parameters as follows: Denote by X the n × d matrix containing n input vectors x i ∈ R d , and let Y be the set consisting of the corresponding multivariate labels y i . In Bayesian linear regression we introduce a set of latent variables that we compute as a linear combination of the input through a set of weights, and we express the likelihood and the prior on the parameters as follows: The posterior of this model is: which implies that the posterior factorizes across the columns of W , with factors Similarly, the marginal likelihood factorizes as the product of the following factors p(Y ·i |X, L, Λ) = N (Y ·i |0, L + XΛX )

C Heteroscedastic Bayesian linear regression
We can extend Bayesian linear regression to the heteroscedastic case where L = diag(σ 2 ) and Λ = I. These yield p(W ·i |Y, X, σ 2 ) = N (W ·i |µ i , Σ i ) with The expression for the marginal likelihood is computationally unconvenient due to the need to deal with an n × n matrix. We can use Woodbury identities 3 to express this calculation using Σ i . In particular, Using Woodbury identites, we can rewrite the algebraic operations as follows: So, wrapping up, we ca express all quantities of interest as: If we factorize Σ −1 i = QQ , we obtain: Predictions follow from the same identities as before -looking at the predicted latent process, we have We can again remove the dependence from the dimensions of W that do not affect the prediction for the ith function as D Full derivation of fully factorized Gaussian posterior approximation to Bayesian linear regression posterior For simplicity of notation, let w be the parameters of interest in Bayesian linear regression for a given output y = Y ·i . We can formulate the problem of obtaining the best approximate factorized posterior of a Bayeian linear model as a minimization of the Kullback-Leibler divergence between q(w) = N (w|m, diag(s 2 )) and the actual posterior p(w|X, y). The expression of the KL divergence between multivariate Gaussians p 0 = N (W |µ 0 , Σ 0 ) and p 1 = N (W |µ 1 , Σ 1 ) is as follows: The KL divergence is not symmetric, so the order in which we take this matters. In case we consider KL[p(w|X, y)||q(w)], the expression becomes: It is a simple matter to show that the optimal mean m is µ as m appears only in the quadratic form which is clearly minimized when m = µ. For the variances s 2 , we need to take the derivative of the KL divergence and set it to zero: Rewriting the trace term as the sum of the Hadamrd product of the matrices in the product ij (diag( This results in s 2 i = Σ ii , which is the simplest way to approximate the correlated posterior over w but it is going to inflate the variance in case of strong correlations. In case we consider KL[q(w)||p(w|X, y)], the expression of the KL becomes: Again, the optimal mean m is µ. For the variances s 2 , we need to take the derivative of the KL divergence and set it to zero: Rewriting the trace term as the sum of the Hadamrd product of the matrices in the product