The exact contents of citations.db main_text.text for this paper — one flattened LaTeX string, title through conclusion, appendix excluded, unmodified except for removing email addresses. This is what our citation measures are computed over.
52,507 characters
Quasi Black-Box Variational Inference with Natural Gradients for Bayesian Learning
\maketitle
\begin{abstract}
We develop an optimization algorithm suitable for Bayesian learning in complex models. Our approach relies on natural gradient updates within a general black-box framework for efficient training with limited model-specific derivations. It applies within the class of exponential-family variational posterior distributions, for which we extensively discuss the Gaussian case for which the updates have a rather simple form. Our Quasi Black-box Variational Inference (QBVI) framework is readily applicable to a wide class of Bayesian inference problems and is of simple implementation as the updates of the variational posterior do not involve gradients with respect to the model parameters, nor the prescription of the Fisher information matrix. We develop QBVI under different hypotheses for the posterior covariance matrix, discuss details about its robust and feasible implementation, and provide a number of real-world applications to demonstrate its effectiveness.
\end{abstract}
\section{Introduction} \label{sec:intro}
Machine Learning (ML) techniques have been proved to be successfully in many prediction and classification tasks across natural-language processing \cite{young_recent_2018}, computer vision \cite{krizhevsky_imagenet_2017}, time-series \cite{langkvist2014review} and finance applications \cite{dixon2020machine}, among the several ones. Recently, Bayesian methods have gained considerable interest in the field as an attractive alternative to point estimation, especially for their ability to address uncertainty via posterior distribution, generalize while reducing overfitting \cite{hoeting1999bayesian}, and for enabling sequential learning \cite{freitas2000hierarchical} while retaining prior and past knowledge. Although Bayesian principles have been proposed in ML decades ago \cite[e.g.][]{mackay1992bayesian,mackay1995probable,lampinen2001bayesian}, it has been only recently that fast and feasible methods boosted a growing use of Bayesian methods in complex models \cite{osawa_practical_2019,khan2018fastscalable,khan2018fastyetsimple}.
The most challenging task is the computation of the posterior \cite{jospin2022hands}. In the typical ML setting characterized by a high number of parameters and a considerable size of data, traditional sampling methods turn unfeasible, yet approximate methods such as Variational Inference (VI) have been shown to be suitable and successful \cite{saul1996mean,wainwright2008graphical,hoffman2013stochastic,blei2017variational}. Furthermore, recent research advocates the use of the natural gradients for boosting the optimum search and the training \cite{wierstra2014natural}, enabling fast and accurate Bayesian learning algorithms that are scalable and versatile. \\
Stochastic Gradient Descent (SGD) methods \cite{robbins1951stochastic} promoted the use of VI for complex high-dimensional DL models \cite{hoffman2013stochastic,salimans2014using}, yet they require significant implementation and tuning \cite{ranganath2014black}. On the other hand, the use of natural gradients \cite{amari1998natural} within the VI framework under exponential family approximations has been shown to be notably efficient and robust \cite{hoffman2013stochastic,wierstra2014natural,khan2018fastscalable}, often leading to simple updates \cite{khan2018fastyetsimple}. In particular, the common choice of a Gaussian approximation to the true posterior \cite{opper2009variational}, results from \cite{bonnet1964transformations} and \cite{price1958useful} leads to simple updates of the variational mean and precision matrix \cite{khan2017conjugate}. Several algorithms have been derived following these results \cite{khan2018fastscalable}, however, they require the model gradients (and perhaps its hessian), and the positive-definiteness constraint on the posterior covariance matrix is a common problem, see. e.g. \cite{tran2021variational,khan2018fastscalable,lin2020handling,ong2018variational}. \\
Traditional optimization algorithms rely on the extensive use of gradients to dynamically adjust the model weights to minimize a given loss, e.g. through backpropagation.
This holds also in Bayesian learning based on VI, as the choice of the likelihood and variational approximation determine a cascade of model-specific derivations, not keen to immediately fit a general setting based on a ready-to-use, plug-and-play optimizer. \\
Black-box methods \cite{ranganath2014black} are of straightforward and general use relying on stochastic sampling without model-specific derivations \cite{salimans2014using,paisley2012variational,kingma2013auto}. Even though black-box methods can benefit from numerous improvements such as variance reduction techniques \cite{paisley2012variational,ranganath2014black}, applying natural gradient updates is challenging, as the specification of the Fisher matrix is required.
We propose Quasi Black-box Variational Inference (QBVI) which introduces natural gradients updates within the black-box VI framework. It combines the flexibility of black-box methods with the SGD theory for exponential family VI in a feasible, scalable and flexible optimizer.
In particular, we rely on VI to approximate the true posterior via parameters' updates that only involve function queries without requiring backpropagation, assumptions on the form of the likelihood, or restrictions on the backbone network. Instead, we employ the typical variational Gaussian assumption \cite{kingma2013auto,khan2018fastyetsimple,osawa_practical_2019}, under which the updates take a rather simple form, embracing both closed-form natural gradient VI elements and black-box elements (thus the word ``quasi'').
We provide results and update rules under both full and diagonal posterior covariance assumptions, and discuss the generalization of the proposed method under the well-established mean-field approximation \cite[e.g.][]{blei2017variational}. We furthermore develop the method of Control variates \cite[e.g.][]{paisley2012variational,salimans2014using} for efficient large-MC sampling, and provide solutions for valid covariance updates adapting developments from existing methods \cite{khan2018fastscalable,tran2021variational,ong2018variational}. Following and discussing the typical practicalities and recommendations in VI applications \cite{tran2021practical}, we provide experiments to validate the proposed optimizer, showing its feasibility in complex learning tasks and promoting its use as a practical and ready-to-use tool for Bayesian optimization.
\section{Variational Inference for Bayesian deep learning}
\subsection{Variational inference}
Let $y$ denote the data and $p\autobracket*{y\vert \theta}$ the likelihood of the data based on a postulated model with $\theta \in \Theta$ a $d$-dimensional vector of model parameters. Let $p\autobracket*{\theta}$ be the prior distribution on $\theta$. The goal of Bayesian inference is the posterior distribution $ p\autobracket*{\theta\vert y} =p\autobracket*{\theta} p\autobracket*{y\vert \theta}/p\autobracket*{y}$.
Bayesian inference is generally difficult due to the fact that the marginal likelihood $p\autobracket*{y}$ is often intractable and of unknown form. In high-dimensional applications, Monte Carlo (MC) methods for sampling the posterior are challenging and unfeasible, and Variational Inference (VI) is an attractive alternative.
VI consists of an approximate method where the posterior distribution is approximated by a probability density $q\autobracket*{\theta}$ (called variational distribution) belonging to some tractable class of distributions $\mathcal{Q}$, such as the exponential family. VI thus turns the Bayesian inference problem into that of finding the best approximation $q^\star \autobracket*{\theta} \in \mathcal{Q}$ to $p\autobracket*{\theta\vert y}$ by minimizing the Kullback-Leibler (KL) divergence from $q\autobracket*{\theta}$ to $p\autobracket*{\theta\vert y}$,
\begin{equation*}
q^\star = \operatorname*{arg\,min}_{q \in \mathcal{Q}} \text{KL}\autobracket*{q \vert \vert p\autobracket*{\theta\vert y}} = \operatorname*{arg\,min}_{q \in \mathcal{Q}} \int q\autobracket*{\theta} \log \frac{q\autobracket*{\theta}}{p\autobracket*{\theta\vert y}} d\theta \text{.}
\end{equation*}
By simple manipulations, it can be shown the KL minimization problem is equivalent to the maximization problem of the so-called Lower Bound (LB) on $\log p\autobracket*{y}$ \cite[e.g.]{tran2021practical},
\begin{equation*}
\mathcal{L}\autobracket*{q} \coloneqq \int q\autobracket*{\theta} \log \frac{p\autobracket*{\theta} p\autobracket*{y\vert \theta}}{q\autobracket*{\theta}} d \theta
= \mathbb{E}_{q}\sparen*{\log \frac{p\autobracket*{\theta} p\autobracket*{y\vert \theta}}{q\autobracket*{\theta}}} \text{.}
\end{equation*}
For any random vector $\theta$ and a function $g\autobracket*{\theta}$ we denote by $\mathbb{E}_{f}[ g\autobracket*{\theta}]$ the expectation of $g\autobracket*{\theta}$ where $\theta$ follows a probability distribution with density $f$, i.e. $\mathbb{E}_{f}[ g\autobracket*{\theta}] =\mathbb{E}_{\theta \sim f}[g\autobracket*{\theta}]$.
To make explicit the dependence of the LB on some vector of parameters $\zeta$ parametrizing the variational posterior we write $\mathcal{L}\autobracket*{\zeta} = \mathcal{L}\autobracket*{q_\zeta} = \mathbb{E}_{q_\zeta}\sparen*{\log p\autobracket*{\theta} - \log q_\zeta\autobracket*{\theta} + \log p\autobracket*{y\vert \theta} }$.
We operate within the fixed-form variational inference (FFVI) framework, where the parametric form of the variational posterior is set.
\subsection{SGD and natural gradients}
A straightforward approach to maximize the LB is that of using a gradient-based method such as SGD, ADAM, RMSprop.
In the FFVI setting with $\mathcal{Q}$ being the exponential family, the LB is often optimized in terms of the natural parameter $\lambda$ \cite{wainwright2008graphical}.
The application of the SGD update based on the standard gradient is problematic because it ignores the information geometry of the distribution $q_\lambda$ \cite{amari1998natural}, as it implicitly relies on the Euclidean norm to capture the dissimilarity between two distributions
which can indeed be a quite poor and misleading measure of dissimilarity \cite{khan2018fastyetsimple}. By replacing the Euclidean norm with the KL divergence, the SGD update results in the following natural gradient update:
\begin{align}\label{eq:sgd_nat_update}
\lambda_{t+1} =\lambda_t+\beta_t\left. \sparen*{\tilde{\nabla}_\lambda\mathcal{L}\autobracket*{\lambda}} \right|_{\lambda = \lambda_t}
\end{align}
The natural gradient update results in better step directions towards the optimum when optimizing the parameter of a distribution. The natural gradient of $\mathcal{L}\autobracket*{\lambda}$ is obtained by rescaling the euclidean gradient $\nabla_\lambda \mathcal{L}\autobracket*{\lambda}$ by the inverse of the Fisher Information Matrix (FIM) $\mathcal{I}_\lambda$,
\begin{equation} \label{eq:def_natgrad}
\tilde{\nabla}_\lambda \mathcal{L}\autobracket*{\lambda} = \mathcal{I}^{-1}_\lambda\nabla_\lambda\mathcal{L}_\lambda \text{.}
\end{equation}
By replacing in the above $\nabla_\lambda \mathcal{L}\autobracket*{\lambda}$ with a stochastic estimate $\hat{\nabla}_\lambda\mathcal{L}\autobracket*{\lambda}$ one obtains a stochastic natural gradient update.
\subsection{Maximization of the Lower Bound}
With respect to an exponential-family prior-posterior pair of the same parametric form, with natural parameters $\eta$ and $\lambda$ respectively, the natural gradient of the LB is given by
\begin{align}
\tilde{\nabla}_\lambda\mathcal{L}\autobracket*{\lambda} &= \tilde{\nabla}_\lambda \mathbb{E}_{q_\lambda}\sparen*{\log \frac{p_\eta\autobracket*{\theta}}{q_\lambda\autobracket*{\theta}}} + \tilde{\nabla}_\lambda\mathbb{E}_{q_\lambda}\sparen*{\log p\autobracket*{y\vert \theta}} \label{eq:LB_gradient_generic}\\
&= \eta - \lambda +\tilde{\nabla}_\lambda\mathbb{E}_{q_\lambda}\sparen*{\log p\autobracket*{y\vert \theta}}. \label{eq:LB_gradient_khan0}
\end{align}
A proof is provided in Appendix \ref{proof:log_trick}. By applying \eqref{eq:LB_gradient_khan0} to \eqref{eq:sgd_nat_update}, we have the compact update \cite{khan2018fastscalable}
\begin{align}
\lambda_{t+1} = \autobracket*{1-\beta}\lambda_t + \beta\autobracket*{\eta + \tilde{\nabla}_\lambda\mathbb{E}_{q_\lambda}\sparen*{\log p\autobracket*{y\vert \theta}}}. \label{eq:update_khan}
\end{align}
Two points are critical. First the $\mathcal{I}^{-1}_\lambda$ requirement, which is implicit in the natural gradient definition. Though impractical, it is perhaps possible to estimate $\mathcal{I}^{-1}_\lambda$ at each iteration with an iterative conjugate gradient method using only matrix-vector products \cite{tran2021practical}, or exploiting a certain (assumed or imposed) structure of the FIM \cite{tran2020bayesian}.
Second, the expectation in \eqref{eq:update_khan} is generally intractable as it depends on the specified log-likelihood. Thus, sampling methods can be used.
\subsection{Exponential family basics}
Assume $q_\lambda\autobracket*{\theta}$ belongs to an exponential family distribution. Its probability density function is parametrized as
\begin{equation}\label{eq:exp_fam_representation}
q_\lambda\autobracket*{\theta} = h\autobracket*{\theta} \exp\autobracket*{\phi\autobracket*{\theta}^\top \lambda - A\autobracket*{\lambda}}{\lambda}}\text{,}
\end{equation}
where $\lambda \in \Omega = \{ \lambda \in \mathbb{R}^d: A\autobracket*{\lambda} < +\infty \}$ is the natural parameter, $\phi\autobracket*{\theta}$ the sufficient statistic, $A\autobracket*{\lambda} = \log \int h\autobracket*{\theta} \exp(\phi\autobracket*{\theta}^\top \lambda) d\nu$ the log-partition function, determined upon the measure $\nu$, $\phi$, and the function $h$.
When $\Omega$ is a non-empty open set, the exponential family is referred to as regular. Furthermore, if there are no linear constraints among the components of $\lambda$ and $\phi\autobracket*{\theta}$, the exponential family \eqref{eq:exp_fam_representation} is of minimal representation. Non-minimal families can always be reduced to minimal families through a suitable transformation and reparametrization, leading to a unique parameter vector $\lambda$ associated with each distribution \cite{wainwright2008graphical}.
The mean (or expectation) parameter $m \in \mathcal{M}$ is defined as a function of $\lambda$, $m\autobracket*{\lambda} = \mathbb{E}_{q_\lambda}\sparen*{\phi\autobracket*{\theta}} = \nabla_\lambda A\autobracket*{\lambda}$. Moreover, for the Fisher Information Matrix $\mathcal{I}_\lambda = -\mathbb{E}_{q_\lambda}\sparen*{\nabla_\lambda^2 \log q_\lambda\autobracket*{\theta}}$ it holds that $\mathcal{I}_\lambda =\nabla^2_\lambda A\autobracket*{\lambda} = \nabla_\lambda m \text{.}$
Under minimal representation, $A\autobracket*{\lambda}$ is convex, thus the mapping $\nabla_\lambda A = m:\Omega \rightarrow \mathcal{M}$ is one-to-one, and $\mathcal{I}_\lambda$ is positive definite and invertible \cite{nielsen2009statistical}.
Therefore, under minimal representation we can express $\lambda$ in terms of $m$ and thus $\mathcal{L}\autobracket*{\lambda}$ in terms of $\mathcal{L}\autobracket*{m}$ and vice versa \cite{khan2018fastyetsimple}. Furthermore, for a generic differentiable function $\mathcal{L}$, by applying the chain rule, $\nabla_\lambda\mathcal{L} = \nabla_\lambda m \nabla_m \mathcal{L} = \nabla_\lambda \autobracket*{\nabla_\lambda A\autobracket*{\lambda}}et*{\lambda}} \mathcal{L}= \nabla^2_\lambda A\autobracket*{\lambda} \mathcal{L} = \mathcal{I}_\lambda \nabla_m \mathcal{L}$, from which
\begin{equation}
\tilde{\nabla}_\lambda \mathcal{L} = \mathcal{I}^{-1}_\lambda \nabla_\lambda\mathcal{L} = \mathcal{I}^{-1}_\lambda \autobracket*{\mathcal{I}_\lambda \nabla_m \mathcal{L}} = \nabla_m \mathcal{L} \text{,} \label{eq:natgrad_lambda_grad_m}
\end{equation}
expressing the natural gradient in the natural parameter space as the gradient in the expectation parameter space. Eq. \eqref{eq:natgrad_lambda_grad_m} enables the computation of natural gradients in the natural parameter space without requiring the FIM \cite{khan2017conjugate,khan2018fastyetsimple}.
This translates into the following equivalence:
\begin{align} \label{eq:natgrad2mgrad}
\tilde{\nabla}_\lambda\mathcal{L} = \eta-\lambda + \tilde{\nabla}_\lambda \mathbb{E}_{q_\lambda} \sparen*{\log p\autobracket*{y\vert \theta}} = \eta-\lambda + \nabla_m \mathbb{E}_{q_\lambda}\sparen*{\log p\autobracket*{y\vert \theta}}\text{.}
\end{align}
\section{Proposed method}
\subsection{Quasi-Black-box natural-gradient VI}
Our approach relies on the following equivalence for the natural gradient of the LB:
\begin{align}
\tilde{\nabla}_\lambda\mathcal{L}\autobracket*{\lambda} &=
\eta-\lambda +\tilde{\nabla}_\lambda\mathbb{E}_{q_\lambda}\sparen*{\log p\autobracket*{y\vert \theta}} \label{eq:LB_gradient_khan}\\
&= \eta-\lambda +\mathbb{E}_{q_\lambda}\sparen*{\tilde{\nabla}_\lambda\sparen*{\log q_\lambda\autobracket*{\theta}} \log p\autobracket*{y\vert \theta}}gpy} \label{eq:LB_gradient_martin}\text{.}
\end{align}
A proof is provided in Appendix \ref{proof:log_trick}. The difference between equations \eqref{eq:LB_gradient_khan} and \eqref{eq:LB_gradient_martin} is substantial. Eq. \eqref{eq:LB_gradient_khan} involves the natural gradient with respect to the log-likelihood, thus implicitly the natural gradient with respect to the underlying function of the covariates over which $p\autobracket*{y\vert \theta}$ is parametrized. Eq. \eqref{eq:LB_gradient_martin} instead requires the gradients of the score function of the variational distribution only. If $\tilde{\nabla}_\lambda \sparen*{\log q\autobracket*{\theta}}$ is available in closed form, \eqref{eq:LB_gradient_martin} prescribes a practical gradient-free method for updating the posterior's natural parameter. Given that $\tilde{\nabla}_\lambda \sparen*{\log q\autobracket*{\theta}}$ has a closed-form, the wording gradient-free is meant in the sense that no further gradients need to be explicitly evaluated, i.e. those of $\log p\autobracket*{y\vert \theta}$, as it typically happens e.g. in backpropagation: the natural gradient of the log-likelihood w.r.t. the network parameters is not required.
This is the case for several exponential-family distributions.
In section \ref{subsec:nat_grad_free_gaussian_vi} we show that $\tilde{\nabla}_\lambda \sparen*{\log q\autobracket*{\theta}_\lambda\autobracket*{\theta}}$ is indeed tractable when $q_\lambda\autobracket*{\theta}$ is Gaussian.
For a function $f_\lambda$ of the natural parameter, \cite{khan2017conjugate} shows that the natural gradient $\tilde{\nabla}_\lambda f_\lambda$ corresponds to the gradient $\nabla_mf_m$ where neither the FIM or its inverse are involved (under minimal representation we can express $f_\lambda$ in terms of $m$, i.e. $f_m := f_\lambda$.
Eqs. \eqref{eq:LB_gradient_martin} and \eqref{eq:natgrad_lambda_grad_m} therefore lead to the following update:
\begin{align}\label{eq:QBVI_update}
\lambda_{t+1} &= \autobracket*{1-\beta}\lambda_t + \beta\autobracket*{\eta + \mathbb{E}_{q_\lambda}\sparen*{\nabla_m \sparen*{\log q_\lambda\autobracket*{\theta}}\log p\autobracket*{y\vert \theta}}gpy}} \text{.}
\end{align}
We refer to our approach based on the above update as Quasi-Black-box (natural-gradient) Variational Inference (QBVI). Indeed the (natural) gradient of the lower bound involves two ingredients: (i) function queries for $\log p\autobracket*{y\vert \theta}$ evaluated at some value $\theta$ drawn from the posterior $q\autobracket*{\lambda_t}$, without requiring model-specific derivations for $\log p\autobracket*{y\vert \theta}$ typical in VI applications, and (ii) $\nabla_m \sparen*{\log q_\lambda\autobracket*{\theta}}$ that needs to be provided. Thus the name quasi-black-box. Algorithm \ref{alg:general_qbvi} summarized the approach in the general case of an exponential-family prior-posterior pair of the same parametric form.
\subsection{QBVI under Gaussian variational posteriors} \label{subsec:nat_grad_free_gaussian_vi}
This section focuses on the implementation of the natural gradient update \eqref{eq:QBVI_update} under the typical multivariate Gaussian variational framework, perhaps the most prominent in VI \cite[][]{kucukelbir2015automatic,kingma2013auto,tran2021variational} and Bayesian learning \cite[][]{graves2011practical,blundell2015weight,khan2018fastscalable,osawa_practical_2019}, among the many others.
The multivariate Gaussian distribution $\mathcal{N}\autobracket*{\mu,S}$ with $d$-dimensional mean vector $\mu$ and covariance matrix $S$ can be seen as a member of the exponential family \eqref{eq:exp_fam_representation}. Its density reads
\begin{equation*}
q_\lambda\autobracket*{\theta} = \autobracket*{2\pi}^{d/2}\exp\{\phi\autobracket*{\theta}^\top \lambda - \frac{1}{2}\mu^\top S^{-1} \mu -\frac{1}{2}\log \vert S\vert\} \text{,}
\end{equation*}
where
\begin{equation*}
\phi\autobracket*{\theta} =\begin{bmatrix}
\theta \\
\theta\theta^\top
\end{bmatrix} \text{,}\quad\quad
\lambda =\begin{bmatrix}
\lambda_1 \\
\lambda_2
\end{bmatrix} =\begin{bmatrix}
S^{-1} \mu \\
-\frac{1}{2}S^{-1}
\end{bmatrix} \text{,}\quad\quad
m = \begin{bmatrix}
m_1 \\
m_2
\end{bmatrix} = \begin{bmatrix}
\mu \\
S + \mu\mu^\top
\end{bmatrix} \text{,}
\end{equation*}
and $A\autobracket*{\lambda} = -\frac{1}{4}\lambda_1^\top\lambda_2^{-1} \lambda_1-\frac{1}{2} \log \autobracket*{-2\lambda_2} $. On the other hand, $\zeta = \sparen*{\zeta_1^\top, \zeta_2^\top}^\top$ with $\zeta_1 = \mu = m_1$ and $\zeta_2 = S = m_2 - \mu \mu^\top$, constitutes the common parametrization of the multivariate Gaussian distribution in terms of its mean and variance-covariance matrix.
\newtheorem{prop_nat_grad_gaussian}{\textbf{Proposition}}
\begin{prop_nat_grad_gaussian} \label{prop:nat_grad_gaussian}
For a differentiable
function $\mathcal{L}$, and for $q_\lambda \sim \mathcal{N}\autobracket*{\mu,S}$,
\begin{align*}
\nabla_{m_1} \mathcal{L} = \nabla_\mu \mathcal{L} -2\sparen*{\nabla_S \mathcal{L}}\mu \text{,} \qquad\qquad
\nabla_{m_2} \mathcal{L} = \nabla_S \mathcal{L}
\end{align*}
where $m_1 = \mu$ and $m_2 = S + \mu \mu^\top$ are the expectation parameters of $q_\lambda$ \cite{khan2017conjugate}.
\end{prop_nat_grad_gaussian}
By applying the above proposition to \eqref{eq:LB_gradient_khan} and \eqref{eq:LB_gradient_martin}, the computation of the gradient of the LB with respect to $m$, reduces to evaluating $\nabla_m \log p\autobracket*{y\vert \theta}$, and thus to its gradients with respect to $\mu$ and $S+\mu \mu^\top$. In the next proposition these gradients are shown to be of a rather simple form.
\begin{prop_nat_grad_gaussian} \label{prop:exp_grad_gaussian_solved}
For $q_\lambda\autobracket*{\theta}$ being a Gaussian distribution with mean $\mu$ and covariance matrix $S$, with natural parameters $\lambda_1$, $\lambda_2$ and expectation parameters $m_1$, $m_2$, be $v = S^{-1}\autobracket*{\theta-\mu}$, then
\begin{align*}
\nabla_{\mu} \log q_\lambda\autobracket*{\theta} = v\text{,} \qquad \qquad \nabla_{S} \log q_\lambda\autobracket*{\theta} = -\frac{1}{2}\autobracket*{S^{-1} - v v^\top}\text{.}
\end{align*}
Furthermore,
\begin{align*}
\nabla_{m_1} \log q_\lambda\autobracket*{\theta} = S^{-1} \theta - v v^\top\mu \text{,} \qquad \qquad
\nabla_{m_2} \log q_\lambda\autobracket*{\theta}= \nabla_{S} \log q_\lambda\autobracket*{\theta} \text{.}
\end{align*}
\end{prop_nat_grad_gaussian}
A proof is provided in Appendix \ref{proof:grad_wrt_m}.
Proposition \ref{prop:exp_grad_gaussian_solved} enables the online computation of the updates for the posteriors' natural parameter, in particular
for a full-covariance Gaussian variational posterior the QBVI translates in the following updates for $\mu$ and $S^{-1}$:
\begin{align}
S^{-1}_{t+1} &= \autobracket*{1-\beta}S^{-1}_t +\beta \sparen*{S^{-1}_0 + \mathbb{E}_{q_\lambda}\sparen*{\autobracket*{S^{-1}_t-v_t v^\top_t}\log p\autobracket*{y\vert \theta}}}py}} \text{,} \label{eq:update_S_martin} \\
\mu_{t+1} &= \mu_t + \betaS_{t+1}[S^{-1}_0\autobracket*{\mu_0-\mu_t} + \mathbb{E}_{q_\lambda}\sparen*{v_t\log p\autobracket*{y\vert \theta}}] \text{,}\label{eq:update_MU_martin}
\end{align}
where $S_0$ and $\mu_0$ respectively denote the mean vector and variance-covariance matrix of the prior distribution. A proof is provided in Appendix \ref{proof:QBVI}. Note that the updates involve a total of $d+d^2$ parameters, and additional $d+d^2$ hyper-parameters are required for the prior specification.
The updates are further simplified under an isotropic Gaussian prior of mean zero $\mu_0=0$ and variance-covariance matrix $S_0 = I/\tau$, with $\tau>0$ a scalar precision parameter. This is a common prior in Bayesian inference, e.g. \cite{ganguly2021introduction,tran2021practical}, and BNN applications, e.g. \cite{graves2011practical,khan2018fastscalable,kingma2013auto}, that furthermore requires the specification of only one prior hyper-parameter.
Furthermore, assuming that the variation posterior is diagonal, the optimization problem reduces to the estimation of $2d$ parameters. Be $s$ and $s^{-1}$ the column vector corresponding to the diagonal of $S$ and $S^{-1}$ respectively, the QBVI update reads:
\begin{align}\label{eq:update_S_diagonal}
s^{-1}_{t+1} &= \autobracket*{1-\beta}s^{-1}_t +\beta \sparen*{\bm{\tau}_d + \mathbb{E}_{q_\lambda}\sparen*{\autobracket*{s^{-1}_t-v_t \odot v_t}\log p\autobracket*{y\vert \theta}}}py}} \\
\mu_{t+1} &= \mu_t+ \beta s_{t+1}\odot [- \tau \mu_t + \mathbb{E}_{q_\lambda}\sparen*{v_t\log p\autobracket*{y\vert \theta}}] \text{,} \nonumber
\end{align}
where $v_t = s^{-1}_t \odot\autobracket*{\theta-\mu_t}$, $\bm{\tau}_d = (\tau,\dots,\tau)^\top \in \mathbb{R}^d$, and $\odot$ is the element-wise product.
The expectations in \eqref{eq:update_S_martin} and \eqref{eq:update_MU_martin} depend on the chosen likelihood and its parametrization in terms of $\theta$, thus cannot be further simplified.
While the typical Bayesian framework is here resembled in terms of prior-posterior assumptions, $\log p\autobracket*{y\vert \theta}$ is unconstrained in its form and complexity, and conjugacy is not a requirement. The updates do not require the gradients of the likelihood, resulting in simple updates (especially under a diagonal posterior and isotropic prior) that exploit natural gradients w.r.t. the variational distribution. Automatic differentiation and backpropagation are here irrelevant potentially enabling the implementation of QBVI in low-level and basic programming languages. The QBVI thus constitutes a generic and ready-to-use solution for Bayesian inference in complex models under a Gaussian variational approximation.
The general QBVI algorithm is summarized in Algorithm \ref{alg:QBVIfull}. The exit function $f_{\text{exit}}$ is discussed in Appendix \ref{app:Z_Appendix_practical_recommendations}.
\begin{figure}
\begin{minipage}[t]{0.48\textwidth}
\begin{algorithm}[H]
\centering
\caption{General QBVI implementation}\label{alg:general_qbvi}
\footnotesize
\begin{algorithmic}[1]
\State \text{Set hyper-parameters: }$ 0<\beta< 1$, $N_s$
\State \text{Set priors and initial values: } $\eta$, $\lambda_0$
\State \text{Set:} $t=1$, $\text{Stop} = \texttt{false}$
\While{$\text{Stop} = \texttt{true}$}
\State \text{Generate: } $\theta_s \sim q_{\lambda_t}$, $s = 1\dots N_s$
\State $\hat{g} = 1/N_s\sum_s \nabla_m \sparen*{\log q_\lambda\autobracket*{\theta}} \log p\autobracket*{y\vert \theta}$
\State $\lambda_{t+1} = \autobracket*{1-\beta}\lambda_t + \beta\autobracket*{\eta + \hat{g}}$
\State $t = t+1$, $\text{Stop} = f_{\text{exit}}\autobracket*{\dots}$
\EndWhile
\end{algorithmic}
\end{algorithm}
\end{minipage}
\hfill
\begin{minipage}[t]{0.48\textwidth}
\begin{algorithm}[H]
\centering
\caption{QBVI, full-covariance Gaussian} \label{alg:QBVIfull}
\footnotesize
\begin{algorithmic}[1]
\State \text{Set hyper-parameters: }$ 0<\beta< 1$, $N_s$
\State \text{Set priors and initial values: } $S_0$, $\mu_0$, $S_1$, $\mu_1$
\State \text{Set:} $t=1$, $\text{Stop} = \texttt{false}$
\While{$\text{Stop} = \texttt{true}$}
\State \text{Generate: } $\theta_s \sim q_{\lambda_t}$, $s = 1\dots N_s$
\State $\hat{g}_S = 1/N_s\sum_s \autobracket*{S^{-1}_t \theta_s - v_t v_t^\top} \log p\autobracket*{y \vert \theta_s}$
\State $\hat{g}_\mu =1/N_s\sum_s v_t \log p\autobracket*{y\vert \theta_s}$
\State $S^{-1}_{t+1} \leftarrow \autobracket*{1-\beta}S^{-1}_t + \beta\sparen*{S^{-1}_0 + \hat{g}_S}$,
\State $\mu_{t+1} \leftarrow \mu_t +\beta S_{t+1}\sparen*{S^{-1}_0\autobracket*{\mu_0-\mu_t} + \hat{g}_\mu}$
\State $t = t+1$, $\text{Stop} = f_{\text{exit}}\autobracket*{\dots}$
\EndWhile
\end{algorithmic}
\end{algorithm}
\end{minipage}
\caption*{ \textit{Algorithm \ref{alg:general_qbvi}}. Line 6: estimates the expectation in \eqref{eq:LB_gradient_martin}, line 7: update \eqref{eq:QBVI_update}, line 9: an exit rule, e.g. $f_{\text{exit}}\autobracket*{\bar{\mathcal{L}_t},P,t}$ (see Appendix \ref{app:Z_Appendix_practical_recommendations}).
\textit{Algorithm \ref{alg:QBVIfull}}. Lines 6-7: estimate the expectations in \eqref{eq:update_S_martin} and \eqref{eq:update_MU_martin}, lines 8-9: updates \eqref{eq:update_MU_martin} and \eqref{eq:update_S_martin}, line 9: see Algorithm \ref{alg:general_qbvi}.}
\end{figure}
\subsection{Contributions, limitations and related methods}\label{subsec:related_works}
QBVI stands as an algorithm that shares several similarities with several existing alternatives while differing from them as well. In particular, it lies in the gap between feasible natural gradient approaches without requiring the FIM (VON \cite{khan2018fastyetsimple} rationale and exponential family properties) and black-box methods (BBVI\cite{ranganath2014black} rationale and the use of the score estimator) that do not require models' gradients. QBVI provides a solution in this direction. Importantly, the complexity of the underlying model or backbone network is of little relevance for the adoption of QBVI over any suitable likelihood. Whatever the complexity of the underlying backbone model, outputs are parsed to the likelihood. Only forward passes are involved, while QBVI does not constrain the form/complexity of the likelihood. Furthermore, the likelihood does not require differentiability.
In addition, in Appendix \ref{app:proofs} we ease the computation of the optimal control variate coefficient $c^\star_i$ providing an analytic solution for the denominator in \eqref{eq:ci_star} under a Gaussian variational posterior.\\
Among the limitations of our approach we recognize that though practical, the adoption of the score estimator leads to higher variances compared to methods that use the model's gradients. The trade-off is between the efficiency of the estimator and the possibility of feasibility obtaining the model's gradients. As for a large number of VI algorithms, QBVI does not handle the positive-definiteness of the covariance matrix: in Appendix \ref{app:Z_Appendix_Positive_definiteness} we discuss some remedies for the diagonal covariance case.
An extensive overview of our approach compared to related words is provided in Appendix ??. Here we mention that as references for empirically testing QBVI we adopt a Monte Carlo Markov Chain (MCMC) sampler the Black-Bock Variational Inference (BBVI) method of \cite{ranganath2014black}, the Cholesky Gaussian Variational Bayes (CGVB) \cite{titsias2014doubly}, and the Manifold Gaussian Variational Bayes (MGVB) of \cite{tran2021variational}. MCMC is indicative of the true non-variational posterior, BBVI is perhaps the close-most related method, using the log-score trick, but not natural parameters, for Gaussian variational CGVB updates the Cholesky factor with euclidean gradients while adoption of the reparametrization trick (thus model gradients) for controlling the variance of the estimated gradients of the LB. Lastly, MGVB is based on manifold optimization. It uses approximate natural gradients and has the advantage of guaranteeing positive-definiteness of the covariance updated throughout the iterations.
\section{Implementation aspects}
\subsection{Gradient estimation and control variates}
For the implementation of the algorithms \ref{alg:general_qbvi} and \ref{alg:QBVIfull} the central aspect to be addressed is the estimation of gradients' expectations and the use of control variates, for which we provide some new analytical results. This is tackled in the following subsection.
With $\lambda = \lambda_t$ and $\theta_s \sim q_\lambda$, $n = 1,\dots,N_s$, in order to make the update applicable in practice, expectations can be approximated via MC sampling by the following naive estimators:
\begin{align*}
\mathbb{E}_{q_\lambda}\sparen*{\autobracket*{S^{-1}_t - v_t v^\top_t} \log p\autobracket*{y\vert \theta} } &\approx \frac{1}{N_s}\sum_{s=1}^{N_s}\sparen*{ (S^{-1}_t - S^{-1}_t\autobracket*{\theta_s-\mu_t}\autobracket*{\theta_s-\mu_t}^\topS^{-1}_t) \log p\autobracket*{y \vert \theta_s}}\text{,}\\
\mathbb{E}_{q_\lambda}\sparen*{v_t \log p\autobracket*{y\vert \theta} } &\approx \frac{1}{N_s} \sum_{s=1}^{N_s}\sparen*{S^{-1}_t \autobracket*{\theta_s - \mu_t}\log p\autobracket*{y \vert \theta_s}} \text{,}
\end{align*}
where the matrix products are replaced with appropriate element-wise products for the diagonal posterior case.
The variance of the above naive MC estimators can be reduced through the use of Control Variates (CV) \cite{paisley2012variational,salimans2014using}, by approximating $\mathbb{E}_{q_\lambda}\sparen*{\nabla_m\sparen*{\log q_\lambda\autobracket*{\theta}}\log p\autobracket*{y\vert \theta}}gpy}$ with
$$
\frac{1}{N_s}\sum_{s=1}^{N_s} \nabla_{m_i}\sparen*{\log q\autobracket*{\theta_s}} \autobracket*{\log p\autobracket*{y\vert \theta_s} - c_i}_s} - c_i} \text{,}
$$
which is an unbiased estimator of the expected gradient but of equal or smaller variance than the naive MC average across samples. For $i=1,2$, the optimal $c_i$ minimizing the variance of the CV estimator is
\begin{equation}\label{eq:ci_star}
c_i^\star = {\,\text{Cov} \autobracket*{\nabla_{m_i}\sparen*{\log q\autobracket*{\theta}} \log p\autobracket*{y\vert \theta},\nabla_{m_i}\log q\autobracket*{\theta}}}/{\mathbb{V}\autobracket*{\nabla_{m_i}\log q\autobracket*{\theta}}} \text{.}
\end{equation}
To ensure the unbiasedness of the estimator, the samples used to estimate $c_i^\star$ must be independent of the ones used to estimate the gradient $\nabla_{m_i}\log q\autobracket*{\theta_s}$. In practice, this is straightforward, as at iteration $t$, $c_i^\star$ is estimated by reusing samples $\theta_s$ previously drawn from $q_{\lambda_{(t-1)}}$, earlier used to estimate the former CV gradient.
Under a variational Gaussian posterior, the variance terms $\mathbb{V}\autobracket*{\nabla_{m_i}\log q\autobracket*{\theta}}$ in \eqref{eq:ci_star} are analytically tractable and correspond to
\begin{equation*}
\mathbb{V}\autobracket*{\nabla_{m_1}\log q\autobracket*{\theta}} = \text{diag}\autobracket*{S^{-1}\autobracket*{S+D}S^{-1}}*{S+D}S^{-1}} \text{,} \qquad
\mathbb{V}\autobracket*{\nabla_{m_2}\log q\autobracket*{\theta}} = \frac{1}{4} \text{vec}^{-1}[\text{diag}\autobracket*{Q}] \text{.}
\end{equation*}
A proof is provided in Appendix \ref{proof:cv_var}, along with the definition of matrices $D$ and $Q$. The above restricts the MC-based estimation of $c_i^\star$ only to the covariance term in \eqref{eq:ci_star}, clearly intractable as it depends on the general form of $\log p\autobracket*{y\vert \theta}$. For the CV estimator, we suggest selecting $N_s \approx 100$, as a compromise between estimates' variance and performance. In BNN frameworks exploiting the reparametrization trick \cite{blundell2015weight} a single MC draw can suffice for providing satisfactory estimates of the gradients, but direct comparison is here difficult as it does require the gradient of the loss with respect to each network's parameter, and the expectation is usually taken with respect to a standard normal distribution. Within Appendix \ref{app:Z_Appendix_practical_recommendations} we empirically illustrate the importance of variance control and address the impact of $N_s$ on the variance of the gradients and on the LB.
\subsection{Further considerations for the implementation}
Besides the gradient estimation, there are a number of practicalities and actual implementation details that nevertheless are important that need to be put in place for a smooth implementation. Such standard provisions are collectively discussed in appendix \ref{app:Z_Appendix_practical_recommendations}. Separately in Appendix \ref{app:Z_MFQBVI}, we discuss the implementation of a mean-field variant of QBVI, while in Appendix \ref{app:Z_Appendix_Positive_definiteness} we tackle the issue related to positive definiteness of the variational covariance matrix.
\label{sec:practical}
\section{Experiments}
The paper introduces a new approach for tacking VI despite the form and complexity of the underlying backbone method. Even for the most elaborate deep-learning applications that may lead to outputs through a sequence of numerous and complex layers, for VB the outputs are parsed into a likelihood function, whose numerical values feed QBVI. In practical applications, the most common forms for the likelihood would those analogous to logistic and linear regression. Indeed a DL network for binary classification would return from its last layer class probabilities for which the very same likelihood of the standard logistic regression applies. Our experiments cover these two common situations and furthermore show that QBVI is applicable (as it should) to non-standard forms of the likelihood function (which is indeed not a part of QBVI, like it is for the class of methods using model's gradients). With this rationale, we decided to suggest experiments that are based on simple models, but relevant (in the form of the likelihood) for more elaborated cases. Extensive results and deeper analyses can be found in Appendix \ref{app:Z_experiments}, along with comparisons between QBVI and the baseline methods introduced in Secection \ref{subsec:related_works}.
We test the QBVI algorithm on different real-world data, in classification and prediction problems.
The most simple scenario for applying our method is that of the logistic regression. In this case, the updates \eqref{eq:update_MU_martin} and \eqref{eq:update_S_martin} immediately apply, and the likelihood has a simple form. We report our results on two datasets. The Labour datasets \cite{mroz1984sensitivity} consists of 753 individuals and 7 variables, with a binary response variable indicating whether the participants are currently in the labour force.\footnote{Publicly available data:\url{//www.key2stats.com/data-set/view/140}} A 75\%-25\% split is applied to extract the training and testing sets, hyperparameters are provided in Appendix \ref{app:Z_experiments}.
We explore and inspect for anomalies in the dynamics of the learning process for QBVI in Fig. \ref{fig:Labour_LB_mu}.
The smoothed LB ($\bar{\mathcal{L}}$ of Appendix \ref{app:Z_Appendix_practical_recommendations}) is increasing reaching a plateau in both training and test data after about 1000 epochs (left plot), accordingly the likelihood of the data given the model parameters and the likelihood of the variational posterior display a similar behavior (rightmost plot), suggesting solidity in the learning phase. Standard performance measures display high learning rates at a very early epoch and reach satisfactory steady levels at about epoch 300 already (middle panel). Figure 2 points out a robust and smooth underlying learning process with a relatively fast convergence of the lower bound, variational likelihood, and performance metrics, where the magnitude of MC sampling noise is furthermore relatively small compared to the growth rate of such curves.
Results corresponding to the posterior parameters at the maximum value of the LB are reported in Table 1:
despite the different optimization objectives (LB as opposed to log-likelihood), the data log-likelihood (LL) under the QBVI and maximum-likelihood (ML) estimates are very close. Performance measures are aligned as well, with rather negligible differences. Note that Table 1 is not meant to address whether the models are adequate and satisfactory in achieving the best performance metrics (LL or others, e.g. accuracy) for the given data and classification task, this is here out of scope. Rather it points out that for the chosen model, QBVI provides a full-Bayesian prescription where the posterior means well-align to the ML estimates leading to negligible differences in the models' LLs evaluated at the ML and QBVI estimates.
Figure \ref{fig:Labour_marginals} depicts the learning of the variational parameters and compares the marginal posterior across different models. From the top row, we observe that variational parameters are readily updated within a few iterations and relatively stable later.
For the marginals (bottom row), we observe a remarkable accordance between the QBVI approximation, the Monte-Carlo-Markov-Chain sampler of the true posterior, and the state-of-the-art Gaussian Manifold VI approach \cite{tran2021variational} taken as reference.
\begin{figure}[h!]
\centering
\includegraphics[width=0.9\textwidth,trim={4cm 0.5cm 4cm 0.8cm}]{Labour_Marginals.pdf}
\caption{Top panel. Variational posterior means (blue) and variances (orange) of the regression coefficients given in the titles. Bottom row. Variational marginals: QBVI (blue), MCMC samples of the true posterior (orange), MGVI \cite{tran2021variational}. Priors (part of) are overlapped (red).}
\label{fig:Labour_marginals}
\end{figure}
The above discussion applies as well to the Statlog (German Credit Data) Data Set\footnote{Publicly available data: \url{archive.ics.uci.edu/ml/datasets/statlog+(german+credit+data)}}, consisting of 24 categorical attributes for 1000 instances for the associated classification task of classifying good or bad credit risk.
Experiments of regression are performed on financial data. We consider the Heterogeneous Autoregressive model (HAR) or realized volatility \cite{corsi2009simple} and the GARCH(1,1) model \cite{bollerslev1986generalized}. For the two, we respectively use daily values of 5-minutes sub-sampled realized volatility \cite{zhang2005tale} and daily close returns for the S\&P 500 index from Jan 1, 2018 to May 15, 2022 (1087 entries).\footnote{Publicly available data: \url{realized.oxford-man.ox.ac.uk}} The regression results in the bottom part of Table \ref{tab:perf_meas} again show a remarkable alignment between QBVI and maximum likelihood performance measures: it is in the availability of (approximate) joint posterior where the actual difference between two approaches lies. Note that regression requires the estimation of the disturbances' variance, for which we use the mean-field approximation in Algorithm \ref{alg:MFQBVI}, see Appendix \ref{app:Z_MFQBVI}. For all the four datasets, further results are reported in Appendix \ref{app:Z_experiments}.
\begin{figure}
\centering
\includegraphics[width=0.6\textwidth,trim={4cm 0.5cm 4cm 1cm}]{Labour_LB.pdf}
\caption{Model learning across iterations. Lower bound (the maximum on the training set marked with a red triangle), performance measures, and loglikelihoods.}
\label{fig:Labour_LB_mu}
\end{figure}
\begin{table}
\centering
\caption{Performance measures for the different datasets. Comparison with the alternative method discussed in Section \ref{subsec:related_works} can be found in Appendix \ref{app:Z_experiments}.}
\scalebox{0.85}{
\begin{tabular}{lrrrrrrrr}
& \multicolumn{4}{c}{Labour data} & \multicolumn{4}{c}{Credit data} \\
\cmidrule(lr){2-5} \cmidrule(lr){6-9} & \multicolumn{2}{c}{Train set} & \multicolumn{2}{c}{Test set} & \multicolumn{2}{c}{Train set} & \multicolumn{2}{c}{Test set} \\
\cmidrule(lr){2-3} \cmidrule(lr){4-5} \cmidrule(lr){6-7} \cmidrule(lr){8-9} & QBVI & ML & QBVI & ML & QBVI & ML & QBVI & ML \\
\midrule
Precision & 0.710 & 0.708 & 0.698 & 0.698 & 0.760 & 0.753 & 0.776 & 0.780 \\
Recall & 0.701 & 0.699 & 0.676 & 0.676 & 0.713 & 0.706 & 0.717 & 0.721 \\
Accuracy & 0.711 & 0.709 & 0.612 & 0.612 & 0.791 & 0.785 & 0.656 & 0.667 \\
F1 & 0.703 & 0.701 & 0.746 & 0.746 & 0.728 & 0.721 & 0.815 & 0.816 \\
LL & -332.99 & -332.98 & -113.89 & -113.80 & -347.68 & -347.48 & -125.80 & -124.50 \\
\midrule
& \multicolumn{4}{c}{HAR model} & \multicolumn{4}{c}{GARCH model} \\
\cmidrule(lr){2-5} \cmidrule(lr){6-9} & \multicolumn{2}{c}{Train set} & \multicolumn{2}{c}{Test set} & \multicolumn{2}{c}{Train set} & \multicolumn{2}{c}{Test set} \\
\cmidrule(lr){2-3} \cmidrule(lr){4-5} \cmidrule(lr){6-7} \cmidrule(lr){8-9} & QBVI & ML & QBVI & ML & QBVI & ML & QBVI & ML \\
MSE $\times 10^5$ & 10946 & 10945 & 10169 & 10338 & 8.75 & 8.73 & 4.96 & 4.94 \\
LL & -249.65 & -249.60 & -73.69 & -74.22 & 2587.27 & 2587.34 & 878.02 & 878.04 \\
\bottomrule
\end{tabular}
}
\label{tab:perf_meas}
\end{table}\label{sec:experiments}
\section{Conclusion}
Whereas Bayesian inference is highly attractive in high-risk domains and applications, several difficulties hinder its use in ML and complex models from a wide audience. Whereas Variational Inference (VI) is an established effective approach for approximating the true posterior with a tractable one, the actual implementation of the parameters' updates is not straightforward. We introduce a Quasi Black-box VI (QBVI) method for the Bayesian learning under a Gaussian variational approximation. Our approach develops on the SGD algorithm with steps in the direction of the natural gradient, to optimize the lower bound on the log-likelihood.
Based on certain properties of the exponential family and of the Gaussian distribution,
we show that the generally complex natural gradient update can be feasibly approximated by sampling terms that require only function queries of models' likelihood, but not their gradients. Our approach extends the scope of Bayesian learning to the wide class of models for which gradients are difficult or costly to compute, and typical VI model-specific derivations are unfeasible. We provide details on the robust and practical implementation of QBVI update and test its performance on well-established datasets and models. Future research might develop in two directions, virtually extending the current research and overcoming its limitations.
On the other hand, from a theoretical perspective, it is relevant to enable the QBVI update for exploiting the information in structured or factor-decomposed covariance matrices, perhaps resulting in faster updates and increased efficiency, and explore extensions over the Gaussian variational framework. Among the limitations of our work are the trade-off between the convenience of the black-bock approach opposed to alternative making use of models' gradients thus achieving a lower approximation error of the stochastic gradients. Future research direction could explore feasible and effective approaches to variance reduction to be applied in combination, or in alternative, to control variates. In addition the positive-definite constraint is not handled within QBVI and future research could develop in this direction with aid of the manifold optimization theory.\label{sec:conclusion}
\section*{Acknowledgments}
The research received funding from the European Union’s Horizon 2020 research and innovation programme under the Marie Sk{\l}odowska-Curie project BNNmetrics (grant agreement No. 890690).
\newpage
\bibliographystyle{abbrvnat}
\bibliography{Biblio}
\newpage
\section*{Checklist}
\begin{enumerate}
\item For all authors...
\begin{enumerate}
\item Do the main claims made in the abstract and introduction accurately reflect the paper's contributions and scope?
\answerYes{}
\item Did you describe the limitations of your work?
\answerYes{Specifically across the text we mention: (i) Actual difficulty in dealing with full-covariance matrices in large-scale problems, leading to diagonal updates (see Sec. \ref{sec:proposed_method} and Sec. \eqref{eq:update_S_diagonal}), (ii) Difficulty in guaranteeing positive-definiteness in covariance updates (Appendix \ref{app:Z_Appendix_Positive_definiteness}), (iii) The simple natural gradient computation is applicable to Gaussian variational approximations only (see \ref{sec:proposed_method} and \ref{sec:conclusion}). The variance of the log-score estimator of the expected gradient may be higher compared to e.g. the reparametrization trick where models' gradients are used}.
\item Did you discuss any potential negative societal impacts of your work?
\answerNA{}
\item Have you read the ethics review guidelines and ensured that your paper conforms to them?
\answerYes{The ethic guidelines are not a concern for the proposed research and the form it is here presented.}
\end{enumerate}
\item If you are including theoretical results...
\begin{enumerate}
\item Did you state the full set of assumptions of all theoretical results?
\answerYes{Proofs in Appendix \ref{app:proofs}}.
\item Did you include complete proofs of all theoretical results?
\answerYes{}
\end{enumerate}
\item If you ran experiments...
\begin{enumerate}
\item Did you include the code, data, and instructions needed to reproduce the main experimental results (either in the supplemental material or as a URL)?
\answerYes{Codes included in the submission as additional material, and planned to be made available online in the future.}
\item Did you specify all the training details (e.g., data splits, hyperparameters, how they were chosen)?
\answerYes{See Sec.\ref{sec:experiments} and Appendix. \ref{app:Z_experiments}}
\item Did you report error bars (e.g., with respect to the random seed after running experiments multiple times)?
\answerYes{Rather than error bars we deal with distributions, and report details e.g. on the marginals.}
\item Did you include the total amount of compute and the type of resources used (e.g., type of GPUs, internal cluster, or cloud provider)?
\answerYes{We report the run-times in Table \ref{tab:par_credit}, for the largest model. Time statistics are deemed to be taken as indicative as MGVB is efficiently implemented in the BVayes-lab package, wheres out current implementation focuses back-testing. We are currently developing an efficient solution for GPU computations in Python. Codes will be made publicly available.}
\end{enumerate}
\item If you are using existing assets (e.g., code, data, models) or curating/releasing new assets...
\begin{enumerate}
\item If your work uses existing assets, did you cite the creators?
\answerYes{}
\item Did you mention the license of the assets?
\answerYes{}
\item Did you include any new assets either in the supplemental material or as a URL?
\answerYes{}
\item Did you discuss whether and how consent was obtained from people whose data you're using/curating?
\answerYes{}
\item Did you discuss whether the data you are using/curating contains personally identifiable information or offensive content?
\answerNA{}
\end{enumerate}
\item If you used crowdsourcing or conducted research with human subjects...
\begin{enumerate}
\item Did you include the full text of instructions given to participants and screenshots, if applicable?
\answerNA{}
\item Did you describe any potential participant risks, with links to Institutional Review Board (IRB) approvals, if applicable?
\answerNA{}
\item Did you include the estimated hourly wage paid to participants and the total amount spent on participant compensation?
\answerNA{}
\end{enumerate}
\end{enumerate}
\newpage