Quantum Expectation-Maximization Algorithm

TL;DR

This paper introduces a quantum expectation-maximization algorithm for Gaussian mixture models, demonstrating its robustness, quantum speedup, and advantages over k-means in complex clustering scenarios.

Contribution

It extends quantum clustering methods by developing a quantum EM algorithm for GMMs, showing improved performance and speedup over existing quantum and classical algorithms.

Findings

The quantum EM algorithm is robust and faster than classical methods.

Numerical results show GMM's advantage over k-means for complex data.

Quantum speedup is demonstrated for the proposed algorithm.

Abstract

Clustering algorithms are a cornerstone of machine learning applications. Recently, a quantum algorithm for clustering based on the k-means algorithm has been proposed by Kerenidis, Landman, Luongo and Prakash. Based on their work, we propose a quantum expectation-maximization (EM) algorithm for Gaussian mixture models (GMMs). The robustness and quantum speedup of the algorithm is demonstrated. We also show numerically the advantage of GMM over k-means for non-trivial cluster data.

Figures

Tables

Table 1. Table 1 : Success rates of the EM algorithm, the δ 𝛿 \delta -EM algorithm, the k 𝑘 k -means algorithm, and the δ 𝛿 \delta - k 𝑘 k -means algorithm. These scores are best ones in one hundred trials with randomized initial inputs.

EM	$\delta$-EM	$k$-means	$\delta$-$k$-means
93.9 %	94.3 %	72.4 %	72.5 %

Table 2. Table 2 : Success rates of the EM algorithm, the δ 𝛿 \delta -EM algorithm, the k 𝑘 k -means algorithm, and the δ 𝛿 \delta - k 𝑘 k -means algorithm. These scores are best ones in one hundred trials with randomized initial inputs.

EM	$\delta$-EM	$k$-means	$\delta$-$k$-means
88.8 %	89.2 %	57.9 %	55.4 %

Equations

$$p(y;\theta )$$

$$r_t^{i,k}$$

$${\pi }_{t+1}^k$$

$${\mu }_{t+1}^k$$

$${\Sigma }_{t+1}^k$$

$$d_{\mathrm{G}}^k(y_i)$$

$$d_{\mathrm{E}}^k(y_i)$$

$$L_{\mathrm{G}}^{\delta }(y_i)$$

$$\displaystyle\frac{1}{\sqrt{N}}\sum_{i=1}^{N}\Ket{i}\Big{(}\otimes_{k\in[K]}\Ket{k}\Ket{0}\Big{)}$$

$$\displaystyle\quad\mapsto\frac{1}{\sqrt{N}}\sum_{i=1}^{N}\Ket{i}\Big{(}\otimes_{k\in[K]}\Ket{k}|\tilde{d}_{\mathrm{G}}^{k}(y_{i})\rangle\Big{)},$$

$$|i⟩|0⟩$$

$$|k⟩|0⟩$$

$$|{\varphi }_{i,k}⟩$$

$$|{\varphi }_{i,k}^C⟩$$

$$\displaystyle|\phi_{i,k}^{H}\rangle=\frac{1}{2}\Ket{i}\Ket{k}\Big{(}$$

$$\displaystyle+\Ket{1}(\hat{G}^{k}\Ket{y_{i}}-\hat{G}^{k}|\mu^{k}\rangle)\Big{)}.$$

$$|{\varphi }_{i,k}^H⟩$$

$${\hat{U}}_1$$

$$\displaystyle\quad\mapsto\Ket{i}\Ket{j}\Big{(}\sqrt{p_{i,k}}\Ket{\mathrm{tar}_{i,k},1}+\sqrt{1-p_{i,k}}\Ket{\mathrm{gar}_{i,k},0}\Big{)}.$$

$$p_{i,k}=\frac{1-⟨y_i|({\hat{G}}^k{)}^2|{\mu }^k⟩}{2}.$$

$${\hat{U}}_2$$

$$\mapsto |i⟩|k⟩(\sqrt{\alpha }|{\tilde{p}}_{i,k},1⟩+\sqrt{1-\alpha }|{\mathrm{gar}}_{i,k},0⟩),$$

$$\Vert |{\Phi }_{i,k}⟩-{|0⟩}^{\otimes L}|{\tilde{p}}_{i,k}⟩{\Vert }_2$$

$${\hat{U}}_3$$

$${\hat{U}}_4$$

$$\mapsto \frac{1}{\sqrt{N}}\sum _{i=1}^N|i⟩|{\mathrm{label}}_t(y_i)⟩,$$

$$|{\psi }_t⟩$$

$$|{\psi }_t⟩$$

$$=\sum _{k=1}^K\sqrt{\frac{N_t^k}{N}}|{\chi }_t^k⟩|k⟩.$$

$$|{\chi }_t^k⟩$$

$$({\chi }_t^k{)}_i$$

$$V_1$$

$$V_2$$

$$V_{0,i}$$

$$|{\mu }_{t+1}^k⟩$$

$$|\mathrm{vec}[{\Sigma }_{t+1}^k]+{\mu }_{t+1}^k\otimes {\mu }_{t+1}^k⟩$$

Full text

††preprint: APS/123-QED

I Introduction

Quantum computing has attracted much attention since the discovery of Shor’s algorithm Shor (1999); Nielsen and Chuang (2002). Recently, with the rapid developments in machine learning, physicists have started to consider utilizing quantum computers for machine learning applications Biamonte et al. (2017); Lloyd et al. (2013); Schuld et al. (2015); Rebentrost et al. (2014); Wiebe et al. (2014); Dunjko et al. (2016). As a result, quantum machine learning has emerged as an interdisciplinary field between quantum computing and machine learning. Furthermore, a quantum algorithm for the $k$-means algorithm Bishop (2007); Murphy (2012) with proven quantum speedup was proposed Kerenidis et al. (2018).

The $k$-means algorithm is an essential tool in many machine learning applications  Bishop (2007); Murphy (2012). However, the $k$-means algorithm is as a special case of the more general Gaussian mixture model (GMM). In the $k$-means algorithm, each Gaussian has the same weight and the covariance matrix of each Gaussian function is the identity. As a result, the $k$-means algorithm may provide poor estimates of the clusters since the assumptions of the $k$-means algorithm are sometimes too strong to capture all properties of complex data sets. The expectation-maximization (EM) algorithm Dempster et al. (1977); Bishop (2007); Murphy (2012) and variational Bayes (VB) inference Bishop (2007); Murphy (2012) with the GMM are often used to improve the clustering, since the general GMM can deal with a wider class of data sets. Recently, one of the authors proposed quantum-inspired algorithms for the EM algorithm Miyahara and Tsumura (2016); Miyahara et al. (2016, 2017) and VB Miyahara and Sughiyama (2018). In Refs. Miyahara et al. (2017); Miyahara and Sughiyama (2018), we have succeeded in improving the performances of the EM algorithm and VB. However, the aim of Refs. Miyahara et al. (2017); Miyahara and Sughiyama (2018) is to make use of quantum fluctuations as a numerical tool and not to provide a quantum speedup over a classical algorithm; as a result, the computational costs are almost the same.

In this paper, we propose a quantum algorithm to estimate the parameters of the GMM, which we call the quantum EM (q-EM) algorithm. To this end, following the spirit of Ref. Kerenidis et al. (2018), we first introduce a randomized variant of the EM algorithm, i.e. the $\delta$-EM algorithm which includes non-deterministic readout of the data. Then, we formulate a quantum algorithm that realizes a speedup of the $\delta$-EM algorithm with respect to the number of data points. The q-EM algorithm may be an important step toward quantum machine learning, since the EM algorithm is an essential algorithm in machine learning.

This paper is organized as follows. In Sec. II, we provide classical preliminaries. In particular, we review the EM algorithm and then introduce the $\delta$-EM algorithm. In Sec. III, we present the detailed procedure of the q-EM algorithm. Then, in Sec. IV, we show its computational cost. In Sec. V, we show numerical simulations of the $\delta$-EM algorithm to confirm that our starting point is valid. In Sec. VI, we discuss the relationship between the EM algorithm and the $k$-means algorithm. Finally, Sec. VII concludes this paper.

Independently of our work, Iordanis Kerenidis, Alessandro Luongo, and Anupam Prakash proposed an extension of the q-means algorithm to Gaussian mixture models similar to this work using soft clustering Kerenidis et al. (2019).

II Classical preliminaries

In this section, we first review the EM algorithm Bishop (2007); Murphy (2012) in detail. We then introduce a randomized variant of the EM algorithm, which we call the $\delta$-EM algorithm. The purpose of introducing the $\delta$-EM algorithm is the need of a robust variant of the EM algorithm against noise which resembles the non-deterministic quantum measurement of a superposed state.

II.1 EM algorithm

The EM algorithm is a generic approach to estimate parameters of probability distributions based on maximum likelihood estimation. For simplicity, we focus on the GMM and review the EM algorithm for the GMM. Let us consider a $d$-dimensional feature space and assume that we have $N$ data points $\{y_{i}\}_{i=1}^{N}$. The GMM for $x\in\mathbb{R}^{d}$ is given by

[TABLE]

where $\mathcal{N}(y;\mu^{k},\Sigma^{k})\coloneqq\frac{1}{(2\pi)^{\frac{d}{2}}|\Sigma^{k}|^{\frac{1}{2}}}e^{-\frac{1}{2}(y-\mu^{k})(\Sigma^{k})^{-1}(y-\mu^{k})}$ is the $d$-dimensional Gaussian function with mean $\mu^{k}$ and covariance $\Sigma^{k}$, and $\sum_{k}\pi^{k}=1$. To simplify the notation, we define $\theta\coloneqq\{\pi^{k},\mu^{k},\Sigma^{k}\}_{k=1}^{K}$.

The EM algorithm, which estimates $\theta$, consists of the following two steps which are iterated until convergence. The first step, which is called the E step, is to compute the responsibilities of cluster $k$ for each datapoint $y_{i}$:

[TABLE]

The second step, which called the M step, is to compute $\theta_{t}$ by using the responsibilities, Eq. (2):

[TABLE]

We iterate the E and M step by substituting Eqs. (3), (4), and (5) until convergence. Note that we can begin either of the E step or the M steps for the first iteration. The EM algorithm is summarized in Algo. 1.

Note that the procedure of the EM algorithm can be generalized for mixture models Bishop (2007); Murphy (2012).

II.2 $\delta$-EM algorithm

As a prerequisite to the q-EM algorithm, we need to modify the original EM algorithm, since we have to take into account randomness associated with quantum measurement.

In Sec. II.1, we explained that the EM algorithm has two steps: the E and M steps. In the $\delta$-EM algorithm, we modify the E step in the spirit of the $\delta$-$k$-means-algorithm in Ref. Kerenidis et al. (2018).

To this end, we first introduce the square GMM distance by

[TABLE]

Note that, when $\pi^{k}=1/K$ and $\Sigma^{k}$ is the identity matrix for $k=1,2,\dots,K$, $d_{\mathrm{G}}^{k}(y_{i})=d_{\mathrm{E}}^{k}(y_{i})$, where $d_{\mathrm{E}}^{k}(\cdot)$ is the square Euclidean distance given by

[TABLE]

We then define the set of labels given by

[TABLE]

where $d_{\mathrm{G}}^{*}(y_{i})\coloneqq\min_{k}d_{\mathrm{G}}^{k}(y_{i})$. In the E step of the $\delta$-EM algorithm, we take random samples from $L_{\mathrm{G}}^{\delta}(y_{i})$ in Eq. (8). We note that for soft clustering, a more precise sampling scheme may be useful, but this simple sampling works well in the approach shown in Ref. Kerenidis et al. (2018).

In the M step of the $\delta$-EM algorithm, we add small noise to the estimated parameters after their estimation. As a result, the $\delta$-EM algorithm becomes robust, and its quantum version will become implementable. In Sec. V, we will show the validity of this algorithm numerically. Remarkably, we find that adding noise can even improve the quality of the studied benchmark examples.

III Quantum algorithm for the EM algorithm with the GMM

In this section, we describe the procedure of the quantum algorithm that realizes a quantum speedup of the EM algorithm for the GMM.

To simplify the notation, we add the tilde for estimates throughout this paper; that is, we denote e.g. $\tilde{a}$ as the estimate of $a$.

III.1 Overview of the q-EM algorithm

We begin with the initialization of the q-EM algorithm. In the EM algorithm, we can begin either with the E step or the M step for the first iteration. For simplicity, in the case of the q-EM algorithm, we consider to start with the E step; then, we set the initial parameter set $\theta_{0}=\{\pi_{0},\mu_{0},\Sigma_{0}\}$ with $\pi_{0}\coloneqq[\pi_{0}^{1},\pi_{0}^{2},\dots,\pi_{0}^{K}]$, $\mu_{0}\coloneqq[\mu_{0}^{1},\mu_{0}^{2},\dots,\mu_{0}^{K}]$, and $\Sigma_{0}\coloneqq[\Sigma_{0}^{1},\Sigma_{0}^{2},\dots,\Sigma_{0}^{K}]$.

The main procedure of the q-EM algorithm is composed of four steps. In step I, we compute the square GMM distance, and in step II, it is minimized for cluster assignment. Then, in step III, we generate quantum states of weight vectors, mean vectors, and covariance matrices, and in step IV, we apply quantum vector state tomography. By using the classical information on mean vectors and covariance matrices, we repeat the whole procedures until convergence. The output of this algorithm is $\theta_{*}=\{\pi_{*},\mu_{*},\Sigma_{*}\}$. The q-EM algorithm is summarized in Algo. 2.

In the rest of this section, we will explain the four steps in detail.

III.2 Step I: Computing the square GMM distance

In this step, we compute the square GMM distance, Eq. (6). Mathematically, we apply the unitary operation:

[TABLE]

where $d_{\mathrm{G}}^{k}(y_{i})$ is the square GMM distance between $y_{i}$ and the $k$-th cluster, and $[K]\coloneqq\{k\}_{k=1}^{K}$. For this computation, we require the precision given by $\|\tilde{d}_{\mathrm{G}}^{k}(y_{i})-d_{\mathrm{G}}^{k}(y_{i})\|\leq\epsilon_{1}$. In the next section, $\epsilon_{1}$ will be used to analyze the runtime.

Equation (9) includes summation, multiplication, and inner products of quantum states. Among them, the computation of summation and multiplication is straightforward with quantum linear algebra while the implementation of inner products is more involved. Let us thus focus on the computation of inner products.

We assume that two unitary operations and their controlled versions are available as follows:

[TABLE]

where $\hat{G}^{k}\coloneqq(\hat{\Sigma}^{k})^{1/2}$ for $k=1,2,\dots,K$. We then begin with the state

[TABLE]

By using controlled versions of Eq. (10) and (11), we create the following state from $\Ket{\phi_{i,k}}$:

[TABLE]

Then we apply the Hadamard gate on the third register of $|\phi_{i,k}^{C}\rangle$ and the resulting state is

[TABLE]

Note that Eq. (14) is also represented as

[TABLE]

where $\Ket{\mathrm{tar}_{i,k},1}\coloneqq\Ket{1}\hat{G}^{k}(\Ket{y_{i}}-|\mu^{k}\rangle)$ and $\Ket{\mathrm{gar}_{i,k},0}$ is a garbage state. That is, we have a unitary operator such that

[TABLE]

We also note that the probability that we get $\Ket{1}$ by measuring the third register is expressed as

[TABLE]

We then apply all the operations except the measurement in amplitude estimation in Ref. Brassard et al. (2002); Kerenidis et al. (2018) on $\hat{U}_{1}$ in Eq. (16). This process realizes the following unitary operation:

[TABLE]

where $\|\tilde{p}_{i,k}-p_{i,k}\|<2\pi\sqrt{p_{i,k}(1-p_{i,k})}/P_{\mathrm{ae}}+\pi^{2}/P_{\mathrm{ae}}^{2}$ and $\alpha>\frac{8}{\pi^{2}}$ Brassard et al. (2002). Here $M$ is a parameter to be determined (see Sec. B.1). Next, applying the mode evaluation method 111This algorithm has no specific name in Ref. Wiebe et al. (2014), and it is called median evaluation in Ref. Kerenidis et al. (2018). But this algorithm realize majority voting; so we call it mode evaluation. in Lemma 8 of Ref. Wiebe et al. (2014) and Thm. 2.2 of Kerenidis et al. (2018) to Eq. (18), we get $\Ket{\Phi_{i,k}}$ such that

[TABLE]

The last step is to estimate the square GMM distance of unnormalized vectors $\|y_{i}\|$ and $\|\mu^{k}\|$ and to multiply the norms of them and adding $\ln\pi^{k}$. A translation operator $\hat{T}(r^{\prime})$ can conduct the adding operation: $\hat{T}(r^{\prime})\Ket{r}=\Ket{r+r^{\prime}}$ for $r,r^{\prime}\in\mathbb{R}^{N}$. Note that we have assumed that we know the norms of $\{\|y_{i}\|\}_{i=1}^{N}$ and the same assumption is used in Ref. Kerenidis et al. (2018).

III.3 Step II: Assignment of clusters

The purpose of step II is cluster assignment. In this step, we utilize the following unitary operation:

[TABLE]

where $|a^{k}\rangle$ is a $(\ln p)$-bit state for $k=1,2,\dots,K$. The computational cost of this operation is $O(K\ln p)$ Kerenidis et al. (2018).

To find cluster assignment, we perform the unitary operation given by

[TABLE]

where $\mathrm{label}_{t}(y_{i})$ is the optimal label of $y_{i}$ at time $t$.

Finally, by uncomputing the square GMM distances, we obtain

[TABLE]

This uncomputation is required to repeat iterations.

III.4 Step III: Generation of the mean and covariance states

In this step, we generate states that store information on the weights, mean vectors, and covariance matrices. Let us recall that Eq. (22) is also expressed as

[TABLE]

Thus, by measuring the label register of $|\psi_{t}\rangle$ in Eq. (22), we obtain, with probability $N_{t}^{k}/N$,

[TABLE]

where $C_{t}^{k}$ is the set of labels that belong to cluster $k$ at time $t$. Then, $\chi_{t}^{k}=[\dots,(\chi_{t}^{k})_{i-1},(\chi_{t}^{k})_{i},(\chi_{t}^{k})_{i+1},\dots]^{\intercal}\in\mathbb{R}^{N}$ is

[TABLE]

for $i=1,2,\dots,N$.

We here define $V_{1}\in\mathbb{R}^{N\times d}$, $V_{2}\in\mathbb{R}^{N\times d\times d}$, and $V_{0,i}$ for $i=1,2,\dots,N$ on a QRAM:

[TABLE]

To obtain mean vectors, we multiply $V_{1}$ to the state $|\chi_{t}^{k}\rangle$ in Eq. (25) by using quantum linear algebra Chakraborty et al. (2018); Kerenidis and Luongo (2018):

[TABLE]

The associated error is $\epsilon_{2}^{\mu}$. Similarly, we compute a state involving information on $\Sigma^{k}$ by using quantum linear algebra Chakraborty et al. (2018); Kerenidis and Luongo (2018):

[TABLE]

The associated error is $\epsilon_{2}^{\Sigma}$. Note that $\epsilon_{2}^{\mu}$ and $\epsilon_{2}^{\Sigma}$ appear only in logarithms; thus, we do not explicitly consider them.

We finally deal with $\{\pi^{k}\}_{k=1}^{K}$; it is relatively easy to compute the weights of the GMM, $\{\pi_{t}^{k}\}_{k=1}^{K}$. We here utilize Eq. (29) as follows:

[TABLE]

Thus, we can estimate $N_{t}^{k}$ similarly. Note that we assume that the sizes of all clusters are $\Omega(N/k)$.

III.5 Step IV: Update of the parameters

At the end of each iteration, we obtain classical information on $\pi_{t+1}$, $\mu_{t+1}$, and $\Sigma_{t+1}$ by performing the quantum state tomography algorithm for $|\chi_{t}^{k}\rangle$, $|\mu_{t+1}^{k}\rangle$, and $|\mathrm{vec}[\Sigma_{t+1}^{k}]\rangle$. Quantum vector state tomography is explained in Ref. Kerenidis et al. (2018). The quantum state tomography algorithm in Ref. Kerenidis et al. (2018) requires a unitary transformation $U:\Ket{0}\mapsto\Ket{x}$; however, the procedure to find $l(y_{i})$ is not deterministic. Then, we have to devise some deterministic methods to find $l(y_{i})$. One solution is to determine $l(y_{i})$ by the rule

[TABLE]

if $d_{\mathrm{G}}^{k}(y_{i})<d_{\mathrm{G}}^{k^{\prime}}(y_{i})-2\delta$ for $k^{\prime}\neq k$, and we discard the points to which no label can be assigned.

By introducing $\epsilon_{4}^{\pi}$, $\epsilon_{3}^{\mu}$, $\epsilon_{4}^{\mu}$, $\epsilon_{3}^{\Sigma}$, and $\epsilon_{4}^{\Sigma}$, we require the following precision in this step: \big{\|}\|\pi^{k}\|-\|\tilde{\pi}^{k}\|\big{\|}\leq\epsilon_{4}^{\pi}, \big{\|}|\mu^{k}\rangle-|\tilde{\mu}^{k}\rangle\big{\|}\leq\epsilon_{3}^{\mu}, \big{\|}\|\mu^{k}\|-\|\tilde{\mu}^{k}\|\big{\|}\leq\epsilon_{4}^{\mu}\|\mu^{k}\|, \big{\|}|\Sigma^{k}\rangle-|\tilde{\Sigma}^{k}\rangle\big{\|}\leq\epsilon_{3}^{\mu}, and \big{\|}\|\Sigma^{k}\|-\|\tilde{\Sigma}^{k}\|\big{\|}\leq\epsilon_{4}^{\mu}\|\Sigma^{k}\|. In the next section, $\epsilon_{4}^{\pi}$, $\epsilon_{3}^{\mu}$, $\epsilon_{4}^{\mu}$, $\epsilon_{3}^{\Sigma}$, and $\epsilon_{4}^{\Sigma}$ will be used to analyze the runtime.

IV Analysis of errors and runtime

This section is dedicated to error and runtime analysis of the q-EM algorithm. We first state the main claim and then explain it.

IV.1 Main result

The runtime of the q-EM algorithm is represented by

[TABLE]

where $\mu(\cdot)$ is given in Eq. (72), $\kappa(\cdot)$ is the condition number, $\eta^{\mu}\coloneqq\max_{i}\|y_{i}\|^{2}$, and $\eta^{\Sigma}\coloneqq\max_{i}\|y_{i}\otimes y_{i}\|^{2}$. The definition of $\tilde{O}(\cdot)$ is given in Appendix A.

This result states that the runtime of each iteration of the q-EM algorithm is exponentially faster than that of the EM algorithm.

IV.2 Error analysis

We first summarize the errors in the q-EM algorithm to analyze the total runtime of the q-EM algorithm in the following subsection. In step I, we compute $d_{\mathrm{G}}^{k}(y_{i})$; the error on this computation is

[TABLE]

For consistency between the EM algorithm and the $\delta$-EM algorithm, we take $\epsilon_{1}<\delta/2$.

In steps III and IV, we compute $\mu$ and $\Sigma$. The errors on $\|\mu^{k}\|$ and $|\mu^{k}\rangle$ are $\sqrt{\eta^{\mu}}\epsilon_{3}^{\mu}$ and $\epsilon_{4}^{\mu}$, respectively. Then, the error on the estimation of $\mu$ takes the form

[TABLE]

See also Appendix C for the above calculation. We need to take $\epsilon_{3}^{\mu}<\frac{\delta}{4\sqrt{\eta^{\mu}}}$ and $\epsilon_{4}^{\mu}<\frac{\delta}{4\sqrt{\eta^{\mu}}}$. Next, we turn our attention to $\Sigma$. The errors on $\sqrt{\eta^{\Sigma}}\epsilon_{3}^{\Sigma}$ and $\sqrt{\eta^{\Sigma}}\epsilon_{4}^{\Sigma}$ are $\sqrt{\eta^{\Sigma}}\epsilon_{3}^{\Sigma}$ and $\epsilon_{4}^{\Sigma}$, respectively. Similarly to the case of $\mu^{k}$, the error on the estimation of $\Sigma$ is shown as

[TABLE]

We also need to take $\epsilon_{3}^{\Sigma}<\frac{\delta}{4\sqrt{\eta^{\Sigma}}}$ and $\epsilon_{4}^{\Sigma}<\frac{\delta}{4\sqrt{\eta^{\Sigma}}}$. Finally, we mention the error associated with the estimation on $\pi^{k}$ for $k=1,2,\dots,K$. The error on $\pi^{k}$ is shown as

[TABLE]

We estimate $\{\pi^{k}\}_{k=1}^{K}$ via the distribution of labels in quantum vector state tomography.

IV.3 Runtime

In the following, Eq. (34), i.e the runtime of each iteration of the q-EM algorithm is derived using Hoeffding’s inequality The required number of quantum vector state tomography of $K$ mean vectors is given as follows Kerenidis et al. (2018):

[TABLE]

Similarly, that for covariance matrices is

[TABLE]

The definition of $\tilde{O}(\cdot)$ is given in Appendix A. Next, we turn our attention to the runtime to prepare single copies of $|\mu^{k}\rangle$ and $|\mathrm{vec}[\Sigma^{k}]\rangle$. The time to prepare a copy of $|\mu_{t}^{k}\rangle$ is

[TABLE]

and that to prepare a copy of $|\mathrm{vec}[\Sigma_{t}^{k}]\rangle$ is

[TABLE]

where $\mu(V_{1})$ is given in Eq. (72), $\kappa(V_{1})$ is the condition number of $V_{1}$, $\eta^{\mu}\coloneqq\max_{i}\|y_{i}\|^{2}$ and $\eta^{\Sigma}\coloneqq\max_{i}\|y_{i}\otimes y_{i}\|^{2}$. Furthermore, $T_{\chi}^{\mu}$, which is the time to prepare $|\chi_{t}^{k}\rangle$ for estimating $\{\mu^{k}\}_{k=1}^{K}$, is given by

[TABLE]

Similarly, $T_{\chi}^{\Sigma}$ is given by

[TABLE]

In addition, $T_{\chi}^{\pi}$ is given by

[TABLE]

We also need to estimate the norms of $|\mu^{k}\rangle$ and $|\mathrm{vec}[\Sigma^{k}]\rangle$. The time for the norm estimation of $|\mu^{k}\rangle$ is

[TABLE]

and that of $|\mathrm{vec}[\Sigma^{k}]\rangle$ is

[TABLE]

We then estimate the runtime for estimating $\{\pi^{k}\}_{k=1}^{K}$. Due to Hoeffding’s inequality Hoeffding (1994), we need to perform sampling $\frac{2K}{(\epsilon_{4}^{\pi})^{2}}\ln\frac{2}{\Delta_{\mathrm{each}}^{\pi}}$ times to realize $\|\tilde{p}_{k}-p_{k}\|\leq\epsilon_{4}^{\pi}$ for $k=1,2,\dots,K$ with probability $(1-\Delta_{\mathrm{each}}^{\pi})^{K}$ since the distribution of $\pi$ is the $K$-state discrete distribution. By setting $1-\Delta^{\pi}\coloneqq(1-\Delta_{\mathrm{each}}^{\pi})^{K}$, we have $\Delta_{\mathrm{each}}^{\pi}=1-(1-\Delta^{\pi})^{1/K}$. Thus, we have

[TABLE]

Furthermore, we have to repeat the estimation process $K$ times for estimation of $\pi$ compared to those of $\mu$ and $\Sigma$, since we have to sample $i$ from $C^{k}$, $K$ times, for $k=1,2,\dots,K$. Thus, the runtime for estimating $\pi$ has an additional multiple of $K$.

Thus, the total runtimes for estimating $\{\pi^{k}\}$, $\{\mu^{k}\}$, and $\{\Sigma^{k}\}$ are, respectively,

[TABLE]

In total, we have obtained Eq. (34).

V Numerical simulation

To devise a quantum version of the EM algorithm, we proposed the $\delta$-EM algorithm in Sec. II. In this section, to see that the EM and $\delta$-EM algorithms are equivalent when $\delta$ is sufficiently small and that the $\delta$-EM algorithm improves upon the $\delta$-$k$-means algorithm, we show numerical simulations of the EM algorithm, the $\delta$-EM algorithm, the $k$-means algorithm, and the $\delta$-$k$-means algorithm.

In Ref. Cornell and Sastry (2015), the comparison between the $k$-means algorithm and the EM algorithm with GMM is shown. Then we use similar synthetic data sets used in Ref. Cornell and Sastry (2015). In the numerical simulations, we set $\delta=0.2$ except the numerical simulation for the $\delta$-dependence of the $\delta$-EM algorithm. For simplicity, we add Gaussian noise to the parameters estimated in the M step of the $\delta$-EM algorithm and the centroids estimated in the $\delta$-$k$-means algorithm 222In the the $\delta$-$k$-means algorithm, we added Gaussian noise whose mean is 0 and variance is 0.01 to each element of all the centroids. In the M step of the $\delta$-EM algorithm, we added Gaussian noise whose mean is 0 and variance is 0.01 to $\pi^{k}$ and each element of $\mu^{k}$, and add Gaussian noise whose mean is 0 and variance is 0.001 to each element of $\Sigma^{k}$. When $\pi$ is not normalized, we normalized $\pi$. When $\Sigma^{k}$ is not symmetric, we symmetrize $\Sigma^{k}$ by $(\Sigma^{k}+(\Sigma^{k})^{\intercal})/2$. When the estimated $\Sigma^{k}$ has a negative eigenvalue $\sigma^{*}$, we make all eigenvalues positive by adding $|\sigma^{*}|I^{k}$ where $I^{k}$ is the $k\times k$ identity matrix.

V.1 Example I

We begin with the explanation of the data set used in this subsection. We generated by drawing 1000 data points from the mixture of two Gaussian functions. The means of the two Gaussian functions are $\mu^{1}=[0.3,0.0]^{\intercal}$ and $\mu^{2}=[-0.3,0.0]^{\intercal}$, respectively, and the covariances are, respectively,

[TABLE]

We also put $\pi^{1}=\pi^{2}=0.5$.

In Figs. 1 and 2, we show the log-likelihood of the $k$-means algorithm and the $\delta$-$k$-means algorithm, and that of the EM algorithm and the $\delta$-EM algorithm, respectively. We plot ten trials for each algorithm. These figures show that the $k$-means and $\delta$-$k$-means algorithms have similar performance and that the EM and $\delta$-EM algorithms have similar performance.

For clarity, we graphically show the parameters estimated by the $\delta$-$k$-means algorithm and the $\delta$-EM algorithm in Fig. 3 333In Fig. 3, we draw the blue ellipses by moving $\theta$ $(0\leq\theta<2\pi)$ in $f(\theta)\coloneqq\mu^{k}+\sqrt{\sigma_{1}^{k}}e_{1}^{k}\cos(\theta)+\sqrt{\sigma_{2}^{k}}e_{2}^{k}\sin\theta$ where $\sigma_{1}^{k}$ and $\sigma_{2}^{k}$ are the eigenvalues of $\Sigma^{k}$, and $e_{1}^{k}$ and $e_{2}^{k}$ are the corresponding eigenvectors of $\Sigma^{k}$, respectively, for $k=1,2$. In Fig. 5, we also draw blue ellipses by the same procedure.. We have chosen the best estimates of the $\delta$-$k$-means algorithm and the $\delta$-EM algorithm in one hundred trials. These figures demonstrate that the $\delta$-EM algorithm outperforms the $\delta$-$k$-means algorithm.

In Table 1, we summarize the success rates of the EM algorithm, the $\delta$-EM algorithm, the $k$-means algorithm, and the $\delta$-$k$-means algorithm. Here, the success rate means that the ratio of the number of the successfully predicted hidden variables 444The hidden variable in the GMM means the index of the Gaussian functions. Note that the GMM has the symmetry on swapping the indices; then, we use the maximum value on swapping. to the number of the total data points. This table shows that the $\delta$-EM algorithm works better than the $\delta$-$k$-means algorithm. Thus, we insist that it is meaningful to devise a quantum version of the $\delta$-EM algorithm.

In Fig. 4, we show the $\delta$-dependence of the best success rates of the $\delta$-EM algorithm in one hundred trials. This figure shows that the $\delta$-EM algorithm is robust for small $\delta$, but the performance decreases rapidly for large values of $\delta$. We need to set $\delta$ small, since the critical value depends on data sets.

V.2 Example II

We again start with the data set used in this subsection. The data points are also generated by the mixture of two Gaussian functions, but the parameters are different. We set $[\pi^{1},\pi^{2}]=[0.7,0.3]$, $\mu^{1}=[0.0,-0.5]^{\intercal}$, $\mu^{2}=[0.0,0.0]^{\intercal}$, and

[TABLE]

Furthermore, we draw 1000 data points from the mixture of two Gaussian functions.

We first show the parameters estimated by the $\delta$-$k$-means algorithm and the $\delta$-EM algorithm in Fig. 5. We have chosen the best estimates of the $\delta$-$k$-means algorithm and the $\delta$-EM algorithm in one hundred trials. These figures represent that the $\delta$-EM algorithm outperforms the $\delta$-$k$-means algorithm. In particular, in the case of the $\delta$-$k$-means algorithm, the covariances are fixed at the identity matrix; then, each cluster tries to exclude each other.

In Table 2, we summarize the success rates of the EM algorithm, the $\delta$-EM algorithm, the $k$-means algorithm, and the $\delta$-$k$-means algorithm. Here, the success rate means that the ratio of the number of the successfully predicted labels to the number of the total data points. This table shows that the $\delta$-EM algorithm works better than the $\delta$-$k$-means algorithm.

VI Discussions

We here discuss the relationship between the EM algorithm with the GMM and the $k$-means algorithm. The EM algorithm with the GMM is an extension of the $k$-means algorithm; thus we explain the two conditions that the EM algorithm with the GMM becomes identical to the $k$-means algorithm.

The first condition is that $r_{t}^{i,k}$ takes $1$ for a certain $k$ and [math] otherwise. This implies that the $k$-means algorithm is an algorithm for hard clustering, while the EM algorithm is one for soft clustering. The second condition is that $\pi^{k}=1/K$ and $\Sigma^{k}=I_{d}$ where $I_{d}$ is the $d$-dimensional identity matrix for $k=1,2,\dots,K$. This is the reason why the $k$-means algorithm does not explicitly deal with weights and covariance matrices. From the viewpoint of a probability distribution, the $k$-means algorithm is an algorithm to estimate $\{\mu_{t}^{k}\}_{k=1}^{K}$ in

[TABLE]

As shown in Sec. V, the weights and the covariances of the GMM play an important role; thus, we also insist that the $\delta$-EM algorithm is a meaningful extension of the $\delta$-$k$-means algorithm.

VII Conclusion

In this paper, we have proposed a quantum algorithm for the EM algorithm and showed that it realize a quantum speedup compared to the classical EM algorithm. The key idea is to generalize the distance that is minimized in the $k$-means algorithm by considering also weights and covariances. Though we have focused on the GMM, we can generalize this condition to other mixture models. In machine learning, the EM algorithm with the GMM is more often used than the $k$-means algorithm; thus, this work is an important step toward quantum machine learning. The algorithm requires a QRAM oracle which has so far not been implemented in experiments yet. As a future direction we will investigate the applicability of novel superposition designs as proposed in Refs. Sieberer and Lechner (2018); Dlaska et al. (2019).

Acknowledgements

H.M. thanks Clemens Dlaska for fruitful discussions. H.M. is supported by JSPS KAKENHI Grant No. JP18J12175. The work is supported by the Austrian Science Fund (FWF) through a START grant under Project No. Y1067-N27 and the SFB BeyondC Project No. F7108-N38, the Hauser-Raspe foundation, and the European Union’s Horizon 2020 research and innovation program under grant agreement No. 817482 PasQuanS.

Appendix A Big $O$ notation

We here introduce the big $O$ notation, which is often used in computer science. For functions $f(x)$ and $g(x)$, one writes

[TABLE]

if and only if

[TABLE]

Similarly, we say

[TABLE]

if and only if

[TABLE]

These definitions will be utilized to describe the q-EM algorithm and to perform error analysis.

Appendix B Quantum preliminaries

We provide some tools that are required for the $q$-means algorithm Kerenidis et al. (2018) in this section. These tools are also utilized in a quantum algorithm for the EM algorithm.

B.1 Amplitude estimation

Here we sumarize the amplitude estimation algorithm that was proposed in Ref. Brassard et al. (2002). Assume that we have $U_{A}$ given by

[TABLE]

Then, there exists an amplitude estimation algorithm that outputs $\tilde{p}$ such that

[TABLE]

with probability at least $8/\pi^{2}$. The algorithm perform $U_{A}$ $P_{\mathrm{ae}}$ times. Note that, if $p=0$, $\tilde{p}=0$, and if $p=0$ and $P_{\mathrm{ae}}$ is even, then $\tilde{p}=0$.

Furthermore, to raise the probability to obtain a good estimate on distances, we utilize a tool in Ref. Wiebe et al. (2014). We make multiple copies of the amplitude estimates, apply the quantum mode evaluation algorithm proposed in Lemma 8 of Ref. Wiebe et al. (2014) in Sec. B.2, and reverse the circuit to remove the garbage state. We note, that very recently an amplitude estimation algorithm without phase estimation was introduced Suzuki et al. (2019).

B.2 Median evaluation

The time complexity of the mode evaluation algorithm is given in Lemma 8 of Ref. Wiebe et al. (2014). Let us summarize the main idea of this Lemma. Let $U$ be a unitary operation given by

[TABLE]

for $1/2<a\leq 1$ in time $T$. Then there exists a quantum algorithm that produces a state $\Ket{\Phi}$ such that

[TABLE]

for $\Delta>0$, $1/2<a_{0}<a$ and integer $L$ in time 2T\Big{\lceil}\frac{\ln(1/\Delta)}{2(|a_{0}|-1/2)^{2}}\Big{\rceil}.

B.3 Quantum random access memory

In the $q$-means and $q$-EM algorithm, it is crucial to prepare data as a quantum state efficiently. To this end, we exploit the quantum random access memory (QRAM) introduced in Refs.  Giovannetti et al. (2008a, b). Here, the authors consider a device that performs the operation

[TABLE]

We follow the application of the QRAM as in Ref. Kerenidis et al. (2018). Let $V_{1}\in\mathbb{R}^{N\times d}$; then, there is a data structure to store the rows of $V_{1}$ such that the time to insert, update, or delete a single entry $v_{i,j}$ is $O(\ln^{2}N)$ and a quantum algorithm on the data structure can be performed in time $O(\ln^{2}N)$ that realizes the following unitaries:

[TABLE]

B.4 Quantum linear algebra

Some useful subroutines that are used in $q$-means and $q$-EM are given as follows:

Theorem 1**.**

Let $M\in\mathbb{R}^{d\times d}$ that satisfies $\|M\|_{2}=1$ and $x\in\mathbb{R}^{d}$. If $M$ is stored in QRAM and the time to prepare $\Ket{x}$ is $T_{x}$, then there exist quantum algorithms that return

• •

a state $\Ket{z}$ such that $\|\Ket{z}-\Ket{Mx}\|\leq\epsilon$ in time \tilde{O}\Big{(}(\kappa(M)\mu(M)+T_{x}\kappa(M))\ln(\epsilon^{-1})\Big{)},

• •

a state $\Ket{z}$ such that $\|\Ket{z}-\Ket{M^{-1}x}\|\leq\epsilon$ in time \tilde{O}\Big{(}(\kappa(M)\mu(M)+T_{x}\kappa(M))\ln(\epsilon^{-1})\Big{)},

• •

the norm $z\in(1+\delta)\|Mx\|$ with relative error $\delta$ in time \tilde{O}\Big{(}T_{x}\kappa(M)\mu(M)\delta^{-1}\ln(\epsilon^{-1})\Big{)}.

where

[TABLE]

and $\kappa(M)$ is the condition number of $M$.

Appendix C Inequality for error analysis

The following holds for general $\vec{a}$ and $\vec{b}$:

[TABLE]

We used the above equation to derive Eq. (36).