A Deterministic Global Optimization Method for Variational Inference

TL;DR

This paper introduces a deterministic optimization method for variational inference that guarantees convergence to the global variational lower bound, addressing the common issue of local optima in standard algorithms.

Contribution

The authors develop a novel iterative algorithm that ensures convergence to the global optimum in variational inference problems, demonstrated on Bayesian Gaussian mixture models.

Findings

Standard variational inference often converges to local optima.

The proposed method always converges to the global variational lower bound.

Choosing the right variational approximation distribution is crucial.

Abstract

Variational inference methods for latent variable statistical models have gained popularity because they are relatively fast, can handle large data sets, and have deterministic convergence guarantees. However, in practice it is unclear whether the fixed point identified by the variational inference algorithm is a local or a global optimum. Here, we propose a method for constructing iterative optimization algorithms for variational inference problems that are guaranteed to converge to the $\epsilon$-global variational lower bound on the log-likelihood. We derive inference algorithms for two variational approximations to a standard Bayesian Gaussian mixture model (BGMM). We present a minimal data set for empirically testing convergence and show that a variational inference algorithm frequently converges to a local optimum while our algorithm always converges to the globally optimal variational lower bound. We characterize the loss incurred by choosing a non-optimal variational approximation distribution suggesting that selection of the approximating variational distribution deserves as much attention as the selection of the original statistical model for a given data set.

Figures

Tables

Table 1. Table 1 : Timing Comparison of GOP, Variational EM, and BARON Algorithms.

Algorithm		Time (s)
GOP
	$ϵ=1.0$	6.93
	$ϵ=0.1$	9.14
	$ϵ=0.01$	10.77
BARON
	$ϵ=1.0$	35
	$ϵ=0.1$	38
	$ϵ=0.01$	49
Variational EM
	$ϵ=1.0$	0.0053
	$ϵ=0.1$	0.0071
	$ϵ=0.01$	0.0073

Equations

$$l_{\mathcal{D}}(\varphi )=\sum _{i=1}^N\log {\int }_{z_i}p(x,z_i\mid \varphi )d{z}_i.$$

$$l_{\mathcal{D}}(\varphi )$$

$$=\sum _{i=1}^N\log {\int }_{z_i}q(z_i\mid \xi )\frac{p(x_i,z_i\mid \varphi )}{q(z_i\mid \xi )}d{z}_i$$

$$\ge \sum _{i=1}^N{E}_q\left[\log p(x_i,z_i|\varphi )\right]-\sum _{i=1}^N{E}_q\left[\log q(z_i\mid \xi )\right]$$

$$\triangleq \mathcal{L}(\xi ,\varphi ).$$

$$\begin{array}{rl}\underset{\alpha ,\beta }{\min }& f(\alpha ,\beta )\\ \text{s.t.}& g(\alpha ,\beta )\le 0\\ & h(\alpha ,\beta )=0\\ & \alpha \in A,\beta \in B,\end{array}$$

$$\begin{array}{rl}\underset{\alpha }{\min }& f(\alpha ,{\beta }^t)\\ \text{s.t.}& g(\alpha ,{\beta }^t)\le 0\\ & h(\alpha ,{\beta }^t)=0\\ & \alpha \in A.\end{array}$$

$$\begin{array}{rl}\underset{\beta }{\min }\underset{\alpha }{\min }& f(\alpha ,\beta )\\ \text{s.t.}& g(\alpha ,\beta )\le 0\\ & h(\alpha ,\beta )=0\\ & \alpha \in A,\beta \in B,\end{array}$$

$$\underset{\beta }{\min }$$

$$v(\beta )=$$

$$\left\{\begin{array}[]{ll}\min\limits_{\alpha\in A}&f(\alpha,\beta^{t})\\ \text{s.t.}&g(\alpha,\beta^{t})\leq 0\\ &h(\alpha,\beta^{t})=0\end{array}\right\}\equiv\sup_{\mu\geq 0,\lambda}\inf_{\alpha\in A}\left\{f(\alpha,\beta^{t})+\mu^{\top}g(\alpha,\beta^{t})+\lambda^{\top}h(\alpha,\beta^{t})\right\}\quad\forall\ \beta^{t}\in B\cap V.$$

$$v(\beta )=\underset{\mu \ge 0,\lambda }{\sup }\underset{\alpha \in A}{\inf }\left\{f(\alpha ,\beta )+{\mu }^{\top }g(\alpha ,\beta )+{\lambda }^{\top }h(\alpha ,\beta )\right\}\forall \beta \in B\cap V.$$

$$v(\beta )\ge \underset{\begin{array}{c}\alpha \in A\\ \mu \ge 0,\lambda \end{array}}{\inf }\left\{f(\alpha ,\beta )+{\mu }^{\top }g(\alpha ,\beta )+{\lambda }^{\top }h(\alpha ,\beta )\right\}\forall \beta \in B\cap V.$$

$$\underset{\beta }{\min }$$

$$v(\beta )\ge$$

$$\begin{array}{rl}\underset{\beta \in B,v_B}{\min }& v_B\\ \text{s.t.}& v_B\ge \underset{\alpha \in A}{\min }L(\alpha ,\beta ,\lambda ,\mu )\mu \ge 0,\lambda ,\end{array}$$

$$L(\alpha ,\beta ,\lambda ,\mu )=f(\alpha ,\beta )+{\mu }^{\top }g(\alpha ,\beta )+{\lambda }^{\top }h(\alpha ,\beta ).$$

$$L^{\ast }({\alpha }^{\ast },\beta ,\lambda ,\mu )=\underset{\alpha \in A}{\min }L(\alpha ,\beta ,\lambda ,\mu ),$$

$$\begin{array}{rl}\underset{\alpha }{\min }& f(\alpha ,{\beta }^t)\\ \text{s.t.}& g(\alpha ,{\beta }^t)\le 0\\ & h(\alpha ,{\beta }^t)=0\\ & \alpha \in A,\end{array}$$

$$\begin{array}{rl}\underset{\beta \in B,v_B}{\min }& v_B\\ \text{s.t.}& v_B\ge \underset{\alpha \in A}{\min }L(\alpha ,\beta ,{\lambda }^t,{\mu }^t),\end{array}$$

$$L(\alpha ,\beta ,{\lambda }^t,{\mu }^t)=f(\alpha ,\beta )+{\mu }^{t^{\top }}g(\alpha ,\beta )+{\lambda }^{t^{\top }}h(\alpha ,\beta ).$$

$$L^{\ast }({\alpha }^{\ast },\beta ,{\lambda }^t,{\mu }^t)=\underset{\alpha \in A}{\min }L(\alpha ,\beta ,{\lambda }^t,{\mu }^t).$$

$$L^{*}(\alpha^{*},\beta,\lambda^{t},\mu^{t})\geq\min_{\alpha\in A}L(\alpha,\beta,\lambda^{t},\mu^{t})\big{|}^{\text{lin}}_{\alpha^{t}}$$

$$L(\alpha,\beta,\lambda^{t},\mu^{t})\big{|}^{\text{lin}}_{\alpha^{t}}\triangleq L(\alpha^{t},\beta,\lambda^{t},\mu^{t})+\left[\nabla_{\alpha}L(\alpha,\beta,\lambda^{t},\mu^{t})\big{|}_{\alpha^{t}}\right]^{T}(\alpha-\alpha^{t}).$$

$$g_{j}^{t}(\beta)\triangleq\frac{\partial L(\alpha,\beta,\lambda^{t},\mu^{t})}{\partial\alpha_{j}}\bigg{|}_{\alpha^{t}},$$

$$\begin{split}\min_{\alpha\in A}L(\alpha,\beta,\lambda^{t},\mu^{t})\big{|}^{\text{lin}}_{\alpha^{t}}&=\min_{\alpha\in A}\left[L(\alpha^{t},\beta,\lambda^{t},\mu^{t})+\sum_{j}g_{j}^{t}(\beta)(\alpha_{j}-\alpha_{j}^{t})\right]\\ &=L(\alpha^{t},\beta,\lambda^{t},\mu^{t})+\min_{\alpha_{j}\in A}\left[\sum_{j}g_{j}^{t}(\beta)(\alpha_{j}-\alpha_{j}^{t})\right].\end{split}$$

$$\min_{\alpha\in A}L(\alpha,\beta^{d},\lambda^{t},\mu^{t})\big{|}^{\text{lin}}_{\alpha^{t}}=L(\alpha^{t},\beta^{d},\lambda^{t},\mu^{t})+\sum_{j}\min_{\alpha_{j}\in A}g_{j}^{t}(\beta^{d})(\alpha_{j}-\alpha_{j}^{t}).$$

$$\{\begin{array}{ll}\text{if}{g}_j^t(\beta )\ge 0& \text{then}\underset{{\alpha }_j}{\min }{g}_j^t(\beta )({\alpha }_j-{\alpha }_j^t)\ge g_j^t(\beta )({\alpha }_j^L-{\alpha }_j^t)\\ \text{if}{g}_j^t(\beta )\le 0& \text{then}\underset{{\alpha }_j}{\min }{g}_j^t(\beta )({\alpha }_j-{\alpha }_j^t)\ge g_j^t(\beta )({\alpha }_j^U-{\alpha }_j^t).\end{array}$$

$$\underset{{\alpha }_j}{\min }{g}_j^t(\beta )({\alpha }_j-{\alpha }_j^t)\ge g_j^t(\beta )({\alpha }_j^B-{\alpha }_j^t),$$

$${\alpha }_j^B=\{\begin{array}{ll}{\alpha }_j^L& \forall \beta :g_j^t(\beta )\ge 0\\ {\alpha }_j^U& \forall \beta :g_j^t(\beta )\le 0.\end{array}$$

$$\begin{split}\min_{\alpha}L(\alpha,\beta^{d},\lambda^{t},\mu^{t})&\geq L(\alpha^{t},\beta^{d},\lambda^{t},\mu^{t})+\sum_{j}g_{j}^{t}(\beta^{d})(\alpha_{j}^{B^{*}}-\alpha_{j}^{t})\\ &\geq L(\alpha^{B^{*}},\beta^{d},\lambda^{t},\mu^{t})\big{|}^{\text{lin}}_{\alpha^{t}}.\end{split}$$

$$g_{j}^{t}(\beta^{t})=\nabla_{\alpha_{j}}L(\alpha,\beta^{t},\lambda^{t},\mu^{t})\big{|}_{\alpha^{t}}=0,\quad\forall j.$$

$$I_c^t\triangleq \left\{j:g_j^t(\beta )\text{is a function of}\beta \right\}.$$

Full text

1 Introduction

Maximum likelihood estimation of latent (hidden) variable models is computationally challenging because one must integrate over the latent variables to compute the likelihood. Often, the integral is high-dimensional and computationally intractable so one must resort to approximation methods such as variational expectation-maximization (Beal and Ghahramani, 2003).

The variational expectation-maximization algorithm alternates between maximizing an evidence lower bound (ELBO) on the log-likelihood with respect to the model parameters and maximizing the lower bound with respect to the variational distribution parameters. Variational expectation-maximization is popular because it is computationally efficient and performs deterministic coordinate ascent on the ELBO surface. However, variational inference has some statistical and computational limitations. First, the ELBO is often multimodal, so the deterministic algorithm may converge to a local rather than a global optimum depending on the initial parameter estimates. Second, the variational distribution is often chosen for computational convenience rather than for accuracy with respect to the posterior distribution, so the lower bound may be far from tight. The magnitude of the cost of using a poor approximation is typically not known. Here, we develop a deterministic global optimization algorithm for variational expectation-maximization to address these the first issue and quantify the optimal ELBO for an example model to address the second issue.

Related work

It is well-known that expectation-maximization algorithms can converge to a local optimum rather than a global optimum due to local modes or saddle-points in evidence lower bound on the log-likelihood function. In a discussion for the Dempster et al. (1977) paper on expectation-maximization, Murray (1977) described a situation he frequently observed where the EM algorithm converged to a fixed point, but not the global optimum of the log-likelihood. Murray provided a simple data set that reproduced the convergence problems and advocated a solution that continues to be the common practice today – restart the algorithm at multiple initial points and hope that it finds a global optimum. Later,  Wu (1983) described regularity conditions that ensure the EM algorithm converges to a global optimum. However, these regularity conditions are difficult to check in practice. Typically, the practical solution remains the one advocated by Murray – restart the algorithm at many initial values and retain the parameter estimates that give the maximum across restarts.

Contributions

In this article, we present a deterministic algorithm for variational inference that converges to the $\epsilon$-global optimum of the evidence lower bound on the log-likelihood. Global optimization methods such as the one we present here are important to study for several reasons. First, in applications where collecting data is expensive or time-consuming (e.g. medical diagnostics and robotic planning in complex environments) we often have large, but not massive, data sets and we need an inference algorithm that maximizes the use of the information present in the sample while also being computationally tractable. Second, global optimization methods provide new avenues for developing computationally-efficient inference algorithms that have not been extensively studied. Lastly, a guaranteed global optimum parameter estimate allows us to benchmark approximate variational inference methods and approach the analysis of the statistical properties of those methods.

The remainder of this article is organized as follows. In this section, we present the general variational inference problem statement for hierarchical exponential family latent variable models. We summarize the properties and algorithmic steps for the GOP algorithm in Section 2. We derive the GOP inference algorithm for two approximations of the Bayesian Gaussian mixture model in Section 3 and Section 4. In those sections, we describe experimental results comparing the GOP algorithm to variational-EM and BARON. In the latter section, we compare the quality of the optimum for two variational models.

Problem statement

Given a statistical model, the observed-data log-likelihood function, or simply the log-likelihood, is $l_{\mathcal{D}}(\phi)=\log p(\mathcal{D}\mid\phi)$, where $\mathcal{D}$ is the observed data set and $\phi$ is the set of model parameters. The maximum likelihood parameter estimate is $\hat{\phi}\triangleq\arg\max_{\phi}l_{\mathcal{D}}(\phi)$. If the model, $p(\mathcal{D}\mid\phi)$, is in the exponential family and all of the random variables are observed, the maximum likelihood optimization problem is convex and there is a unique maximum. However, if the model has unobserved random variables, we must integrate over those random variables and the log-likelihood may be multi-modal. When the integral is computationally intractable, we must resort to approximation methods to compute the maximum likelihood parameter estimate.

In a latent variable exponential family model, we introduce the unobserved or latent variable $z_{i}$ with the conditional distribution, $p(x_{i}\mid z_{i},\eta)$, and marginal distribution, $p(z_{i}|\theta)$. Our observed data set is i.i.d. $\mathcal{D}=\{(x_{i}),i=1,\ldots,N\}$. The complete data set is $\mathcal{D}_{C}=\{(x_{i},z_{i}),i=1,\ldots,N\}$. Taking $\phi\triangleq\{\eta,\theta\}$, the log-likelihood of the model is

[TABLE]

Unfortunately, the integral over $z_{i}$ in the log-likelihood function is often computationally intractable. By introducing an averaging distribution, $q(z\mid\xi)$, we can form a lower bound on the observed-data log-likelihood using Jensen’s inequality,

[TABLE]

Then, we can obtain an approximate maximum likelihood estimate by maximizing this evidence lower bound (ELBO) over both the averaging distribution and the model parameters, $\max_{\phi,\xi}\mathcal{L}(\xi,\phi)$. The first term in the ELBO is the expected complete-data log-likelihood and the second term is the entropy of the variational distribution.

In the expected complete log-likelihood, there is a coupling between the variational distribution parameters, $\xi$, and the model parameters, $\phi$. This coupling makes optimizing the ELBO difficult in practice because the resulting optimization problem is nonconvex.

2 GOP Variational Inference

The general Global OPtimization (GOP) algorithm was first developed by Floudas and Visweswaran (1990). They extended the decomposition ideas of Benders (1962) and Geoffrion (1972) in the context of solving nonconvex optimization problems to global optimality (Floudas and Visweswaran, 1993; Floudas, 2000). While the framework has been successfully used for problems in process design, control, and computational chemistry, it has not, to our knowledge, been applied to variational statistical inference.

[Add reference to linear regression here.]

Here, we extend the GOP algorithm for variational inference in hierarchical exponential family models. In Section 2.1, we state the problem conditions necessary for the GOP algorithm. Then, we briefly review the mathematical theory for the GOP algorithm in Section 2.2. Finally, in Section 2.3, we outline the general GOP algorithmic steps. The review of the GOP algorithm in this section is not novel and necessarily brief. A more complete presentation of the general algorithm can be found in Floudas (2000).

2.1 Problem Statement

GOP addresses biconvex optimization problems that can be formulated as

[TABLE]

where $A$ and $B$ are convex compact sets, and $g(\alpha,\beta)$ and $h(\alpha,\beta)$ are vectors of inequality and equality constraints respectively. We require the following conditions on (2):

$f(\alpha,\beta)$ is convex in $\alpha$ for every fixed $\beta$, and convex in $\beta$ for every fixed $\alpha$; 2. 2.

$g(\alpha,\beta)$ is convex in $\alpha$ for every fixed $\beta$, and convex in $\beta$ for every fixed $\alpha$; 3. 3.

$h(\alpha,\beta)$ is affine in $\alpha$ for every fixed $\beta$, and affine in $\beta$ for every fixed $\alpha$; 4. 4.

first-order constraints qualifications (e.g. Slater’s conditions) are satisfied for every fixed $\beta$.

If these conditions are satisfied, we have a biconvex optimization problem and we can use the GOP algorithm to find the $\epsilon$-global optimum.

Partitioning and transformation of decision variables

Variational inference is typically formulated as an alternating coordinate ascent algorithm where the evidence lower bound is iteratively maximized with respect to the model parameters and the variational distribution parameters. However, the objective function construed under that partitioning of the decision variables is not necessarily biconvex. Instead, we consider all of the decision variables (model parameters and variational parameters) in the variational inference optimization problem jointly. Then, we are free to partition variables as we choose to ensure the GOP problem conditions are satisfied.

One may be able to transform variables in the original variational optimization problem, adding equality constraints as necessary, to satisfy the GOP problem conditions. For larger problems and more complex models, the partitioning of the decision variables may not be apparent. Hansen et al. (1993) proposed an algorithm to bilinearize quadratic and polynomial function problems, rational polynomials, and problems involving hyperbolic functions.

Biconvexity for exponential family hierarchical latent variable models

The variational EM problem is a biconvex optimization problem under certain conditions. We cast the variational EM problem into a biconvex convex optimization form by partitioning the model parameters, $\phi$, and the variational parameters $\xi$ into $\alpha$ and $\beta$ such that the GOP conditions are satisfied. All of the functions are analytical and differentiable. For two choices of approximating distribution for a Bayesian Gaussian mixture model, we find a suitable partition of the decision variables. However, what conditions guaranetee such a partition for general exponential family hierarchical models remains an open question.

2.2 GOP Theory

Primal problem

Because the original problem (2) satisfies the biconvexity conditions, fixing $\beta$ (say to $\beta^{t}\in B$) results in a convex subproblem,

[TABLE]

called the primal problem. The restriction of $\beta$ to a fixed value is equivalent to solving the original problem with additional constraints. Therefore, the solution of the primal problem provides an upper bound on the original problem (2). Since the problem is convex, it can be solved with any convex optimization solver.

In variational inference, this step would typically be considered the E-step if $\alpha$ was the variational parameters and the M-step if $\alpha$ was the model parameters. However, we have partitioned the decision variables to formulate the problem as a biconvex optimization problem, therefore the direct EM-step interpretation is lost.

Relaxed dual problem

We can write the original optimization problem (2) as

[TABLE]

in order to cast the problem as inner and outer optimization problems using a projection of the problem onto the space of $\beta$ variables:

[TABLE]

For a given $\beta^{t}$, the inner minimization problem in (2.2) is simply the primal problem. Therefore, the function $v(\beta)$ can be defined for a set of solutions of the primal problem (3) for different values of $\beta$. However, since $v(\beta)$ is defined implicitly (2.2) can be as difficult to solve as the original problem (2).

Because the original problem satisfies Slater’s conditions, by nonlinear duality theory, the optimal objective function value of the primal problem (3) for any $\beta=\beta^{t}$ is identical to the solution of its corresponding dual problem ,

[TABLE]

Substituting into the primal problem (3) gives

[TABLE]

Using the definition of the supremum we can form the bound

[TABLE]

Then, we have an equivalent problem to (2.2),

[TABLE]

The last two conditions in the dual problem (2.2) define an implicit set of constraints making the solution of the problem difficult.

Dropping the last two conditions gives a relaxed dual problem,

[TABLE]

where $v_{B}\in\mathbb{R}$ and the Lagrange function for the primal problem (3) is

[TABLE]

We call the minimization problem in (6) the inner relaxed dual problem,

[TABLE]

where $\alpha^{*}$ denotes the value of the $\alpha$ variables at the optimal solution of the inner relaxed dual problem.

In summary, the primal problem (3) contains more constraints than the original problem because $\beta$ is fixed and so provides an upper bound on the original problem. The relaxed dual problem (6) contains fewer constraints than the original problem and so provides a lower bound on the original problem. Because the primal and relaxed dual problems are formed from an inner-outer optimization problem, they can be used to form an iterative alternating algorithm to determine the global solution of the original problem. Making the iteration $t$ explicit, the primal problem is

[TABLE]

and the relaxed dual problem is

[TABLE]

where

[TABLE]

The key here is that information is passed from the primal problem to the dual problem through the Lagrange multipliers, $\{\lambda^{t},\mu^{t}\}$, and information is passed from the dual problem to the primal problem through the optimal dual variables, $\beta^{t}$.

Relaxed dual problem decomposition

The form of the relaxed dual problem (6) is difficult to solve because it contains an inner relaxed dual problem that is parametric in $\beta$. Here, we decompose the relaxed dual problem into a set of independent relaxed dual subproblems that are combined to find the solution of the relaxed dual problem. The algorithm iterates between solving the primal problem and the relaxed dual problem in a way that provides guaranteed convergence to the global optimum.

Recall, the inner relaxed dual problem (8) is

[TABLE]

Because the inner relaxed dual problem is convex in $\alpha$, we can form a lower bound by taking the first-order Taylor series approximation about some fixed $\alpha^{t}$,

[TABLE]

where

[TABLE]

For notational convenience we define

[TABLE]

the $j$-th component of the gradient of the original problem Lagrange function with respect to $\alpha$ evaluated at $\alpha^{t}$. Note that we are using $g$ to denote the inequality constraints in the original problem; this different function has a superscript $t$ that indicates the iteration. Then, we can write the lower bound on the inner relaxed dual problem as

[TABLE]

For any fixed $\beta=\beta^{d}$ the summation and minimization can be exchanged,

[TABLE]

Because the $j$-th component of the second term on the right hand side is linear in $\alpha_{j}$, the minimum will be at a bound of $\alpha_{j}\in A$.

The specific nature of the bound (lower or upper) is determined by the sign of g_{j}^{t}(\beta)=\nabla_{\alpha_{j}}L(\alpha,\beta,\lambda^{t},\mu^{t})\big{|}_{\alpha^{t}}. There are two cases:

[TABLE]

where $\alpha_{j}^{L}$ and $\alpha_{j}^{U}$ are the lower and upper bounds of $\alpha_{j}\in A$ respectively. We can write the two cases compactly as

[TABLE]

where

[TABLE]

So, for any fixed $\beta^{d}$, there exists a combination of bounds $B^{*}$ for the $\alpha$ variables such that

[TABLE]

The vector $B$ is the same size as $\alpha$ and each element $j$ indicates whether $\alpha_{j}$ is set at its upper or lower bound. The variable $\alpha_{j}^{B^{*}}$ is the value of $\alpha_{j}$ at the upper/lower bound such that (16) is valid. For any discrete $\beta^{d}$, we can obtain a lower bound on the inner relaxed dual problem $L^{*}(\alpha^{*},\beta^{d},\lambda^{t},\mu^{t})$ by taking the minimum of the linearized Lagrange function over all combinations of bounds $B\in CB$. Since this is true for every $\beta^{d}$, it must be true for all $\beta$.

Computing the linearized Lagrange over all combinations of bounds requires that we compute $2^{|\alpha|}$ linearized Lagrange functions, which can be computationally expensive. The Karush-Kuhn-Tucker (KKT) optimality conditions for the primal problem require that

[TABLE]

So, if the partial derivative of the primal Lagrange function with respect to $\alpha_{j}$ is found not to be a function of $\beta$, we do not need to consider the bounds of $\alpha_{j}$. We reduce the computational burden by a factor of $2$ for every such variable.

Every $\alpha_{j}$ for which $g_{j}^{t}(\beta)$ is a function of $\beta$ is called a connected variable. The set of all such connected variables is

[TABLE]

Therefore, the solution of the inner relaxed dual problem (8) with Lagrange function replaced by its linearization about $\alpha^{t}$ depends only on the connected $\alpha$ variables, and we can write the inner relaxed dual minimization problem as

[TABLE]

Now, we can solve for a lower bound of the inner relaxed dual problem by setting the connected variables at their bounds and minimizing the linearized primal Lagrange function. Setting the connected variables at their bounds is equivalent to selecting a half-space defined by the cut (hyper-plane) $g_{j}^{t}(\beta)=0$ in the domain of $\beta$. Therefore, we are solving for a lower bound of the linearized Lagrange function in a unique region in the domain of $\beta$ for each $B\in CB$. Since we solve for such a lower bound for all combinations of the connected variables, we cover the entire domain of $\beta$.

Since $g_{i}^{t}(\beta)$ can enter as either $\leq 0$ or $\geq 0$, we require that $g_{j}^{t}(\beta)$ be linear in $\beta$ to ensure all of the regions defined by $B\in CB$ are convex. However, in general $g_{j}(\beta)$ is convex. When $g_{j}(\beta)$ is convex, but not linear, we can linearize the Lagrange function with respect to $\beta$ by taking the first-order Taylor series approximation with respect to $\beta$ at $\beta^{t}$. This linearization provides a lower bound on the Lagrange function and the resulting $g_{j}^{t}(\beta)$ cuts are linear in $\beta$.

Relaxed dual problem as a sequence of convex problems

We have shown that we can approximate the inner relaxed dual problem by a minimization over a finite combination of settings of $\alpha$ at its bounds. Here, we show that we can solve the relaxed dual problem as a sequence of convex subproblems. Finally, we show we achieve $\epsilon$-global convergence by iteratively alternating between solving the primal problem and one step in the sequence of relaxed dual subproblems. These optimality results are based on proofs by Geoffrion (1972).

Suppose we have a sequence of relaxed dual subproblems such that $v_{B}^{*^{T}}$ is the optimal value of the $T$-th relaxed dual problem.

[TABLE]

Recall that for any $T$, the relaxed dual problem contains fewer constraints than the original dual problem and is a valid under-estimator of the original problem. Further, since the relaxed dual problem at iteration $T$ contains more constraints than the relaxed dual problem at iteration $t<T$, the under-estimator is nondecreasing in $T$. Now, we use the approach for solving the inner relaxed dual problem by setting $\alpha$ at its bounds to derive an iterative algorithm for (19).

Iteration 1

For $t=1$ we have

[TABLE]

As this holds for all $\beta$, it holds for the minimum over $\beta$,

[TABLE]

Because the right hand side is only a function of $\beta$, the minimization operators can be interchanged,

[TABLE]

or equivalently

[TABLE]

So, for iteration 1, we have an optimization problem that provides a valid lower bound on the inner relaxed dual problem. For each setting of $\alpha$ at its bounds $B\in CB$, we solve a convex optimization problem – the relaxed dual subproblem.

Iteration T

Recall that $g_{j}^{t}(\beta)$ defines a cut in the domain of $\beta$ and all $j\in I_{C}^{t}$ cuts intersect at $\beta^{t}$ such that any point $\beta^{d}\in B$ lies in one and only one region defined by the cuts (with the exception of boundary points). Further, recall that there is a one-to-one correspondence between $\alpha^{B}$, the primal problem variables set at their boundaries, and the region defined by the cuts, $g^{t}(\beta)$. Now, let $UL(t,T)$ be the set of Lagrange functions and associated qualifying constraints from the $t$-th iteration whose qualifying constraints are satisfied at $\beta^{T}$, the current value of the relaxed dual problem decision variables. The set of active qualifying constraints uniquely identifies $\alpha^{B_{t}}$, the primal problem variable bounds, where we denote the combination of bounds that uniquely identify the active qualifying constraints from iteration $t$ as $B_{t}$.

We can write the inner relaxed dual problems from (19) for iteration $t<T$ as

[TABLE]

where $B_{t}$ is identified by $UL(t,T)$. Since the $\alpha$ variables are set for all iterations $t<T$, the $\min_{B\in CB}$ operator only applies to the $T$-th iteration and we have

[TABLE]

Therefore, the relaxed dual problem can be solved as a series of convex problems. Each iteration adds constraints to the relaxed dual problem ensuring that the lower bound is non-decreasing. Since the relaxed dual problem is a valid lower bound on the global optimum for every iteration, it is a valid lower bound for the $T$-th iteration. The proof that the upper bound provided by the primal and the lower bound provided by the relaxed dual converge can be found in (Floudas, 2000)[Theorem 3.6.1].

2.3 GOP Algorithm Steps

Step 0 – Initialize parameters.

Initialize the upper bound $P^{UBD}=\infty$ and the lower bound $R^{LBD}=-\infty$. Define lower and upper bounds for the primal decision variables $\alpha_{j}^{L}$ and $\alpha_{j}^{U}$ respectively. Select a feasible initial value of the dual decision variables $\beta^{1}$. Set the iteration counter to $T=1$. Finally, select a convergence tolerance parameter $\epsilon$.

Step 1 – Solve primal problem

Solve the primal problem with the relaxed dual variables fixed at $\beta^{T}$ and store the optimal Lagrange multipliers $\lambda^{T}$ and $\mu^{T}$. Update the upper bound $P^{UBD}=\min(P^{UBD},P^{T}(\beta^{T}))$, where $P^{T}(\beta^{T})$ is the solution of the primal problem.

Step 2 – Select Lagrange functions from previous iterations

While $T>1$, for $t=1,\ldots,T-1$ identify the settings of the primal problem variables at their bounds $B_{t}$ such that all of the qualifying constraints from iteration $t$ are active at the current value of the relaxed dual variables $\beta^{T}$. Select those qualifying constraints and their corresponding Lagrange function to be in the set of constraints $UL(t,T)$ of the current relaxed dual problem.

Step 3 – Solve relaxed dual problem

Determine the set of connected variables at iteration $T$, $I_{C}^{T}$. For each setting of the connected variables at their bounds $\alpha^{B}$ where $B\in CB$, solve the relaxed dual problem:

[TABLE]

where $t=1,\ldots,T-1$.

If feasible, store the solution in $v_{B}^{\text{stor}}(T,B)$ and $\beta^{\text{stor}}(T,B)$. If not feasible, then store a null value.

Step 4 – Update lower bound and $\beta^{(T+1)}$

Select the minimum from the entire set $v_{B}^{\text{stor}}$, $v_{B}^{\text{min}}$. Note that $v_{B}^{\text{stor}}$ contains the $2^{|I_{C}^{T}|}$ solutions of the relaxed dual problem from the current iteration and the solutions from all of the previous iterations. Update the value of the lower bound $R^{LBD}\leftarrow v_{B}^{\text{min}}$. Update the value of $\beta$ with the corresponding optimal relaxed dual decision variable $\beta^{(T+1)}\leftarrow\beta^{\text{min}}$. Remove $v_{B}^{\text{min}}$ and $\beta^{\text{min}}$ from the stored sets so that they are not selected again in a subsequent iteration.

Step 5 – Check for convergence

Check if $R^{LBD}>P^{UBD}-\epsilon$. If the convergence criterion is satisfied, return the optimal parameter values, else increment the iteration counter $T\leftarrow T+1$ and go to step 1.

3 BGMM – Point Mass Approximation

Here, we consider a variant of the Gaussian mixture model where cluster mean is endowed with a Gaussian prior. We consider the variational model for the posterior distribution to be a point mass function for discrete latent random variables and a Dirac delta function for continuous random variables. Variational expectation-maximization with a point-mass approximating distribution is exactly classical expectation-maximization (Gelman et al., 2013, p 337).

3.1 GOP Algorithm Derivation

Model structure

A Gaussian mixture model with a prior distribution on the cluster means is

[TABLE]

with model parameters $\phi\triangleq\{\pi,\Gamma\}=\{\pi_{1},\ldots,\pi_{K},\Gamma\}$, unobserved/latent variables $\{Z,M_{1},\dots,M_{K}\}$, and observed variable $Y=y$.

Our inferential aim is the joint posterior distribution $Z_{i},M\mid Y_{i},\hat{\phi}$, for $i=1,\ldots,N$, where $\hat{\phi}$ is the maximum likelihood parameter estimate. The log likelihood for the model is $l_{\mathcal{D}}(\phi)\triangleq\log f(y|\phi)$.

Variational distribution

We take the variational distribution to be the following fully factorized (mean-field) distribution

[TABLE]

where

[TABLE]

We define $\xi\triangleq\{\tau,\nu\}$ be the set of variational distribution parameters.

Variational lower bound for the log-likelihood

Given the variational distribution (23), the evidence lower bound (ELBO) is found by applying Jensen’s inequality. The final ELBO is construed as a function of the model parameters, $\phi=\{\pi,\eta_{m}\}$, and the variational parameters $\xi=\{\tau,\nu\}$, where we have replaced $\Gamma$ by its natural parameter $\eta_{m}=-1/(2\Gamma)$ to make the ELBO linear in the parameter. The ELBO can be expanded to (see Appendix B.2 for the full derivation)

[TABLE]

Partitioning and transforming decision variables

The variational inference problem in standard form is a minimization problem, where the objective function is the negative ELBO.

The objective function is not biconvex in the separation of variables into model parameters, $\phi=\{\pi,\eta_{m}\}$, and variational parameters, $\xi=\{\tau,\nu\}$ because of the first term involving $\tau_{ik}(y_{i}-\nu_{k})^{2}$. This term is cubic in the variational parameters and thus nonconvex. However, if we separate $\tau$ and $\nu$, the first term is bi-convex. The second term is biconvex in $\pi$ and $\tau$. The third term is biconvex in $\eta_{m}$ and $\nu$. The remaining two terms are convex in $\eta_{m}$ and $\tau$ respectively. So, we must separate $\tau$ from $\nu$, $\tau$ from $\pi$, and $\eta_{m}$ from $\nu$. A partition that satisfies these conditions is $\alpha\triangleq\{\nu,\pi\}$ and $\beta\triangleq\{\tau,\eta_{m}\}$. A proof that all four GOP algorithm conditions are satisfied is provided in Appendix A.3.

Our selection of which variable to optimize over in the primal problem affects computational scalability. If we choose to optimize over $\beta$ in the primal problem, we must solve at most $2^{|\beta|}=2^{NK+1}$ relaxed dual subproblems for each iteration where each subproblem has $|\alpha|=2K$ decision variables. This is disadvantageous because the number of relaxed dual subproblems grows exponentially in the sample size. Instead, we will optimize over $\alpha$ in the primal problem, which means we only need to solve at most $2^{|\alpha|}=2^{2K}$ relaxed dual subproblems where each subproblem has $|\beta|=NK+1$ decision variables. While each subproblem may be large, it will be a linear program which can be solved efficiently for millions of decision variables.

The variational inference problem is now

[TABLE]

where $\alpha\triangleq\{\nu,\pi\}$ and $\beta\triangleq\{\tau,\eta_{m}\}$.

3.2 Experimental Results

First, we examine the convergence of the upper and lower bounds on the ELBO for a small data set. Then, we compare the accuracy and computational efficiency of GOP and variational inference algorithms. Last, we compare the optimal solution of GOP and variational inference under many random initializations.

It is well-known that the expectation-maximization algorithm is only guaranteed to converge to a fixed-point which may be arbitrarily far from a global optimum. Others have offered anecdotes and evidence that it is actually rather common for expectation-maximization to converge to a local optimum or saddle point rather than a global optimum. Murray (1977) described a simple data set that exhibits this local convergence issue for a Gaussian model. We sought to identify a minimal data set for the (Bayesian) Gaussian mixture model that exhibits the local convergence problem. Archambeau (2003) found that repeated observations and outliers are particularly problematic for expectation-maximization and the Gaussian mixture model. We consider these properties to construct the data set $y=[-10,-10,5,25]$.

We find that this minimal data set is useful for stress-testing any inference algorithm for the (Bayesian) Gassian mixture model.

Convergence experiment

We examined the optimum ELBO values for GOP and variational EM (VEM) by randomly drawing initial model and variational parameters and running each algorithm to convergence. The initial $\pi$ and $\tau_{i}$ are generated from a Dirichlet distribution with $\alpha=[1,1]$. The initial $\Gamma$ is generated from a Gamma with shape parameter equal to the range of $y$. The initial $\nu_{k}$ is drawn from uniform distribution over the range of $y$. We drew 100 random initial values and we observed that VEM converges to a local optimum 87 times out of 100 restarts, while GOP converges to the global optimum for all 100 restarts.

We verified the global optimal value using BARON, an independent global optimization algorithm (Sahinidis, 1996; Tawarmalani and Sahinidis, 2005). These results show that GOP, empirically, always converges to the global optimum. The VEM algorithm failed to converge to the global optimum more often than it succeeded. These results quantify and agree with anecdotal evidence in the literature (Murray, 1977).

Computational time experiment

We explore two aspects of computational efficiency: “What is the behavior of the upper and lower bounds from the GOP algorithm across iterations?” and “How does the timing of the GOP algorithm compare to other local and global inference methods?”.

Figure 1 shows the convergence iterations for the GOP algorithm for two initial seeds. For one of the seeds the VEM algorithm converges to a suboptimal local optimum Figure 1(a) and for the other seed the VEM algorithm converges to a suboptimal local optimum Figure 1(b). The optimality of the solution is verified by BARON. The globally optimal ELBO is $-84.04$ and the locally optimal ELBO is $-108.8$.

In early iterations, the GOP bounds include the local optimum, but at approximately 3 sec, the bound interval no longer covers the local optimum. So, even though the GOP algorithm is slower than other local methods, one way it can be used is as a running oracle on an algorithm with only local guarantees. When the GOP bound interval no longer covers the VEM solution, we can restart VEM until the solution found by VEM is covered by the GOP interval.

We compared the time for convergence to VEM and a standard global optimization algorithm, BARON for our minimal data set. The results are summarized in Table 1.

Our implementation of the GOP algorithm is faster than the global optimization algorithm, BARON and slower than VEM. Indeed, VEM is several orders of magnitude faster. However, the key point of this work is to demonstrate a method to obtain a global optimum. Our implementation of GOP can be substanially enhanced by (1) using the reduction method that shrinks the decision variable search space at each iteration (Visweswaran and Floudas, 1996), (2) using mixed-integer linear programming to identify the optimal relaxed dual sub-problem at each iteration (Visweswaran and Floudas, 1996), and (3) using multi-threading to compute each master problem. The multi-threading option is most appealing because each relaxed dual problem is completely independent, thereby suggesting how to drastically reduce the most expensive computation step.

4 BGMM – Gaussian Approximation

In this section we explore the quality of the approximation provided by two different variational models for the Bayesian Gaussian mixture model. The original model structure is the same as in Section 3, so we start with the variational model structure here.

4.1 GOP Algorithm Derivation

Variational distribution

We take the variational distribution to be the following fully factorized (mean-field) distribution

[TABLE]

where

[TABLE]

We define $\xi\triangleq\{\tau,\nu,\gamma\}$ to be the set of variational distribution parameters.

Variational lower bound for the log-likelihood

Given the variational distribution (26), the evidence lower bound (ELBO) is found by applying Jensen’s inequality. The final ELBO is construed as a function of the model parameters, $\phi=\{\pi,\eta_{m}\}$, and the variational parameters $\xi=\{\tau,\nu,\gamma\}$, where we have replaced $\Gamma$ by its natural parameter $\eta_{m}=-1/(2\Gamma)$ to make the ELBO linear in the parameter.

[TABLE]

The ELBO can be expanded (see Appendix C.1) to

[TABLE]

Partitioning and transforming decision variables

The variational inference problem in standard form is a minimization problem, where the objective function is the negative ELBO. Using the same argument as in Section 3, we find a partition that satisfies the GOP conditions – $\alpha\triangleq\{\nu,\pi,\gamma\}$ and $\beta\triangleq\{\tau,\eta_{m}\}$. A full proof that all four GOP algorithm conditions are satisfied is provided in the Appendix C.3.

The variational inference problem is now

[TABLE]

where $\alpha\triangleq\{\pi,\nu,\gamma\}$ and $\beta\triangleq\{\tau,\eta_{m}\}$.

4.2 Experimental Results

In variational inference, the choice of approximating distribution is typically motivated by computational convenience rather than accuracy with respect to the posterior distribution. In addition, the loss incurred by choosing a particular approximating distribution is generally not known. The task of characterizing this loss is daunting because local methods such as VEM do not provide a performance guarantee.

Quality of optimal ELBO

We investigate the quality of the optimal ELBO provided by (1) the point mass approximation detailed in Section 3 and (2) the Gaussian approximation presented earlier in this section. We run GOP, adapted to both approximations, using the minimal data set described in 3. The globally optimal ELBO for point mass approximation and for Gaussian approximation is -84.04 and -82.75, respectively. We observe that the globally optimal ELBO for Gaussian approximation is higher than that of point mass approximation meaning that the Gaussian variational distribution is a better approximating distribution for our model compared to the point mass variational distribution. This result shows that we do indeed incur a loss in accuracy when selecting an inferior approximating distribution suggesting that careful consideration must be made when choosing an approximating distribution.

5 Conclusion

There are several contributions of this article. We describe an adaptation of a general-purpose biconvex global optimization method (GOP) for variational expectation-maximization problems. We derive GOP inference algorithms for two variational approximations of the Bayesian Gaussian mixture model. We provide a simple data set for which the variational evidence lower bound is multimodel that is useful for testing approximation methods. We compare our GOP inference algorithm to one locally optimal algorithm (variational-EM) and one state-of-the-art global optimization algorithm (BARON) in terms of accuracy and computational time and show that GOP always converges to the global optimum and is not prohibitive in terms of computational time for the small data set tested. We compare the globally optimal variational evidence lower bounds for two variational approximate distributions and show that the global variational evidence lower bound is better for a more flexible variational model. The focus in this article is on the properties of the global optimization algorithm. We describe several future improvements that can significantly improve the computational efficiency.

Appendix A Appendix for GMM

A.1 Derivation of Log-likelihood Lower Bound

The variational lower bound is

[TABLE]

Writing out each term of the bound gives

[TABLE]

The ELBO requires expected values under the variational distribution for $E_{q}[z_{ik}]$ and $E_{q}[z_{ik}^{2}]$. We can simplify these expected values:

[TABLE]

Plugging these expectations back into the full ELBO gives us the result

[TABLE]

A.2 Variational Inference Algorithm

The standard variational EM algorithm for the model is an alternating coordinate ascent algorithm on the ELBO. The M-step maximizes the ELBO with respect to the model parameters, $\phi$. The E-step maximizes the ELBO with respect to the variational distribution parameters, $\xi$. Since the only variational quantity of interest is the expected value of the latent variables,$E_{q}[z_{ik}]$, variational inference is equivalent to the EM algorithm for this model. Recall the ELBO is

[TABLE]

and the ELBO Lagrange function is

[TABLE]

M-step

We can maximize the ELBO with respect to $\pi$ by setting the derivative of the ELBO with respect to $\pi$ and solving for $\pi$. Solving for the maximum of the resulting Lagrange function gives the update formula,

[TABLE]

We see that the MLE for $\pi$ is the usual sample mean of $\tau$.

Setting the derivative of the ELBO Lagrange function with respect to $\mu_{k}$ and solving for $\mu_{k}$ gives the update formula,

[TABLE]

E-step

Setting the derivative of the ELBO Lagrange function with respect to $\tau_{ik}$ and solving for $\tau_{ik}$ gives the update formula,

[TABLE]

where we have solved for $\lambda_{qi}$ by plugging $\hat{\tau}_{ik}$ into its constraint.

Algorithm

Putting the two update equations for the M-step together with the two update equations for the E-step gives the variational inference algorithm.

A.3 Proof that Gaussian Mixture Model ELBO satisfies GOP conditions

The objective function of interest is

[TABLE]

where $\phi=\{\pi,\mu\}$ and $\xi=\{\tau\}$.

The Hessian with respect to $\phi$ for fixed $\hat{\xi}$ is

[TABLE]

where

[TABLE]

and

[TABLE]

The matrix is diagonal and all of the diagonal elements are non-negative, so the Hessian is positive semidefinite and $f(\phi,\hat{\xi})$ is convex.

The Hessian with respect to $\xi$ for fixed $\hat{\phi}$ is

[TABLE]

where $\tilde{\tau}$ is the matrix $\tau$ construed as a vector taken column-wise. Again, the matrix is diagonal and all of the diagonal elements are non-negative, so the Hessian is positive semidefinite and $f(\hat{\phi},\xi)$ is convex.

Therefore, the first condition of the GOP algorithm is satisfied. The second condition is satisfied because we have no inequality constraints that are functions of both $\phi$ and $\xi$ – all constraints can be absorbed into the sets $A,B$. The third condition is satisfied because the equality constraints are each convex combinations and therefore affine. The fourth condition (Slater’s condition) is satisfied because a point exists within the interior of the feasible region. Interestingly, for the Gaussian mixture model, the problem is biconvex in the original separation of decision variables into model parameters and variational parameters.

Appendix B Appendix for BGMM – Point Mass Approximation

B.1 Derivation of Log-Likelihood

The likelihood function for the Bayesian Gaussian mixture model is

[TABLE]

The log-likelihood function is then,

[TABLE]

Since the model has a hierarchical structure, the joint density function can be expanded as a product of conditional distributions,

[TABLE]

The integral with respect to $z_{i}$ can be replaced with a summation since it is finite,

[TABLE]

The conditional density functions are:

[TABLE]

The conditional density functions can be substituted in to the log-likelihood to give

[TABLE]

Now, if $z_{i}=k$ the integrals for $m_{k^{\prime}}$ where $k^{\prime}\neq k$ do not involve the density for $y_{i}$ or $z_{i}$ and are constant. The kernels for those integrals are then $p(m_{k^{\prime}}\mid\Gamma)$ and the integrals are equal to one. Therefore, the log-likelihood can be simplified as

[TABLE]

where

[TABLE]

After completing the square, the integral can be substituted by the normalizing factor for a Gaussian distribution,

[TABLE]

Simplifying and substituting $g(z_{i})$ back into the log-likelihood gives

[TABLE]

The only term left involving $z_{i}$ is $\pi_{z_{i}}$, and $\sum_{z_{i}=1}^{K}\pi_{z_{i}}=1$, so the parameter $pi$ drops out of the log-likelihood,

[TABLE]

Taking the derivative, setting it equal to zero and solving for $\Gamma$ gives the maximum likelihood estimate

[TABLE]

and substituting back into the log-likelihood gives the value of the maximum log-likelihood.

B.2 Derivation of Log-likelihood Lower Bound

The variational lower bound is

[TABLE]

We can write out each term in the ELBO:

[TABLE]

where we have substituted $\eta_{m}\triangleq-(2\Gamma)^{-1}$. The ELBO requires expected values under the variational distribution for $E_{q}\left[\sum_{k=1}^{K}z_{ik}m_{k}\right]$, $E_{q}\left[\left(\sum_{k=1}^{K}z_{ik}m_{k}\right)^{2}\right]$, $E_{q}[z_{ik}]$, and $E_{q}[m_{k}^{2}]$. We can simplify these expected values:

[TABLE]

Substituting these expectations back into the full ELBO gives us the result

[TABLE]

B.3 Variational Inference Algorithm

The standard variational EM algorithm for the model is an alternating coordinate ascent algorithm on the ELBO. The M-step maximizes the ELBO with respect to the model parameters, $\phi$. The E-step maximizes the ELBO with respect to the variational distribution parameters, $\xi$.

M-step

We can maximize the ELBO with respect to $\pi$ by setting the derivative of the ELBO function with respect to $\pi$ and solving for $\pi$. However, we also require that $\sum_{k=1}^{K}\pi_{k}=1$, so we add a Lagrange multiplier to enforce this constraint. Solving for the maximum of the resulting Lagrange function gives the update formula,

[TABLE]

We see that the MLE for $\pi$ is the usual sample mean of $\tau$.

Setting the derivative of the ELBO Lagrange function with respect to $\eta_{m}$ and solving for $\eta_{m}$ gives the update formula,

[TABLE]

This natural parameter estimate gives the usual MLE estimate for the corresponding moment parameter $\hat{\Gamma}=\frac{1}{K}\sum_{k=1}^{K}\nu_{k}^{2}$.

E-step

Setting the derivative of the ELBO Lagrange function with respect to $\tau_{ik}$ and solving for $\tau_{ik}$ gives the update formula,

[TABLE]

where we have solved for $\lambda_{qi}$ by plugging $\hat{\tau}_{ik}$ into its constraint. The update for $\tau_{ik}$ can be rewritten using the moment parameters as

[TABLE]

This form shows that $\hat{\tau}_{ik}$ is proportional to the product of the prior for belonging to the $k$-th cluster and the likelihood of sampling $y_{i}$ from a Gaussian with mean $\nu_{k}$.

Setting the derivative of the ELBO Lagrange function with respect to $\nu_{k}$ and solving for $\nu_{k}$ gives the update formula,

[TABLE]

The update for $\nu_{k}$ can be rewritten using the moment parameters as

[TABLE]

This form shows that $\hat{\nu}_{k}$ is the sample mean of the weighted data assigned to the $k$-th cluster when the prior precision $\Gamma^{-1}=0$. This value of the precision corresponds to an improper prior on the cluster means and reduces the Bayesian Gaussian mixture model to a standard Gaussian mixture model. When $\Gamma^{-1}$ is greater than zero $\hat{\nu}_{k}$ is regularized (shrunken) towards the prior mean, which in this model is zero.

Algorithm

Putting the two update equations for the M-step together with the two update equations for the E-step gives the variational inference algorithm.

B.4 Proof of GOP conditions

The objective function of interest is

[TABLE]

where $\alpha\triangleq\{\nu,\pi\}$ and $\beta\triangleq\{\tau,\eta_{m}\}$.

The Hessian with respect to $\beta$ for fixed $\hat{\alpha}$ is

[TABLE]

where

[TABLE]

and

[TABLE]

We define $\tilde{\tau}$ to be the matrix $\tau$ construed as a vector taken column-wise. Since the matrix is diagonal and all of the diagonal elements are non-negative, the Hessian is positive semidefinite and $f(\hat{\alpha},\beta)$ is convex.

The Hessian with respect to $\alpha$ for fixed $\hat{\beta}$

[TABLE]

where

[TABLE]

and

[TABLE]

Again, the matrix is diagonal and all of the diagonal elements are non-negative, so the Hessian is positive semidefinite and $f(\alpha,\hat{\beta})$ is convex.

Therefore, the first condition of the GOP algorithm is satisfied. The second condition is satisfied because we have no inequality constraints that are functions of both $\alpha$ and $\beta$ – all constraints can be absorbed into the sets $A,B$. The third condition is satisfied because the equality constraints are each convex combinations and therefore affine. The fourth condition (Slater’s condition) is satisfied because the interior of the feasible region has an interior point.

Appendix C Appendix for BGMM – Gaussian Approximation

C.1 Derivation of Log-likelihood Lower Bound

The variational lower bound is

[TABLE]

We can write out each term in the ELBO:

[TABLE]

where we define $\eta_{m}=-(2\Gamma)^{-1}$. The ELBO requires expected values under the variational distribution for $E_{q}\left[\sum_{k=1}^{K}z_{ik}m_{k}\right]$, $E_{q}\left[\left(\sum_{k=1}^{K}z_{ik}m_{k}\right)^{2}\right]$, $E_{q}[z_{ik}]$, and $E_{q}[m_{k}^{2}]$. We can simplify these expected values:

[TABLE]

Substituting these expectations back into the full ELBO gives us the result

[TABLE]

C.2 Variational Inference Algorithm

The standard variational EM algorithm for the model is an alternating coordinate ascent algorithm on the ELBO. The M-step maximizes the ELBO with respect to the model parameters, $\phi$. The E-step maximizes the ELBO with respect to the variational distribution parameters, $\xi$.

M-step

We can maximize the ELBO with respect to $\pi$ by setting the derivative of the ELBO function with respect to $\pi$ and solving for $\pi$. However, we also require that $\sum_{k=1}^{K}\pi_{k}=1$, so we add a Lagrange multiplier to enforce this constraint. Solving for the maximum of the resulting Lagrange function gives the update formula,

[TABLE]

We see that the MLE for $\pi$ is the usual sample mean of $\tau$.

Setting the derivative of the ELBO Lagrange function with respect to $\eta_{m}$ and solving for $\eta_{m}$ gives the update formula,

[TABLE]

This natural parameter estimate gives the usual MLE estimate for the corresponding moment parameter $\hat{\Gamma}=\frac{1}{K}\sum_{k=1}^{K}(\nu_{k}^{2}+\gamma_{k})$.

E-step

Setting the derivative of the ELBO Lagrange function with respect to $\tau_{ik}$ and solving for $\tau_{ik}$ gives the update formula,

[TABLE]

where we have solved for $\lambda_{qi}$ by plugging $\hat{\tau}_{ik}$ into its constraint. The update for $\tau_{ik}$ can be rewritten using the moment parameters as

[TABLE]

This form shows that $\hat{\tau}_{ik}$ is proportional to the product of the prior for belonging to the $k$-th cluster and the likelihood of sampling $y_{i}$ from a Gaussian with mean $\nu_{k}$.

Setting the derivative of the ELBO Lagrange function with respect to $\nu_{k}$ and solving for $\nu_{k}$ gives the update formula,

[TABLE]

The update for $\nu_{k}$ can be rewritten using the moment parameters as

[TABLE]

This form shows that $\hat{\nu}_{k}$ is the sample mean of the weighted data assigned to the $k$-th cluster when the prior precision $\Gamma^{-1}=0$. This value of the precision corresponds to an improper prior on the cluster means and reduces the Bayesian Gaussian mixture model to a standard Gaussian mixture model. When $\Gamma^{-1}$ is greater than zero $\hat{\nu}_{k}$ is regularized (shrunken) towards the prior mean, which in this model is zero.

Setting the derivative of the ELBO Lagrange function with respect to $\gamma_{k}$ and solving for $\gamma_{k}$ gives the update formula,

[TABLE]

C.3 Proof of GOP conditions

[TABLE]

where $\alpha\triangleq\{\pi,\nu,\gamma\}$ and $\beta\triangleq\{\tau,\eta_{m}\}$.

The Hessian with respect to $\beta$ for fixed $\hat{\alpha}$ is

[TABLE]

where

[TABLE]

and

[TABLE]

We define $\tilde{\tau}$ to be the matrix $\tau$ construed as a vector taken column-wise. Since the matrix is diagonal and all of the diagonal elements are non-negative, the Hessian is positive semidefinite and $f(\hat{\alpha},\beta)$ is convex.

The Hessian with respect to $\alpha$ for fixed $\hat{\beta}$

[TABLE]

where

[TABLE]

[TABLE]

and

[TABLE]

Again, the matrix is diagonal and all of the diagonal elements are non-negative, so the Hessian is positive semidefinite and $f(\alpha,\hat{\beta})$ is convex.

Therefore, the first condition of the GOP algorithm is satisfied. The second condition is satisfied because we have no inequality constraints that are functions of both $\alpha$ and $\beta$ – all constraints can be absorbed into the sets $A,B$. The third condition is satisfied because the equality constraints are each convex combinations and therefore affine. The fourth condition (Slater’s condition) is satisfied because the interior of the feasible region has an interior point.