Message-passing algorithm of quantum annealing with nonstoquastic Hamiltonian

TL;DR

This paper develops a message-passing algorithm to simulate quantum annealing with nonstoquastic Hamiltonians, enabling performance assessment without quantum Monte Carlo methods, and validates it against the replica method.

Contribution

It introduces a novel message-passing algorithm for classical simulation of nonstoquastic Hamiltonians in quantum annealing, expanding the tools for analyzing complex quantum systems.

Findings

The algorithm accurately assesses QA performance with nonstoquastic Hamiltonians.

Classical simulation of certain nonstoquastic Hamiltonians is feasible at low temperatures.

Validation shows consistency with the replica method results.

Abstract

Quantum annealing (QA) is a generic method for solving optimization problems using fictitious quantum fluctuation. The current device performing QA involves controlling the transverse field; it is classically simulatable by using the standard technique for mapping the quantum spin systems to the classical ones. In this sense, the current system for QA is not powerful despite utilizing quantum fluctuation. Hence, we developed a system with a time-dependent Hamiltonian consisting of a combination of the formulated Ising model and the "driver" Hamiltonian with only quantum fluctuation. In the previous study, for a fully connected spin model, quantum fluctuation can be addressed in a relatively simple way. We proved that the fully connected antiferromagnetic interaction can be transformed into a fluctuating transverse field and is thus classically simulatable at sufficiently low temperatures. Using the fluctuating transverse field, we established several ways to simulate part of the nonstoquastic Hamiltonian on classical computers. We formulated a message-passing algorithm in the present study. This algorithm is capable of assessing the performance of QA with part of the nonstoquastic Hamiltonian having a large number of spins. In other words, we developed a different approach for simulating the nonstoquastic Hamiltonian without using the quantum Monte Carlo technique. Our results were validated by comparison to the results obtained by the replica method.

Equations

$${\hat{H}}_0=f({\hat{\mathbf{\sigma }}}^z)$$

$${\hat{H}}_1=-Ng\left(\frac{1}{N}\sum _{i=1}^N{\hat{\sigma }}_i^x\right).$$

$$Z=\mathrm{Tr}\left\{\exp \left(-\beta {\hat{H}}_0-\beta {\hat{H}}_1\right)\right\}$$

$$Z=\underset{M\to \infty }{\lim }\mathrm{Tr}\left\{\prod _{k=1}^M\exp \left(-\frac{\beta }{M}f({\hat{\mathbf{\sigma }}}_k^z)\right)\exp \left(\frac{N\beta }{M}g\left(\frac{1}{N}\sum _{i=1}^N{\hat{\sigma }}_{i,k}^x\right)\right)\right\}$$

$$1=\int d{m}_{x,k}\delta (N{m}_{x,k}-\sum _{i=1}^N{\sigma }_{i,k}^x)=\int d{\tilde{m}}_{x,k}\int d{m}_{x,k}\exp \left\{-\frac{\beta {\tilde{m}}_{x,k}}{M}\left(N{m}_{x,k}-\sum _{i=1}^N{\sigma }_{i,k}^x\right)\right\}$$

$$Z\approx \underset{M\to \infty }{\lim }\mathrm{Tr}\left\{\prod _{k=1}^M\int d{\tilde{m}}_{x,k}\int d{m}_{x,k}\exp \left(\frac{N\beta }{M}g\left(m_{x,k}\right)-\frac{N\beta }{M}{\tilde{m}}_{x,k}{m}_{x,k}\right)\exp \left(-\frac{\beta }{M}f({\hat{\mathbf{\sigma }}}_k^z)\right)\exp \left(\frac{\beta {\tilde{m}}_{x,k}}{M}\sum _{i=1}^N{\hat{\sigma }}_{i,k}^x\right)\right\}.$$

$$Z\approx \mathrm{Tr}\left\{\int d{\tilde{m}}_x\int d{m}_x\exp \left(N\beta g\left(m_x\right)-N\beta {\tilde{m}}_x{m}_x\right)\exp \left(-\beta f({\hat{\mathbf{\sigma }}}^z)+\beta {\tilde{m}}_x\sum _{i=1}^N{\hat{\sigma }}_i^x\right)\right\}.$$

$${\hat{H}}_{\mathrm{eff}}=f({\hat{\mathbf{\sigma }}}^z)-{\tilde{m}}_x\sum _{i=1}^N{\hat{\sigma }}_i^x.$$

$$P({\tilde{m}}_x)\propto \mathrm{Tr}\left\{\int d{m}_x\exp \left(N\beta g\left(m_x\right)-N\beta {\tilde{m}}_x{m}_x\right)\exp \left(-\beta f({\hat{\mathbf{\sigma }}}^z)+\beta {\tilde{m}}_x\sum _{i=1}^N{\hat{\sigma }}_i^x\right)\right\}.$$

$$m_x$$

$${\tilde{m}}_x$$

$$⟨{\mathbf{\sigma }}_{k+1}^z|\exp \left(\frac{\beta {\tilde{m}}_{x,k}}{M}{\hat{\sigma }}_{i,k}^x\right)|{\mathbf{\sigma }}_k^z⟩\propto ⟨{\mathbf{\sigma }}_{k+1}^z|\exp \left(\beta J_k^{\ast }\sum _{i=1}^N{\sigma }_{i,k}^z{\sigma }_{i,k+1}^z\right)|{\mathbf{\sigma }}_k^z⟩$$

$$m_{x,k}$$

$$f(\mathbf{\sigma })=-\sum _{\mu }{J}_{\mu }\prod _{l\in \partial \mu }{\sigma }_l-\sum _{i=1}{h}_i{\sigma }_i.$$

$$M_{\mu \to i}({\sigma }_i)\propto \sum _{\mathbf{\sigma }/{\sigma }_i}\exp \left(\beta J_{\mu }\prod _{l\in \partial \mu }{\sigma }_l\right)\prod _{l\in \partial \mu /i}{M}_{l\to \mu }({\sigma }_l)$$

$$M_{i\to \mu }({\sigma }_i)\propto \exp \left(\beta h_i{\sigma }_i\right)\prod _{\nu \in \partial i/\mu }{M}_{\nu \to i}({\sigma }_i),$$

$$Q_{\mu }({\mathbf{\sigma }}_{\partial \mu })\propto \sum _{\mathbf{\sigma }}\exp \left(\beta J_{\mu }\prod _{i\in \partial \mu }{\sigma }_{\mu }\right)\prod _{i\in \partial \mu }{M}_{i\to \mu }({\sigma }_i)$$

$$Q_i({\sigma }_i)\propto \exp \left(\beta h_i{\sigma }_i\right)\prod _{\mu \in \partial i}{M}_{\mu \to i}({\sigma }_i).$$

$$P(\mathbf{\sigma })\approx \prod _{\mu }\frac{Q_{\mu }({\mathbf{\sigma }}_{\partial \mu })}{{\prod }_{i\in \partial \mu }{Q}_i({\sigma }_i)}\prod _{i=1}^N{Q}_i({\sigma }_i).$$

$$M_{\mu \to i}({\sigma }_i)$$

$$M_{i\to \mu }({\sigma }_i)$$

$${\tilde{m}}_{i\to \mu }=h_i+\sum _{\nu \in \partial i/\mu }{\tilde{m}}_{\nu \to i}$$

$$M_{\mu \to i}({\sigma }_i)\propto \cosh \left(\beta J_{\mu }\right)+{\sigma }_i\prod _{l\in \partial \mu }\left(\sum _{{\sigma }_l}{M}_{l\to \mu }({\sigma }_l){\sigma }_l\right)\sinh \left(\beta J_{\mu }\right).$$

$$\tanh \beta {\tilde{m}}_{\mu \to i}=\tanh (\beta J_{\mu })\prod _{l\in \partial \mu /i}{m}_{l\to \mu }.$$

$$m_{l\to \mu }\equiv \frac{1}{Z_{l\to \mu }}\sum _{{\sigma }_l}{\sigma }_l\exp \left(\beta {\tilde{m}}_{l\to \mu }{\sigma }_l\right)=\tanh \left(\beta {\tilde{m}}_{l\to \mu }\right).$$

$$m_{i\to \mu }$$

$${\tilde{m}}_{i\to \mu }$$

$$m_i$$

$${\tilde{m}}_i$$

$$M_{i\to \mu }({\sigma }_i)\propto \exp \left(\beta {\tilde{m}}_{i\to \mu }{\sigma }_i\right)\to {\hat{M}}_{i\to \mu }\propto \exp \left(\beta {\tilde{m}}_{i\to \mu }{\hat{\sigma }}_i^z+\beta {\tilde{m}}_i^x{\hat{\sigma }}_i^x\right).$$

$$M_{i\to \mu }(\mathbf{\sigma })\propto \exp \left(\beta {\tilde{m}}_{i\to \mu }\sum _{k=1}^M{\sigma }_{i,k}+\frac{\beta }{2}{\tilde{r}}_{i\to \mu }{\left(\sum _{k=1}^M{\sigma }_{i,k}\right)}^2+\beta J_i^{\ast }{\sigma }_{i,k}{\sigma }_{i,k+1}\right).$$

$$M_{i\to \mu }(\mathbf{\sigma })=\frac{1}{Y_{i\to \mu }\left({\tilde{m}}_{i\to \mu },{\tilde{r}}_{i\to \mu },{\tilde{m}}_i^x\right)}\int Dy\prod _{k=1}^M\exp \left(\beta \left({\tilde{m}}_{i\to \mu }+\sqrt{{\tilde{r}}_{i\to \mu }/\beta }y\right){\sigma }_{i,k}+\beta J_i^{\ast }{\sigma }_{i,k}{\sigma }_{i,k+1}\right),$$

$$Y_{i\to \mu }\left({\tilde{m}}_{i\to \mu },{\tilde{r}}_{i\to \mu },{\tilde{m}}_i^x\right)=\sum _{\mathbf{\sigma }}\int Dy\prod _{k=1}^M\exp \left(\beta \left({\tilde{m}}_{i\to \mu }+\sqrt{{\tilde{r}}_{i\to \mu }/\beta }y\right){\sigma }_{i,k}+\beta J_i^{\ast }{\sigma }_{i,k}{\sigma }_{i,k+1}\right),$$

Full text

1 Introduction

Quantum annealing (QA) is a generic algorithm aimed at solving optimization problems by exploiting the quantum tunneling effect. The scheme was originally proposed as an algorithm for numerical computation [1] inspired by simulated annealing (SA) [2] and exchange Monte Carlo simulation [3]. However, its experimental realization was accomplished recently, attracting considerable attention. QA is advantageous for solving an optimization problem formulated with discrete variables using a platform, the Ising model. In QA, we developed a system with a time-dependent Hamiltonian consisting of a combination of the formulated Ising model and the “driver” Hamiltonian with only quantum fluctuation. The often-used driver Hamiltonian is the transverse field, which generates the superposition of the up and down spins. The first stage of QA is initialized in the trivial ground state of the driver Hamiltonian. The quantum effect will be gradually turned off and eventually ends so that only the formulated Ising model with the nontrivial ground state remains. When the transverse field changes sufficiently slowly, the quantum adiabatic theorem ensures that we find the nontrivial ground state at the end of QA [4, 5, 6]. Numerous reports state that QA outperforms SA [7, 8, 9]. Possibly, the performance of QA stems from the quantum tunneling effect penetrating the valley of the potential energy. The protocol of QA is realized in an actual quantum device using current technology, namely, quantum annealer [10, 11, 12, 13]. The output solution from the current system comprising quantum annealer is not always optimal due to device limitations and environmental effects [14]. Therefore, several protocols based on QA do not follow several conditions in adiabatic quantum computation or maintain a system in the ground state to reach the optimal solution in the final stage of QA; rather, they employ a nonadiabatic counterpart [15, 16, 17, 18]. Nevertheless, quantum annealer has been tested for numerous applications such as portfolio optimization [19], protein folding [20], the molecular similarity problem [21], computational biology [22], job-shop scheduling [23], traffic optimization [24], election forecasting [25], machine learning [26, 27, 28], and automated guided vehicles in plants [29]. In addition, studies on solving various problems by using quantum annealer have been performed [30, 31, 32, 33]. The potential of QA might be boosted by the nontrivial quantum fluctuation, referred to as the nonstoquastic Hamiltonian, for which efficient classical simulation is intractable [34, 35, 36, 37].

The efficiency of QA is characterized by the energy gap between the ground state and the first excited state in the intermediate Hamiltonian. The adiabatic theorem formulates the relationship between the energy gap $\Delta$ and the necessary computation time $\tau_{\rm QA}$ for QA as $\tau_{\rm QA}\sim 1/\Delta^{2}$ [4]. The necessary computation time can be longer when the size of the problems, which is characterized by the number of spins $N$, increases because the energy gap decreases. When the system is involved in phase transition, the energy gap collapses exponentially or polynomially. The former is a first-order phase transition, and the latter is a second-order one. Correspondingly, the order of the phase transition describes the performance of QA. In particular, the existence of the first-order phase transition during hampers efficient manipulation of QA to find the nontrivial ground state at the end. In this sense, avoiding the first-order phase transition is one of the central issues in QA.

One of the possible ways to avoid the first-order phase transition is to introduce the antiferromagnetic XX interaction into the conventional time-dependent Hamiltonian in QA [34, 35, 38]. In particular, the quantum fluctuation of the antiferromagnetic XX interaction is useful for eliminating the sudden change in the ferromagnetic ordering of the Ising spins. Thus, the first-order phase transition can be modified into the second-order one. One of the nontrivial properties of the antiferromagnetic XX interaction is the emergence of the negative sign in the standard way to map the quantum spin systems with the effective classical spin systems. For instance, the Suzuki–Trotter decomposition [39] yields a negative sign for the antiferromagnetic XX interaction in a straightforward way. In this sense, the Ising spin models with antiferromagnetic XX interactions are nonstoquastic [40]. The nonstoquastic Hamiltonian is expected to enhance the power of quantum computation because it is capable of addressing intractable quantum fluctuation on a classical computer. Thus, the antiferromagnetic XX interaction is also expected to be with certain potential for quantum speed up. However, beyond modifying ferromagnetic ordering, the advantage of antiferromagnetic XX interactions is not yet clarified.

In the previous study [36], they proved that the fully connected antiferromagnetic interaction can be transformed into a fluctuating transverse field and is thus classically simulatable at sufficiently low temperatures; this can be implemented in a system with a large number of spins using the quantum Monte Carlo simulation. In the present study, we formulated the message-passing algorithm and its approximated version. The resultant algorithm is capable of assessing the performance of QA with part of the nonstoquastic Hamiltonian having a large number of spins.

The remainder of the paper is organized as follows: In the next section, we discuss the transformation of the antiferromagnetic XX interactions into a fluctuating transverse field and building of the quantum Monte Carlo simulation for part of the nonstoquastic Hamiltonian. In the following section, we describe the message-passing algorithm for the effective model including the antiferromagnetic XX interactions. We then discuss the asymptotic behavior for a large number of spins. In the conclusion section, we summarize our findings and discuss future research directions.

2 Transformation of Non-stoquastic Hamiltonian

Let us consider the following quantum system for solving the optimization problem using QA. The target Hamiltonian, which represents the optimization problem to be solved, can be written as

[TABLE]

where $\hat{{\boldsymbol{\sigma}}}^{z}=(\hat{\sigma}_{1}^{z},\hat{\sigma}_{2}^{z},\cdots,\hat{\sigma}_{N}^{z})$ and $\hat{\sigma}_{i}^{z}$ is the $z$-component of the Pauli matrices. The driver Hamiltonian is chosen with the following special form as

[TABLE]

The function $g$ includes various types of quantum fluctuations. For instance, the transverse field is the case with $g(x)=\Gamma x$ and the antiferromagnetic XX interaction, which is a typical example of the nonstoquastic Hamiltonian; that is, $g(x)=-\gamma x^{2}/2$. Using the following decomposition, we can simulate part of the nonstoquastic Hamiltonian without facing the sign problem. Let us consider the composite system comprising the target and driver Hamiltonians by computing the partition function.

[TABLE]

Here, we employed the Suzuki–Trotter decomposition to separate the noncommutable operators in the exponential.

[TABLE]

We may regard each operator in the Hamiltonian as the $c$ numbers after inserting the completeness relation $1=\sum_{{\boldsymbol{\sigma}}^{z}_{k}}|{\boldsymbol{\sigma}}^{z}_{k}\rangle\langle{\boldsymbol{\sigma}}^{z}_{k}|$ and $1=\sum_{{\boldsymbol{\sigma}}^{x}_{k}}|{\boldsymbol{\sigma}}^{x}_{k}\rangle\langle{\boldsymbol{\sigma}}^{x}_{k}|$. In addition, we used the following identity and changed the delta function in terms of the Fourier integration as

[TABLE]

After again considering the c numbers as the operators, the resultant expression of the partition function can be written as

[TABLE]

Here, we utilized the inverse relation of the Suzuki–Trotter decomposition and obtained

[TABLE]

This expression led to the interpretation of the original composite system comprising the target and driver Hamiltonians as the following simple system with a fluctuating transverse field:

[TABLE]

The strength of the transverse field is governed by the probability distribution as follows:

[TABLE]

The resulting expression states that the nonstoquastic term $g(\sum_{i=1}^{N}\hat{\sigma}^{x}_{i}/N)$ can be reduced to the transverse-field term. This fact implies that the present form of the nonstoquastic term does not represent strong quantum effects as the interference of the wave function. However, the same type of nonstoquastic term changes the first-order phase transition, which is harmful to the second-order one and can be easily solved using QA [34].

In QA, we were interested in extremely low-temperature regions as $\beta\to\infty$. Therefore, we considered the saddle point. The saddle point equations are given as

[TABLE]

where $\hat{\rho}_{\rm eff}=\exp(-\beta\hat{H}_{\rm eff})/Z(\tilde{m}_{x})$ and $Z(\tilde{m}_{x})={\rm Tr}\exp\left(-\beta\hat{H}_{\rm eff}\right)$. Instead of the saddle-point equation, one may utilize the Langevin dynamics to generate the probability distribution. However, it might take a long time to equilibrate the dynamics. One would need to accelerate the equilibration as demonstrated in some studies [41, 42, 43]. After determining the transverse magnetization $m_{x}$ for a given value of $\tilde{m}_{x}$, we updated the value of $\tilde{m}_{x}$, which is the strength of the effective transverse magnetic field. Iterative manipulation of the above procedure yields the equilibrium state of the system. We can then compute the physical quantity in the equilibrium state. In the current system implementing QA in the physical device using D-Wave 2000Q, the transverse field can be tuned by following the predetermined schedule. The schedule allows for pausing the control of the transverse field in the intermediate period of the QA protocol. We can estimate the physical quantities in the equilibrium state with the transverse field, as explained in a previous study [44]. Nevertheless, unfortunately, the current system is not capable of estimating the transverse magnetization directly. Therefore, we did not directly perform the above procedure using D-Wave 2000Q. Instead of the quantum device, we may employ a classical simulation to estimate the transverse magnetization, similar to a previous study [36]. They utilized the Suzuki–Trotter decomposition and employed the following relationship:

[TABLE]

where $\exp\left(-2\beta J^{*}\right)=\tanh\beta\tilde{m}_{x,k}/M$. In other words, the original quantum system can be mapped onto a classical effective Hamiltonian with spin variables $\sigma_{i,k}^{z}=\pm 1$ and continuous variables $\tilde{m}_{x,k}$ and $m_{x,k}$. Instead of using Eq. (10), we estimated the transverse magnetization using the following relationship, which was straightforwardly obtained by taking the derivative with respect to $\tilde{m}_{x,k}$,

[TABLE]

We should emphasize that the target Hamiltonian is not restricted to the case with all-to-all connections. The target Hamiltonian can only be used with short-range interactions such as those in the finite-dimensional Edwards–Anderson model.

In the present study, we formulated the belief-propagation method to reduce the computational cost [45, 46]. In addition, we derived an approximate message-passing algorithm, which led to results consistent with the analytical solution at the level of replica symmetry.

3 Message-passing algorithm

We hereafter assume that the typical spin-glass model takes the following form:

[TABLE]

where $\mu$ denotes the factor nodes to which the interactions between spins are assigned, and $k$ denotes variable nodes to which the spin variables are allocated. This assumption is reasonable for the standard definition of spin-glass models.

In the classical case without the transverse field, we iteratively computed the following update equation:

[TABLE]

where $M_{\mu\to i}(\sigma_{i})$ and $M_{i\to\mu}(\sigma_{i})$ are the messages passing through the factor and variable nodes, respectively. The notation ${\boldsymbol{\sigma}}/\sigma_{i}$ denotes that the summation over the variables except for $\sigma_{i}$. These messages form the approximation of the distribution function through

[TABLE]

These ingredients approximate the distribution function according to the following relationship:

[TABLE]

This approximation can be obtained by minimization of the Kullback–Liebler divergence under constraints such as $\sum_{{\boldsymbol{\sigma}}_{\partial\mu/i}}Q_{\mu}({\boldsymbol{\sigma}}_{\partial\mu})=Q_{i}(\sigma_{i})$ and normalization. The optimal solution satisfies the previous equations for messages (15) and (16).

Let the message be parametrized as

[TABLE]

These parametrizations are validated by the fact that the message is a function of the local Ising spin variable. We obtained the first update equation as

[TABLE]

The second update equation was determined by taking the summation over the Ising variable. The message was assessed as follows:

[TABLE]

Thus, we obtained

[TABLE]

We define the following quantity for convenience:

[TABLE]

where $Z_{l\to\mu}=\sum_{\sigma_{l}}\exp\left(\beta\tilde{m}_{l\to\mu}\sigma_{l}\right)$. In other words, $m_{l\to\mu}$ is the local expectation of the Ising spin variable. Using these messages, we obtained the following relationships:

[TABLE]

These equations are well-known results for the belief propagation based on the Ising spin-glass model. For the local magnetization, we obtained

[TABLE]

Considering the above calculations, we found that the essential calculation was the expectation of the Ising spin variable following the local distribution characterized by $\tilde{m}_{l\to\mu}$. We considered generalization of the local distribution to include the quantum fluctuations. One of the simplest approaches is changing the local distribution $M_{i\to\mu}(\sigma_{i})$ into the density matrix following the prescription of quantum mechanics as shown below:

[TABLE]

We performed the Suzuki–Trotter decomposition of this density matrix and further generalized the message as long as we could construct the closed update equation as follows:

[TABLE]

This is rewritten using the Hubbard–Stratonovich transformation as follows:

[TABLE]

where

[TABLE]

and $\int Dy=\int^{\infty}_{-\infty}dy\exp\left(-y^{2}/2\right)/\sqrt{2\pi}$. Using the above form of the message, we did not always obtain the closed update equations. When we omit higher-order terms other than $\beta J_{\mu}^{2}$, we obtain the following closed update equations. For the longitudinal magnetization, we have Eq. (27) and

[TABLE]

In addition, we determined another set of update equations as follows:

[TABLE]

We then estimated both the longitudinal and transverse local magnetization to be

[TABLE]

where

[TABLE]

We iterated the update equations for the messages in Eqs. (27), (34), (35), and (36) and let them converge to fixed-point solutions. For $\tilde{m}_{i}^{x}$, we utilized the saddle-point equations (11). Then, we assessed $Y_{i\to\mu}\left(\tilde{m}_{i\to\mu},\tilde{r}_{i\to\mu},\tilde{m}_{i}^{x}\right)$. It is possible to straightforwardly manipulate the calculation as

[TABLE]

The derived update equations allow us to elucidate the nontrivial results at the level of static approximation as shown below. Notice that it is just approximation in terms of the belief propagation at the level of the Bethe approximation in terms of statistical mechanics. We may use the above update equation for an arbitrary shape of the graph as in the loopy belief propagation, which is sometimes applied in signal processing and image restoration. Similar to the case of the loopy belief propagation, we may encounter some nonconvergent behaviors depending on the problem.

We might generalize the message in a different form. For instance, one of the candidates is

[TABLE]

This is the same computation for assessing the normalization constant as that for the partition function of the one-dimensional Ising model under the “random field” defined by $\tilde{m}_{i\to\mu,k}$. Using the transfer matrix method, we can compute the normalization constant [47]. This approach corresponds to the path-integral representation in the quantum cavity method according to literature. In this sense, the form of the message is essentially important to characterize the level of the approximation of the approach. To obtain the nontrivial results beyond the static approximation, we have to tackle a generic form of the normalization constant [48]. In the present study, we analyzed the effect of the nonstoquastic Hamiltonian at the level of the static approximation as shown below. Several terms, $\tilde{m}_{i\to\mu}$ and $\tilde{r}_{i\to\mu}$, consist of the summation of the random variables. When we consider $|\partial i|\sim N\to\infty$, namely the fully connected models, according to the central limit theorem, $\tilde{m}_{i\to\mu}$ and $\tilde{r}_{i\to\mu}$ can be characterized by their mean and covariance. As discussed below, the above message-passing algorithm for the fully connected models can be validated at the level of the replica symmetric and static approximation.

4 Validation of our results

4.1 p-body spin-glass model

We assumed that the interaction term covering all the combinations of the Ising variables and their strength scales as $J_{\mu}\sim O(1/N^{(p-1)/2})$, where $p$ is the number of the Ising spin variables involved in the factor $\mu$. Then, $\tanh(\beta J_{\mu})\approx\beta J_{\mu}$ and $\tanh^{-1}(\beta J_{\mu}\prod_{l\in\partial\nu/i}m_{l\to\nu})\approx\beta J_{\mu}\prod_{l\in\partial\nu/i}m_{l\to\nu}$. In a large-sized limit such as $N\to\infty$, the self-averaging property assures that $J_{\mu}$ and $h_{i}$ are random variables following their distribution function. The update equation of $\tilde{m}_{i\to\mu}$ for the case without the transverse field is reduced to

[TABLE]

Similarly, $\tilde{m}_{i}$ is given by

[TABLE]

The local magnetization and message have the following relationship:

[TABLE]

The local magnetization $m_{i}=\tanh(\beta m_{i})$ is thus given by, though Eq. (28),

[TABLE]

This is the Thouless–Anderson—Palmer (TAP) equation [49, 50]. For the quenched interaction and field, the local magnetization can fluctuate over the site $i$. We then estimated the local magnetization through the TAP equation. On the other hand, the dynamics of the expectation value of the local magnetization over the site dependence can be assessed using the central limit theorem. We computed the expectation and variance of $\tilde{m}_{i}$ as

[TABLE]

where ${\mathbb{E}}\left[\cdots\right]$ is the expectation on the probability distribution of $J_{\mu}$ with a mean equal to $J_{0}/N^{p-1}$ and variance equal to $\sqrt{J^{2}/N^{p-1}}$; we define the magnetization as $m=\sum_{i=1}^{N}m_{i}/N$ and the spin-glass order parameter as $q=\sum_{i=1}^{N}m_{i}^{2}/N$. Here, for simplicity, we set the magnetic field as constant. Then, $\tilde{m}_{i}$ becomes

[TABLE]

where $z_{i}$ is the random variable with zero mean and unit variance. Therefore, the magnetization can be evaluated as follows:

[TABLE]

Similar to the magnetization, we obtained the spin-glass order parameter as follows:

[TABLE]

These results are consistent with the result of the replica method for the $p$-body spin-glass model under the replica symmetric assumption.

Next, we discuss the quantum system. We again computed the expectations and variances of $\tilde{m}_{i}$ and $\tilde{r}_{i}$ as

[TABLE]

where we define the magnetization as $m=\sum_{i=1}^{N}m_{i}/N$, the spin-glass order parameter as $q=\sum_{i=1}^{N}m_{i}^{2}/N$, and an overlap between different imaginary times as $R=\sum_{i=1}^{N}r_{i}/N$. Then, we obtain

[TABLE]

We set $a(y,z)\equiv h+J_{0}pm^{p-1}+\sqrt{pJ^{2}q^{p-1}}z+\sqrt{pJ^{2}(R^{p-1}-q^{p-1})}y$. Using these results, we can reproduce the update equations, which are consistent with the saddle-point equations for the $p$-body spin-glass model under a uniform transverse field $\tilde{m}_{i}^{x}=\tilde{m}_{x}$ [51] as follows:

[TABLE]

where we omit the arguments of the function $Y_{i\to\mu}$ to lighten the writing. Instead of the quantum Monte Carlo simulation, we employed the above update equations for estimating the order parameters conditioned on $\tilde{m}^{x}$. The effective magnetic field changes according to the saddle-point equation (11). Iterative manipulation of these procedures yields the results for the equilibrium state of the fully connected $p$-body Ising model with antiferromagnetic XX interactions. The results are generalizations of those obtained in previous studies [34, 35]. In other words, we developed a different approach for the nonstoquastic Hamiltonian without using the quantum Monte Carlo simulations and the replica method.

4.2 Application to restricted Boltzmann machine

We introduce an application of our message passing algorithm for a part of the non-stoquastic Hamiltonian. We take the restricted Boltzmann machine, whose Hamiltonian is given as

[TABLE]

where $v_{i}\in{-1,1}$ and $h_{j}\in{-1,1}$ are the visible and hidden variables, respectively. We regard these variables as the $z$-components of the Pauli matrices. We then apply the transverse field and antiferromagnetic XX interactions for each variable. We construct the message passing algorithm by defining two sets of the messages for visible and hidden variables.

[TABLE]

where $\mu=(ij)$, $M_{i\to\mu}$, and $M_{j\to\mu}$ similarly to Eq. (16). The direct manipulation of the above calculation yields the parameters of the messages as

[TABLE]

for the visible variables and similarly for the hidden variables by replacement of $i\leftrightarrow j$ and $b^{v}_{i}\leftrightarrow b^{h}_{j}$. We find that $m_{i\to\mu}$, $m_{j\to\mu}$, $r_{i\to\mu}$, and $r_{j\to\mu}$ follow the same forms as Eqs. (34) and (36). Both the longitudinal and transverse local magnetization can be obtained by Eqs. (37), (38) and (39). The numbers of the visible and hidden variables are $N_{v}$ and $N_{h}$. We assume $N_{v}$ and $N_{v}$ go infinity while $\alpha_{v}=N_{v}/(N_{v}+N_{h})$ and $\alpha_{h}=N_{h}/(N_{v}+N_{h})$ are finite. Then we obtain, by the similar way to the previous case,

[TABLE]

for the visible variables and similarly for the hidden variables by replacement of $i\leftrightarrow j$ and $b^{v}_{i}\leftrightarrow b^{h}_{j}$. The restricted Boltzmann machine generates many samples in order to calculate the expectation value and variance in learning. However our message-passing algorithm make sampling easier and faster because we omit relaxation time to generate the Gibbs-Boltzmann sampling, which needs relatively longer time. In order to validate our result, we assume that $\mathbb{E}[W_{\mu}]=w/\sqrt{N_{v}N_{h}}$ and $\mathbb{E}[W_{\mu}^{2}]=W^{2}/\sqrt{N_{v}N_{h}}$ and no biases exist for simplicity. The expectations and variances of $\tilde{m}_{i}$ and $\tilde{m}_{j}$ are

[TABLE]

for the visible variables and similarly for the hidden variables by replacement of $i\leftrightarrow j$ and $(b_{v},B_{v},m_{v},r_{v},q_{v})\leftrightarrow(b_{h},B_{h},m_{h},r_{h},q_{h})$. Then $\tilde{m}_{i}$ is expressed by, according to the central limit theorem,

[TABLE]

Replacement of $a(y,z)=\sqrt{\alpha_{v}/\alpha_{h}}wm_{h}+\left(\alpha_{v}/\alpha_{h}\right)^{1/4}W\sqrt{q_{h}}z+\left(\alpha_{v}/\alpha_{h}\right)^{1/4}W\sqrt{R_{h}-q_{h}}y$ in Eqs. (56), (57), (58) and (59) leads to a natural generalization of the results of the previous study [52]. We emphasize that our results include the quantum fluctuation expressed by the transverse field and beyond. In the restricted Boltzmann machine, we need many samples to estimate the expectation and variance in learning. When we add the quantum fluctuation to the system, it is a harmful computation to generate many samples through the quantum Monte-Carlo simulation. Instead, we here propose the message-passing algorithm in a deterministic way and it mitigate the difficulty in sampling. The results are validated by comparison with the natural generalization into the case with the quantum fluctuation. As far as our knowledge, this is the first investigation of the quantum effect in the restricted Boltzmann machine in theory while several algorithms for learning are proposed [53, 54]. The detailed analysis on the restricted Boltzmann machine in the transverse field and beyond will be reported elsewhere.

We again emphasize that the proposed method is not restricted to the case of a fully connected Ising model. At the level of belief propagation, we estimated both the longitudinal and transverse local magnetizations. For obtaining more accurate results, we may employ the cluster variational method [55] and its variants [56, 57, 58, 59].

5 Summary

We propose a message-passing algorithm for part of the nonstoquastic Hamiltonian such as the Ising model with antiferromagnetic XX interactions. Previous studies focused on the adaptive quantum Monte Carlo simulation, which generates the spin configuration in the equilibrium state of the limited range of the nonstoquastic Hamiltonian. However, it is often the case with long-time equilibration. Thus, belief propagation, despite being just an approximation, is valuable for investigating the equilibrium state of the nonstoquastic Hamiltonian. One of the nontrivial features in the nonstoquastic Hamiltonian is mitigating the computational complexity in QA by avoiding the first-order phase transition. This phenomenon is just a special case of several spin models showing a remarkable improvement in the performance of QA. Although such a nontrivial feature is not universal when using the nonstoquastic Hamiltonian, it is a significant step toward understanding the potential of QA. However, avoiding the first-order phase transition involved in the antiferromagnetic interaction is revealed at the level of static approximation [60]. In our approach, we propose belief propagation for part of the nonstoquastic Hamiltonian at the level of static approximation. Building on the present analysis, we intend to employ the nontrivial form of the messages because of which we may have to deal with multipoint correlations. Future studies beyond static approximation would reveal nontrivial aspects of the nonstoquastic Hamiltonian and bring out the potential of QA.

Acknowledgements.

The authors would like to thank Masamichi J. Miyama, Shuntaro Okada, Shunta Arai and Shu Tanaka for fruitful discussions. The present work is financially supported by JSPS KAKENHI Grant No. 15H03699 and 16H04382, and Next Generation High-Performance Computing Infrastructures and Applications R&D Program by MEXT, JST CREST (No. JPMJCR1402), ImPACT, and JST START.