Permutation Matrix Representation Quantum Monte Carlo

TL;DR

This paper introduces a versatile quantum Monte Carlo algorithm based on permutation matrix representations that unifies the simulation of various quantum and classical models without Trotter errors.

Contribution

The novel algorithm employs permutation group elements to represent quantum dimensions, enabling parameter-free, Trotter error-free simulations across diverse models.

Findings

Clarifies the origin of the sign problem in non-stoquastic models

Demonstrates improved performance over existing algorithms on transverse-field Ising models

Shows flexibility in simulating models with random two-body interactions

Abstract

We present a quantum Monte Carlo algorithm for the simulation of general quantum and classical many-body models within a single unifying framework. The algorithm builds on a power series expansion of the quantum partition function in its off-diagonal terms and is both parameter-free and Trotter error-free. In our approach, the quantum dimension consists of products of elements of a permutation group. As such, it allows for the study of a very wide variety of models on an equal footing. To demonstrate the utility of our technique, we use it to clarify the emergence of the sign problem in the simulations of non-stoquastic physical models. We showcase the flexibility of our algorithm and the advantages it offers over existing state-of-the-art by simulating transverse-field Ising model Hamiltonians and comparing the performance of our technique against that of the stochastic series expansion algorithm. We also study a transverse-field Ising model augmented with randomly chosen two-body transverse-field interactions.

Tables

Table 1. Table 1: Cycle completion moves for the transverse field Ising model with two-body X 𝑋 X interactions. The moves include insertions or removals of pairs of identical one-body and two-body X 𝑋 X terms [i) and and ii)], the breaking up of a two-body term to its one-body constituents and the inverse operation [iii)], swapping [iv)] and the contraction of two operators into one [v)].

		Move		Change in $q$
i)	$\mathrm{𝟙}$	$\leftrightarrow$	($X_i,X_i$)	$\pm 2$
ii)	$\mathrm{𝟙}$	$\leftrightarrow$	($X_i{X}_j,X_i{X}_j$)	$\pm 2$
iii)	$X_i{X}_j$	$\leftrightarrow$	($X_i,X_j$)	$\pm 1$
iv)	($X_i,X_j$)	$\leftrightarrow$	($X_j,X_i$)	no change
v)	($X_i{X}_j,X_j{X}_k$)	$\leftrightarrow$	$X_i{X}_k$	$\pm 1$

Equations

$$H=\sum _{j=0}^M{\tilde{P}}_j=\sum _{j=0}^M{D}_j{P}_j=H_c+\sum _{j=1}^M{D}_j{P}_j,$$

$$H=\sum _j{R}_j\left(e^{i{\Phi }_j}{P}_j+e^{-i{\Phi }_j}{P}_j^{-1}\right),$$

$$H={\alpha }_{0}1+{\alpha }_{1}X+{\alpha }_{2}Y+{\alpha }_{3}Z,$$

$$H=D_0{P}_0+D_1{P}_1$$

$$H=\sum _{i<j}\sum _{\begin{array}{c}{K}_i\in \left\{1_i,X_i,Y_i,Z_i\right\}\\ K_j\in \left\{1_j,X_j,Y_j,Z_j\right\}\end{array}}{\alpha }_{ij,K_i,K_j}{K}_i{K}_j.$$

$$P_0=1=\left[\begin{array}{ccc}1& 0& 0\\ 0& 1& 0\\ 0& 0& 1\end{array}\right],P_1=\left[\begin{array}{ccc}0& 0& 1\\ 1& 0& 0\\ 0& 1& 0\end{array}\right],P_2=\left[\begin{array}{ccc}0& 1& 0\\ 0& 0& 1\\ 1& 0& 0\end{array}\right],$$

$$H=-t\sum _{⟨i,j⟩}{\hat{b}}_i^{†}{\hat{b}}_j+\frac{U}{2}\sum _i{\hat{n}}_i\left({\hat{n}}_i-1\right)-\mu \sum _i{\hat{n}}_i.$$

$$Z$$

$$d_{z_j}^{(i_j)}=⟨z_j|D_{i_j}|z_j⟩,$$

$$Z$$

$$\sum _{\{k_i\}}\frac{(-\beta {)}^q}{(q+\sum k_i)!}\prod _{j=0}^q(-\beta E_{z_j}{)}^{k_j}=e^{-\beta [E_{z_0},\dots ,E_{z_q}]}$$

$$,$$

$$F[E_{z_0},\dots ,E_{z_q}]\equiv \sum _{j=0}^q\frac{F(E_{z_j})}{{\prod }_{k\ne j}(E_{z_j}-E_{z_k})}$$

$$F[E_{z_0},\dots ,E_{z_q}]=e^{-\beta [E_{z_0},\dots ,E_{z_q}]}.$$

$$Z=\sum _z\sum _{q=0}^{\infty }\sum _{\{S_q\}}⟨z|S_{{\mathbf{i}}_q}|z⟩D_{(z,S_{{\mathbf{i}}_q})}{e}^{-\beta [E_{z_0},\dots ,E_{z_q}]},$$

$$D_{(z,S_{{\mathbf{i}}_q})}=\prod _{j=1}^q{d}_{z_j}^{(i_j)}.$$

$$e^{-\beta [E_0,\dots ,E_q]}={\int }_{\Omega }\text{d}{t}_0\cdots \text{d}{t}_q{e}^{-\beta \left(E_0{t}_0+E_1{t}_1+\dots +E_q{t}_q\right)},$$

$$W_{\mathcal{C}}=D_{(z,S_{{\mathbf{i}}_q})}{e}^{-\beta [E_{z_0},\dots ,E_{z_q}]},$$

$$Z=\sum _z\sum _{S_{{\mathbf{i}}_q}=1}{D}_{(z,S_{{\mathbf{i}}_q})}{e}^{-\beta [E_{z_0},\dots ,E_{z_q}]}.$$

$$(-1{)}^q{D}_{(z,S_{{\mathbf{i}}_q})}=\mathrm{Re}\left[\prod _{j=1}^q(-d_{z_j}^{(i_j)})\right]$$

$$H=\sum _{i,j}{J}_{ij}{Z}_i{Z}_j+\sum _j{h}_j{Z}_j+\Gamma \sum _j{X}_j.$$

$$H=\sum _{i,j}{J}_{ij}{Z}_i{Z}_j+\Gamma \sum _{⟨i,j⟩}{X}_i{X}_j,$$

$$W_{{\mathcal{C}}_0}=e^{-\beta [E_z]}=e^{-\beta E_z},$$

$$p=\min \left(1,\frac{W_{{\mathcal{C}}^{\prime }}}{W_{\mathcal{C}}}\right)=\min \left(1,\frac{e^{-\beta [E_{{\mathcal{C}}^{\prime }}]}}{e^{-\beta [E_{\mathcal{C}}]}}\right),$$

$$P_{i_j}\cdots P_{i_{j+C}}=1.$$

$$P_{i_j}\cdots P_{i_k}{P}_{i_{k+1}}\cdots P_{i_{j+C}}\to P_{i_{k+1}}\cdots P_{i_{j+C}}{P}_{i_j}\cdots P_{i_k}.$$

$$|z^{\prime }⟩=S_1|z⟩=P_{i_k}\cdots P_{i_1}|z⟩,$$

$$P_{i_j}\cdots P_{i_{j+C}}=1$$

$${\left(P_{i_{k+1}}\cdots P_{i_{j+C}}\right)}^{-1}=P_{i_{j+C}}^{-1}\cdots P_{i_{k+1}}^{-1}.$$

Full text

I Introduction

Quantum Monte Carlo (QMC) algorithms Landau and Binder (2005); Barkema (1999) are extremely useful for studying equilibrium properties of large quantum many-body systems, with applications ranging from superconductivity and novel quantum materials Aspuru-Guzik et al. (2005); Kassal et al. (2008); Lanyon et al. (2010) through the physics of neutron stars Lonardoni et al. (2014) and quantum chromodynamics Chandrasekharan and Wiese (1999); Kieu and Griffin (1994). The algorithmic development of QMC remains an active area of research, with the dual goal of extending the scope of QMC applicability and improving convergence rates of existing algorithms in order to facilitate the discovery of new phenomena Sandvik (1992, 1999); Prokof’ev et al. (1998).

While QMC algorithms have been adapted to the simulation of a wide variety of physical systems, different models typically require the development of distinct model-specific update rules and measurement schemes. A notable recent example is the transverse-field Ising model (TFIM), which traditionally includes only single-body $X$ terms (the transverse field), supplemented with two-body $X$ terms. While the updates associated with single-body $X$ terms can be implemented by local (in space) updates, the inherently non-local nature of the two-body $X$ terms requires novel cluster updates Mazzola and Troyer (2017). Thus, a proper treatment of the Hamiltonian with single-body and two-body $X$ terms requires updates that are different than if the Hamiltonian included only single-body or only two-body $X$ terms.

In this paper, we provide a QMC scheme that has the flexibility to simulate a broad range of quantum many-body models. The technique we propose here builds on a power-series expansion of the canonical quantum partition function about the classical partition function first introduced in Refs. Albash et al. (2017); Hen (2018)). While strongly inspired by earlier methods that expand the partition function in powers of the inverse-temperature Handscomb (1962, 1964); Sandvik and Kurkijärvi (1991); Sandvik (1992), our expansion is more accurately described as an expansion in off-diagonal operators of the Hamiltonian.

Our formalism enables a very general treatment of Hamiltonians, allowing us to develop a QMC scheme that is applicable to a wide variety of models, ranging from highly interacting models with multi-body terms to non-interacting ones and from strongly quantum models to purely classical ones, using the same updating formalism. In order to demonstrate the advantage of the new method, we study the TFIM with random two-body longitudinal interactions, comparing the performance of our technique against that of stochastic series expansion QMC Sandvik (2003). In addition, we study in detail a variant of the model with added random $XX$ interactions. This model poses a challenge for traditional QMC approaches Barker (1979); Mazzola and Troyer (2017) due to the random connectivity of the Ising couplings and $XX$ interactions, which vary between instances. While we focus most of our attention in what follows on finite-dimensional Hamiltonians, the technique we present here should apply with equal rigor to infinite-dimensional systems.

The paper is organized as follows. In Sec. II, we describe the ‘permutation matrix representation’ (PMR) of Hamiltonians on which the partition function expansion detailed in Sec. III is founded. We discuss the emergence of the sign problem within the formulation and its sometimes-intricate relation with the concept of non-stoquasticity in Sec. IV and in Sec. V we present the QMC algorithm we have devised based on the expansion. Results of TFIM simulations are presented in Sec. VI. We conclude in Sec. VII with additional discussions and some caveats.

II The permutation matrix representation

We consider many-body systems whose Hamiltonians we cast as the sum

[TABLE]

where $\{\tilde{P}_{j}\}$ is a set of $M+1$ distinct generalized permutation matrices Joyner (2008), i.e., matrices with precisely one nonzero element in each row and each column (this condition can be relaxed to allow for rows and columns with only zero elements). Each operator $\tilde{P}_{j}$ can be written, without loss of generality, as $\tilde{P}_{j}=D_{j}P_{j}$ where $D_{j}$ is a diagonal matrix111The diagonal matrix $D_{j}$ will be invertible, i.e., will not contain zero elements along the diagonal, if $\tilde{P}_{j}$ is a bonafide generalized permutation matrix. and $P_{j}$ is a permutation matrix with no fixed points (equivalently, no nonzero diagonal elements) except for the identity matrix $P_{0}=\mathbb{1}$. We will refer to the basis in which the operators $\{D_{j}\}$ are diagonal as the computational basis and denote its states by $\{|z\rangle\}$. We will call the diagonal matrix $D_{0}$ the ‘classical Hamiltonian’ and will sometimes denote it by $H_{\text{c}}$. The permutation matrices appearing in $H$ will be treated as a subset of a permutation group, wherein $P_{0}$ is the identity element.

The $\{D_{j}P_{j}\}$ off-diagonal operators (in the computational basis) give the system its ‘quantum dimension’. Each term $D_{j}P_{j}$ obeys $D_{j}P_{j}|z\rangle=d_{j}(z^{\prime})|z^{\prime}\rangle$ where $d_{j}(z^{\prime})$ is a possibly complex-valued coefficient and $|z^{\prime}\rangle\neq|z\rangle$ is a basis state. While the above formulation may appear restrictive, we show in Appendix A that any finite-dimensional matrix can be written in the form of Eq. (1).

We also note that $H=\sum_{j}D_{j}P_{j}$ is hermitian if and only if for every index $j$ there is an associated index $j^{\prime}$ such that $P_{j}=P_{j^{\prime}}^{-1}$ and $D_{j}=D_{j^{\prime}}^{*}$ where the indices $j$ and $j^{\prime}$ can be the same (see Appendix B). This in turn implies that any Hamiltonian $H$ can be written as

[TABLE]

where $R_{j},\Phi_{j}$ are real-valued diagonal matrices. In the case where a permutation matrix $P_{j}$ is its own inverse, the corresponding $\Phi_{j}$ will necessarily be the zero matrix.

To further elucidate PMR, we now provide several examples.

II.1 Example I: A single spin-$1/2$ particle

The Hamiltonian of a single spin-$1/2$ particle can most generally be written as

[TABLE]

where $X,Y$ and $Z$ are the matrix representations of the usual Pauli operators in the basis that diagonalizes the Pauli-$Z$ operator. In PMR, the Hamiltonian becomes

[TABLE]

with $P_{0}=\mathbb{1}$, $P_{1}=X$, $D_{0}=H_{c}=\alpha_{0}\mathbb{1}+\alpha_{3}Z$ and $D_{1}=\alpha_{1}\mathbb{1}-i\alpha_{2}Z$.

II.2 Example II: Two-local spin-$1/2$ models

A general two-local $n$-particle spin-$1/2$ Hamiltonian has similarly the following form

[TABLE]

Here, the basis states are tensor products of the single spin states. We can cast the Hamiltonian in the form of Eq. (1) by grouping together elements that change a given basis state $|z\rangle$ to the same basis state $|z^{\prime}\rangle$. For example, the terms $X_{i}$, $Y_{i}$, $X_{i}Z_{j}$, $Y_{i}Z_{j}$ are grouped together as the action of the combined term $V_{i}=\alpha_{ij10}X_{i}+\alpha_{ij20}Y_{i}+\sum_{j}(\alpha_{ij13}X_{i}Z_{j}+\alpha_{ij23}Y_{i}Z_{j})$ can be written as $D_{i}X_{i}$, where $D_{i}$ is a (generally complex-valued) diagonal matrix. This approach can be straightforwardly generalized to three- and higher-local Hamiltonians. The permutation matrices $P_{i}$ for general spin-1/2 Hamiltonians are by extension $\{\mathbb{1},X_{i},\ldots,X_{i}X_{j},\ldots,X_{i}X_{j}X_{k},\ldots\}$.

II.3 Example III: Spin-one particles (qutrits)

The permutation matrix representation generalizes straightforwardly to higher dimensional systems. The Hamiltonian for a single qutrit can be written as $H=D_{0}P_{0}+D_{1}P_{1}+D_{2}P_{2}$ where

[TABLE]

and $D_{2}=D_{1}^{\ast}$, a condition imposed by the hermiticity of the Hamiltonian.

II.4 Example IV: The Bose-Hubbard model

Another model that can just as easily be represented in permutation matrix form is the Bose-Hubbard model. This discretely infinite dimensional model captures the physics of interacting spinless bosons on a lattice Lewenstein et al. (2012) and is commonly used to describe superfluid-insulator transitions Fisher et al. (1989), bosonic atoms in an optical lattice Jaksch and Zoller (2005) and certain magnetic insulators Giamarchi et al. (2008).

The Bose-Hubbard Hamiltonian is given by

[TABLE]

Here, $\left\langle i,j\right\rangle$ denotes summation over all neighboring lattice sites $i$ and $j$, while $\hat{b}^{\dagger}_{i}$ and $\hat{b}_{i}$ are regular bosonic creation and annihilation operators such that $\hat{n}_{i}=\hat{b}^{\dagger}_{i}\hat{b}_{i}$ gives the number of particles at the $i$-th site. The model is parametrized by the hopping amplitude $t$ and the on-site interaction $U$.

In the bosonic number basis where states are described by the number of bosons in each site $|n_{1}\rangle\ldots|n_{L}\rangle$ (here $L$ is the number of lattice sites) we identify the diagonal part to be $D_{0}=\frac{U}{2}\sum_{i}\hat{n}_{i}\left(\hat{n}_{i}-1\right)-\mu\sum_{i}\hat{n}_{i}$ and the off-diagonal (infinite dimensional) permutation operators as $P_{\left\langle i,j\right\rangle}=\hat{b}^{\dagger}_{i}\hat{b}_{j}$. The diagonal operators associated with $P_{\left\langle i,j\right\rangle}$ are $D_{\left\langle i,j\right\rangle}$ whose entries are $-t$ for states whose $n_{j}$ is positive (the $j$-th boson can be annihilated) and zero otherwise.

III Off-diagonal partition function expansion

We are now in a position to discuss the off-diagonal series expansion of the partition function $Z={\mbox{Tr\,}}\left[{{e}}^{-\beta H}\right]$ as it applies to Hamiltonians cast in the form given in Eq. (1).

We begin by replacing the trace operation ${\mathrm{Tr}}[\cdot]$ with the explicit sum $\sum_{z}\langle z|\cdot|z\rangle$ and then expanding the exponent in the partition function in a Taylor series in $\beta$:

[TABLE]

In the last step we have expressed $(-H)^{n}$ in terms of all sequences of length $n$ composed of products of basic operators $H_{\text{c}}$ and $D_{j}P_{j}$, which we have denoted by the set $\{S_{{\bf i}_{n}}\}$. Here ${\bf i}_{n}=(i_{1},i_{2},\ldots,i_{n})$ is a set of indices, each of which runs from [math] to $M$, that denotes which of the $M+1$ operators in $H$ appear in $S_{{\bf i}_{n}}$.

We proceed by stripping away all the diagonal Hamiltonian terms from the sequence $\langle z|{S}_{{\bf i}_{n}}|z\rangle$. We do so by evaluating the action of these terms on the relevant basis states, leaving only the off-diagonal operators unevaluated inside the sequence (see Refs. Albash et al. (2017); Hen (2018) for a more detailed derivation).

The partition function may then be written as

[TABLE]

where $E_{z_{i}}=\langle z_{i}|H_{\text{c}}|z_{i}\rangle$ and ${\{S_{{\bf{i}}_{q}}\}}$ denotes the set of all products of length $q$ of ‘bare’ off-diagonal operators $P_{j}$. Also

[TABLE]

which can be considered as the ‘hopping strength’ of $P_{i_{j}}$ with respect to $|z_{j}\rangle$. Note that while the partition function is positive and real-valued, the $d^{(i_{j})}_{z_{j}}$ elements do not necessarily have to be so.

The term in parentheses in Eq. (8) sums over the diagonal contribution of all $\langle z|S_{{\bf i}_{n}}|z\rangle$ terms that correspond to the same $\langle z|S_{{\bf{i}}_{q}}|z\rangle$ term. The various $\{|z_{i}\rangle\}$ states are the states obtained from the action of the ordered $P_{j}$ operators in the product $S_{{\bf{i}}_{q}}$ on $|z_{0}\rangle$, then on $|z_{1}\rangle$, and so forth. For example, for $S_{{\bf{i}}_{q}}=P_{i_{q}}\ldots P_{i_{2}}P_{i_{1}}$, we obtain $|z_{0}\rangle=|z\rangle,P_{i_{1}}|z_{0}\rangle=|z_{1}\rangle,P_{i_{2}}|z_{1}\rangle=|z_{2}\rangle$, etc. The proper indexing of the states $|z_{j}\rangle$ along the path is $|z_{(i_{1},i_{2},\ldots,i_{j})}\rangle$ to indicate that the state at the $j$-th step depends on all $P_{i_{1}}\ldots P_{i_{j}}$. We will use the shorthand $|z_{j}\rangle$. The sequence of basis states $\{|z_{i}\rangle\}$ may be viewed as a ‘walk’ on the hypercube of basis states Albash et al. (2017); Hen (2018, 2019) (see Fig. 1).

After a change of variables, $n\to n+q$, we arrive at:

[TABLE]

Noting that the various $\{E_{z_{i}}\}$ are the classical energies of the states $|z_{i}\rangle$ states created by the operator product $S_{{\bf{i}}_{q}}$, the partition function is now given by:

[TABLE]

A feature of the above infinite sum is that the term in parentheses can be further simplified to give the exponent of divided differences of the $E_{z_{i}}$’s (a short description of divided differences and an accompanying proof of the above assertion can be found in Ref. Albash et al. (2017)); it can therefore be succinctly rewritten as:

[TABLE]

where $[E_{z_{0}},\ldots,E_{z_{q}}]$ is a multiset of energies and where a function $F[\cdot]$ of a multiset of input values is defined by

[TABLE]

and is called the divided differences Whittaker and Robinson (1967); de Boor (2005) of the function $F[\cdot]$ with respect to the list of real-valued input variables $[E_{z_{0}},\ldots,E_{z_{q}}]$. In our case, $F[\cdot]$ is the function

[TABLE]

Therefore, the infinite sum over energies in Eq. (III) may be simplified to

[TABLE]

where we have denoted

[TABLE]

We stress that the partition function expanded as such is not an expansion in $\beta$. Similar to the techniques introduced by Handcomb in the 1960s Handscomb (1962, 1964) and further developed in the stoquastic series expansion (SSE) scheme pioneered by Sandvik Sandvik and Kurkijärvi (1991); Sandvik (1992), the off-diagonal series expansion begins with a Taylor series expansion of the exponential function in the inverse temperature $\beta$, but the regrouping of terms into the exponent of divided-differences means that it is no longer a high-temperature expansion. SSE writes the trace of products of the Hamiltonian as a sum of products of matrix elements written in a suitably chosen basis which are then sampled (thereby overcoming the limitation of Handscomb’s scheme which required the evaluation of traces of products of the Hamiltonian). This is made possible by breaking up the Hamiltonian into a sum of local bonds. In PMR, this is not the case — all terms may remain non-local but are grouped according to their action on basis states.

Specially, the diagonal portion of the Hamiltonian remains ‘intact’ which then allows us to regroup a large portion of all terms. A single divided-difference term thus corresponds to a sum of an infinite number of SSE terms, meaning that a single PMR configuration represents very many standard SSE configurations and a single PMR weight sums up very many standard SSE weights. This can be immediately seen in the derivation above, particularly Eq. (8), which relates the standard SSE weight, which involves sequences of diagonal as well as off-diagonal bonds, to the weights of the current approach that only involve off-diagonal bonds.

It is worth noting that the cost of the massive grouping of the off-diagonal series expansion is manifested in the computational cost associated with calculating PMR terms (or ratios thereof). As we discuss in Sec. V in more detail, these can be calculated (or more precisely, updated) using $O(q)$ basic arithmetic operations (where $q$ is the number of operators in a sequence or equivalently the size of the imaginary-time dimension), as opposed to SSE’s $O(1)$. We refer the reader to Ref. Albash et al. (2017) for a more detailed comparison between the off-diagonal expansion and SSE. In Sec. VI.1 we provide a runtime comparison of PMR vs SSE for simulations of TFIM instances.

Aside from SSE, it is also interesting to observe that the exponent of divided differences also has close relations to continuous-time QMC (e.g., Ref. Prokof’ev et al. (1998)), via the Hermite-Genocchi formula de Boor (2005):

[TABLE]

where $t_{i}\geq 0$ and the area of integration $\Omega$ is bounded by $t_{0}+t_{1}+\ldots+t_{q}$ from above.

Having derived the expansion Eq. (15) for any Hamiltonian cast in the form Eq. (1), we are now in a position to interpret the partition function expansion as a sum of weights, i.e., $Z=\sum_{\{{\cal C}\}}W_{{\cal C}}$, where the set of configurations $\{{\cal C}\}$ is all the distinct pairs $\{|z\rangle,S_{{\bf{i}}_{q}}\}$. Because of the form of $W_{{\cal C}}$,

[TABLE]

we refer to it as a ‘generalized Boltzmann weight’ (or, a GBW). It can be shown Albash et al. (2017) that $(-1)^{q}{{e}}^{-\beta[E_{z_{0}},\ldots,E_{z_{q}}]}$ is strictly positive. Another feature of divided differences is that they are invariant under rearrangement of the input values.

We note that as written, the weights $W_{\cal C}$ are complex-valued, despite the partition function being real (and positive). Since for every configuration $\mathcal{C}=\{|z\rangle,S_{{\bf{i}}_{q}}\}$ there is a conjugate configuration $\bar{\mathcal{C}}=\{|z\rangle,S^{\dagger}_{{\bf i}_{q}}\}$222For $S_{{\bf{i}}_{q}}=P_{i_{q}}\ldots P_{i_{2}}P_{i_{1}}$, the conjugate sequence is simply $S^{\dagger}_{{\bf i}_{q}}=P_{i_{1}}^{-1}P_{i_{2}}^{-1}\ldots P_{i_{q}}^{-1}$. that produces the conjugate weight $W_{\bar{\mathcal{C}}}=\bar{W}_{\mathcal{C}}$, the imaginary contributions cancel out. Expressed differently, the imaginary portions of complex-valued weights do not contribute to the partition function and may be disregarded altogether. We may therefore redefine $D_{(z,S_{{\bf{i}}_{q}})}=\operatorname{Re}\left[\prod_{j=1}^{q}d^{(i_{j})}_{z_{j}}\right]$, obtaining strictly real-valued weights.

Before we move on, we note that $\langle z|S_{{\bf{i}}_{q}}|z\rangle$ evaluates either to 1 or to zero. Moreover, since the permutation matrices with the exception of $P_{0}$ have no fixed points, the condition $\langle z|S_{{\bf{i}}_{q}}|z\rangle=1$ implies $S_{{\bf{i}}_{q}}=\mathbb{1}$, i.e., $S_{{\bf{i}}_{q}}$ must evaluate to the identity element $P_{0}$ (note that the identity element does not appear in the sequences $S_{{\bf{i}}_{q}}$). The expansion can thus be more succinctly rewritten as

[TABLE]

IV Non-stoquasticity and emergence of the sign problem

An attractive property of the formalism introduced above is that it allows us to identify the emergence of the sign problem in QMC via inspection of the weights $W_{\cal C}$, thereby making more apparent the connection between the notion of non-stoquasticity — the existence of positive or complex-valued off-diagonal Hamiltonian matrix entries — which has garnered increasing attention with the advent of quantum computers in recent years Bravyi et al. (2008); Bravyi and Hastings (2014); Marvian et al. (2019) and the onset of the sign problem.

To interpret the real-valued weight terms $W_{\mathcal{C}}$ as actual weights (equivalently, un-normalized probabilities), they must be nonnegative. The occurrence of negative weights marks the onset of the infamous sign problem. A weight is positive iff

[TABLE]

is positive, that is, a QMC algorithm will encounter a sign problem, equivalently a negative weight, during a simulation if and only if there exists a closed walk on the hypercube of basis states along which $\operatorname{Re}\left[\prod_{j=1}^{q}(-d^{(i_{j})}_{z_{j}})\right]<0$. It is thus clear that it is not mere non-stoquasticity (equivalently, the sign of off-diagonal entries) that creates the sign problem, but rather the sign of closed walks on the hypercube of basis states that determines its occurrence.

A special class of models where the sign problem does not emerge, i.e., where $\operatorname{Re}\left[\prod_{j=1}^{q}(-d^{(i_{j})}_{z_{j}})\right]\geq 0$ for all configurations, is that of ‘stoquastic’ Hamiltonians Bravyi et al. (2008); Bravyi and Hastings (2014) for which all $d^{(i_{j})}_{z_{j}}$ are negative, which is equivalent to having only nonpositive off-diagonal elements in the matrix representation of the Hamiltonian. In this case, all products trivially yield positive-valued walks.

The existence of positive off-diagonal terms does not however immediately imply a sign problem for QMC (the reader is referred to Ref. Gupta and Hen (2019) for a more detailed discussion). Another example of a sign-problem-free family of models is one where all $d^{(i_{j})}_{z_{j}}$ elements are positive but closed walks are all of even length. One such model is the transverse-field Ising Hamiltonian

[TABLE]

for $\Gamma>0$. A slightly less trivial example is the two-body model

[TABLE]

provided that the underlying connectivity $\langle i,j\rangle$ of the two-body $X$ terms is bi-partite (allowing only even cycles).

It is also interesting to note that any single-qubit Hamiltonian is necessarily also sign-problem-free. In this case, the Hamiltonian is $H=D_{0}P_{0}+D_{1}P_{1}$ as described in Sec. II.1. Since here the $S_{{\bf{i}}_{q}}$ are sequences consisting of only one type of non-identity permutation matrices, namely $P_{1}=X$, the expansion order $q$ must be even for $S_{{\bf{i}}_{q}}$ to evaluate to the identity element. This in turn results in $\left[\prod_{j=1}^{q}(-d^{(i_{j})}_{z_{j}})\right]={(\alpha^{2}_{1}+\alpha^{2}_{2})}^{q/2}$ being strictly nonnegative. The same is however not true for a single qutrit in which case a sign problem may arise.

V The QMC algorithm

Having derived the series expansion of the partition function for permutation-represented Hamiltonians, we are now in a position to discuss a QMC algorithm that can be associated with the above expansion.

V.1 QMC configurations and GBW calculation

As was discussed above, a configuration ${\cal C}=\{|z\rangle,S_{{\bf{i}}_{q}}\}$ is a pair of a classical state and a product $S_{{\bf{i}}_{q}}$ of permutation operators that must evaluate to the identity element $P_{0}=\mathbb{1}$. To take full advantage of our partition function decomposition above, we treat the off-diagonal permutation terms $\{P_{j}\}$ in our QMC algorithm as elements in a permutation group $G$ (with matrix product as the group operation). Since the elements $\{P_{j}\}$ appearing in the Hamiltonian may not form a complete group, we shall treat any additional element $P_{j^{\prime}}$ required to complete the set to form a group as appearing in the Hamiltonian with an associated diagonal matrix $D_{j^{\prime}}=0$ [see Eq. (1)].

The pair ${\cal C}$ induces a list of states $\{|z_{0}\rangle=|z\rangle,|z_{1}\rangle,\ldots,|z_{q}\rangle=|z\rangle\}$, which in turn also generates a corresponding multiset of diagonal energies $E_{{\cal C}}=\{E_{z_{0}},E_{z_{1}},\ldots,E_{z_{q}}\}$ of not-necessarily-distinct values (recall that $E_{z_{i}}=\langle z_{i}|H_{c}|z_{i}\rangle$). For systems with discrete energy values, the multiset can be stored efficiently in a ‘multiplicity table’ $M_{{\cal C}}=\{m_{0},m_{1},\ldots,m_{j},\ldots\}$, where $m_{j}$ is the multiplicity of the energy $E_{z_{j}}$ in the multiset.

Given the multiset $E_{{\cal C}}$, the evaluation of the GBW $W_{{\cal C}}$ follows from its definition as a function of divided differences (the reader is referred to Ref. Albash et al. (2017) for a more detailed description). The calculation of a GBW consisting of $q$ permutation operators requires the evaluation of a divided-differences exponential with $(q+1)$ energies. This calculation can be accomplished with at most $O(q)$ operations Zivcovich (2019); Gupta et al. (2020).

V.2 Initial state

At this point we can consider a QMC algorithm based on the partition function expansion generating the weights $W_{{\cal C}}$, Eq. (18). The Markov process would start with the initial configuration $\mathcal{C}_{0}=\{|z\rangle,S_{0}=\mathbb{1}\}$ where $|z\rangle$ is a randomly generated initial classical state. The weight of this initial configuration is

[TABLE]

i.e., the classical Boltzmann weight of the initial random state $|z\rangle$.

V.3 Updates

We next describe the basic update moves for the algorithm. These are also succinctly summarized in Fig. 2.

V.3.1 Classical moves

Classical moves are any moves that involve a manipulation of the classical state $|z\rangle$ while leaving $S_{{\bf{i}}_{q}}$ unchanged [see Fig. 2(a)]. In a single bit-flip classical move, a spin from the classical bit-string state $|z\rangle$ of ${\cal C}$ is picked randomly and is flipped, generating a state $|z^{\prime}\rangle$ and hence a new configuration ${\cal C}^{\prime}$. Calculating the weight of ${\cal C}^{\prime}$ requires recalculating the energies associated with the product $S_{{\bf{i}}_{q}}$ leading to a new energy multiset $E_{{\cal C}^{\prime}}$ and can become computationally intensive if $q$ is large. Classical moves should therefore be attempted with low probabilities if $q$ is large. Simply enough, the acceptance probability for a classical move is

[TABLE]

where ${{e}}^{-\beta[E_{{\cal C}}]}$ is a shorthand for ${{e}}^{-\beta[E_{z_{0}},E_{z_{1}},\ldots,E_{z_{q}}]}$ of configuration ${\cal C}$ and likewise for ${\cal C}^{\prime}$.

In the absence of a quantum part to the Hamiltonian ($D_{j}=0$ for all $j>0$), not only are classical moves the only moves necessary, they are also the only moves that have nonzero acceptance probabilities. Since the initial configuration of the QMC algorithm is a random classical configuration $|z\rangle$ and an empty operator sequence $S_{0}=\mathbb{1}$, for a purely classical Hamiltonian, the algorithm automatically reduces to a classical thermal algorithm keeping the size of the imaginary-time dimension at zero ($q=0$) for the duration of the simulation.

V.3.2 Cyclic rotations

The ‘cyclic rotation’ move, [Fig. 2(b)], consists of identifying short sub-sequences, or cycles, of consecutive operators in the sequence $S_{{\bf{i}}_{q}}$, whose product is the identity element, i.e., sub-sequences that obey

[TABLE]

Depending on the nature of the operators, preparing a lookup table of short cycles that evaluate to the identity may prove useful. Once a cycle is identified, a random cycle rotation is attempted. Here, a random internal insertion point within the sub-sequence is picked and a rotation is attempted:

[TABLE]

The rotated sequence also evaluates to the identity. Since the internal classical states between the elements in the cycle may change by the rotation, the rotation involves adding new energies $\{E_{z^{\prime}}\ldots\}$ and removing old ones $\{E_{z^{\prime}}\ldots\}-\{E_{z}\ldots\}$ from the energy multiset. Short cycles should therefore be preferred. The acceptance probability for the move is as in Eq. (23) with $E_{{\cal C}^{\prime}}=E_{{\cal C}}+\{E_{z^{\prime}}\ldots\}-\{E_{z}\ldots\}$.

V.3.3 Block-swap

A block swap [Fig. 2(c)] is an update that involves a change of the classical state $z$. Here, a random position $k$ in the product $S_{{\bf{i}}_{q}}$ is picked such that the product is split into two (non-empty) sub-sequences, $S_{{\bf{i}}_{q}}=S_{2}S_{1}$, with $S_{1}=P_{i_{k}}\cdots P_{i_{1}}$ and $S_{2}=P_{i_{q}}\cdots P_{i_{k+1}}$. The classical state $|z^{\prime}\rangle$ at position $k$ in the product is given by

[TABLE]

where $|z\rangle$ is the classical state of the current configuration. The state $|z^{\prime}\rangle$ has energy $E_{z^{\prime}}$, and the state $|z\rangle$ has energy $E_{z}$. The new block-swapped configuration is ${\cal C}^{\prime}=\{|z^{\prime}\rangle,S_{1}S_{2}\}$. The multiplicity table of this configuration differs from that of the current configuration by having one fewer $E_{z}$ state and one additional $E_{z^{\prime}}$ state. The weight of the new configuration is then proportional to $e^{-\beta[E_{{\cal C}^{\prime}}]}$ where the multiset $E_{{\cal C}^{\prime}}=E_{{\cal C}}+\{E_{z^{\prime}}\}-\{E_{z}\}$. The acceptance probability is as in Eq. (23) with the aforementioned $E_{{\cal C}^{\prime}}$.

V.3.4 Cycle completion

The moves presented so far have left the number of group elements in the sequence, or expansion order, namely $q$, unchanged. The cycle completion move has the effect of changing the value of $q$. A lookup table of short cycles obeying

[TABLE]

will be helpful in this case. The cycle completion move identifies a sub-cycle in the sequence $S_{{\bf{i}}_{q}}$, e.g., $P_{i_{j}}\cdots P_{i_{k}}$ and replaces it with its complement

[TABLE]

Note that the inverses of permutation matrices are also permutation matrices and are therefore also present in $G$.

For concreteness, let us consider the case of subsequences of length two. We randomly pick a point $k\in[0,q]$ in the sequence. With probability $1/4$, the subsequence is taken to be $P_{0}P_{0}$, $P_{0}P_{i_{k}}$, $P_{i_{k-1}}P_{0}$, $P_{i_{k-1}}P_{i_{k}}$.333Note that if $k=0$ or $q$, namely, the edges of the sequence, then two of the choices correspond to non-starters as there is an operator only to one side of the insertion point and so either $P_{i_{k-1}}$ or $P_{i_{k}}$ are not defined. The identified subsequence is replaced by its complement, resulting in a new configuration $\mathcal{C}^{\prime}$. Because we can interpret $P_{0}^{-1}=P_{0}=P_{j}P_{j}^{-1}$ (for any arbitrary index $j$) and so on, the cycle completion move can grow and shrink the sequence. The acceptance probability is as in Eq. (23) with the new configuration.

V.4 Measurements

Having reviewed the various update moves we next turn to discuss measurements within the algorithm.

V.4.1 Diagonal measurements

A diagonal operator $\Lambda$ obeys $\Lambda|z\rangle=\lambda(z)|z\rangle$ where $\lambda(z)$ is a number that depends both on the operator and the state it acts on. Since $\langle z|\Lambda S_{{\bf{i}}_{q}}|z\rangle=\lambda(z)\langle z|S_{{\bf{i}}_{q}}|z\rangle$, for any given configuration ${\cal C}=(|z\rangle,S_{{\bf{i}}_{q}})$, there is a contribution $\lambda=\lambda(z)$ to the diagonal operator thermal average $\langle\Lambda\rangle$. To improve statistics, one may also consider rotations in (the periodic) imaginary time. To do that, we may consider ‘virtual’ block-swap moves (see Sec. V.3.3) that rotate $S_{{\bf{i}}_{q}}$ and as a result also change the classical configuration from $|z\rangle$ to $|z_{i}\rangle$. The contribution to the expectation value of a diagonal operator $\Lambda$ thus becomes:

[TABLE]

where $E_{{\cal C}_{i}}$ is the energy multiset associated with configuration ${\cal C}_{i}$ whose multiset is $E_{{\cal C}_{i}}=E_{{\cal C}}+\{E_{z_{i}}\}-\{E_{z}\}$ (recall that $z_{0}\equiv z$, so $E_{{\cal C}_{0}}=E_{{\cal C}}$). The normalization factor ${\cal Z}$ above is the sum

[TABLE]

over all nonzero multiplicities $m_{j}$. In the case where $\Lambda=H_{\text{c}}$ the above expression simplifies to:

[TABLE]

V.4.2 Off-diagonal measurements

We next consider the case of measuring the expectation value of an off-diagonal operator $P_{k}$, namely, $\langle P_{k}\rangle$. To do this, we interpret the instantaneous configuration as follows

[TABLE]

where ${\cal C}^{\prime}$ is the configuration associated with the multiset $E_{{\cal C}^{\prime}}=E_{{\cal C}}-\{E_{z}\}$. In the above form, we can reinterpret the weight $W_{{\cal C}}$ as contributing

[TABLE]

to $\langle P_{k}\rangle$ where $d^{(i_{q})}_{z_{q}}$ is the ’hopping strength’ of $P_{k}$ Eq. (9).

As in the case of the diagonal measurements, one can take advantage of the periodicity in the imaginary time direction to improve statistics by rotating the sequence such that any of the elements of $S_{{\bf{i}}_{q}}$ becomes the last element of the sequence (see Sec. V.3.3), weighted accordingly by the block-swap probability. By doing so, $P_{k}$ becomes

[TABLE]

where $E_{{\cal C}_{i}}=E_{{\cal C}}+\{E_{z_{i}}\}-\{E_{z}\}$, the sum $\sum_{j}$ is over all rotated configurations ${\cal C}^{\prime}$.

V.4.3 Products of off-diagonal measurements

The sampling of expectation values of the form $\langle P_{k_{1}}P_{k_{2}}\rangle$ proceeds very similarly to the single operator case except that now both operators must appear at the end of the sequence. The argument proceeds similarly to the single off-diagonal measurement, and we have that the contribution to the expectation value of $\langle P_{k_{1}}P_{k_{2}}\rangle$ is

[TABLE]

with $E_{{\cal C}^{\prime}}=E_{{\cal C}}-\{E_{z},E_{z_{q-1}}\}$. As in the single off-diagonal operator case, we can use the block-swap move to alter the elements at the end of the sequence, and for each pair of adjacent operators in the sequence obtain an improved contribution. By doing so, $\langle P_{k_{1}}P_{k_{2}}\rangle$ becomes

[TABLE]

where $E_{{\cal C}_{k}}=E_{{\cal C}}+\{E_{z_{k}}\}-\{E_{z}\}$, $E_{{\cal C}^{\prime}_{i}}=E_{{\cal C}}-\{E_{z},E_{z_{q-1}}\}$ with $|z^{\prime\prime}\rangle=P_{k_{2}}|z^{\prime}\rangle$ and $|z^{\prime}\rangle$ is the classical state after the block swap. Similar to the single off-diagonal operator case, the sum $\sum_{j}$ is over all rotated configurations ${\cal C}^{\prime}$ whose $S_{{\bf{i}}_{q}}$ ends with $P_{k_{1}}P_{k_{2}}$.

Measurements of thermal averages of products of more than two off-diagonal operators can also be derived in a straightforward manner.

V.4.4 Improved measurements

As will often happen, certain physical operators will have more than one representation as group element. E.g., if $P_{3}=P_{1}P_{2}$, one could measure both the single operator $\langle P_{3}\rangle$ and the operator product $\langle P_{1}P_{2}\rangle$ and combine the results.

VI Results

In this section we present some results that highlight some of the advantages that PMR has to offer over existing methods. Specifically, we compare the performance of PMR over SSE on random 3-regular MAX2SAT instances augmented with a transverse field. This class of instances corresponds to a particular choice of the Ising Hamiltonian given in Eq. (20) whereby each spin is coupled antiferromagnetically (with strength $J_{ij}=1$) with exactly three other spins picked at random. This class of instances is known to exhibit a quantum spin-glass phase transition and is notoriously difficult to simulate by standard QMC techniques, making it suitable to illustrate the strengths of the PMR algorithm (see Ref. Farhi et al. (2012) for more details).

In a subsequent section, we demonstrate the versatility of our formalism by considering the performance of PMR on a transverse-field Ising model augmented with random two-body $XX$ connections, which to the authors’ knowledge cannot be readily implemented using existing methods Mazzola and Troyer (2017).

VI.1 PMR vs SSE: Transverse-field Ising model simulations

To demonstrate the advantages of PMR over existing state-of-the-art, we study random 3-regular MAX2SAT instances augmented with a transverse field. We study the thermal properties by utilizing a parallel tempering scheme Hukushima and Nemoto (1996); Marinari et al. (1998); Hen and Young (2011); Farhi et al. (2012) for both PMR and SSE with 11 replicas at inverse temperatures, $\beta\in\left\{0.1,0.2,0.5,1,2,5,10,20,30,40,50\right\}$. We carry out simulations for 50 random MAX2SAT instances at sizes $N=96,128$ and at transverse field strengths $\Gamma=0.1,0.4$.

To quantify the performance of the algorithms, we fix the total number of measurements of our observables, and we vary the total number of updates between measurements. Thus, we are able to measure the dependence of the estimate of the thermal expectation value with the total time spent de-correlating measurement samples. The performance comparisons of the PMR algorithm against SSE are summarized in Figs. 3 and 4. Both figures depict the thermal average of $x$-magnetization as a function of simulation runtime (other observables exhibit similar behavior).

As is evident from the figures, PMR is four or more orders of magnitude faster than SSE, converging on average after only a few seconds in all cases. On the other hand, for a large fraction of the instances, the SSE simulations did not finish running over the 24 hour window allocated for each run. As expected, the difference in performance is even more pronounced in the more ‘classical’ $\Gamma=0.1$ case (Fig. 4).

VI.2 Ising model with random $XX$ interactions

We next illustrate the utility of our technique by studying a model that likely requires highly non-trivial implementations if studied by other QMC algorithms. We consider a transverse-field Ising model with random $XX$ interactions whose Hamiltonian is given by

[TABLE]

Here, $s$ is a parameter in the range $(0,1)$, and $b\in\{0,1\}$ determines whether the two-body $X$ terms are absent ($b=0$) or present ($b=1$). We examine underlying connectivity graphs that are Erdős–Rényi random, meaning we randomly pick a pair of spins to connect. The total number of edges for each instance is taken to be $nm/2$, where $m\in\left\{3,4,5\right\}$ is the average degree of the graph (we focus only on single component graphs for simplicity).

Hamiltonians of the above general form appear widely in the context of quantum annealing processes Kadowaki and Nishimori (1998), where the system is evolved according to the above Hamiltonian by varying the parameter $s$ slowly in time from $s=0$ to $s=1$. The goal in quantum annealing is for the system to reach a state at the end of the anneal that has considerable overlap with the ground state manifold of the $Z$-dependent ‘problem’ Hamiltonian, which in this case is a MaxCut instance (or a random antiferromagnetic) Farhi et al. (2012). While in standard quantum annealing the two-body $X$ terms are normally absent (i.e., $b=0$), one is often interested in understanding the effects of augmenting the Hamiltonian with a ‘catalyst’ — an extra term that is hoped to reduce the amount of time required for the annealing process to take place (see, e.g., Refs. Durkin (2019); Albash (2019); Crosson et al. (2014); Hormozi et al. (2017); Nishimori and Takada (2017)). Setting $b=1$ can be viewed as an example of such a situation.

For $H$ in Eq. (37), we have $D_{0}=H_{\text{c}}=s\sum_{\langle ij\rangle}Z_{i}Z_{j}$ as well as one-body and two-body (in the $b=1$ case) off-diagonal $P_{j}$ operators: $\{X_{i}\}\cup\{X_{i}X_{j}\}_{\langle ij\rangle}$. The $D_{j}$ operators are all of the form $D_{j}=d_{j}\cdot\mathbb{1}$ where for the one body operators $d_{j}=-(1-s)$ and for the two-body operators $d_{j}=-s(1-s)$. That $d_{j}\leq 0$ implies that the model is stoquastic and hence sign-problem-free. Given a configuration $\{|z\rangle,S_{q}=P_{i_{1}}P_{i_{2}}\ldots P_{i_{q}}\}$, the QMC algorithm generated new configuration by either changing basis state $|z\rangle$ or subsequence of off-diagonal operator $S_{q}$. For our updates, we restrict to subsequences of length two which is enough to ensure ergodicity. The possible completion moves are summarized in Table 1. For illustration, say $P_{i_{1}}=X_{i}X_{k}$ then new configuration $\{|z\rangle,S^{\prime}_{q}=P^{\prime}_{i_{0}}P^{\prime}_{i_{1}}P_{i_{2}}\ldots P_{i_{q}}\}$ where $P^{\prime}_{i_{0}}=X_{i}X_{j}$, $P^{\prime}_{i_{1}}=X_{j}X_{k}$ can be generated by replacing $X_{i}X_{k}\rightarrow X_{i}X_{j},X_{j}X_{k}$ and accepting the move with a probability satisfying detailed balance condition.

By inspecting the results of QMC simulations of the above Hamiltonian we are able to answer a number of questions that are relevant to quantum annealing. We first examine the variance of $H$ (denoted $\sigma^{2}_{H}$), which when close to 0 indicates that the thermal state is close to being purely in the ground state of the system (technically, any energy eigenstate of the system will give $\sigma^{2}_{H}=0$). For sufficiently low temperatures, this will always be the case, but the energy gap and the density of states determine how low the temperature needs to be.

We therefore study the dependence of $\sigma^{2}_{H}$ on the instances’ tree-widths. This is shown in Fig. 5. We see that while instances with different $m$ values may have the same tree-width, in the presence of $XX$ interactions there can be a significant difference in their $\sigma^{2}_{H}$ values. We make several observations. First, we find that the Hamiltonian with $XX$ interactions requires more sweeps in order to de-correlate, i.e. thermalize. Second, the differences for different $m$ are more substantial with $XX$ interactions than without them, indicating that the $XX$ interaction makes the spectrum much more susceptible to $m$. Third, larger $m$ values with $XX$ interactions tend to correspond to lower $\sigma^{2}_{H}$ values. This suggests that in the presence of $XX$ interactions, larger $m$ values are effectively ‘colder’. Finally, we find that the tree-width makes little difference to the $\sigma^{2}_{H}$ values, with or without $XX$ interactions.

Figure 6 shows the average diagonal energy as a function of the annealing parameter $s$ for the $b=0$ (no $XX$) and the $b=1$ case and for two different values of average graph degree, namely, $m=3$ and $m=5$. We find that the presence of the $XX$ catalyst has two important consequences: it minimizes the effect of the graph degree and significantly raises the average value. The former effect suggests that the presence of an $XX$ catalyst will minimize performance differences in solving random MaxCut problems with different graph degrees. The latter effect is not surprising since the presence of both $X$ and $XX$ in the Hamiltonian means that we can expect the eigenstates to remain disordered for a larger region of $s$.

VII Summary and discussion

We presented a parameter-free, Trotter-error free, universal quantum Monte Carlo scheme for the simulation of a broad range of physical models under a single unifying framework. Our technique enables the study of essentially any model on an equal footing. In our approach, the quantum dimension consists of products of elements of permutation groups, allowing us to formulate update rules and measurement schemes independently of the model being studied.

We used our approach to clarify the emergence of the sign problem in the simulation of non-stoquastic physical models. In addition, studying the thermal properties of transverse-field Ising models, we illustrated the advantages of our technique over existing state-of-the-art, specifically the stochastic series expansion (SSE) algorithm. We showed that one of the features of the permutation matrix representation QMC that distinguishes it from SSE is that it bundles infinitely many SSE weights into a single weight. We demonstrated that this translates to orders-of-magnitude runtime advantages in most cases.

We further illustrated the flexibility of our method by studying models with variable localities of interactions and underlying connectivities, namely, a transverse-field Ising model augmented with randomly placed two-body $X$ interactions, which to the authors’ knowledge cannot be readily implemented using existing methods. For these, we found that the presence of the $XX$ interactions can mitigate differences associated with connectivity graphs of different degree. We expect that our algorithm can be further improved by implementing updates that utilize more ‘global’ moves that, e.g., update longer cycles of operators. The extent to which this can translate to further runtime advantages remains an open question that we leave for future work.

We believe that the generality and flexibility of our algorithm will make it a useful tool in the study of physical models that have so far been inaccessible, cumbersome or too large to implement with existing techniques.

Acknowledgements.

Computation for the work described here was supported by the University of Southern California’s Center for High-Performance Computing (http://hpcc.usc.edu) and by ARO grant number W911NF1810227. IH is supported by the U.S. Department of Energy, Office of Science, Office of Advanced Scientific Computing Research (ASCR) Quantum Computing Application Teams (QCATS) program, under field work proposal number ERKJ347. Work by LG was supported by the U.S. Department of Energy (DOE), Office of Science, Basic Energy Sciences (BES) under Award No. DE-SC0020280. Work by TA is supported by the Office of the Director of National Intelligence (ODNI), Intelligence Advanced Research Projects Activity (IARPA), via the U.S. Army Research Office contract W911NF-17-C-0050. The U.S. Government is authorized to reproduce and distribute reprints for Governmental purposes notwithstanding any copyright annotation thereon. The views and conclusions contained herein are those of the authors and should not be interpreted as necessarily representing the official policies or endorsements, either expressed or implied, of the ODNI, IARPA, or the U.S. Government.

Appendix A Finite dimensional permutation matrix representations

To show that any finite-dimensional matrix can be written in the form of Eq. (1), we will make use of ‘cycle notation’ Rotman (1995) — a compact representation of permutations — to represent the permutation matrices $P_{j}$. We start with some terminology.

A cycle is a string of integers that represents an element of the symmetric permutation group $S_{n}$, which cyclically permutes these integers and fixes all other integers. For example cycle $(a_{1},a_{2},\ldots,a_{m})$ is the permutation that sends $a_{i}$ to $a_{i+1}$, $1\leq i\leq m-1$ and sends $a_{m}$ to $a_{1}$. The cycle given in the above example is an $m$-cycle. In general, any element $\sigma\in S_{n}$ can be written as a product of $k$ cycles as $(a_{1}\;a_{2}\ldots a_{m_{1}})(a_{m_{1}+1}\;a_{m_{1}+2}\ldots a_{m_{2}})\ldots$ $(a_{m_{k-1}+1}\;a_{m_{k-1}+2}\ldots a_{m_{k}})$. The order of a permutation $\sigma$ is defined as the smallest positive integer $p$ such that $\sigma^{p}$ is the identity element. In this notation, the action of $\sigma$ on any number from $1$ to $n$ can be determined as follows. If $a$ appears at the right end of one of the $k$ cycles, then $\sigma(a)$ is the integer at the start of the cycle to which $a$ belongs. If an integer $a$ does not appear at the right end of one of the $k$ cycles, then $\sigma(a)$ is the integer to the right of $a$ in the cycle to which $a$ belongs.

For concreteness, let us write the $3\times 3$ permutation matrices used in Sec. II.3 in cycle notation.

[TABLE]

[TABLE]

The identity operation can thus be written as

[TABLE]

We illustrate the evaluation of a product of two cycles by computing $P^{2}$ in cycle notation. Since we defined $P=(1,2,3)$ we have $P_{2}=P^{2}=(1,2,3)(1,2,3)$. By the above definition

$P_{2}(1)=(1,2,3)(1,2,3)(1)=(1,2,3)(2)=(3)$,

$P_{2}(2)=(1,2,3)(1,2,3)(2)=(1,2,3)(3)=(1)$,

$P_{2}(3)=(1,2,3)(1,2,3)(3)=(1,2,3)(1)=(2)$.

Therefore $P_{2}=(1,3,2)$. If we enumerate the basis states as

[TABLE]

then with the action described above one can see that $|k\rangle=P_{i}^{\text{matrix notation}}|j\rangle$ corresponds to $k=P_{i}^{\text{cycle notation}}(j)$.

In the above notation, the set of $n\times n$ permutation matrices can be seen as groups generated by an $n$-cycle. Any given $n$-cycle $\sigma=(a_{0},a_{1},a_{2},a_{3},...,a_{n-1})\in S_{n}$, where $S_{n}$ is the symmetric permutation group, has order $n$. To see why this is so, observe that $\sigma^{k}(a_{0})=a_{k}$ for $0<k<n$ so its order cannot be less than $n$ and $\sigma^{n}$ is the identity as $\sigma^{n}(a_{i})=a_{i}$ for $i\in\{0,n-1\}$. Therefore the group generated by $\sigma$ has $n$ elements. More generally we have: $\sigma^{k}(a_{j})=a_{(j+k)\;mod\;n}$. This implies $\sigma^{k_{1}}(a_{j})\neq\sigma^{k_{2}}(a_{j})$ for $k_{1}\neq k_{2}$.

Let $P$ be the permutation matrix corresponding to $\sigma$. Then, $P^{k_{1}}|a_{j}\rangle\neq P^{k_{2}}|a_{j}\rangle$ for $k_{1}\neq k_{2}$ and basis vector $|a_{j}\rangle$. Since any row (or column) of a permutation matrix has value 0 at all positions but one, where it has the value 1, this implies that no two permutation matrices generated from $P$ have the same row otherwise that would mean $P^{k_{1}}|a_{j}\rangle=P^{k_{2}}|a_{j}\rangle$ for some $k_{1},k_{2},a_{j}$.

Next, we will show that for any given matrix entry $(i,j)$ there is at least one permutation matrix $P^{k}$ having value 1 at $(i,j)$, that is $P^{k}_{i,j}=1$ for some $k$. Let $|i\rangle$ denote the basis vector which has value 1 at the $i$-th index and 0 at all others. Since these permutation matrices form a group there must be some matrix $P^{k}$ such that $P^{k}|j\rangle=|i\rangle$. This however means that $P^{k}$ has value 1 at $(i,j)$. Combining the two statements above, we find that for any entry $(i,j)$ there is exactly one permutation matrix generated from $P$ that has value 1 at that entry. The diagonal matrices $D_{j}$ may be used to convert these 1’s to any desired value. We have thus shown that by choosing our $P_{j}$ permutation matrices as $n$-cycles, one can construct arbitrary Hamiltonians $H=\sum_{j}D_{j}P_{j}$. This proof also provides a prescription as to how to explicitly choose the permutations $P_{j}$.

Appendix B Hermiticity of $H$

Here we show that $H=\sum D_{j}P_{j}$ is hermitian if and only if for every index $j$ there is an associated index $j^{\prime}$ such that $P_{j}=P_{j^{\prime}}^{-1}$ and $D_{j}=D_{j^{\prime}}^{*}$ (in general, $j$ and $j^{\prime}$ may correspond to the same index).

We first prove the if direction. Let $H$ above be a hermitian matrix. We show that this implies that for every index $j$ there is an associated index $j^{\prime}$ such that $P_{j}=P_{j^{\prime}}^{-1}$. Let $P_{j}$ be the permutation sending a basis vector $p$ to another basis vector $q$. Thus in ‘cycle notation’ (see Appendix A) $P_{j}=\ldots(\ldots p,q\ldots)\ldots$. Let us assume that $D_{j}$ is not the zero matrix (otherwise $D_{j}P_{j}$ is trivially zero). Case I: Let ${D_{i}}_{(q,q)}\neq 0$. Since $H$ is hermitian thus there exist $P_{j^{\prime}}=\ldots(\ldots q,p\ldots)\ldots$ in the decomposition of $H$. But then $P_{j}P_{j^{\prime}}=\ldots(\ldots q,q\ldots)\ldots$ which has to be equal to $\mathbb{1}$ (as otherwise this would imply the existence of a permutation matrix that has a fixed point contradictory to our initial setup). Case II: Let ${D_{j}}_{(q,q)}=0$. Since $D_{j}$ is not identically zero, there exists an element $s$ such that ${D_{j}}_{(s,s)}\neq 0$. Let us denote by $r$ the element that is sent to $s$ in $P_{j}$, that is $P_{j}=\ldots(\ldots p,q\ldots r,s\ldots)\ldots$. Again since $H$ is hermitian there must exist $P_{j^{\prime}}=\ldots(\ldots s,r\ldots)\ldots$ in the decomposition of $H$. But then $P_{j}P_{j^{\prime}}=\ldots(\ldots s,s\ldots)\ldots$ which has to be equal to $\mathbb{1}$. Thus $P_{j}$ has an inverse $P_{j^{\prime}}$ in the decomposition of $H$.

Next, we prove that there can be no two permutation matrices in $H$ with the same nonzero element. In cycle notation, this assertion translates to the assertion that there can be no two distinct permutations that send a basis state to the same basis state. We prove this by contradiction. Let $P_{j}$ and $P_{k}$ be two distinct permutations both of which send a basis vector $p$ to $q$. In ’cycle notation’ this mean $P_{j}=\ldots(\ldots p,q\ldots)\ldots$ and $P_{k}=\ldots(\ldots p,q\ldots)\ldots$ with $P_{j}\neq P_{k}$. We showed above that $P_{j}$ has an inverse $P_{j^{\prime}}$, that is $P_{j}P_{j^{\prime}}=\mathbb{1}$. But $P_{j^{\prime}}P_{k}=\ldots(\ldots p,p\ldots)\ldots$ has fixed point and cannot be the identity due to the uniqueness of the inverse. Thus we have reached a contradiction.

We now prove the full if part. By definition, that $H$ is hermitian implies $\sum D_{j}P_{j}=\sum D_{j}^{*}{P_{j}}^{-1}$. Let $P_{j^{\prime}}$ be the inverse of $P_{j}$ which as we proved should exist in the decomposition with a nonzero weight $D_{j}$. Equating the right-hand side and the left-hand side of the equality gives $D_{j}=D_{j^{\prime}}^{*}$.

Proving the other direction is simpler. We assume that in the summation $H=\sum D_{j}P_{j}$ there is for every index $j$ an associated index $j^{\prime}$ such that $P_{j}={P_{j^{\prime}}}^{-1}$ and $D_{j}={D_{j^{\prime}}}^{*}$. Thus $H^{\dagger}=\sum{D_{j^{\prime}}}^{*}{P_{j^{\prime}}}^{-1}=\sum{D_{j}}{P_{j}}=H$.