On the Convergence of Extended Variational Inference for Non-Gaussian Statistical Models

TL;DR

This paper analyzes the convergence properties of extended variational inference (EVI) in non-Gaussian models, comparing approximation strategies and conditions, and demonstrating the advantages of the single lower-bound approach through theoretical and experimental results.

Contribution

It clarifies the convergence conditions of EVI, compares approximation strategies, and highlights the benefits of the single lower-bound method in non-Gaussian Bayesian models.

Findings

SLB approximation shows superior convergence properties.

Weak and strong conditions have distinct theoretical implications.

Extensive experiments validate the advantages of the SLB approach.

Abstract

Variational inference (VI) is a widely used framework in Bayesian estimation. For most of the non-Gaussian statistical models, it is infeasible to find an analytically tractable solution to estimate the posterior distributions of the parameters. Recently, an improved framework, namely the extended variational inference (EVI), has been introduced and applied to derive analytically tractable solution by employing lower-bound approximation to the variational objective function. Two conditions required for EVI implementation, namely the weak condition and the strong condition, are discussed and compared in this paper. In practical implementation, the convergence of the EVI depends on the selection of the lower-bound approximation, no matter with the weak condition or the strong condition. In general, two approximation strategies, the single lower-bound (SLB) approximation and the multiple lower-bounds (MLB) approximation, can be applied to carry out the lower-bound approximation. To clarify the differences between the SLB and the MLB, we will also discuss the convergence properties of the aforementioned two approximations. Extensive comparisons are made based on some existing EVI-based non-Gaussian statistical models. Theoretical analysis are conducted to demonstrate the differences between the weak and the strong conditions. Qualitative and quantitative experimental results are presented to show the advantages of the SLB approximation.

Tables

Table 1. TABLE I: Required conditions for EVI.

	Auxiliary function	Form of the Auxiliary Function	Systematic Gap
Strong condition	$p(𝐗,𝐙)\ge {\tilde{p}}_s(𝐗,𝐙)$	${𝐄}_{𝐙{\backslash }_{{𝐙}_i}}\left[\ln \tilde{p}(𝐗,𝐙)\right]\approxeq \ln p_i\left({𝐙}_i\right)$ †	${𝒢}_{\text{s}}>{𝒢}_{\text{w}}$
Weak condition	${𝐄}_{𝐙}\left[\ln p(𝐗,𝐙)\right]\ge {𝐄}_{𝐙}\left[\ln {\tilde{p}}_w(𝐗,𝐙)\right]$

Table 2. TABLE II: Comparisons of the objective functions for Bayesian BMM.

Model	Parameters	${ℒ}_{\text{SLB}}-{ℒ}_{\text{MLB}}$	${\text{KL}}_{\text{SLB}}-{\text{KL}}_{\text{MLB}}$
A	${\pi }_1=0.3,u_1=2,v_1=8$	$3.6\times {10}^{-3}$	$-2.8\times {10}^{-3}$
	${\pi }_2=0.7,u_2=15,v_2=4$
B	${\pi }_1=0.3,u_1=10,v_1=2$	$1.3\times {10}^{-3}$	$-0.58\times {10}^{-3}$
	${\pi }_2=0.4,u_2=2,v_2=12$
	${\pi }_3=0.3,u_3=10,v_3=10$

Equations

$$\mathcal{L}={\mathbf{E}}_{\mathbf{Z}}\left[\ln p(\mathbf{X},\mathbf{Z})-\ln q(\mathbf{Z})\right],$$

$$p(\mathbf{Z}|\mathbf{X})\approx q(\mathbf{Z})=\prod _{i=1}^M{q}_i({\mathbf{Z}}_i),$$

$$\ln q_i^{\ast }({\mathbf{Z}}_i)={\mathbf{E}}_{\mathbf{Z}{\backslash }_{{\mathbf{Z}}_i}}\left[\ln p(\mathbf{X},\mathbf{Z})\right]+\text{const}.$$

$${\mathbf{E}}_{\mathbf{Z}}\left[\ln p(\mathbf{X},\mathbf{Z})\right]\ge {\mathbf{E}}_{\mathbf{Z}}\left[\ln \tilde{p}(\mathbf{X},\mathbf{Z})\right],$$

$$\begin{array}{rl}\mathcal{L}=& {\mathbf{E}}_{\mathbf{Z}}\left[\ln p(\mathbf{X},\mathbf{Z})\right]-{\mathbf{E}}_{\mathbf{Z}}\left[\ln q(\mathbf{Z})\right]\\ \ge & {\mathbf{E}}_{\mathbf{Z}}\left[\ln \tilde{p}(\mathbf{X},\mathbf{Z})\right]-{\mathbf{E}}_{\mathbf{Z}}\left[\ln q(\mathbf{Z})\right]\\ \triangleq & \tilde{\mathcal{L}}.\end{array}$$

$$\ln {\tilde{q}}_i^{\ast }({\mathbf{Z}}_i)={\mathbf{E}}_{\mathbf{Z}{\backslash }_{{\mathbf{Z}}_i}}\left[\ln \tilde{p}(\mathbf{X},\mathbf{Z})\right]+\text{const}.$$

$$p(\mathbf{X},\mathbf{Z})\ge \tilde{p}(\mathbf{X},\mathbf{Z}).$$

$$\begin{array}{rl}p(\mathbf{X},\mathbf{Z})\ge & {\tilde{p}}_1(\mathbf{X},{\mathbf{Z}}_1)\\ ⋮& \\ p(\mathbf{X},\mathbf{Z})\ge & {\tilde{p}}_i(\mathbf{X},{\mathbf{Z}}_i)\\ ⋮& \\ p(\mathbf{X},\mathbf{Z})\ge & {\tilde{p}}_M(\mathbf{X},{\mathbf{Z}}_M).\end{array}$$

$$\ln {\tilde{q}}_i^{\ast }({\mathbf{Z}}_i)={\mathbf{E}}_{\mathbf{Z}{\backslash }_{{\mathbf{Z}}_i}}\left[\ln {\tilde{p}}_i(\mathbf{X},{\mathbf{Z}}_i)\right]+\text{const}.$$

$$\begin{array}{rl}{\tilde{\mathcal{L}}}_1=& {\mathbf{E}}_{\mathbf{Z}}\left[\ln {\tilde{p}}_1(\mathbf{X},{\mathbf{Z}}_1)-\ln q(\mathbf{Z})\right]\\ {\tilde{\mathcal{L}}}_2=& {\mathbf{E}}_{\mathbf{Z}}\left[\ln {\tilde{p}}_2(\mathbf{X},{\mathbf{Z}}_2)-\ln q(\mathbf{Z})\right].\end{array}$$

$$\ln {\tilde{q}}_1^{\ast }({\mathbf{Z}}_1)=$$

$$\ln {\tilde{q}}_2^{\ast }({\mathbf{Z}}_2)=$$

$$2\times \mathcal{L}\ge$$

$$=$$

$$+$$

$$\begin{array}{rl}f(\mathbf{X};\mathbf{\Pi },\mathbf{U},\mathbf{V})=& \prod _{n=1}^N\sum _{i=1}^I{\pi }_i\mathrm{Beta}(x_n;u_i,v_i),\end{array}$$

$$\mathrm{Beta}(x;u,v)=\frac{\Gamma (u+v)}{\Gamma (u)\Gamma (v)}{x}^{u-1}(1-x{)}^{v-1},u,v>0.$$

$$\begin{array}{rl}f(\mathbf{X},\mathbf{Z}|\mathbf{U},\mathbf{V},\mathbf{\Pi })=& f(\mathbf{X}|\mathbf{U},\mathbf{V},\mathbf{\Pi },\mathbf{Z})f(\mathbf{Z}|\mathbf{\Pi })\\ =& f(\mathbf{X}|\mathbf{U},\mathbf{V},\mathbf{Z})f(\mathbf{Z}|\mathbf{\Pi })\\ =& \prod _{n=1}^N\prod _{i=1}^I{\left[{\pi }_i\mathrm{Beta}(x_n|u_i,v_i)\right]}^{z_{ni}}.\end{array}$$

$${\mathbf{E}}_{u_i,v_i}\left[\mathrm{LIB}(u_i,v_i)\right]={\mathbf{E}}_{u_i,v_i}\left[\frac{\Gamma (u_i+v_i)}{\Gamma (u_i)\Gamma (v_i)}\right].$$

$$\begin{array}{rl}& {\mathbf{E}}_{\mathbf{Z}}\left[\ln {\tilde{p}}_w(\mathbf{X},\mathbf{Z})\right]\\ =& {\tilde{\mathcal{L}}}_{\text{SLB}}\\ =& \ln \frac{\Gamma ({\overline{u}}_i+{\overline{v}}_i)}{\Gamma ({\overline{u}}_i)\Gamma ({\overline{v}}_i)}\\ & +{\overline{u}}_i\left[\psi ({\overline{u}}_i+{\overline{v}}_i)-\psi ({\overline{u}}_i)\right](\mathbf{E}\left[\ln u_i\right]-\ln {\overline{u}}_i)\\ & +{\overline{v}}_i\left[\psi ({\overline{u}}_i+{\overline{v}}_i)-\psi ({\overline{v}}_i)\right](\mathbf{E}\left[\ln v_i\right]-\ln {\overline{v}}_i),\end{array}$$

$$\begin{array}{rl}{\tilde{p}}_{s_{u_i}}(\mathbf{X},\mathbf{Z})=& \ln \frac{\Gamma ({\overline{u}}_i+{\overline{v}}_i)}{\Gamma ({\overline{u}}_i)\Gamma ({\overline{v}}_i)}\\ & +{\overline{u}}_i\left[\psi ({\overline{u}}_i+{\overline{v}}_i)-\psi ({\overline{u}}_i)\right](\ln u_i-\ln {\overline{u}}_i)\\ & +{\overline{v}}_i\left[\psi ({\overline{u}}_i+{\overline{v}}_i)-\psi ({\overline{v}}_i)\right](\ln v_i-\ln {\overline{v}}_i)\\ & +{\overline{u}}_i{\overline{v}}_i{\psi }^{{}^{\prime }}({\overline{u}}_i+{\overline{v}}_i)(\ln u_i-\ln {\overline{u}}_i),\end{array}$$

$$\begin{array}{rl}{\tilde{\mathcal{L}}}_{{\text{MLB}}_{u_i}}=& {\mathbf{E}}_{\mathbf{Z}}\left[{\tilde{p}}_{s_{u_i}}(\mathbf{X},\mathbf{Z})\right]\\ =& \ln \frac{\Gamma ({\overline{u}}_i+{\overline{v}}_i)}{\Gamma ({\overline{u}}_i)\Gamma ({\overline{v}}_i)}\\ & +{\overline{u}}_i\left[\psi ({\overline{u}}_i+{\overline{v}}_i)-\psi ({\overline{u}}_i)\right](\mathbf{E}\left[\ln u_i\right]-\ln {\overline{u}}_i)\\ & +{\overline{v}}_i\left[\psi ({\overline{u}}_i+{\overline{v}}_i)-\psi ({\overline{v}}_i)\right](\mathbf{E}\left[\ln v_i\right]-\ln {\overline{v}}_i)\\ & +{\overline{u}}_i\cdot {\overline{v}}_i\cdot {\psi }^{{}^{\prime }}({\overline{u}}_i+{\overline{v}}_i)(\mathbf{E}\left[\ln u_i\right]-\ln {\overline{u}}_i).\end{array}$$

$$\begin{array}{rl}{\tilde{\mathcal{L}}}_{{\text{MLB}}_{v_i}}=& \ln \frac{\Gamma ({\overline{u}}_i+{\overline{v}}_i)}{\Gamma ({\overline{u}}_i)\Gamma ({\overline{v}}_i)}\\ & +{\overline{u}}_i\left[\psi ({\overline{u}}_i+{\overline{v}}_i)-\psi ({\overline{u}}_i)\right](\mathbf{E}\left[\ln u_i\right]-\ln {\overline{u}}_i)\\ & +{\overline{v}}_i\left[\psi ({\overline{u}}_i+{\overline{v}}_i)-\psi ({\overline{v}}_i)\right](\mathbf{E}\left[\ln v_i\right]-\ln {\overline{v}}_i)\\ & +{\overline{u}}_i\cdot {\overline{v}}_i\cdot {\psi }^{{}^{\prime }}({\overline{u}}_i+{\overline{v}}_i)(\mathbf{E}\left[\ln v_i\right]-\ln {\overline{v}}_i).\end{array}$$

$$\begin{array}{rl}& {\tilde{p}}_{s_{z_{ni}}}(\mathbf{X},\mathbf{Z})\\ =& \ln \frac{\Gamma ({\overline{u}}_i+{\overline{v}}_i)}{\Gamma ({\overline{u}}_i)\Gamma ({\overline{v}}_i)}\\ & +{\overline{u}}_i\left[\psi ({\overline{u}}_i+{\overline{v}}_i)-\psi ({\overline{u}}_i)\right](\ln u_i-\ln {\overline{u}}_i)\\ & +{\overline{v}}_i\left[\psi ({\overline{u}}_i+{\overline{v}}_i)-\psi ({\overline{v}}_i)\right](\ln v_i-\ln {\overline{v}}_i)\\ & +0.5\cdot {\overline{u}}_i^2\left[{\psi }^{{}^{\prime }}({\overline{u}}_i+{\overline{v}}_i)-{\psi }^{{}^{\prime }}({\overline{u}}_i)\right](\ln u_i-\ln {\overline{u}}_i{)}^2\\ & +0.5\cdot {\overline{v}}_i^2\left[{\psi }^{{}^{\prime }}({\overline{u}}_i+{\overline{v}}_i)-{\psi }^{{}^{\prime }}({\overline{v}}_i)\right](\ln v_i-\ln {\overline{v}}_i{)}^2\\ & +{\overline{u}}_i\cdot {\overline{v}}_i\cdot {\psi }^{{}^{\prime }}({\overline{u}}_i+{\overline{v}}_i)(\ln u_i-\ln {\overline{u}}_i)(\ln v_i-\ln {\overline{v}}_i).\end{array}$$

$$\begin{array}{rl}& {\tilde{\mathcal{L}}}_{{\text{MLB}}_{z_{ni}}}\\ =& {\mathbf{E}}_{\mathbf{Z}}\left[{\tilde{p}}_{s_{z_{ni}}}(\mathbf{X},\mathbf{Z})\right]\\ =& \ln \frac{\Gamma ({\overline{u}}_i+{\overline{v}}_i)}{\Gamma ({\overline{u}}_i)\Gamma ({\overline{v}}_i)}\\ & +{\overline{u}}_i\left[\psi ({\overline{u}}_i+{\overline{v}}_i)-\psi ({\overline{u}}_i)\right](\mathbf{E}\left[\ln u_i\right]-\ln {\overline{u}}_i)\\ & +{\overline{v}}_i\left[\psi ({\overline{u}}_i+{\overline{v}}_i)-\psi ({\overline{v}}_i)\right](\mathbf{E}\left[\ln v_i\right]-\ln {\overline{v}}_i)\\ & +0.5\cdot {\overline{u}}_i^2\left[{\psi }^{{}^{\prime }}({\overline{u}}_i+{\overline{v}}_i)-{\psi }^{{}^{\prime }}({\overline{u}}_i)\right]\mathbf{E}\left[(\ln u_i-\ln {\overline{u}}_i{)}^2\right]\\ & +0.5\cdot {\overline{v}}_i^2\left[{\psi }^{{}^{\prime }}({\overline{u}}_i+{\overline{v}}_i)-{\psi }^{{}^{\prime }}({\overline{v}}_i)\right]\mathbf{E}\left[(\ln v_i-\ln {\overline{v}}_i{)}^2\right]\\ & +{\overline{u}}_i\cdot {\overline{v}}_i\cdot {\psi }^{{}^{\prime }}({\overline{u}}_i+{\overline{v}}_i)(\mathbf{E}\left[\ln u_i\right]-\ln {\overline{u}}_i)(\mathbf{E}\left[\ln v_i\right]-\ln {\overline{v}}_i).\end{array}$$

$$\begin{array}{rl}\Delta {\tilde{\mathcal{L}}}_{\text{SLB}\text{vs.}{\text{MLB}}_{u_i}}=& {\tilde{\mathcal{L}}}_{\text{SLB}}-{\tilde{\mathcal{L}}}_{{\text{MLB}}_{u_i}}\\ =& -{\bar{u}}_i{\bar{v}}_i{\psi }^{\prime }({\bar{u}}_i+{\bar{v}}_i)(\mathbf{E}\left[\ln v_i\right]-\ln {\bar{v}}_i)\\ \ge & 0,\end{array}$$

$$\begin{array}{rl}& \Delta {\tilde{\mathcal{L}}}_{\text{SLB}\text{vs.}{\text{MLB}}_{z_{ni}}}\\ =& {\tilde{\mathcal{L}}}_{\text{SLB}}-{\tilde{\mathcal{L}}}_{{\text{MLB}}_{z_{ni}}}\\ =& -\{0.5\cdot {\overline{u}}_i^2\left[{\psi }^{{}^{\prime }}({\overline{u}}_i+{\overline{v}}_i)-{\psi }^{{}^{\prime }}({\overline{u}}_i)\right]\mathbf{E}\left[(\ln u_i-\ln {\overline{u}}_i{)}^2\right]\\ & +0.5\cdot {\overline{v}}_i^2\left[{\psi }^{{}^{\prime }}({\overline{u}}_i+{\overline{v}}_i)-{\psi }^{{}^{\prime }}({\overline{v}}_i)\right]\mathbf{E}\left[(\ln v_i-\ln {\overline{v}}_i{)}^2\right]\\ & +{\overline{u}}_i\cdot {\overline{v}}_i\cdot {\psi }^{{}^{\prime }}({\overline{u}}_i+{\overline{v}}_i)(\mathbf{E}\left[\ln u_i\right]-\ln {\overline{u}}_i)(\mathbf{E}\left[\ln v_i\right]-\ln {\overline{v}}_i)\}.\end{array}$$

Full text

1 Introduction

Gaussian distribution is the ubiquitous probability distribution used in statistics, signal processing, and pattern recognition areas [1]. However, not all the data being processed are Gaussian distributed [2]. In many real-life applications, the distribution of data is asymmetric and, therefore, is not Gaussian distributed [3]. For example, the image pixel values [4, 5], the reviewer’s rating to an item in a recommendation system [6, 7, 8], and the DNA methylation level data [9] are distributed in a range with bounded support. The diversity gain over the $K_{G}$ fading [10] and the periodogram coefficients in speech enhancement [11, 12] are semi-bounded (nonnegative). The spatial fading correlation [13] and the yeast gene expressions [14] have directional property so that the $l_{1}$ norm equals one. In signal processing, the acoustic noise with colored spectra [15] and the measurement noise in state-space model [16] are heavy-tailed. In the stock market, the asymptotic behavior of the first-order autoregressive (AR) process is clearly non-Gaussian [17] and the underlying Bayesian copula model for the stock index series are non-Gaussian as well [18]. Although the above mentioned data represent diverse characteristics, a common property is that, these data not only have specific support range, but also have non-bell distribution shape. The natural properties of Gaussian distribution (the definition domain is unbounded and the distribution shape is symmetric) do not fit such data well. Hence, these data are non-Gaussian distributed. It has been found in recent studies that explicitly utilizing the non-Gaussian characteristics can significantly improve the practical performance [2, 4, 5, 9, 10, 11, 13, 12, 14, 15, 16]. Hence, it is of particular importance and interest to make thorough studies of the non-Gaussian data and non-Gaussian statistical models.

Bayesian analysis plays an essential role in parameter estimation of statistical models [19, 20, 21, 22]. Unlike the conventionally used maximum-likelihood (ML) estimation [23], Bayesian estimation assumes that the parameters potentially follow underlying distributions and derives the posterior distributions of the parameters by applying the Bayes’ theorem [24] through combining the prior distributions with the likelihood function obtained from the observed data [21, 25]. Estimation of the posterior distribution via Bayesian estimation has several advantages over the ML estimation. Firstly, it gives statistical description to the parameters, rather than the simple point estimate that yield by the ML estimation. This makes Bayesian estimation more robust and reliable, by including the resulting uncertainty into the estimation [22]. Secondly, it can potentially prevent the overfitting problem, which is one of the drawbacks the ML estimation suffers. This is mainly due to the advantage of Occam’s razor effect in Bayesian estimation. Last but not the least, Bayesian estimation can estimate the model complexity automatically from the data. In ML estimation, model complexity decision usually requires cross validations and, therefore, is computationally costly [21, 23].

Varitional inference (VI) framework, among others, is a widely used strategy to infer the posterior distribution of the parameters in Bayesian analysis [21, 26, 27, 28]. In a full Bayesian model where all the parameters are assigned with prior distributions, we minimize the Kullback-Leibler (KL) divergence of the true posterior distribution from the approximating one to obtain an optimal approximation to the posterior distribution [21, Ch. 10]. This procedure is equivalent to maximizing the lower-bound to the marginal likelihood (model evidence). The optimal posterior distribution can be obtained by iteratively updating one variable (or one variable group) while fixing the rest. However, unlike the famous Gaussian distribution [21, 29, 1], most of the non-Gaussian statistical models (e.g., beta mixture model (BMM) [5, 4], Dirichlet mixture model (DMM) [31], von-Mises Fisher mixture model (VMM) [14], beta-Gamma nonnegative matrix factorization (BG-NMF) [6]) do not have analytically tractable solution to estimate the posterior distribution of the parameters. Numerical methods, e.g., Newton-Raphson algorithm, Gibbs sampling, Markov Chain Monte Carlo, are usually employed to sample from the posterior distribution [27, 5]. Numerical method often depends on Markov chain convergence and is in general computationally costly, especially in the high-dimensional space [32].

Recently, an improved framework, namely the extended variational inference (EVI) [33, 34, 35, 4, 31, 14, 6], has become popular in solving the above mentioned problem. Similar as the VI framework, EVI also seeks an optimal approximation to the posterior distribution. The difference is that EVI relaxes the objective function (the evidence lower-bound to the marginal likelihood) by constructing lower-bound approximation to the objective function. This lower-bound relaxation, which uses the convexity or relative convexity [36] of the objective function, can yield analytically tractable solution so that the parameter estimation is facilitated. Although systematic bias has been introduced due to the lower-bound approximation, several works have demonstrated the advantages of EVI in Bayesian estimation of statistical models [34, 4, 14, 6]. In Bayesian estimation of BMM, Ma et al. [4] derived an analytically tractable solution which outperforms the numerical Gibbs sampling based method [5]. As an extension work, Bayesian estimation of DMM via EVI have been proposed in [31], respectively. For directional data, von-Mises Fisher distribution is an important model in several applications. Analytically tractable solution to Bayesian estimation of VMM has been proposed by using EVI to provide lower-bound approximation [14]. For non-negative matrix factorization (NMF), EVI was also applied in deriving analytically tractable solutions for Poisson process (discrete) NMF [38], Gamma process NMF in music recording [34], and beta-Gamma NMF for bounded support data [6].

Convergence is an important issue in parameter estimation algorithm. For VI-based method, the objective function maximized during each iteration is convex or relatively convex in terms of the target variable’s posterior distribution [21, 36]. Hence, convergence is theoretically guaranteed. In EVI, the introduced lower-bound approximation to the objective function can be obtained either via a single extension over the whole variable group or multiple extensions, one for a subset of the whole variable group. Based on this, two lower-bound approximation strategies are obtained, one is the single lower-bound (SLB) approximation [34, 14, 6] and the other is the multiple lower-bounds (MLB) approximation [4, 31]. For EVI with SLB approximation, convergence is also guaranteed because the original objective function is replaced by one single lower-bound and this new objective function (i.e., the single lower-bound to the original objective function in VI) is convex (or relatively convex) and maximized during each iteration. However, when applying EVI with MLB approximation, the variable group is divided into different disjoint subsets and there exists different lower-bound approximations to the objective function. During each iteration, different lower-bounds, one for each variable subset, are maximized iteratively. Since the new objective function is not unique, convergence can not be theoretically guaranteed.

In order to clarify the convergence property of the EVI framework, we will discuss and summarize the conditions that required in EVI implementation. The SLB and MLB approximations will also be analyzed and compared qualitatively and quantitatively. Experimental results based on the recently proposed EVI-based BMM and DMM estimation algorithms will be presented to demonstrate the advantages of the SLB approximations. We draw some conclusions in the end.

2 Variational Inference and Extended Variational Inference

2.1 Variational Inference

In Bayesian estimation, a universal solution to the variational inference (VI) framework [26] is to approximate the posterior distribution by a product of several factor distributions and then update each factor distribution individually [21]. This method is the so-called factorized approximation (FA) which was developed from the mean field theory in physics [40]. With the FA method, the variational objective function that we want to maximize can be represented as the negative KL divergence as

[TABLE]

where $\mathbf{X}$ is the observed data, $\mathbf{Z}$ denotes all the random variables. If $\mathbf{Z}$ can be (approximately) factorized into $M$ disjoint groups as $\mathbf{Z}=\left\{\mathbf{Z}_{1},\ldots,\mathbf{Z}_{i},\ldots,\mathbf{Z}_{M}\right\}$ and we approximate the true posterior distribution $p(\mathbf{Z}|\mathbf{X})$ as

[TABLE]

the optimal solution can be written as

[TABLE]

The operator $\mathbf{E}_{\mathbf{Z}\backslash_{\mathbf{Z}_{i}}}$ means expectation with respect all the variables in $\mathbf{Z}$, except for $\mathbf{Z}_{i}$. If the optimal solution to the posterior distribution of $\mathbf{Z}_{i}$, which is $\ln q_{i}^{*}(\mathbf{Z}_{i})$ in (3), has the same logarithmical form as the prior distribution, the conjugate match between the prior and the posterior distributions are satisfied. Then we have obtained an analytically tractable solution. However, this conjugate match is not satisfied in most of the practical problems [4, 31, 6]. This is due to the fact that the optimal solution depends on the expectation computed with respect to the factor distribution [21].

2.2 Extended Variational Inference

In order to satisfy the conjugate match requirement, some approximations can be applied to get a nearly optimally analytically tractable solution. Braun et al. [35] considered the zeroth-order and first-order delta method for moments [39] to derive an alternative for the objective function to simplify the calculation. Blei et al. [33] proposed a correlated topic model (CTM) and used a first-order Taylor expansion to preserve a bound such that an intractable expectation was avoided. Similar idea was also applied in [4, 31, 14] for approximating the posterior distributions in BMM, DMM, and VMM, respectively. Using Jensen’s inequality has become commonplace in variational inference. In [34], the concavity of the function $-x^{-1}$ and the convexity of $-\log x$ were studied and the Jensen’s inequality and the first-order Taylor expansion were applied to approximately calculated the posterior distribution. Moreover, the EVI strategy was also applied in low rank matrix approximation area [6], where the Taylor expansion and Jensen’s inequality were both applied for the purpose of deriving analytically tractable solution.

All the aforementioned works utilized the following property. Given an auxiliary function $\widetilde{p}(\mathbf{X},\mathbf{Z})$ which satisfies

[TABLE]

the variational objective function (see [21], pp. 465 for more details) can be lower-bounded as

[TABLE]

Then we can maximize $\mathcal{\widetilde{L}}$, which is an lower-bound to the original objective function $\mathcal{\widetilde{L}}$, to asymptotically reach the maximum value of $\mathcal{{L}}$ [34, 4, 31]. The approximated optimal solution in this case is written as

[TABLE]

This method is the so-called EVI framework [33, 34, 35, 4, 31, 14, 6]. Although it introduces systematic gap when involving the lower-bound approximation, the EVI allows more flexibility when calculating intractable integrations in non-Gaussian statistical models and provides a convenient way to obtain an analytically tractable solution.

3 Convergence of EVI

3.1 Weak Condition and Strong Condition

As mentioned in Sec. 2.2, finding an auxiliary function $\widetilde{p}(\mathbf{X},\mathbf{Z})$ is an essential yet difficult part in EVI implementation. Generally speaking, this auxiliary function should satisfy the relation presented in (4) or it should satisfy

[TABLE]

It is obvious that an auxiliary function satisfies (7) should also satisfy (4). Hence, the condition in (4) is named as the weak condition and the one in (7) is referred to as the strong condition. When using an auxiliary function to lower-bound the original objective function, the EVI will introduce a systematic gap. Generally speaking, the gap incurred by the applying weak condition is relatively smaller than that introduced by using the strong condition. Fig. 1 illustrates the different gaps introduced by the weak and strong conditions, respectively.

It is worthwhile to note that the auxiliary function $\widetilde{p}(\mathbf{X},\mathbf{Z})$ is not necessary to be a normalized probability density function (PDF)111Actually, an auxiliary function that satisfies the strong condition cannot be a normalized PDF, as ${p}(\mathbf{X},\mathbf{Z})$ itself is a normalized PDF.. This will not affect the final solution since either VI or EVI will re-normalize the obtained optimal posterior distribution in the end.

In practice, in addition to the above mentioned weak or strong condition, an auxiliary function should also has a specific mathematical form so that the optimal solution in (6) has the same logarithmic form as the prior distribution and the conjugate match between the prior and the posterior distributions is satisfied. This is another required condition for choosing the auxiliary function. Table I lists the required conditions when implementing EVI.

Generally speaking, it is usually not feasible to find an auxiliary function that satisfies the strong condition, except that the original function $p(\mathbf{X},\mathbf{Z})$ is globally concave in terms of $\mathbf{Z}$ 222According to our experience, globally concavity holds only for Gaussian distribution. For (most of) the non-Gaussian statistical models, the original function is not globally concave.. Compared to the strong condition, it is easy to find an auxiliary function to fulfill the weak condition, although the ordinal function $p(\mathbf{X},\mathbf{Z})$ might be partially concave with respect to part of $\mathbf{Z}$ [6]. For example, the multivariate log-inverse-beta (MLIB) function in the Dirichlet distribution is not globally concave in terms of all of its variables. It is only relatively concave w.r.t. one of its variable when fixing the rest. Iteratively taking this property, an auxiliary function that satisfies the weak condition and the requirement of the mathematical form can be found so that an analytically tractable solution was derived. Moreover, the weak conditions yields smaller systematic gap. Therefore, the weak condition is more preferable in practice.

In summary, in order to apply the EVI to derive an analytically tractable solution for the Bayseisn estimation of non-Gaussian statistical models, an auxiliary function should 1) satisfies either the weak or the strong condition and 2) have the same mathematical form as the prior distribution (up to a constant difference).

3.2 SLB Approximation and MLB Approximation

If we can find an auxiliary function $\widetilde{p}(\mathbf{X},\mathbf{Z})$ that contains all the variables $\mathbf{Z}$ and satisfies the aforementioned required conditions, the convergence of EVI is naturally guaranteed as this new objective function is convex or relatively convex in terms of $q_{i}(\mathbf{Z}_{i})$ [21]. Since only one lower-bound approximation is applied to the original objective function, this approach is referred to as the single lower-bound (SLB) approximation and has been applied in, e.g., [14, 6].

When dividing $\mathbf{Z}$ into $M$ disjoint groups as $\mathbf{Z}=\left\{\mathbf{Z}_{1},\ldots,\mathbf{Z}_{i},\ldots,\mathbf{Z}_{M}\right\}$, there might exist several auxiliary functions. For example, we could have $M$ auxiliary functions as

[TABLE]

This approach is referred to as the multiple lower-bound (MLB) approximation. As each of the above mentioned auxiliary functions satisfies the required conditions in Sec. 3.1, the optimal solution in (6) is

[TABLE]

In this case, the new objective function that maximized during each iteration is not unique. Hence, there is no globally objective function that is maximized during each iteration. The convergence cannot be theoretically guaranteed. Such procedure has been applied in [4] and [31]. Although it is not guaranteed theoretically, the convergence was observed empirically.

Let’s study a simple case with two disjoint groups in the MLB approximation. Assuming that $\mathbf{Z}=\{\mathbf{Z}_{1},\mathbf{Z}_{2}\}$ and we have two auxiliary functions $\widetilde{p}_{1}(\mathbf{X},\mathbf{Z}_{1})$ and $\widetilde{p}_{2}(\mathbf{X},\mathbf{Z}_{2})$ for $\mathbf{Z}_{1}$ and $\mathbf{Z}_{2}$, respectively. As mentioned above, two different lower-bounds are obtained as

[TABLE]

If we maximize each lower-bound separately, the optimal solutions to these two disjoint groups are

[TABLE]

With these solutions, it looks like what we are maximizing is just two times of the original lower-bound as

[TABLE]

When performing the update strategy (11a), we get (12b) to be maximized. This maximization makes the distribution of $\mathbf{Z}_{1}$ to be less uncertain. As $-\mathbf{E}_{\mathbf{Z}}\left[\ln q(\mathbf{Z})\right]$ in (12c) is the differential entropy of $\mathbf{Z}$, (12c) is decreasing while (12b) is maximizing. It is hard to evaluate if (12b) changes more than (12c) or not. Thus, the overall lower-bound, i.e., $\mathcal{\widehat{L}}_{1}+\mathcal{\widehat{L}}_{2}$ in (12a), might decrease during some iterations. On the one hand, as the lower-bound (i.e., $\mathcal{\widetilde{L}}_{1}+\mathcal{\widetilde{L}}_{2}$) to the original objective function can not be guaranteed to be maximized all the time, this strategy may not promise convergence. On the other hand, if the change to (12b) is larger than that to (12c), the convergency is still guaranteed. There is no general judgement for the convergence. It should be studied case by case. Similar arguments can be applied to the case with more than two auxiliary functions. Thus, the convergency of MLB approximation is underdetermined.

In summary, SLB approximation can theoretically guarantee the convergence while MLB approximation, in general, cannot promise convergence.333In practice (e.g., [4, 31]), the EVI-based algorithm may also converge with MLB approximation. However, it is empirical result without proof.

4 Experimental Results and Discussions

Recently, several EVI-based parameter estimation algorithms for non-Gaussian statistical models have been proposed. Among others, the EVI-based Bayesian BMM [4] and the EVI-based Bayesian DMM [31] took the strong condition and the related analytically tractable solutions were derived with MLB approximation. With SLB approximation, the improved work about EVI-based Dirichlet mixture model was proposed. Regarding the non-negative matrix factorization, Hoffman et al. and Ma et al. proposed the EVI-based strategies for musical signal [34] and bounded support data [6], respectively. The EVI-based von Mises-Fisher mixture model was proposed in [14] where a structural factorization was considered. For all the aforementioned SLB approximation-based method, the weak condition is fulfilled.

In this section, we first compare the weak and strong conditions quantitatively. Secondly, we intensively compare the performance of the MLB approximation-based methods with the SLB approximation-based methods.

4.1 *Comparisons of Weak and Strong Condition *

Since Dirichlet distribution is a multivariate case of beta distribution, the EVI-based Bayesian BMM that constructs a auxiliary function with weak condition can be obtained based on the work in [37] by simply setting the dimension $K=2$. The EVI-based Bayesian BMM proposed in [4] utilized the strong condition to choose the auxiliary function. We compare these two different methods to demonstrate the differences between the strong and weak conditions.

Following the same notation in [4], we denote a multivariate Bayesian BMM with observation data $\mathbf{X}$ as

[TABLE]

where $\pi_{i}$ is the mixture weigh for the $i$th mixture component and $\mathrm{Beta}(x;u,v)$ is the beta distribution, which can be denoted as

[TABLE]

We consider the observation ${x}_{n}$ and the unobserved indication vector $\mathbf{z}_{n}$ as the complete data. The conditional distribution of $\mathbf{X}=\{{x}_{1},\ldots,{x}_{N}\}$ and $\mathbf{Z}=\{\mathbf{z}_{1},\ldots,\mathbf{z}_{N}\}$ given the latent variables $\{\mathbf{U},\mathbf{V},\mathbf{\Pi}\}$ is

[TABLE]

The ultimate goal is to estimate the posterior distributions of $u_{i}$, $v_{i}$, and $z_{ni}$, respectively.

In order to derive an analytically tractable solution for the posterior distributions, the most challengeable part with the EVI framework is to calculate the expectation of the bivariate log-inverse-beta (LIB) function

[TABLE]

4.1.1 EVI-based Bayesian BMM with Weak Condition [37]

In the Bayesian BMM with SLB approximation 444A Bayesian BMM with SLB approximation can be derived from the Bayesian DMM with SLB approximation [37] by setting the dimensions of the Dirichlet variable equal to two., the new objective function that we are maximizing is

[TABLE]

where $\overline{x}$ is the expected value of $x$ and $\psi(x)$ is the digamma function defined as $\psi(x)=\frac{\partial\ln\Gamma(x)}{\partial x}$. This lower-bound satisfies the weak condition such that $\mathbf{E}_{\mathbf{Z}}\left[\ln p(\mathbf{X},\mathbf{Z})\right]\geq\mathbf{E}_{\mathbf{Z}}\left[\ln\widetilde{p}_{w}(\mathbf{X},\mathbf{Z})\right]$. Moreover, this lower-bound is identical for all the variables $u_{i}$, $v_{i}$, and $z_{ni}$

4.1.2 EVI-based Bayesian BMM with Strong Condition [4]

For the case with strong condition, an auxiliary function $\widetilde{p}_{s}(\mathbf{X},\mathbf{Z})$ is required. In [4], three different auxiliary functions were derived for the variables $u_{i}$, $v_{i}$, and $z_{ni}$, respectively. To specify, for $u_{i}$, the auxiliary function is

[TABLE]

where $\psi^{{}^{\prime}}(x)=\frac{\partial\psi(x)}{\partial x}$. Hence, when considering $u_{i}$ as the variable, the objective function that was maximized is [4]

[TABLE]

Similarly, due to the symmetry of $u_{i}$ and $v_{i}$, the objective function, when treating $v_{i}$ as the variable, is [4]

[TABLE]

When taking $z_{ni}$ as the only variable, the auxiliary function that proposed in [4] is

[TABLE]

Correspondingly, the objective function for updating the posterior distribution of $z_{ni}$ can be represented as

[TABLE]

It has been analyzed in Sec. 3.1 that both the strong condition and the weak condition incur systematic gaps. We now quantitatively compare the gaps. It is worth to note that the EVI-based Bayesian BMM with strong condition is also a MLB approximation. We focus only on the comparisons of weak and strong conditions in thie section. The comparisons about the SLB approximation with the MLB approximation will be presented in the next section.

When taking $u_{i}$ as the variable, the difference between the objective functions obtained via weak and strong conditions, respectively, can be calculated as

[TABLE]

where we used the fact that $\psi^{{}^{\prime}}(x)>0$ and $\ln x$ is a convex function in terms of $x$. For $v_{i}$, it is straightforward to show the difference is also positive by using the symmetric properties.

When comparing $\mathcal{\widetilde{L}}_{\text{SLB}}$ with $\mathcal{\widetilde{L}}_{\text{MLB}_{z_{ni}}}$, the difference is

[TABLE]

It can be proved that the difference $\Delta\mathcal{\widetilde{L}}_{\text{SLB}\ \text{vs.}\ \text{MLB}_{z_{ni}}}$ is also greater than or equal to [math]. More details for this proof can be found in Appendix LABEL:Appendix-1.

The aforementioned three positive differences indicate that the new objective function with weak condition [37] is tighter (i.e., closer to the original objective function) than that with strong condition [4]. Thus, for the EVI-based Bayesian BMM, the systematic gap incurred by the weak condition is smaller than that incurred by the strong condition. This makes the weak condition more favorable in practice [37, 14, 6, 34].

Similar analysis can be applied to the Bayesian DMM with MLB [31] and the Bayesian DMM with SLB [37], as Dirichlet distribution is an multivariate extension of beta distribution.

4.2 Comparisons of MLB and SLB Approximations

In the previous section, we analyzed and compared the weak and the strong conditions for the EVI framework. Another important issue in EVI implementation is to distinguish the MLB and the SLB approximations, as the latter one can guarantee convergence but the first one may not. To this end, we compare the MLB approximation-based algorithm with the SLB approximation-based algorithm in this section.

4.2.1 Observations of Non-convergence

As discussed in Sec. 3.2, the convergence of the MLB method is not guaranteed. We ran the MLB approximation-based Bayesian BMM algorithm [4] and Bayesian DMM algorithm [31], respectively, and monitored the value of the objective function during each iteration. It can be observed that, for some rounds of simulations 555Here, one simulation round means that we ran the estimation algorithm until it stops according to some criterion., the objective function is decreasing during some iterations. This phenomenon has been observed for several times, both for BMM and DMM. Figure 3 lists the decreasing objective function values and the corresponding iterations. For the SLB approximation-based Bayesian BMM and Bayesian DMM [37], the monitored objective function was always increasing until converging. The observation of non-convergence demonstrates that the convergence with MLB approximation is underdetermined.

4.2.2 Comparisons of Estimation Accuracy

In this section, we compare the MLB approximation with the SLB approximation quantitatively. With a known BMM or DMM, $2,000$ samples were generated, respectively. The above-mentioned Bayesian estimation algorithms were applied to estimate the posterior distributions, respectively. We calculated the original variational objective function in (1) to examine which approximation is better. With the obtained posterior distribution $q^{*}(\mathbf{Z})$, the original variational objective function is calculated numerically by sampling method. Hence,we got two different values, ${\mathcal{{L}}_{\text{SLB}}}$ and $\mathcal{{L}}_{\text{MLB}}$, from the SLB approximation and the MLB approximation, respectively. Larger value means closer lower-bound approximation. In addition to this, we also measure the estimation accuracy by the KL divergence of the estimated PDF from the true one as $\text{KL}(p(\mathbf{X}|\boldsymbol{\Theta})\|p(\mathbf{X}|\boldsymbol{\widehat{\Theta}}))$, where $\boldsymbol{\Theta}$ is the true parameter vector and $\boldsymbol{\widehat{\Theta}}$ is the estimated one. Similarly, we numerically calculated $\text{KL}_{\text{SLB}}$ and $\text{KL}_{\text{MLB}}$ from the SLB and MLB approximations 666For the MLB approximation, we only take those simulation rounds that always converge into consideration., respectively. The smaller the KL divergence is, the more accurate the estimation is.

For Bayesian BMM, the comparisons are presented in Table II and Figure 4. The simulations were run $20$ rounds and the mean values are reported. The comparisons of the Bayesian DMM via SLB [37] and MLB [31] approximations are illustrated in Table LABEL:Tab:_DMM and Figure LABEL:Fig:_DMM. It can be observed that, for both Bayesian BMM and Bayesian DMM, the SLB approximation yields higher objective function value than the MLB approximation. Meanwhile, the KL divergences obtained by the SLB approximation are all smaller than those obtained by the MLB. These facts demonstrate that the SLB approximation is superior to the MLB approximation.

5 Conclusions

The extended variational inference (EVI) framework can be applied in efficiently estimation of non-Gaussian statistical models. We discussed and summarized the required conditions for selection of the auxiliary functions in the EVI framework. Moreover, we also analyzed and compared the single lower-bound (SLB) approximation and the multiple lower-bounds (MLB) approximation. Theoretical analysis showed that the weak condition, in general, can incur smaller systematic gap than the strong condition. Hence, the weak condition is more preferable in practice. Furthermore, quantitative evaluations based on Bayesian beta mixture model and Bayesian Dirichlet mixture model demonstrated that the SLB approximation can theoretically guarantee convergence and is superior to the MLB approximation.