Variational Bayesian Approximation (VBA): A Comparison between Three Optimization Algorithms ^†

Abstract

In many Bayesian computations, we first obtain the expression of the joint distribution of all the unknown variables given the observed data. In general, this expression is not separable in those variables. Thus, obtaining the marginals for each variable and computing the expectations is difficult and costly. This problem becomes even more difficult in high dimensional quandaries, which is an important issue in inverse problems. We may then try to propose a surrogate expression with which we can carry out approximate computations. Often, a separable expression approximation can be useful enough. The variational Bayesian approximation (VBA) is a technique that approximates the joint distribution p with an easier, for example separable, distribution q by minimizing the Kullback–Leibler divergence $KL(q|p)$. When q is separable in all the variables, the approximation is also called the mean field approximation (MFA), and so q is the product of the approximated marginals. A first standard and general algorithm is the alternate optimization of $KL(q|p)$ with respect to $q_i$. A second general approach is its optimization in the Riemannian manifold. However, in this paper, for practical reasons, we consider the case where p is in the exponential family and so is q. For this case, $KL(q|p)$ becomes a function of the parameters $\theta$ of the exponential family. Then, we can use any other optimization algorithm to obtain those parameters. In this paper, we compare three optimization algorithms, namely a standard alternate optimization, a gradient-based algorithm and a natural gradient algorithm, and study their relative performances in three examples.

1. Introduction

In many applications, with direct or indirect observations, the use of the Bayesian computations starts with obtaining the expression of the joint distribution of all the unknown variables given the observed data. Then, we must use it for inference. In general, this expression is not separable in all the variables of the problem. So, the computations become hard and costly. For example, obtaining the marginals for each variable and computing the expectations are difficult and costly. This problem becomes even more crucial in high dimensional quandaries, which is an important issue in inverse problems. We may then need to propose a surrogate expression with which we can carry out approximate computations.

The variational Bayesian approximation (VBA) is a technique that approximates the joint distribution p with an easier, for example a separable one, q, by minimizing the Kullback–Leibler divergence $KL(q|p)$, which makes the marginal computations much easier. For example, in the case of two variables, $p(x,y)$ is approximated by $q(x,y)=q_1(x)q_2(y)$ via minimizing $KL(q_1{q}_2|p)$. When q is separable in all the variables of p, the approximation is also called mean field approximation (MFA).

To obtain the approximate marginals $q_1$ and $q_2$, we have to minimize $KL(q_1{q}_2|p)$. The first standard and general algorithm is the alternate optimization of $KL(q_1{q}_2|p)$ with respect to $q_1$ and $q_2$. By finding the expression of the functional derivatives of $KL(q_1{q}_2|p)$ with respect to $q_1$ and $q_2$ and then equating them to zero alternatively, we obtain an iterative optimization algorithm. A second general approach is its optimization in the Riemannian manifold. However, in this paper, for practical reasons, we consider the case where p is in the exponential family and so are $q_1$ and $q_2$. For this case, $KL(q_1{q}_2|p)$ becomes a function of the parameters $\theta$ of the exponential family. Then, we can use any other optimization algorithm to obtain those parameters.

In this paper, we compare three optimization algorithms: a standard alternate optimization (Algorithm 1), a gradient-based algorithm [1,2] (Algorithm 2) and a natural gradient algorithm [3,4,5] (Algorithm 3). The aim of this paper is to consider the first algorithm as the VBA method and compare it with the two other algorithms.

Of the main advantages of the VBA for inference problems, such as inverse problems and machine learning, we can mention the following:

-

First, VBA builds a sufficient model according to prior information and the final posterior distribution. Especially in the mean field approximation (MFA), the result ends in an explicit form for each unknown component using conjugate priors and works well for small sample sizes [6,7,8].

-

The second benefit is, for example in machine learning, that it is a robust way for classification based on the predictive posterior distribution and diminishes over-trained parameters [7].

-

The third privilege is that the target structure has uncertainty in the VBA recursive processes. This feature prevents further error propagation and increases the robustness of VBA [9].

Besides all these preponderances, the VBA has some weaknesses, such as difficulty regarding the solution of integrals and expectations in terms of obtaining a posterior distribution, and there is no evidence of finding an exact posterior [6]. Its most significant drawback arises when there are strong dependencies between unknown parameters, and the VBA ignores them. Then estimates, computed based on this approximation, may be very far from the exact values. However, it works well when the number of dependencies are low [8].

In this article, we examine three different estimating algorithms of the unknown parameters in a model concerning prior information. The first iterative algorithm is a standard alternate optimization based on VBA, which begins a certain initial points. Sometimes, the points are estimated from an available dataset, but most of the time, we do not have enough data on the parameters to make certain pre-estimations of them. To solve this obstacle, we can start the algorithm with certain desired points, and then by repeating the process, they approach the true values using the posterior distribution. The second two algorithms are gradient-based and natural gradient algorithms, whose base function is Kullback–Leibler divergence. First, the gradient of Kullback–Leibler for all unknown parameters must be found, then must start from points either estimated from data or desired choices. Then, we repeat the iterative algorithm until it converges to certain points. If we denote the unknown parameter space with $\theta$, then the recursive formula is ${\tilde{\theta }}^{(k+1)}={\tilde{\theta }}^{(k)}-\gamma \nabla KL({\tilde{\theta }}^{(k)})$ for gradient-based and natural gradient algorithms with different values of $\gamma$.

Additionally, we consider three examples, normal-inverse-gamma, multivariate normal and linear inverse problem for checking the performance and convergence speed of the algorithms.

We propose the following organization of this paper: In Section 2, we present a brief explanation of the basic analytical aspect of VBA. In Section 3, we explain our first example related to normal-inverse-gamma distribution analytically and, in practice, explain the outcomes of three algorithms to estimate the unknown parameters. In Section 5, we study a more complex example of a multivariate normal distribution whose means and variance–covariance matrix are unknown and have normal-inverse-Wishart distribution. The aim of this section is to demonstrate the marginal distributions of $\tilde{\mu }$ and $\tilde{\mathsf{\Sigma }}$ using a set of multivariate normal observations using the mean and variance. In Section 6, the example is closer to realistic situations and is a linear inverse problem. In Section 7, we present a summary of the work carried out in the article and compare the three recursive algorithms through three different examples.

2. Variational Bayesian Approach (VBA)

As we mentioned previously, VBA uses Kullback–Leibler divergence. Kullback–Leibler divergence [10] $KL(q|p)$ is an information measure of discrepancy between two probability functions defined as follows. Let $p(x)$ and $q(x)$ be two density functions of a continuous random variable x with respect to support set ${\mathbb{S}}_X$. $KL(q|p)$ function is introduced as:

$$KL(q|p)={\int }_{x\in {\mathbb{S}}_X}q(x)ln\frac{q(x)}{p(x)}\mathrm{d}x.$$

(1)

For simplicity, we assume a bivariate case of distribution $p(x,y)$ and want to assess it via VBA; therefore, we have:

$$KL(q|p)=-H(q_1)-H(q_2)-{\langle lnp(x,y)\rangle }_{q_1{q}_2},$$

(2)

where

$$H(q_1)=-{\int }_{x\in {\mathbb{S}}_X}{q}_1(x)ln{q}_1(x)\mathrm{d}x\mathrm{and}H(q_2)=-{\int }_{y\in {\mathbb{S}}_Y}{q}_2(y)ln{q}_2(y)\mathrm{d}y$$

are, respectively, the Shannon entropies of x and of y, and

$${\langle lnp(x,y)\rangle }_{q_1{q}_2}=\int {\int }_{(x,y)\in {\mathbb{S}}_{XY}}{q}_1(x)q_2(y)lnp(x,y)\mathrm{d}x\mathrm{d}y.$$

Now, differentiating the Equation (2) with respect to $q_1$ and then with respect to $q_2$ and equating them to zero, we obtain:

$$q_1(x)\propto exp\left\{{\langle lnp(x,y)\rangle }_{q_2(y)}\right\}\mathrm{and}{q}_2(y)\propto exp\left\{{\langle lnp(x,y)\rangle }_{q_1(x)}\right\}$$

(3)

These results can be easily extended to more dimensions [11]. They do not have any closed form because they depend on the expression of $p(x,y)$ and that of $q_1$ and $q_2$. An interesting case is that of exponential families and conjugate priors, where writing

$$p(x,y)=p(x|y)p(y),\mathrm{and}p(y|x)=\frac{p(x,y)}{p(x)}=\frac{p(x|y)p(y)}{p(x)},$$

(4)

we can consider $p(y)$ as prior, $p(x|y)$ as the likelihood, and $p(y|x)$ as the posterior distributions. Then, if $p(y)$ is a conjugate prior for the likelihood $p(x|y)$, then the posterior $p(y|x)$ will be in the same family as the prior $p(y)$. To illustrate all these properties, we provide details of these expressions for a first simple example of normal-inverse-gamma $p(x,y)=\mathcal{N}(x|\mu ,y)\mathcal{I}\mathcal{G}(y|\alpha ,\beta )$ with $q_1(x)=\mathcal{N}(x|\mu ,v)$ and $q_2=\mathcal{I}\mathcal{G}(y|\alpha ,\beta )$. For this simple case, first we give the expression of $KL(q|p)$ with $q_1(x)=\mathcal{N}(x|\tilde{\mu },\tilde{v})$ and $q_2(y)=\mathcal{I}\mathcal{G}(y|\tilde{\alpha },\tilde{\beta })$ as a function of the parameters $\theta =(\tilde{\mu },\tilde{v},\tilde{\alpha },\tilde{\beta })$ and then the expressions of the three above-mentioned algorithms; after which, we study their convergence.

3. Normal-Inverse-Gamma Distribution Example

The purpose of this section is to explain in detail the process of performing calculations in VBA. For this we consider a simple case for which we have all the necessary expressions. The objective here is to compare the three different algorithms mentioned above. Additionally, its practical application can be explained as follows:

We have a sensor which measures a quantity X, N times $x_1,\dots ,x_N$. We want to model these data. In a first step, we model them as $N(x|\mu ,v)$ with fixed $\mu$ and v. Then, it is easy to estimate the parameters $(\mu ,v)$ either by maximal likelihood or Bayesian strategy. If we assume that the model is Gaussian with unknown variance and call this variance y and assign an $\mathcal{I}\mathcal{G}$ prior to it, then we have a model $\mathcal{N}\mathcal{I}\mathcal{G}$ for $p(x,y)$. The $\mathcal{N}\mathcal{I}\mathcal{G}$ priors were applied to the wavelet context with correlated structures because they were conjugated with normal priors [12]. We chose the normal-inverse-gamma distribution because of this conjugated property and ease of handling.

We showed that the margins are $St$ and $\mathcal{I}\mathcal{G}$. Working directly with $St$ is difficult. So, we want to approximate it with a Gaussian $q_1(x)$. This is equivalent to approximating $p(x,y)$ with $q_1(x)q_2(y)$. Now, we want to find the parameters $\mu$, v, $\alpha$, and $\beta$, which minimize $KL(q_1{q}_2|p)$. This process is called VBA. Then, we want to compare three algorithms to obtain the parameters which minimize $KL(·|·)$. $KL(·|·)$ is convex with respect to $q_1$ if $q_2$ is fixed and is convex with respect to $q_2$ if $q_1$ is fixed. So, we hope that the iterative algorithm converges. However, $KL(·|·)$ may not be convex in the space of parameters. So, we have to study the shape of this criterion concerning the parameters $\tilde{v}$, $\tilde{\alpha }$ and $\tilde{\beta }$.

The practical problem considered here is the following: A sensors delivers a few samples $x=\{x_1,x_2,\cdots ,x_N\}$ of a physical quantity X. We want to find $p(x)$. For this process, we assume a simple Gaussian model but with unknown variance y. Thus, the forward model can be written as $p(x,y)=\mathcal{N}(x|\mu ,y)\mathcal{I}\mathcal{G}(y|\alpha ,\beta )$. In this simple example, we know that $p(x)$ is Student’s t-distribution obtained by:

$$S(x|\mu ,\alpha ,\beta )=\int \mathcal{N}(x|\mu ,y)\mathcal{I}\mathcal{G}(y|\alpha ,\beta )\mathrm{d}y$$

(5)

Our objective is to find the three parameters $\theta =(\mu ,\alpha ,\beta )$ from the data $x$ and an approximate marginal $q(x)$ for $p(x)$.

The main idea is to find such $q_1(x)q_2(y)$ as an approximation of $p(x,y)$. Here, we show the VBA, step by step. For this, we start by choosing the conjugate families $q_1(x)=\mathcal{N}(x|\tilde{\mu },\tilde{v})$ and $q_2(y)=\mathcal{I}\mathcal{G}(y|\tilde{\alpha },\tilde{\beta })$.

In the first step, we have to calculate $lnp(x,y)$

$$lnp(x,y)=c-\frac{1}{2}lny-\frac{1}{2y}{(x-\tilde{\mu })}^2-(\tilde{\alpha }+\frac{1}{2})lny-\frac{\tilde{\beta }}{y}.$$

(6)

where c is a constant value term independent of x and y. First of all, to use the iterative algorithm given in (3), starting by $q_1=N(x|{\mu }^{\prime },v^{\prime })$ we have to find $q_2(y)$, so we have to start by finding $q_2(y)$. The integration of $lnp(x,y)$ concerns $q_1(x)$

$${\langle lnp(x,y)\rangle }_{q_1}=c-\frac{1}{2y}{\langle {(x-\tilde{\mu })}^2\rangle }_{q_1}-(\tilde{\alpha }+1)lny-\frac{\tilde{\beta }}{y}.$$

(7)

Since the mean of x is the same in prior and posterior distribution, $\tilde{\mu }={\tilde{\mu }}^{\prime }$ then ${\langle {(x-\tilde{\mu })}^2\rangle }_{q_1}=\tilde{v}$, otherwise ${\langle {(x-\tilde{\mu })}^2\rangle }_{q_1}=\tilde{v}+{(\tilde{\mu }-{\tilde{\mu }}^{\prime })}^2$. Thus when $\tilde{\mu }={\tilde{\mu }}^{\prime }$

$$q_2(y)\propto exp\left[-(\tilde{\alpha }+1)lny-(\frac{\tilde{v}}{2}+\tilde{\beta })\frac{1}{y}\right].$$

(8)

Thus, the function $q_2(y)$ is equivalent to an inverse gamma distribution $\mathcal{I}\mathcal{G}(\tilde{\alpha },\frac{\tilde{v}}{2}+\tilde{\beta })$. Similarly, $q_2(y)$ is $\mathcal{I}\mathcal{G}(\tilde{\alpha },\frac{\tilde{v}+{(\tilde{\mu }-{\tilde{\mu }}^{\prime })}^2}{2}+\tilde{\beta })$ when $\tilde{\mu }\ne {\tilde{\mu }}^{\prime }$. We have to take integral of $lnp(x,y)$ over $q_2$ to find $q_1$

$${\langle lnp(x,y)\rangle }_{q_2}=c-(\tilde{\alpha }+1){\langle lny\rangle }_{q_2}-(\tilde{\beta }+\frac{1}{2}{(x-\tilde{\mu })}^2){\langle \frac{1}{y}\rangle }_{q_2}$$

(9)

Note that the first term does not depend on x and ${\langle \frac{1}{y}\rangle }_{q_2}=\frac{2\tilde{\alpha }}{2\tilde{\beta }+\tilde{v}}$, so

$$q_1(x)\propto exp[-\frac{2\tilde{\alpha }}{2\tilde{\beta }+\tilde{v}}(\tilde{\beta }+\frac{1}{2}{(x-\tilde{\mu })}^2)]\propto exp[-\frac{{(x-\tilde{\mu })}^2}{2\frac{2\tilde{\beta }+\tilde{v}}{2\tilde{\alpha }}}].$$

(10)

We see that $q_1$ is, again, a normal distribution but with updated parameters $\mathcal{N}(\tilde{\mu },\frac{2\tilde{\beta }+\tilde{v}}{2\tilde{\alpha }})$, so $\tilde{v}=\frac{2\tilde{\beta }+\tilde{v}}{2\tilde{\alpha }}$. Note that we obtained the conjugacy property: if $p(x|y)=\mathcal{N}(x|\mu ,y)$ and $p(y)=\mathcal{I}\mathcal{G}(y|\alpha ,\beta )$, then $p(y|x)=\mathcal{I}\mathcal{G}(y|{\alpha }^{\prime },{\beta }^{\prime })$ where ${\mu }^{\prime }$, ${\alpha }^{\prime }$ and ${\beta }^{\prime }$ are ${\mu }^{\prime }=\mu$, ${\alpha }^{\prime }=\alpha$, ${\beta }^{\prime }=\beta +\frac{2\beta +v}{2\alpha }$. In this case, we also know that $p(x|\alpha ,\beta )=St(x|{\mu }^{\prime },\alpha ,\beta )$.

In standard alternate optimization based on VBA (Algorithm 1), there is no need for an iterative process for $\tilde{\mu }$ and $\tilde{\alpha }$, which are approximated by $\tilde{\mu }={\mu }_0$ and $\tilde{\alpha }={\alpha }_0$, respectively. The situation for $\tilde{\beta }$ and $\tilde{v}$ is different because there are circular dependencies among them. So, the approximation needs an iterative process, staring from ${\tilde{\mu }}^{(1)}={\mu }_0$, ${\tilde{v}}^{(1)}=v_0$, ${\alpha }^{(1)}={\alpha }_0$ and ${\beta }^{(1)}={\beta }_0$. As a conclusion for this case is that the values of $\alpha$ and $\mu$ do not change during the iterations, and so depend on the initial values. However, the values of $\beta$ and v are interdependent and change during the iterations. This algorithm is summarized below.

Alternate optimization algorithm based on VBA (Algorithm 1):

$$\begin{array}{ccc}{\tilde{\alpha }}^{(k+1)}& =& {\tilde{\alpha }}^{(k)},\hfill \\ {\tilde{\beta }}^{(k+1)}& =& {\tilde{\beta }}^{(k)}+\frac{{\tilde{v}}^{(k)}}{2},\hfill \\ {\tilde{\mu }}^{(k+1)}& =& {\tilde{\mu }}^{(k)},\hfill \\ {\tilde{v}}^{(k+1)}& =& \frac{2{\tilde{\beta }}^{(k)}+{\tilde{v}}^{(k)}}{2{\tilde{\alpha }}^{(k)}}.\hfill \end{array}$$

This algorithm converges to $v=(2\tilde{\beta }+\tilde{v})/(2\tilde{\alpha })$, which gives $\tilde{v}=(2\tilde{\beta })/(2\tilde{\alpha }-1)$ and $\tilde{\beta }=0$, so $\tilde{v}=0$, which is a degenrate solution.

The two other algorithms, gradient- and natural gradient-based, require to find the expression of $KL(q_1{q}_2:p)$ as a function of the parameters $\theta =(\alpha ,\beta ,\mu ,v)$:

$$\begin{array}{cc}\hfill KL(\tilde{\theta })=& 2ln\mathsf{\Gamma }(\tilde{\alpha }-\frac{1}{2})-(2\tilde{\alpha }+\frac{3}{2}){\psi }_0(\tilde{\alpha }-\frac{1}{2})+\frac{(2\tilde{\alpha }-1)(\tilde{v}+2\tilde{\beta })}{4\tilde{\beta }}+\tilde{\alpha }+\frac{5}{2}ln\tilde{\beta }-\frac{1}{2}ln\tilde{v}-1.\hfill \end{array}$$

(11)

We also need the gradient expression of $\nabla KL(\tilde{\theta })$ for $\tilde{\theta }$:

$$\begin{array}{cc}\hfill \nabla KL(\tilde{\theta })=& \left(\frac{\tilde{v}}{2\tilde{\beta }}-(2\tilde{\alpha }+\frac{3}{2}){\psi }_1(\tilde{\alpha }-\frac{1}{2})+2,-(2\tilde{\alpha }-1)\frac{\tilde{v}}{4{\tilde{\beta }}^2}+\frac{5}{2\tilde{\beta }},0,\frac{2\tilde{\alpha }-1}{4\tilde{\beta }}-\frac{1}{2\tilde{v}}\right).\hfill \end{array}$$

(12)

As we can see, these expressions do not depend on $\tilde{\mu }$, so their derivatives with respect to $\tilde{\mu }$ are zero. The means are preserved.

Gradient and natural gradient algorithms (Algorithms 2 and 3):

$$\begin{array}{ccc}{\tilde{\alpha }}^{(k+1)}& =& {\tilde{\alpha }}^{(k)}-\gamma \left(\frac{{\tilde{v}}^{(k)}}{2{\tilde{\beta }}^{(k)}}-(2{\tilde{\alpha }}^{(k)}+\frac{3}{2}){\psi }_1({\tilde{\alpha }}^{(k)}-\frac{1}{2})+2\right),\hfill \\ {\tilde{\beta }}^{(k+1)}& =& {\tilde{\beta }}^{(k)}-\gamma \left(-(2{\tilde{\alpha }}^{(k)}-1)\frac{{\tilde{v}}^{(k)}}{4{\left[{\tilde{\beta }}^{(k)}\right]}^2}+\frac{5}{2{\tilde{\beta }}^{(k)}}\right),\hfill \\ {\tilde{\mu }}^{(k+1)}& =& {\tilde{\mu }}^{(k)},\hfill \\ {\tilde{v}}^{(k+1)}& =& {\tilde{v}}^{(k+1)}-\gamma \left(\frac{2{\tilde{\alpha }}^{(k)}-1}{4{\tilde{\beta }}^{(k)}}-\frac{1}{2{\tilde{v}}^{(k)}}\right).\hfill \end{array}$$

Here, $\gamma$ is fixed for the gradient algorithm and is proportional to $1/\parallel \nabla KL\parallel$ for the natural gradient algorithm.

4. Numerical Experimentations

To show the relative performances of these algorithms, we generate $n=100$ samples from the model $p(x,y)=\mathcal{N}(x|1,y)\mathcal{I}\mathcal{G}(y|3,1)$ for the numerical computations. Thus, it should be noted that we know the exact values of the unknown parameters which can be used to show the performances of the proposed algorithms. The following results: ${\tilde{\theta }}_1$, ${\tilde{\theta }}_2$ and ${\tilde{\theta }}_3$ are the estimated parameters using, respectively, Algorithm 1, Algorithm 2 and Algorithm 3. The contour plots of the corresponding probability density functions are shown in Figure 1 compared with original model.

$${\tilde{\theta }}_1=\left\{\begin{array}{c}\tilde{\mu }=0.044724,\hfill \\ \tilde{v}=0.590894,\hfill \\ \tilde{\alpha }=3.792337,\hfill \\ \tilde{\beta }=1.765735\hfill \end{array}\right.{\tilde{\theta }}_2=\left\{\begin{array}{c}\tilde{\mu }=0.044724,\hfill \\ \tilde{v}=7.910504,\hfill \\ \tilde{\alpha }=4.991621,\hfill \\ \tilde{\beta }=3.002594\hfill \end{array}\right.{\tilde{\theta }}_3=\left\{\begin{array}{c}\tilde{\mu }=0.044724,\hfill \\ \tilde{v}=0.423761,\hfill \\ \tilde{\alpha }=4.415239,\hfill \\ \tilde{\beta }=0.706049.\hfill \end{array}\right.$$

All three algorithms try to minimize the same criterion. So, the objectives are all the same, but the number of steps may differ. The requirements must reach the minimum $KL(·)$. In this simple example of the normal-inverse-gamma distribution, the convergence step numbers of VBA, gradient-based and natural gradient are 1, 2 and 1 using moment initializations. The overall performance of the standard alternate optimization (VBA) is more precise than any other. The poorest estimation is from a gradient-based algorithm. So, the algorithms are able to approximate the joint density function with a separable one but with different accuracy. In the following section, we will tackle a more complex model.

5. Multivariate Normal-Inverse-Wishart Example

In previous section, we explained how to preform VBA in order to approximate a complicated joint distribution function by tractable marginal factorials using a simple case study. In this section, a multivariate normal case $p(x)=\mathcal{N}(x|\tilde{\mu },\tilde{\mathsf{\Sigma }})$ is considered, which is approximated by $q(x)={\prod }_i\mathcal{N}(x_i|{\tilde{\mu }}_i,{\tilde{v}}_i)$ for different shapes of the covariance matrix $\tilde{\mathsf{\Sigma }}$.

We assume that the basic structure of an available dataset is multivariate normal with unknown mean vector $\tilde{\mu }$ and variance–covariance matrix $\tilde{\mathsf{\Sigma }}$. Their joint prior distribution is a normal-inverse-Wishart distribution of $\mathcal{N}\mathcal{I}\mathcal{W}(\tilde{\mu },\tilde{\mathsf{\Sigma }}|{\tilde{\mu }}_0,\tilde{\kappa },\tilde{\mathsf{\Psi }},\tilde{\nu })$, which is the generalized form of the classical $\mathcal{N}\mathcal{I}\mathcal{G}$. The posteriors are multivariate normal for mean vector and inverse-Wishart in the variance–covariance matrix. Since the normal-inverse-Wishart distribution is a conjugate prior distribution for multivariate normal, the posterior distribution of $\tilde{\mu }$ and $\tilde{\mathsf{\Sigma }}$ again belongs to the same family, and their corresponding margins are

$$\mathcal{M}\mathcal{N}\left(\frac{\tilde{\kappa }{\tilde{\mu }}_0+n\overline{x}}{\tilde{\kappa }+n},\frac{1}{\tilde{\kappa }+n}\tilde{\mathsf{\Lambda }}\right),\mathcal{I}\mathcal{W}\left(\tilde{\mathsf{\Lambda }}+\tilde{\mathsf{\Psi }}+\sum _{i=1}^n(x_i-\overline{x}){(x_i-\overline{x})}^T,\tilde{\nu }+n\right),$$

(13)

where n is the sample size. To present the performance of the three algorithms, we examined on a dataset based on $x\sim \mathcal{N}\mathcal{I}\mathcal{W}(x|\mu ,\mathsf{\Sigma })$, whose parameters have the following low-density structure:

$$\mu \sim \mathcal{M}\mathcal{N}\left(\mu |\left[\begin{array}{c}2\\ 1\end{array}\right],\frac{1}{2}\left[\begin{array}{cc}3& -1\\ -1& 1\end{array}\right]\right),\mathsf{\Sigma }\sim \mathcal{I}\mathcal{W}\left(\mathsf{\Sigma }|\left[\begin{array}{cc}3& -1\\ -1& 1\end{array}\right],6\right).$$

We used only the data of $x$ in the estimation processes. The results of algorithms are presented in Figure 2 along with the true contour plot of the model. The VBA estimation is the most separable distribution compared with gradient and natural gradient methods. The next best case is the natural gradient algorithm, but its weakness is transferring the dependency slightly to the approximation. The results for the gradient-based algorithm show the dependency completely, along with its inability to obtain a separable model.

6. Simple Linear Inverse Problem

Finally the third example is the case of linear inverse problems with $g=Hf+ϵ$ with priors $p(ϵ)=\mathcal{N}(ϵ|0,v_{ϵ}I)$, $p(f)=\mathcal{N}(f|0,\mathrm{diag}\left[v\right])$ and $p(v|\alpha ,\beta )={\prod }_j\mathcal{I}\mathcal{G}(v_j|\alpha ,\beta )$, where $f=[f_1,f_2,\cdots ,f_N]$ and $v=[v_{f_1},v_{f_2},\cdots ,v_{f_N}]$. Using these priors, we get $p(f,v|g)\propto p(g|f,v_{ϵ})p(f|v)p(v)$ with $p(g|f,v_{ϵ})=\mathcal{N}(g|Hf,v_{ϵ}I)$. See [13] for details.

Thus, the joint distribution of $g$, $f$, and $v$:

$$p(g,f,v)\propto p(g|f,v_{ϵ})p(f|v)p(v).$$

(14)

is approximated by $q(g,\tilde{f},\tilde{v})=q_1(g|\tilde{f})q_2(\tilde{f})q_3(\tilde{v})$ using the VBA method. Even if the main interest is the estimation of $q_2(\tilde{f})$, but in the recursive process, $q_1(g|\tilde{f})$ and $q_3(\tilde{v})$ are also updated. For simplicity, we suppose that the transfer matrix $H$ is an identical matrix $I$. The final outputs are as follows:

$$\begin{array}{c}\hfill \tilde{f}\sim \mathcal{M}\mathcal{N}\left(\frac{{\tilde{\mu }}_{\tilde{f}}}{1+\frac{2\widetilde{v_{ϵ}}\tilde{\alpha }}{n(\tilde{v}+{\tilde{\mu }}_{\tilde{f}}^2+2\tilde{\beta })}},\mathrm{diag}[\frac{\widetilde{v_{ϵ}}(\tilde{v}+{\tilde{\mu }}_{\tilde{f}}^2+2\tilde{\beta })}{n(\tilde{v}+{\tilde{\mu }}_{\tilde{f}}^2)+2n\tilde{\beta }+2\widetilde{v_{ϵ}}\tilde{\alpha }}]\right),\\ \hfill g\sim \mathcal{N}({\tilde{\mu }}_{\tilde{f}},\widetilde{v_{ϵ}}I),\mathrm{and}{\tilde{v}}_j\sim \mathcal{I}\mathcal{G}\left({\tilde{\alpha }}_j,\frac{{\tilde{v}}_j+{\tilde{\mu }}_{{\tilde{f}}_j}^2}{2}+{\tilde{\beta }}_j\right),j=1,\cdots ,N.\end{array}$$

(15)

For simulations, we chose a model to examine the performance of these margins and compare them with gradient and natural gradient algorithms. The selected model is $g=Hf+ϵ$ with the following assumptions:

$$H=I,f\sim \mathcal{M}\mathcal{N}(f|0,\mathrm{diag}\left[v_1,v_2\right]),v_1\sim \mathcal{I}\mathcal{G}(v_1|3,2),v_2\sim \mathcal{I}\mathcal{G}(v_2|4,3),ϵ\sim \mathcal{M}\mathcal{N}(ϵ|0,I).$$

(16)

In the assessment procedure, we do not apply the above information. The output of algorithms are shown in Figure 3, as well as the actual contour plot. In this example, the best diagnosis is from the natural gradient algorithm. The VBA by construction is separable and cannot be the same as the original.

7. Conclusions

This paper presents an approximation method of the unknown density functions for hidden variables called VBA. It is compared with gradient and natural gradient algorithms. We also consider three examples normal-inverse-gamma, normal-inverse-Wishart and linear inverse problem. We provided details of the first model and showed examples of two other examples throughout the whole paper. In all three models, the parameters are unexplored and need to be estimated by recursive algorithms. We attempted to approximate the joint complex distribution with a simpler version of the margin factorials that appeared to be independent cases. The VBA and natural gradient converged fairly early. The major discrepancy in algorithms comes from the accuracy of the results. They estimate the intricate joint distribution with separable ones. Here, the best overall performance of VBA is demonstrated.