Exact sampling hardness of Ising spin models

TL;DR

This paper demonstrates that classically sampling from the output distribution of certain Ising spin models is computationally hard, linking quantum dynamics to complex computational problems and extending the understanding of quantum sampling hardness.

Contribution

The paper constructs a specific Ising Hamiltonian whose exact output sampling is classically hard, advancing the classification of quantum sampling problems and connecting to complexity theory.

Findings

Exact sampling from the Ising model is computationally hard.

The results extend the classification of two-qubit commuting Hamiltonians.

Physical implementations may soon reach problem sizes where classical sampling is infeasible.

Abstract

We study the complexity of classically sampling from the output distribution of an Ising spin model, which can be implemented naturally in a variety of atomic, molecular, and optical systems. In particular, we construct a specific example of an Ising Hamiltonian that, after time evolution starting from a trivial initial state, produces a particular output configuration with probability very nearly proportional to the square of the permanent of a matrix with arbitrary integer entries. In a similar spirit to BosonSampling, the ability to sample classically from the probability distribution induced by time evolution under this Hamiltonian would imply unlikely complexity theoretic consequences, suggesting that the dynamics of such a spin model cannot be efficiently simulated with a classical computer. Physical Ising spin systems capable of achieving problem-size instances (i.e. qubit numbers) large enough so that classical sampling of the output distribution is classically difficult in practice may be achievable in the near future. Unlike BosonSampling, our current results only imply hardness of exact classical sampling, leaving open the important question of whether a much stronger approximate-sampling hardness result holds in this context. As referenced in a recent paper of Bouland, Mancinska, and Zhang (A. Bouland, L. Mancinska, and X. Zhang, CCC 2016, pp. 28:1-28:33), our result completes the sampling hardness classification of two-qubit commuting Hamiltonians.

Equations

$$\mathcal{H}=\sum _{i,j}{J}_{i,j}{\hat{\sigma }}_i^x{\hat{\tau }}_j^x.$$

$$|\psi (0)⟩=\underset{i\in \mathcal{A}}{⨂}{|\downarrow ⟩}_i\underset{j\in \mathcal{B}}{⨂}{|\downarrow ⟩}_j.$$

$$P_t({\sigma }_1$$

$$=$$

$$P_t\equiv P_t(\uparrow ,\dots ,\uparrow )=|⟨\uparrow ,\dots ,\uparrow |\exp (-it\mathcal{H})|\downarrow ,\dots ,\downarrow ⟩{|}^2\equiv |M_t{|}^2,$$

$${\hat{\sigma }}_i^x{\hat{\tau }}_j^x={\hat{\sigma }}_i^{+}{\hat{\tau }}_j^{+}+{\hat{\sigma }}_i^{+}{\hat{\tau }}_j^{-}+{\hat{\sigma }}_i^{-}{\hat{\tau }}_j^{+}+{\hat{\sigma }}_i^{-}{\hat{\tau }}_j^{-},$$

$$M_t$$

$$=(-it{)}^N\mathrm{Per}(J)+O(t^{N+2}),$$

$$\displaystyle P_{t}=t^{2N}\big{(}\mathscr{P}+O(t^{2})\big{)}.$$

$$p(m)$$

$$=p(N_{++}+N_{--})$$

$$=p(N_{++}-N_{--})$$

$$=P(N),$$

$$M_t=\sum _{\alpha =0}^{\infty }⟨\uparrow ,\dots ,\uparrow |\frac{(-it\mathcal{H}{)}^{N+2\alpha }}{(N+2\alpha )!}|\downarrow ,\dots ,\downarrow ⟩\equiv \sum _{\alpha =0}^{\infty }{M}_t^{(\alpha )},$$

$$M_t^{(0)}=(-it{)}^N\mathrm{Per}(J).$$

$$P_t$$

$$\displaystyle=t^{2N}\big{(}\mathscr{P}+\eta_{t}\big{)},$$

$${\eta }_t$$

$$\le |\delta M_t|(2|M_t^{(0)}|+|\delta M_t|)/t^{2N}.$$

$$|M_t^{(0)}|$$

$$|\delta M_t|$$

$${\eta }_t$$

$$\le 6N^4{t}^2\frac{N^{4N}}{(N!{)}^2}\le t^2\mathrm{poly}(\mathrm{N})e^{2N(\ln N+1)},$$

$$t=o(e^{-2N\ln N}).$$

Full text

I Introduction

It is often taken for granted that quantum computers can efficiently perform certain computational tasks that classical computers cannot. But finding a quantum task that, on the one hand, admits compelling complexity-theoretic arguments against efficient classical simulation, and on the other hand admits experimental demonstration with technology that is feasible in the near future, remains an important and challenging task in the field of quantum information science. An extremely exciting line of work, starting with results of Terhal and DiVincenzo and Bremner, Jozsa, and Shepherd, has shown that quantum computers are capable of sampling from distributions that cannot be sampled exactly by randomized classical algorithms Terhal and DiVincenzo (2004); Bremner et al. (2010). The BosonSampling protocol Aaronson and Arkhipov (2013), proposed by Aaronson and Arkhipov, gives a hardness of sampling result that may be within reach for near-term quantum experiments. The basic idea is to send photons through a network of linear optical devices, arranged in such a way that the probabilities of typical output configurations of the photons are proportional to the squares of permanents of matrices with independent and Gaussian-distributed random entries. Given reasonable assumptions about the hardness of computing permanents of such matrices, the ability to efficiently classically sample from any distribution even close (in total variation distance) to this distribution would imply extremely unlikely complexity theoretic consequences.

A number of proof-of-principle experiments implementing BosonSampling have already been carried out Broome et al. (2013); Spring et al. (2013); Tillmann et al. (2013); Crespi et al. (2013). However, a remaining bottleneck to producing an experimentally convincing demonstration of BosonSampling is the technical difficulty of building linear-optical systems that are large enough and clean enough to realize BosonSampling instances for which classical sampling is actually difficult. By comparison, state preparation and readout of individual spins can be done with high fidelity and relative ease, and the ability to massively parallelize spin-spin interactions between large numbers of qubits is reasonably sophisticated; experiments have successfully implemented some simple instances of the Ising model with system sizes ranging from tens Richerme et al. (2014) to many hundreds of spins Bohnet et al. (2016). Moreover, recent developments in ion-trapping experiments raise the exciting prospect of implementing arbitrary Ising interaction graphs in systems of (potentially) many tens of trapped ions Debnath et al. (2016). For this reason, finding results analogous to BosonSampling for simple spin models is highly desirable, and potentially affords a simpler route towards the experimental demonstration of an efficient quantum task that, under extremely plausible assumptions about classical complexity theory, cannot be efficiently performed by a classical system Bremner et al. (2010); Peropadre et al. (2015).

Our goal in this manuscript is to show that the dynamics of an experimentally implementable commuting spin model—the Ising model with no transverse field—can induce an output distribution over the spin states that is hard to sample from classically. The general strategy, which will be elaborated on below, is to divide a set of Ising spins into two mutually interacting registers, each having $N$ spins (see Fig. 1). The $N$ spins in the first and second register can be placed in correspondence with the $N$ row and column labels, respectively, of an $N\times N$ matrix $J$; each of the $N^{2}$ pairwise Ising couplings $J_{i,j}$ between a spin ($i$) in one register and a spin ($j$) in the other is a matrix element of $J$. By initializing the system in a spatially homogeneous product state and then letting it evolve under Ising interactions for a short time, it can be shown that a single probability of the output distribution induced by measurement is proportional to the square of the permanent of $J$, plus an $o(1)$ correction. This is enough, using a tool known as “Stockmeyer counting” Stockmeyer (1985), to imply a hardness of “exact sampling” result: no efficient classical randomized algorithm can sample from exactly this distribution, under a ubiquitous hardness assumption (namely, that the Polynomial-time Hierarchy does not collapse). Note that in a recent paper Peropadre et al. (2015), BosonSampling was directly generalized to the context of spin Hamiltonians. However, our work encounters the permanent in a fundamentally different way; an important difference is that our results do not rely on a “diluteness criterion”, and thus $N$ is set by—as opposed to much less than—the number of physical qubits. Much like other “exact sampling” results, our result also demonstrates hardness to classically sample from any distribution in which all probabilities are within a constant multiplicative factor of the ideal quantum distribution. However, unlike BosonSampling, a recent proposal of Bremner, Montanaro and Shepherd (sometimes called “IQP” sampling), and Quantum Fourier Sampling, it is not yet clear whether the distributions we consider can be used to show an “approximate-sampling” hardness result Aaronson and Arkhipov (2013); Bremner et al. (2010); Fefferman and Umans (2016). This would show something far stronger: there is no classical algorithm that can sample from any distribution inverse polynomial in total variation distance from the ideal quantum distribution.

II The model

The model we consider consists of $2N$ spin-1/2 particles, which we divide into two sublattices of $N$ spins each, denoted $\mathcal{A}$ and $\mathcal{B}$ (blue and red spins in Fig. 1). We consider quench dynamics under an Ising Hamiltonian with exclusively two-body inter-sublattice interactions (but no interactions within either sublattice), which can take arbitrary integer values,

[TABLE]

Here, Pauli operators $\hat{\sigma}$ act on the spins of sublattice $\mathcal{A}$, while Pauli operators $\hat{\tau}$ act on the spins of sublattice $\mathcal{B}$. These spins could be, for example, two subsets of ions in a Paul trap, where the $\ket{\downarrow}$ and $\ket{\uparrow}$ are, respectively, the electronic ground state and some long-lived metastable state (in general either an excited hyperfine level of the electronic ground-state manifold or a dipole-forbidden optical excitation). The Ising interactions can then be implemented via a spatially-structured Mølmer-Sørensen interaction Mølmer and Sørensen (1999); Korenblit et al. (2012); Debnath et al. (2016).

We consider a quantum quench in which the system is initialized at time $t=0$ with all of the spins (in both registers) in the spin-down state along the $z$-direction,

[TABLE]

We then allow the system to evolve under the Hamiltonian in Eq. (1) for a time $t$.

III Output distribution

After evolution for a time $t$ under the action of $\mathcal{H}$, measurement in the $z$ basis samples from the induced probability distribution

[TABLE]

where $\sigma_{j},\tau_{j}=\,\downarrow\,,\uparrow$. We are interested in just one such probability,

[TABLE]

to end in the state with all spins in both registers pointing up. By writing an individual term in the Hamiltonian as

[TABLE]

it is straightforward to see that repeated applications of $\mathcal{H}$, and thus time evolution, generates population in all possible spin states in the $z$ basis. Expanding $e^{-i\mathcal{H}t}$ as a power series in time, the lowest-order in time non-vanishing contribution to the matrix element $M_{t}=\bra{\uparrow,\dots,\uparrow}\exp(-it\mathcal{H})\ket{\downarrow,\dots,\downarrow}$ arises at order $t^{N}$, because every spin needs to be flipped at least once. The contributing terms contain exactly $N$ powers of operators $\hat{\sigma}^{+}_{i}\hat{\tau}_{j}^{+}$, with no repetitions of the indices $i$ and $j$, so that each qubit gets flipped from $\ket{\downarrow}$ to $\ket{\uparrow}$ exactly one time; see Fig. 2 for an illustration of such a term for $N=3$. It is straightforward to show that, to order $t^{N}$, the matrix element $M_{t}$ is given by

[TABLE]

where the summation is over all permutations $\sigma$ of the integers $i=1,\dots,N$. As a result, and defining $\mathscr{P}=|{\rm Per}(J)|^{2}$, we have

[TABLE]

We next aim to place a constraint on how $t$ must scale with $N$ in order to ensure that the $O(t^{2})$ additive error to the permanent is $o(1)$ with respect to the system size $N$.

IV Higher orders in time

As discussed above, the lowest-order in time contribution to the matrix element $M_{t}$ comes at order $N$. It is not hard to see that all other contributing terms occur at order $m$ such that $m-N$ is a positive even integer. In particular, take $N_{+-}$ to be the number of times an operator $\hat{\sigma}^{+}_{i}\hat{\tau}^{-}_{j}$ occurs inside the matrix element, and similarly for $N_{-+}$, $N_{++}$, and $N_{--}$, such that $N_{++}+N_{--}+N_{+-}+N_{-+}=m$. Since we need to flip the same number of qubits in both registers, we must have $N_{+-}=N_{-+}$. Also, the total number of flipped qubits is equal to $2(N_{++}-N_{--})$, and since all qubits need to be flipped, we have $N_{++}-N_{--}=N$. Now, defining $p(n)$ to be the parity of the integer $n$, we have

[TABLE]

which shows that $m-N$ is an even integer. The matrix element in question can therefore be expanded as

[TABLE]

and from above we have

[TABLE]

Defining $\delta{M}_{t}=\sum_{\alpha=1}^{\infty}M_{t}^{(\alpha)}$, such that $M_{t}=M_{t}^{(0)}+\delta M_{t}$, we can write

[TABLE]

where

[TABLE]

For notational simplicity, here we will assume that the entries of $J$ are drawn from the set $\{-1,0,1\}$; note that nothing about our argument would change if arbitrary integers were used, except that the time $t$ would be rescaled in the bounds below by $\max(J_{i,j})$. Using $\bra{\uparrow,\dots,\uparrow}\mathcal{H}^{m}\ket{\downarrow,\dots,\downarrow}\leq N^{2m}\lVert\hat{\sigma}^{x}\rVert^{2m}=N^{2m}$, $M_{t}^{(\alpha)}$ can be bounded as $|M_{t}^{(\alpha)}|\leq(N^{2}t)^{N+2\alpha}/(N+2\alpha)!$. Therefore,

[TABLE]

The final inequality in Eq. (13) is valid for $t^{2}\leq 1/(2N^{4})$, because $0\leq\sum_{\alpha=1}^{\infty}x^{\alpha}\leq 2x$ whenever $0\leq x\leq 1/2$. Plugging Eqs. (12,13) into Eq. (11) leads to

[TABLE]

with the final inequality obtained by Stirling’s approximation. It follows immediately that $\eta_{t}=o(1)$ is guaranteed as long as

[TABLE]

V Hardness of sampling

Here we prove our main theorem, establishing a very unlikely complexity theoretic consequence which would arise naturally from the presumed existence of a classical algorithm that samples exactly from the output distribution described in the prior sections. Similar arguments to the one sketched here are implicit in other works on quantum hardness of sampling results starting with the BosonSampling proposal Aaronson and Arkhipov (2013).

We first begin with a very brief overview of the computational complexity theoretic components necessary to understand this hardness of sampling result. Computing exactly the permanent of an $N\times N$ matrix $X$ with integer entries is as hard as computing the number of satisfying assignments to a boolean formula. We therefore say it is a ${\sf{\#{\sf{P}}}}$-hard problem, as established by Valiant Valiant (1979). When $X$ has nonnegative integer entries this problem is also in ${\sf{\#{\sf{P}}}}$.

For our purposes, we will be interested in the complexity of computing multiplicative estimates to the permanent. We say an algorithm $\mathscr{A}$ efficiently computes a multiplicative estimate to a function $f$ if, given input $x$, the output of $\mathscr{A}$ is within a $1\pm\epsilon$ multiplicative factor of $f(x)$ in time polynomial in $N$ and $1/\epsilon$. A famous result of Jerrum, Sinclair and Vigoda gives an algorithm for efficiently computing a multiplicative estimate to the permanent of a matrix with nonnegative entries Jerrum et al. (2004). On the other hand, it can be shown using a binary search and padding argument that computing such an estimate to the permanent (or even the square of the permanent) of a matrix with general integer entries is in fact ${\sf{\#{\sf{P}}}}$-hard (see e.g., Aaronson and Arkhipov (2013); Aaronson (2011)). Therefore computing these estimates are as hard as computing the permanent exactly. How powerful is ${\sf{\#{\sf{P}}}}$? We know from Toda’s Theorem that any problem in the Polynomial-time hierarchy, or ${\sf{PH}}$, can be solved using the ability to solve a ${\sf{\#{\sf{P}}}}$-hard problem Toda (1991). Being a bit more formal, Toda’s theorem tells us that ${\sf{PH}}\subseteq{\sf{P}}^{{\sf{\#{\sf{P}}}}}$.

Now, for any $N\times N$ matrix $X$ define $\mathcal{D}_{X}$ to be the outcome distribution from Section III that arises from starting in the $\ket{\downarrow,\dots,\downarrow}$ state, evolving for a particular time $t$ under the action of the Hamiltonian from Eq. (1) with coupling constants $J_{i,j}$ set to the entries of $X$, and measuring in the $z$ basis. As shown in Sections III and IV, the probability of observing the $\ket{\uparrow,\dots,\uparrow}$ outcome at time $t$ is proportional to the square of the permanent of $X$ plus an $o(1)$ correction, provided that $t$ is chosen to be $o(e^{-2N\ln N})$. Notice that this probability is exponentially small. Therefore, to get any reasonable estimate by repeated sampling we would need an exponential number of samples. Indeed, this does not imply an efficient quantum algorithm for computing the permanent. Nonetheless, we can use the fact that a single exponentially small amplitude is proportional to the permanent to argue about the classical intractability of sampling from this distribution.

Suppose we have an efficient classical sampler which samples from the same distribution. We define this to be an efficient randomized algorithm that takes as input an $N\times N$ integer matrix $X$ and outputs a sample from the distribution $\mathcal{D}_{X}$. A classic result of Stockmeyer gives an algorithm for computing a multiplicative estimate to the probability of any given outcome of an efficient classical sampler in the third level of the ${\sf{PH}}$, or $\sf\Sigma_{3}$ Stockmeyer (1985). Using this result, together with the presumed existence of an efficient classical sampler for our quantum distribution, we can compute a multiplicative estimate to the square of the permanent of an arbitrary integer matrix in the third level of the ${\sf{PH}}$. As mentioned above, this is a ${\sf{\#{\sf{P}}}}$-hard problem. This tells us we can solve any problem in ${\sf{\#{\sf{P}}}}$ in the third level of the Polynomial-time hierarchy, or formally, that ${\sf{P}}^{{\sf{\#{\sf{P}}}}}\subseteq\sf\Sigma_{3}$. Combining this with Toda’s theorem, we have that ${\sf{PH}}\subseteq{\sf{P}}^{{\sf{\#{\sf{P}}}}}\subseteq\Sigma_{3}$, and so the entire Polynomial-time Hierarchy collapses to the third level, as claimed. Therefore, it is very unlikely that an efficient classical sampler for the distribution with probabilities given by Eq. (3) exists.

VI Discussion and implications

These results extend several key ideas of BosonSampling to the context of spin dynamics under Ising spin Hamiltonians. Just like non-interacting bosons, the Ising model without a transverse field is often viewed—from the perspective of many-body quantum physics—to be trivial, since it can be trivially diagonalized. However, just as with non-interacting bosons, this point of view stems from a restricted notion of what it means to “simulate” a quantum system. As in the case of non-interacting bosons, it is indeed classically efficient to compute low-order correlation functions of operators in the model we study van den Worm et al. (2013); Foss-Feig et al. (2013), but sampling from the output distribution is simply a more general (and less trivial) task.

Another interesting motivation for our result comes from the desire to classify all two-qubit commuting Hamiltonians. Suppose we start in a computational basis state of $n$ qubits, and can apply a fixed two-qubit Hamiltonian to any pair of qubits. A recent result of Bouland, Mancinska, and Zhang gave a hardness of sampling classification for this model Bouland et al. (2016). They prove, in all cases except the one we consider (in which the two qubit Hamiltonian is $X\otimes X$) that the corresponding sampling task is classically hard, as long as the commuting Hamiltonian is capable of generating entanglement from a computational basis state. Otherwise, the output is in a product state and clearly classically simulable. Thus our hardness result completes the sampling hardness classification of the complete class of two-qubit commuting Hamiltonians (see their paper for additional details Bouland et al. (2016)).

VII Acknowledgments

We thank A. Deshpande and A. Bouland for helpful discussions. We also thank A. Bouland, Laura Mancinska and Xue Zhang for sharing an early version of their results. A.V.G. acknowledges support by ARL CDQI, ARO MURI, NSF QIS, ARO, NSF PFC at JQI, and AFOSR. This material is based upon work supported by, or in part by, the U. S. Army Research Laboratory and the U. S. Army Research Office under contract/grant number 025989-001.