Relaxation of the EM Algorithm via Quantum Annealing for Gaussian Mixture Models

TL;DR

This paper introduces DQAEM, a novel algorithm combining quantum annealing with EM to improve Gaussian mixture model optimization, reducing local optima trapping and enhancing stability.

Contribution

It presents the DQAEM algorithm that integrates quantum annealing into EM, providing theoretical stability and demonstrating improved performance in Gaussian mixture models.

Findings

DQAEM outperforms traditional EM in avoiding local optima.

Theoretical proof of stability for DQAEM.

Numerical simulations confirm efficiency improvements.

Abstract

We propose a modified expectation-maximization algorithm by introducing the concept of quantum annealing, which we call the deterministic quantum annealing expectation-maximization (DQAEM) algorithm. The expectation-maximization (EM) algorithm is an established algorithm to compute maximum likelihood estimates and applied to many practical applications. However, it is known that EM heavily depends on initial values and its estimates are sometimes trapped by local optima. To solve such a problem, quantum annealing (QA) was proposed as a novel optimization approach motivated by quantum mechanics. By employing QA, we then formulate DQAEM and present a theorem that supports its stability. Finally, we demonstrate numerical simulations to confirm its efficiency.

Tables

Table 1. TABLE I : Ratios of success and failure for DQAEM and EM.

		DQAEM
		Success	Fail	Total
	Success	55.9 %	0.7 %	56.6 %
EM	Fail	41.5 %	1.9 %	43.4 %
	Total	97.4 %	2.6 %	100.0 %

Table 2. TABLE II : Ratios of success for DQAEM, EM and DSAEM.

DQAEM	EM	DSAEM
97.4 %	56.6 %	77.8 %

Equations

$$\mathcal{L}(Y_{\mathrm{obs}};\theta )$$

$$=\sum _{i=1}^N\log \sum _{{\sigma }^{(i)}\in {\Omega }^{(i)}}p(y^{(i)},{\sigma }^{(i)};\theta ).$$

$$\mathcal{L}(Y_{\mathrm{obs}};\theta )$$

$$\ge \sum _{i=1}^N\sum _{{\sigma }^{(i)}\in {\Omega }^{(i)}}P({\sigma }^{(i)}|y^{(i)};{\theta }^{\prime })\log \frac{p(y^{(i)},{\sigma }^{(i)};\theta )}{P({\sigma }^{(i)}|y^{(i)};{\theta }^{\prime })}$$

$$\ge Q(\theta ;{\theta }^{\prime })$$

$$-\sum _{i=1}^N\sum _{{\sigma }^{(i)}\in {\Omega }^{(i)}}P({\sigma }^{(i)}|y^{(i)};{\theta }^{\prime })\log P({\sigma }^{(i)}|y^{(i)};{\theta }^{\prime }),$$

$$Q(\theta ;{\theta }^{\prime })$$

$$P({\sigma }^{(i)}|y^{(i)};{\theta }^{\prime })$$

$$p(y^{(i)};{\theta }^{\prime })$$

$${\theta }^{(t+1)}=\underset{\theta }{\mathrm{argmax}}Q(\theta ;{\theta }^{(t)}).$$

$$1_k$$

$$p(y;\theta )$$

$$p(y|\sigma =1_k;\theta )$$

$$P(\sigma =1_k;\theta )$$

$$p(y,\sigma ;\theta )$$

$$=\prod _k{\left[{\pi }_{k}g(y;{\mu }_k,{\Sigma }_k)\right]}^{{\sigma }_k},$$

$$H(y,\sigma ;\theta )$$

$${\hat{\sigma }}_k$$

$$=\mathrm{diag}(\underset{k-1}{\underbrace{0,\dots ,0}},1,\underset{K-k}{\underbrace{0,\dots ,0}}),$$

$$H(y,\hat{\sigma };\theta )$$

$$=\mathrm{diag}(h_1,h_2,\dots ,h_K),$$

$$⟨\sigma =1_i|H(y,\hat{\sigma };\theta )|\sigma =1_j⟩$$

$$\mathcal{L}(Y_{\mathrm{obs}};\theta )$$

$$Q$$

$$p(y^{(i)},{\hat{\sigma }}^{(i)};\theta )$$

$$P({\hat{\sigma }}^{(i)}|y^{(i)};\theta )$$

$${\mathcal{Z}}^{(i)}(\theta )$$

$$p_{\Gamma }(y^{(i)},{\hat{\sigma }}^{(i)};\theta )$$

$$F_{\Gamma }(\theta )=-\log {\mathcal{Z}}_{\Gamma }(\theta ),$$

$${\mathcal{Z}}_{\Gamma }(\theta )$$

$${\mathcal{Z}}_{\Gamma }^{(i)}(\theta )$$

$$F_{\Gamma =0}(\theta )=-\mathcal{L}(Y_{\mathrm{obs}};\theta ).$$

$$U_{\Gamma }$$

$$P_{\Gamma }({\hat{\sigma }}^{(i)}|y^{(i)};\theta )$$

$${\theta }^{(t+1)}=\underset{\theta }{\mathrm{argmin}}{U}_{\Gamma }(\theta ,{\theta }^{(t)}),$$

$${\pi }_k^{(t+1)}$$

Full text

I Introduction

Combinatorial optimization is a fundamental issue in both science and engineering. Although some problems in such optimization can be efficiently solved by well-known algorithms [1, 2], other problems in a class of NP-hard, e.g. the traveling salesman problem, are essentially difficult to solve.

One of the effective approaches for NP-hard problems is simulated annealing (SA), which was proposed by Kirkpatrick et al. [3, 4]. SA is a generic approach for optimization, in which random numbers that mimic thermal fluctuations are used to go over potential barriers in objective functions. Furthermore, its global convergence is in some sense guaranteed by Geman and Geman et al. [5]. After that, a quantum extension of SA, which is called quantum annealing (QA), was proposed in physics [6, 7, 8], and has been intensively studied [9, 10, 11, 12, 13, 14, 15, 16, 17]. In QA, instead of thermal fluctuations, quantum fluctuations are used to overcome potential barriers in objective functions, and it has been reported that QA is more effective than SA for some problems [12]. Especially, due to quantum fluctuations, QA exhibits better performance than SA when objective functions have steep multimodality.

Such combinatorial optimization also appears in machine learning, which has attracted much interest recently [18, 19]. For example, some class of data clustering is known to be NP-hard problems [20]. One of common methods for data clustering is as follows. Assuming data points are generated by Gaussian mixture models (GMMs), we estimate the parameters in GMMs by the expectation-maximization (EM) algorithm [21]. However, parameter estimation sometimes fails since EM depends on initial values and suffers from the problem of local optima. To relax the problem, Ueda and Nakano proposed a deterministic simulated annealing expectation-maximization (DSAEM) algorithm 111This algorithm is called the deterministic annealing expectation-maximization algorithm in Ref. [22]., and it succeeds to relax the difficulty of the multimodality in EM. This algorithm is based on deterministic simulated annealing (DSA) 222This algorithm is called deterministic annealing in Ref. [23]., which was proposed by Rose et al. [23, 24]. The essence of these approaches is to make objective functions smooth by introducing thermal fluctuations without random numbers, and the non-convex problem in optimization is considerably managed without increase of numerical cost.

As we have explained, QA is considered to be effective than SA in some conditions [12], and thus the quantum version of DSA is expected to be superior to it. In this paper, we propose a deterministic quantum annealing expectation-maximization (DQAEM) algorithm for Gaussian mixture models because it is expected that quantum fluctuations can relax the problem of local optima in parameter estimation. In our previous paper [25], we proposed DQAEM for continuous latent variables, and obtained the result that DQAEM outperformed EM. However, its applicability is limited because the latent variables are assumed to be continuous and most difficulties in parameter estimation come from optimization of discrete latent variables, such as Gaussian mixture models. Thus, in this paper, we develop DQAEM for discrete latent variables and apply it to GMMs. After the formulation of the algorithm, we present a theorem that guarantees its stability. Finally, to illustrate its efficiency compared to EM, we show numerical simulations, in which DQAEM is applied to GMMs for data clustering.

This paper is organized as follows. In Sec. II, we review GMMs and EM to prepare for DQAEM. In Sec. III, which is the main section of this paper, we describe the formulation of DQAEM in detail and present a theorem on its convergence. In Sec. IV, we demonstrate numerical simulations and discuss its efficiency. In Sec. V, we conclude this paper.

II Review of the expectation-maximization (EM) algorithm and Gaussian mixture models (GMMs)

In this section, we review EM to prepare for introducing our DQAEM, and consider an estimation problem of GMMs to formulate DQAEM because it is one of the simplest models with discrete variables.

II-A Maximum likelihood estimation (MLE) and the expectation-maximization (EM) algorithm

The aim of this subsection is to describe EM because DQAEM is based on it. First, we review maximum likelihood estimation (MLE) briefly. Suppose we have $N$ data points $Y_{\mathrm{obs}}=\{y^{(1)},y^{(2)},\dots,y^{(N)}\}$ and they are independent and identically distributed obeying $p(y^{(i)};\theta)$ where $\theta$ is a parameter. Moreover we define $p(y^{(i)},\sigma^{(i)};\theta)$ as the probability density functions for complete data with the unobservable variables $\{\sigma^{(1)},\sigma^{(2)},\dots,\sigma^{(N)}\}$. Namely, $p(y^{(i)};\theta)=\sum_{\sigma^{(i)}\in\Omega^{(i)}}p(y^{(i)},\sigma^{(i)};\theta)$, where $\Omega^{(i)}$ represents the domain of $\sigma^{(i)}$. Then the log likelihood function is given by

[TABLE]

Note that $i$ in $y^{(i)}$ and $\sigma^{(i)}$ is the index for each observed data point. MLE is a technique to estimate the parameter $\theta$ in model distributions that maximize the log likelihood function $\mathcal{L}(Y_{\mathrm{obs}};\theta)$.

In general, maximizing the log likelihood function $\mathcal{L}(Y_{\mathrm{obs}};\theta)$ with respect to $\theta$ is difficult because it is sometimes a non-convex optimization, and then we replace it with its lower bound. Using Jensen’s inequality, we have the following inequality

[TABLE]

where $\theta^{\prime}$ is an arbitrary parameter and $P(\sigma^{(i)}|y^{(i)};\theta^{\prime})$ is the conditional probability. Then, the procedure of EM consists of the following two steps. The first one, which is called the E step, is to compute the conditional probability $P(\sigma^{(i)}|y^{(i)};\theta)$ by

[TABLE]

Here we have used Bayes’ rule. The second one, which is called the M step, is to maximize the $Q$ function (2) with respect to $\theta$ instead of $\mathcal{L}(Y_{\mathrm{obs}};\theta)$. Denoting the tentative estimated parameter at the $t$-th iteration by $\theta^{(t)}$, the estimated parameter is updated by

[TABLE]

At the end of this subsection, we summarize EM in Algo. 1.

II-B Gaussian mixture models (GMMs)

Here we introduce GMMs and its quantum mechanical representation. We follow the notations in Refs. [18, 19]. Let $y$ and $\sigma$ denote continuous observable and discrete unobservable variables. Here, we assume that $\Omega$, which is the domain of $\sigma$, is given by $\{1_{k}\}_{k=1}^{K}$, where

[TABLE]

for $k=1,\dots,K$, and then the number of elements in $\Omega$ is $K$. Specifically, $\sigma=1_{k}$ when $\sigma$ denotes the $k$-th element in $\Omega$.

Using the above notation, the probability density function of GMMs is given by

[TABLE]

where

[TABLE]

$\{\pi_{k}\}_{k=1}^{K}$ satisfies $\sum_{k}\pi_{k}=1$, $g(y;\mu_{k},\Sigma_{k})$ is a Gaussian function with mean $\mu_{k}$ and covariance $\Sigma_{k}$ for $k=1,\dots,K$, and $\theta=\{\pi_{k},\mu_{k},\Sigma_{k}\}_{k=1}^{K}$. The joint probability density function for GMMs is therefore given by

[TABLE]

where $\sigma_{k}$ is the $k$-th element of $\sigma$.

To introduce quantum fluctuations, we need to rewrite the above equations in the Hamiltonian formulation. Taking the logarithm of (4), the Hamiltonian for GMMs can then be written as

[TABLE]

where $h_{k}=-\log\{\pi_{k}g(y;\mu_{k},\Sigma_{k})\}$ for $k=1,\dots,K$. Here, we introduce ket vectors, bra vectors and “spin” operators to rewrite (5) in the manner of quantum mechanics. First, we define the ket vector $|\sigma=1_{k}\rangle$ by $1_{k}$ and the “spin” operator $\hat{\sigma}_{k}$ by

[TABLE]

respectively, where the bra vector $\langle\sigma=1_{j}|$ satisfies the orthonormal condition $\langle\sigma=1_{j}|\sigma=1_{i}\rangle=\delta_{ij}$. Replacing $\sigma$ with $\hat{\sigma}$, we have the Hamiltonian operator

[TABLE]

and this satisfies

[TABLE]

where $\delta_{ij}$ is the Kronecker delta. We use this formulation to describe DQAEM in the following section. Note that a similar expression is presented in Ref. [26].

III Deterministic quantum annealing expectation-maximization algorithm (DQAEM)

First, we formulate DQAEM by using the quantum representation described in the previous section. Then we discuss its stability by showing the monotonicity of the free energy during the algorithm.

III-A Formulation

In this subsection, we formulate DQAEM by employing the concept of quantum annealing [8] (also see App. -A). First, we rewrite EM in the quantum representation. The log likelihood function (1) is rewritten as

[TABLE]

Note that $\mathrm{Tr}\left[\cdot\right]=\sum_{k}\langle\sigma^{(i)}=1_{k}|\left[\cdot\right]|\sigma^{(i)}=1_{k}\rangle$. As we have explained in Sec. II-A, the $Q$ function (2) is maximized in the M step of EM. Similarly to (7), the quantum representation of the $Q$ function (2) is given by

[TABLE]

where

[TABLE]

and $H(y^{(i)},\hat{\sigma}^{(i)};\theta)$ is in Eq. (6). Furthermore, the conditional probability $P(\hat{\sigma}^{(i)}|y^{(i)};\theta)$ is computed using Bayes’ rule. That is,

[TABLE]

Here, the normalization factor, which is called the partition function in physics, has the form

[TABLE]

Now we begin to formulate DQAEM. To introduce quantum fluctuations, we add $H^{\prime}(\Gamma)=\Gamma\hat{\sigma}^{\prime}$ whose $\hat{\sigma}^{\prime}$ satisfies $[\hat{\sigma}_{k},\hat{\sigma}^{\prime}]\neq 0$ for $k=1,\dots,K$ to the original Hamiltonian $H$, and then (8) is converted to

[TABLE]

In MLE, the log likelihood function (7) is optimized. On the other hand, the objective function in DQAEM, which is called the free energy, is given by

[TABLE]

where

[TABLE]

By taking in into account $H^{\prime}(0)=0$ and comparing to Eq. (7), we obtain the relation between the free energy and the log likelihood function as

[TABLE]

Thus we can say that the negative free energy at $\Gamma=0$ is the log likelihood function.

Next, we define the function $U_{\Gamma}(\theta;\theta^{\prime})$ to formulate DQAEM, which corresponds to the $Q$ function in EM. Using (9), the function $U_{\Gamma}(\theta;\theta^{\prime})$ has the form

[TABLE]

where

[TABLE]

Then DQAEM is composed of the following two steps. The first one is to compute the conditional probability (12), and this is called the E step of DQAEM. The second one is to update the parameter $\theta^{(t)}$ by minimizing the function $U_{\Gamma}(\theta;\theta^{\prime})$ (11). That is,

[TABLE]

and this is called the M step in DQAEM. Furthermore, we decrease $\Gamma$ during the iterations. We summarize DQAEM in Algo. 2.

III-B Convergence theorem

We have proposed DQAEM in the previous subsection. Here, we present the theorem that guarantees its stability via iterations.

Theorem 1

Let $\theta^{(t+1)}=\operatorname*{arg\,min}_{\theta}U_{\Gamma}(\theta;\theta^{(t)})$. Then $F_{\Gamma}(\theta^{(t+1)})\leq F_{\Gamma}(\theta^{(t)})$ holds. Moreover, the equality holds if and only if $U_{\Gamma}(\theta^{(t+1)};\theta^{(t)})=U_{\Gamma}(\theta^{(t)};\theta^{(t)})$ and $S_{\Gamma}(\theta^{(t+1)};\theta^{(t)})=S_{\Gamma}(\theta^{(t)};\theta^{(t)})$, where $S_{\Gamma}(\theta;\theta^{\prime})=\sum_{i=1}^{N}\mathrm{Tr}\left[P_{\Gamma}(\hat{\sigma}^{(i)}|y^{(i)};\theta^{\prime})\log P_{\Gamma}(\hat{\sigma}^{(i)}|y^{(i)};\theta)\right]$.

This theorem insists that DQAEM converges at least the global optimum or a local optimum. We mention that the global convergence of EM is discussed by Dempster et al. [21] and Wu [27], and their discussions apply to DQAEM.

IV Numerical simulations

In this section, we carry out numerical simulations to confirm the performance of DQAEM. In the first subsection, we present the setup of numerical simulations, and, in the following subsection, we provide numerical results.

IV-A Mathematical setup

We estimate the parameters of GMMs by using both DQAEM and EM. Suppose $N$ data points $Y_{\mathrm{obs}}=\{y^{(1)},y^{(2)},\dots,y^{(N)}\}$ are identically sampled by GMMs with $K=3$. Here, a GMM is given by (4). In EM, the updating equations for $\theta=\{\pi_{k},\mu_{k},\Sigma_{k}\}_{k=1}^{K}$ are determined by the derivative of the $Q$ function (2) with respect to $\theta$. The parameter $\theta^{(t+1)}$ of GMMs at the $t+1$-th iteration is then given by

[TABLE]

where $\theta^{(t)}$ is the tentative estimated parameter at the $t$-th iteration.

In DQAEM, the updating equations for $\theta$ are determined by the derivative of the function $U_{\Gamma}(\theta,\theta^{\prime})$ in (11) with respect to $\theta$, and then $P(\sigma^{(i)}=1_{k}|y^{(i)};\theta^{(t)})$ in (13), (14) and (15) are replaced by $P_{\mathrm{QA}}(\sigma^{(i)}=1_{k}|y^{(i)};\theta^{(t)})=\langle\sigma^{(i)}=1_{k}|P_{\Gamma}(\hat{\sigma}^{(i)}|y^{(i)};\theta^{(t)})|\sigma^{(i)}=1_{k}\rangle$. That is, the updating equations for DQAEM are given by

[TABLE]

Note that the quantum effects for parameter estimation comes from $P_{\mathrm{QA}}(\sigma^{(i)}=1_{k}|y^{(i)};\theta^{(t)})$. The annealing parameter $\Gamma$ are varied from initial values to [math] via iterations.

In this section, assume that, in matrix notation, $\hat{\sigma}^{\prime}$ is given by

[TABLE]

Obviously $[\hat{\sigma}_{k},\hat{\sigma}^{\prime}]\neq 0$ is satisfied. Note that the size of the Hamiltonian is determined by assumed number of mixtures.

IV-B Numerical results

In this subsection, using the data set shown in Fig. 1(a), we compare DQAEM, EM, and DSAEM, which was proposed in Ref. [22]. This data set is generated by the GMM that consists of three two-dimensional Gaussian functions whose means are $(X,Y)=(-3,0)$, $(0,0)$ and $(3,0)$. Here we set $\Gamma_{\mathrm{init}}=1.0$ in DQAEM to discuss the effect of quantum fluctuations simply. We also choose the annealing parameter in DSAEM as $\beta_{\mathrm{init}}=0.7$. Note that, in DSAEM, the annealing parameter is given by temperature. Furthermore, we exponentially vary $\beta$ and $\Gamma$ to 1 and 0, respectively. We plot transitions of the log likelihood functions of EM and the negative free energies of DSAEM and DQAEM in Fig. 1(b) by red lines, orange lines, and blue lines, respectively. The value of $-712.1$ depicted by the green line in Fig. 1(b) is the optimal value in these numerical simulations. DQAEM, EM, and DSAEM give the optimal estimate or suboptimal estimates depending on initial optimization values.

To understand visually how DQAEM and EM behave in parameter estimation, we illustrate estimated Gaussian functions in the case where the log likelihood function is $-712.1$ in Fig. 2(a) and in one of the cases where the log likelihood function is lower than the optimal value in Fig. 2(b). The case demonstrated in Fig. 2(b) clearly fails in data clustering.

However, the ratios of success for DQAEM, EM, and DSAEM are much different. To see the ratios of success and failure for DQAEM and EM, we performed DQAEM and EM with same initial optimization values $1000$ times, respectively, and summarized the results in Table I.

Here, we have defined the “success” of DQAEM and EM when square errors between the estimated means of three Gaussian functions and the true means are less than $0.3$ times the covariances of three Gaussian functions. Table I shows that DQAEM succeeds with the ratio of $97.4$ % while EM succeeds with the ratio of $56.6$ %, and that DQAEM is superior to EM. In Table II, we show the ratios of success for DQAEM, EM, and DSAEM in parameter estimation. This table also shows that DQAEM is superior to DSAEM.

V Conclusion

In this paper, we have proposed the deterministic quantum annealing expectation-maximization (DQAEM) algorithm for Gaussian mixture models (GMMs) to relax the problem of local optima of the expectation-maximization (EM) algorithm by introducing the mechanism of quantum fluctuations into EM. Although we have limited our attention to GMMs in this paper to simplify the discussion, the derivation presented in this paper can be straightforwardly applied to any models which have discrete latent variables. After formulating DQAEM, we have presented the theorem that guarantees its convergence. We then have given numerical simulations to show its efficiency compared to EM and DSAEM. It is expect that the combination of DQAEM and DSAEM gives better performance than DQAEM. Finally, one of our future works is a Bayesian extension of this work. In other words, we are going to propose a deterministic quantum annealing variational Bayes inference.

-A Quantum annealing

Here, we briefly introduce “quantum” annealing (QA) to prepare for DQAEM. First we consider the minimization problem of the Ising model. That is,

[TABLE]

where

[TABLE]

$\sigma_{z}^{(i)}=\pm 1$ for each $i$, and $J_{ij}$ is the coupling constant between spins at site $i$ and site $j$. Note that this problem can describe many combinatorial problems such as the traveling salesman problem and the max-cut problem [28].

In QA, we quantize the Ising model (16) by applying magnetic fields along the $x$ axis to the model and solve the Schrödinger equation on this system while decreasing the magnetic fields. Then the Hamiltonian of this system is given by

[TABLE]

where

[TABLE]

using

[TABLE]

and $\Gamma$ represents the strength of the magnetic fields. This is called the Transverse Ising model. Thus the Schrödinger equation that we solve in QA is given by

[TABLE]

where $i$ is the imaginary unit, $\hbar$ is the Dirac constant, and $|\psi\rangle$ is the ket vector. The magnetic field $\Gamma$ is set to be large at the beginning of QA, and then $|\psi\rangle$ is initially equal or close to the eigenstate of $\sum_{i}\hat{\sigma}_{x}^{(i)}$. During solving (17), we gradually decrease $\Gamma$ and finally make $\Gamma$ go to zero. Therefore, $|\psi\rangle$ gives a solution for the original Hamiltonian (16). The efficiency of QA is discussed in Refs. [8, 11, 12].