Extracted main text — title through conclusion, appendix excluded. This is what our citation measures are computed over, published so the extraction can be checked by eye.
52,507 characters · 17 sections · 61 citation commands
Quasi Black-Box Variational Inference with Natural Gradients for Bayesian Learning
Machine Learning (ML) techniques have been proved to be successfully in many prediction and classification tasks across natural-language processing young_recent_2018, computer vision krizhevsky_imagenet_2017, time-series langkvist2014review and finance applications 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 hoeting1999bayesian, and for enabling sequential learning freitas2000hierarchical while retaining prior and past knowledge. Although Bayesian principles have been proposed in ML decades ago mackay1992bayesian,mackay1995probable,lampinen2001bayesian, it has been only recently that fast and feasible methods boosted a growing use of Bayesian methods in complex models osawa_practical_2019,khan2018fastscalable,khan2018fastyetsimple.
The most challenging task is the computation of the posterior 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 saul1996mean,wainwright2008graphical,hoffman2013stochastic,blei2017variational. Furthermore, recent research advocates the use of the natural gradients for boosting the optimum search and the training wierstra2014natural, enabling fast and accurate Bayesian learning algorithms that are scalable and versatile. \\ Stochastic Gradient Descent (SGD) methods robbins1951stochastic promoted the use of VI for complex high-dimensional DL models hoffman2013stochastic,salimans2014using, yet they require significant implementation and tuning ranganath2014black. On the other hand, the use of natural gradients amari1998natural within the VI framework under exponential family approximations has been shown to be notably efficient and robust hoffman2013stochastic,wierstra2014natural,khan2018fastscalable, often leading to simple updates khan2018fastyetsimple. In particular, the common choice of a Gaussian approximation to the true posterior opper2009variational, results from bonnet1964transformations and price1958useful leads to simple updates of the variational mean and precision matrix khan2017conjugate. Several algorithms have been derived following these results 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. 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 ranganath2014black are of straightforward and general use relying on stochastic sampling without model-specific derivations salimans2014using,paisley2012variational,kingma2013auto. Even though black-box methods can benefit from numerous improvements such as variance reduction techniques 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 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 blei2017variational. We furthermore develop the method of Control variates paisley2012variational,salimans2014using for efficient large-MC sampling, and provide solutions for valid covariance updates adapting developments from existing methods khan2018fastscalable,tran2021variational,ong2018variational. Following and discussing the typical practicalities and recommendations in VI applications 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.
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}$,
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}$ tran2021practical,
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.
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$ 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$ 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 khan2018fastyetsimple. By replacing the Euclidean norm with the KL divergence, the SGD update results in the following natural gradient update:
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$,
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.
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
A proof is provided in Appendix (ref). By applying (ref) to (ref), we have the compact update khan2018fastscalable
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 tran2021practical, or exploiting a certain (assumed or imposed) structure of the FIM tran2020bayesian. Second, the expectation in (ref) is generally intractable as it depends on the specified log-likelihood. Thus, sampling methods can be used.
Assume $q_\lambda\autobracket*{\theta}$ belongs to an exponential family distribution. Its probability density function is parametrized as
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 (ref) 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 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 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 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
expressing the natural gradient in the natural parameter space as the gradient in the expectation parameter space. Eq. (ref) enables the computation of natural gradients in the natural parameter space without requiring the FIM khan2017conjugate,khan2018fastyetsimple. This translates into the following equivalence:
Our approach relies on the following equivalence for the natural gradient of the LB:
A proof is provided in Appendix (ref). The difference between equations (ref) and (ref) is substantial. Eq. (ref) 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. (ref) 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, (ref) 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) 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, 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. (ref) and (ref) therefore lead to the following update:
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) summarized the approach in the general case of an exponential-family prior-posterior pair of the same parametric form.
This section focuses on the implementation of the natural gradient update (ref) under the typical multivariate Gaussian variational framework, perhaps the most prominent in VI kucukelbir2015automatic,kingma2013auto,tran2021variational and Bayesian learning 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 (ref). Its density reads
where
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}{Proposition} \begin{prop_nat_grad_gaussian} For a differentiable function $\mathcal{L}$, and for $q_\lambda \sim \mathcal{N}\autobracket*{\mu,S}$,
where $m_1 = \mu$ and $m_2 = S + \mu \mu^\top$ are the expectation parameters of $q_\lambda$ khan2017conjugate. \end{prop_nat_grad_gaussian}
By applying the above proposition to (ref) and (ref), 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} 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
Furthermore,
\end{prop_nat_grad_gaussian} A proof is provided in Appendix (ref). Proposition (ref) 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}$:
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). 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. ganguly2021introduction,tran2021practical, and BNN applications, e.g. 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:
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 (ref) and (ref) 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). The exit function $f_{\text{exit}}$ is discussed in Appendix (ref).
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 khan2018fastyetsimple rationale and exponential family properties) and black-box methods (BBVIranganath2014black 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) we ease the computation of the optimal control variate coefficient $c^\star_i$ providing an analytic solution for the denominator in (ref) 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) 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 ranganath2014black, the Cholesky Gaussian Variational Bayes (CGVB) titsias2014doubly, and the Manifold Gaussian Variational Bayes (MGVB) of 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.
For the implementation of the algorithms (ref) and (ref) 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:
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) 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
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 (ref) are analytically tractable and correspond to
A proof is provided in Appendix (ref), 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 (ref), 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 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) 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.
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). Separately in Appendix (ref), we discuss the implementation of a mean-field variant of QBVI, while in Appendix (ref) we tackle the issue related to positive definiteness of the variational covariance matrix.
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), along with comparisons between QBVI and the baseline methods introduced in Secection (ref).
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 (ref) and (ref) immediately apply, and the likelihood has a simple form. We report our results on two datasets. The Labour datasets 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).
We explore and inspect for anomalies in the dynamics of the learning process for QBVI in Fig. (ref). The smoothed LB ($\bar{\mathcal{L}}$ of Appendix (ref)) 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) 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 tran2021variational taken as reference.
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 corsi2009simple and the GARCH(1,1) model bollerslev1986generalized. For the two, we respectively use daily values of 5-minutes sub-sampled realized volatility 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) 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), see Appendix (ref). For all the four datasets, further results are reported in Appendix (ref).
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.
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).