# A rapidly mixing Markov chain from any gapped quantum many-body system

Sergey Bravyi 111IBM Quantum, IBM T.J. Watson Research Center, Yorktown
Heights, USA Giuseppe Carleo 222École Polytechnique Fédérale de Lausanne
(EPFL), Institute of Physics, CH-1015 Lausanne, Switzerland David Gosset
333Institute for Quantum Computing, University of Waterloo, Canada
444Department of Combinatorics and Optimization, University of Waterloo,
Canada 555Perimeter Institute for Theoretical Physics, Canada Yinchen Liu
††footnotemark:  ††footnotemark:  ††footnotemark:

###### Abstract

We consider the computational task of sampling a bit string $x$ from a
distribution $\pi(x)=|\langle x|\psi\rangle|^{2}$, where $\psi$ is the unique
ground state of a local Hamiltonian $H$. Our main result describes a direct
link between the inverse spectral gap of $H$ and the mixing time of an
associated continuous-time Markov Chain with steady state $\pi$. The Markov
Chain can be implemented efficiently whenever ratios of ground state
amplitudes $\langle y|\psi\rangle/\langle x|\psi\rangle$ are efficiently
computable and the starting state of the chain satisfies a mild technical
condition that can be efficiently checked. This extends a previously known
relationship between sign-problem free Hamiltonians and Markov chains. The
tool which enables this generalization is the so-called fixed-node Hamiltonian
construction, previously used in Quantum Monte Carlo simulations to address
the fermionic sign problem. We implement the proposed sampling algorithm
numerically and use it to sample from the ground state of Haldane-Shastry
Hamiltonian with up to $56$ qubits. We observe empirically that our Markov
chain based on the fixed-node Hamiltonian mixes more rapidly than the standard
Metropolis-Hastings Markov chain.

The task of generating samples from a given probability distribution underlies
almost all randomized algorithms used in computational physics, machine
learning, and optimization. In many applications the target distribution is
efficiently computable in the sense that there is a polynomial-time subroutine
that computes the relative probabilities of any two given elements. A
paradigmatic example is the (classical) Boltzmann distribution $\pi(x)\sim
e^{-E(x)/T}$ associated with an efficiently computable energy function $E(x)$
and a temperature $T$. In this case the ratio of Boltzmann probabilities
$\pi(y)/\pi(x)$ is simply related to the energy difference $E(y)-E(x)$. The
fundamental obstacle in such cases is that distributions that are efficiently
computable may still be challenging to sample from. In particular, there are
many examples when generating a sample from an efficiently computable
distribution is an NP-hard problem. This includes low-temperature Boltzmann
distributions with energy function $E(x)$ that is given by a 3-SAT formula or
an Ising spin glass [1] and certain distributions described by neural networks
of RBM type [2].

It is therefore interesting to ask:

Q1:Which efficiently computable distributions can also be efficiently sampled?

and, more broadly,

Q2:Which distributions admit an efficient reduction from sampling to computing
probabilities?

Here we note that a polynomial-time reduction between the two tasks may exist
for a broader family of distributions, some of which may not be efficiently
computable.

A natural way to address these questions is to use Markov Chain Monte Carlo
(MCMC), an empirically successful and versatile algorithmic tool for sampling
probability distributions. MCMC methods work by constructing a Markov chain
$M$ such that the target distribution $\pi$ is the unique steady distribution
of $M$. A distribution $\pi_{t}$ generated after implementing $t$ steps of the
chain $M$ approximates the steady distribution $\pi$ provided that the number
of steps $t$ is large compared with the mixing time of $M$. A well known
sampling algorithm in this category is the Metropolis-Hastings Markov chain
[3] and its variations. Although MCMC methods are widely used in practice,
their main limitation is the difficulty of obtaining rigorous upper bounds on
the mixing time of Markov chains. Such bounds can be established only in
certain special cases using techniques such as the canonical paths method,
coupling of Markov chains, or the conductance bound, see e.g. [4].

In this work we are interested in variants of the questions Q1 and Q2 for
probability distributions originating from ground states of quantum many-body
systems. Can we design specialized MCMC sampling algorithms that exploit their
structure? Are efficiently computable distributions that arise from ground
states efficiently samplable?

In particular, we consider a system of $n$ qubits with few-qubit interactions
described by a $k$-local Hamiltonian $H=\sum_{a=1}^{m}H_{a}$, where $H_{a}$
are Hermitian operators acting non-trivially on subsets of at most $k=O(1)$
qubits. We choose the energy scale such that $\|H_{a}\|\leq 1$ for all $a$.
Let $\psi$ be the ground state of $H$. We assume that the ground state of $H$
is non-degenerate and separated from excited states by an energy gap
$\gamma>0$. Our goal is to sample a bit string $x\in\\{0,1\\}^{n}$ from the
ground state distribution

$\pi(x)=|\langle x|\psi\rangle|^{2}.$ (1)

It can be viewed as the zero-temperature quantum analogue of a classical
Boltzmann distribution describing classical spins with few-spin interactions.

To accomplish this goal one may construct a suitable quantum-to-classical
mapping that converts a quantum Hamiltonian with the ground state $\psi$ to a
classical Markov chain $M$ with the steady distribution $\pi(x)=|\langle
x|\psi\rangle|^{2}$. A possible choice for $M$ is a Metropolis-Hastings (MH)
Markov chain with local updates. In a simple, typical setting, each step of
the MH chain flips a randomly chosen bit (or a subset of bits) of $x$ to
propose a candidate next state $y$. The proposed state $y$ is accepted with
the probability $\min\\{1,\pi(y)/\pi(x)\\}$ to ensure the detailed balance
condition. As was shown in Ref. [5], the mixing time of the MH chain can be
upper bounded as

$T\leq\tilde{O}(n^{k}s/\gamma),$ (2)

where $\tilde{O}$ notation hides certain logarithmic factors and $s$ is a
sensitivity parameter defined as $s=\max_{x\neq y}|\langle y|H|x\rangle\langle
y|\psi\rangle/\langle x|\psi\rangle|$. For the widely studied family of sign-
problem free666A Hamiltonian $H$ with real matrix elements is sign-problem
free if $\langle x|H|y\rangle\leq 0$ for all $x\neq y$. Hamiltonians the
sensitivity parameter can be bounded as [5]

$s\leq 2\|H\|\quad\qquad\text{(if $H$ is sign-problem free)}.$ (3)

For a local Hamiltonian of the type we consider we have
$\|H\|=O(\mathrm{poly}(n))$ and so Eqs. (2,3) give a polynomial upper bound on
the mixing time of the MH chain. Since each step of the chain makes use of a
ratio $\pi(y)/\pi(x)$, this constitutes an efficient reduction from sampling
to computing (ratios of) probabilities. In other words, we obtain a positive
answer to Q2 for sign-problem free Hamiltonians 777We note that a different
efficient reduction from sampling to computing probabilities was previously
known for this family of Hamiltonians [6, 7].. If we further specialize to
frustration-free and sign-problem free Hamiltonians then the requisite ratios
of probabilities can be computed efficiently [7] and we get a partial answer
to Q1 as well.

Unfortunately, for more general local Hamiltonians—those that may not be sign-
problem free— it is unknown whether the sensitivity $s$ admits a polynomial
upper bound in terms of $n$ and $1/\gamma$, and significant differences may
thwart this approach altogether. In the sign-problem free case the Perron
Frobenius theorem implies that the ground state amplitudes $\langle
x|\psi\rangle$ are nonnegative in the computational basis. For more general
local Hamiltonians the entries of the ground state wavefunction may have
nontrivial relative phases, but the MH chain ignores this information. Indeed,
this chain only depends on the ratios of probabilities $\pi(y)/\pi(x)$. For
these reasons it appears unlikely that the MH chain can be shown to mix
rapidly for general gapped local Hamiltonians.

To proceed, we use a quantum-to-classical mapping based on a proposal from
Ref. [8] which was originally introduced to circumvent the fermionic sign
problem in Quantum Monte Carlo simulations [9, 10]. This construction maps any
$k$-local $n$-qubit Hamiltonian $H$ with real matrix elements, unique ground
state $\psi$, and a spectral gap $\gamma$ to a new $n$-qubit Hamiltonian
$F$—defined in Eq. (12)—such that: (i) $F$ and $H$ have the same ground state
$\psi$, (ii) the spectral gap of $F$ is at least $\gamma$, and (iii) the
Hamiltonian $F$ is sign-problem free, modulo a simple basis change, see
Section 1.1 for details. Following Ref. [8] we refer to $F$ as a fixed-node
Hamiltonian. Importantly, matrix elements of $F$ are efficiently computable
given an efficient subroutine for computing the ratio of amplitudes

$\qquad\qquad\qquad\qquad\qquad x,y\longrightarrow\langle
y|\psi\rangle/\langle x|\psi\rangle.\qquad(\textbf{Amplitude computation
subroutine})$

We are thus led to consider natural variants of the questions
$\textbf{Q1},\textbf{Q2}$ where probabilities are replaced by amplitudes. This
allows us to exploit the information encoded in the amplitudes’ relative
phase, which is an important feature of the ground state. We note that the
model where a quantum state can be accessed via the amplitude computation
subroutine has been recently studied in [11].

Since the fixed node Hamiltonian is sign-problem free, we can sample the
ground state distribution $\pi(x)$ using the MH chain and upper bound the
mixing time of the chain using Eqs. (2,3). Unfortunately this is not useful
because, as we shall see below, the norm of the fixed-node Hamiltonian can be
unbounded.

Instead, we introduce a quantum-to-classical mapping that yields a continuous-
time Markov chain rather than a regular discrete-time chain. Recall that a
continuous-time Markov chain (CTMC) with a state space ${\cal S}$ defines a
family of probability distributions $\pi_{t}(x)$ where $t\geq 0$ is the
evolution time and $x\in{\cal S}$ is the state reached at time $t$. The time
evolution of $\pi_{t}(x)$ is governed by a differential equation

$\frac{d}{dt}\pi_{t}(x)=\sum_{y\in{\cal S}}\langle x|G|y\rangle\pi_{t}(y),$
(4)

where $G$ is a generator matrix. Rows and columns of $G$ are labeled by
elements of ${\cal S}$. A matrix element $\langle x|G|y\rangle$ with $x\neq y$
can be viewed as the rate of transitions from $y$ to $x$. Accordingly, all
off-diagonal elements of $G$ must be non-negative. The normalization condition
$\sum_{x\in{\cal S}}\pi_{t}(x)=1$ is satisfied for all $t\geq 0$ as long as
each column of $G$ sums to zero. A solution of Eq. (4) has the form
$\pi_{t}(x)=\langle x|e^{Gt}|x_{in}\rangle$, where $x_{in}\in{\cal S}$ is the
starting state at time $t=0$ and $e^{Gt}$ denotes the matrix exponential.

Our CTMC is defined by a generator $G$ which is a suitably rescaled version of
the fixed-node Hamiltonian $F$ associated with $H$. By design, it has a steady
distribution $\pi(x)=\lim_{t\to\infty}\pi_{t}(x)$. It is given by

$\langle x|G|y\rangle=\max\\{0,-\langle x|H|y\rangle\langle
x|\psi\rangle/\langle y|\psi\rangle\\}$ (5)

for $x\neq y$. Here and below we assume for simplicity that

$\pi(x)>0\qquad\text{for all}\qquad x\in\\{0,1\\}^{n}.$ (6)

Note that Eq. (5) also determines the diagonal matrix elements of $G$, since
each column of $G$ sums to zero (due to the normalization condition). Our main
results are as follows.

######  Theorem 1 (Rapid mixing).

Let $H$ be a $k$-local $n$-qubit Hamiltonian with real matrix elements in the
standard basis, unique ground state $\psi$, and a spectral gap $\gamma$. Then
a continuous-time Markov chain with the state space $\\{0,1\\}^{n}$ and a
generator matrix $G$ defined in Eq. (5) has a unique steady distribution
$\pi(x)=|\langle x|\psi\rangle|^{2}$ and obeys

$\|\pi_{t}-\pi\|_{1}\leq\frac{e^{-\gamma t}}{\sqrt{\pi(x_{in})}}$ (7)

for any $t\geq 0$ and any starting state $x_{in}\in\\{0,1\\}^{n}$. Here
$\pi_{t}(x)=\langle x|e^{Gt}|x_{in}\rangle$ is the distribution achieved by
the Markov chain at a time $t$.

As we show below (Lemma 1), the restriction to Hamiltonians $H$ with real
matrix elements is not essential and can be avoided by adding one ancillary
qubit. Theorem 1 shows that we may approximately sample from $\pi$ by running
the continuous-time Markov chain for a total time
$T\sim\gamma^{-1}\log(\pi(x_{in})^{-1})$. However it is not immediately clear
how to simulate this process using resources polynomial in $T$, because of the
significant caveat that the norm of $G$ may be large. This may lead to many
transitions of the chain occuring in a very short time, and it prevents us
from approximating the continuous-time chain by a discrete-time one obtained
by naively discretizing the interval $[0,T]$.

Our saving grace is that we are able to establish a mild upper bound on the
mean number of transitions of the chain $G$ within a given interval, when the
starting state is sampled from the steady distribution. This allows us to
directly simulate the Markov chain of Theorem 1 using a truncated version of
the well-known Gillespie’s algorithm [12] in which we impose an upper limit on
the total number of transitions of the chain. In this way we obtain the
following result.

######  Theorem 2 (Ground state sampling).

Let $\pi^{*}=\min_{x}\pi(x)$. There exists a classical randomized algorithm
that takes as input a precision $\epsilon>0$, a starting state
$x_{in}\in\\{0,1\\}^{n}$, makes at most

$T=\frac{O(\epsilon^{-1}n^{2k}\|H\|)}{\gamma}\log{\left(\frac{1}{\epsilon\pi^{*}}\right)}$

calls to the amplitude computation subroutine, and either outputs a bit string
$y\in\\{0,1\\}^{n}$ or declares an error. Let
$S_{\epsilon}\subseteq\\{0,1\\}^{n}$ be the set of starting states $x_{in}$
for which the algorithm declares an error with probability at most
$\epsilon/4$. The set $S_{\epsilon}$ is nonempty. Moreover, if the algorithm
is run with starting state $x_{in}\in S_{\epsilon}$ and does not declare an
error, then its output $y$ is sampled from a distribution $\epsilon$-close to
$\pi$.

The aforementioned caveat that $G$ has large matrix elements is also direcly
related to the additional requirement above that we are provided with a good
starting state $x_{in}$. Strictly speaking, Theorem 2 falls short of giving an
efficient reduction from sampling to computing amplitudes of the ground state,
since it requires this extra input. However, a good starting state $x_{in}$
can at least be verified using polynomial resources (and the amplitude
computation subroutine): given $x\in\\{0,1\\}^{n}$ and $\epsilon>0$ we can
decide with high probability whether or not $x\in S_{\epsilon}$ by running the
above algorithm $O(\epsilon^{-2})$ times with starting state $x$, and using
the results to compute an estimate of the probability that the algorithm
declares an error.

In Section 1 we describe the quantum-to-classical mapping in detail and we
prove Theorems 1 and 2. Then in Section 2 we demonstrate our algorithm for a
concrete example (the Haldane-Shastry spin chain) where the amplitudes of the
ground state are efficiently computable, and we compare our approach with the
Metropolis-Hastings algorithm.

##  1 Quantum-to-classical mapping

In this section we prove Theorems 1 and 2.

In the following, for any matrix $M$ with real eigenvalues we write
$\lambda_{i}(M)$ for the $i$-th smallest eigenvalue of $M$. If $M$ is an
$n$-qubit operator then $i=1,2,\ldots,2^{n}$. As described in the
introduction, we shall consider an $n$-qubit, $k$-local Hamiltonian with a
unique ground state $|\psi\rangle$ and we are interested in sampling from the
distribution $\pi(x)=|\langle x|\psi\rangle|^{2}$. We shall assume Eq. (6)
holds, i.e., $|\langle x|\psi\rangle|>0$ for all $x\in\\{0,1\\}^{n}$.

###  1.1 The continuous-time Markov chain

We first establish the following Lemma that shows we may restrict our
attention to Hamiltonians with real matrix elements in the standard basis.

######  Lemma 1 (Reduction to real Hamiltonians).

Let $H$ be a local Hamiltonian with unique ground state $\psi$ and spectral
gap $\gamma>0$, satisfying Eq. (6). There is a $O(n^{k})$-sparse $(n+1)$-qubit
Hamiltonian $H_{\mathrm{R}}$ with unique ground state
$|\phi\rangle=\mathrm{Re}(|\psi\rangle)|0\rangle+\mathrm{Im}(|\psi\rangle)|1\rangle$
and spectral gap at least $\mathrm{min}\\{1,\gamma\\}$. The $j$th nonzero
entry of $H_{\mathrm{R}}$ in a given row can be computed using one call to the
amplitude computation subroutine and efficient classical computation.

###### Proof.

Suppose $H$ is an $n$-qubit, $k$-local Hamiltonian with a unique ground state
$\psi$ such that $\langle x|\psi\rangle\neq 0$ for all $x$. Let us fix the
global phase of $\psi$ such that $\langle 0^{n}|\psi\rangle\in\mathbb{R}$.

We may write

$H=A+i(B-B^{T})$

for real matrices $A,B$ which have at most $O(n^{k})$ entries in each row. Now
adjoin one ancilla qubit and consider the Hamiltonian

$H^{\prime}=A\otimes I+(B-B^{T})\otimes\left(\begin{array}[]{cc}0&-1\\\
1&0\end{array}\right).$

Note that $H^{\prime}$ is a real symmetric matrix. Let us suppose that
$|\phi\rangle$ is an eigenvector of $H$ with eigenvalue $\lambda$. Write
$|\phi\rangle=|a\rangle+i|b\rangle$ where $a,b$ are real vectors, and consider
the $n+1$-qubit states

$\displaystyle|\phi_{0}\rangle$
$\displaystyle=|a\rangle|0\rangle+|b\rangle|1\rangle$ (8)
$\displaystyle|\phi_{1}\rangle$
$\displaystyle=-|b\rangle|0\rangle+|a\rangle|1\rangle.$ (9)

Observe that $|\phi_{0}\rangle$ is orthogonal to $|\phi_{1}\rangle$ and that
each of these states is an eigenstate of $H^{\prime}$ with eigenvalue
$\lambda$. By letting $|\phi\rangle$ range over all $2^{n}$ eigenvectors of
$H$ we obtain a complete set of $2\cdot 2^{n}$ eigenvectors for $H^{\prime}$.
In particular, the Hamiltonian $H^{\prime}$ has the same spectral gap as $H$
but has a groundspace spanned by two eigenstates $\psi_{0},\psi_{1}$.

In order to get rid of this degeneracy while preserving the spectral gap, we
add a suitable positive semidefinite term to $H^{\prime}$. For each
$x\in\\{0,1\\}^{n}$ write the complex phase of $\langle x|\psi\rangle$ as

$\frac{\langle x|\psi\rangle}{|\langle
x|\psi\rangle|}=\cos(\theta_{x})+i\sin(\theta_{x})$

where $\theta_{x}\in\mathbb{R}$. Then consider a Hamiltonian

$H^{\prime\prime}=H^{\prime}+\sum_{x\in\\{0,1\\}^{n}}|x\rangle\langle
x|\otimes|v_{x}\rangle\langle v_{x}|\quad\text{ where
}\quad|v_{x}\rangle=-\sin(\theta_{x})|0\rangle+\cos(\theta_{x})|1\rangle.$

Note that $H^{\prime\prime}$ has $O(n^{k})$ nonzero entries in each row and
each one can be computed efficiently using one call to the amplitude
computation subroutine888In particular, the angle $\theta_{x}$ can be obtained
from the complex phase of the ratio $\langle x|\psi\rangle/\langle
0^{n}|\psi\rangle$..

Now observe that

$\langle x|\langle v_{x}|\psi_{0}\rangle=-\sin(\theta_{x})\mathrm{Re}(\langle
x|\psi\rangle)+\cos(\theta_{x})\mathrm{Im}(\langle x|\psi\rangle)=0$

where we used the definition of $\theta_{x}$. Therefore
$H^{\prime\prime}|\psi_{0}\rangle=H^{\prime}|\psi_{0}\rangle=\lambda_{1}(H)|\psi_{0}\rangle$.
Moreover,

$\langle x|\langle v_{x}|\psi_{1}\rangle=+\sin(\theta_{x})\mathrm{Im}(\langle
x|\psi\rangle)+\cos(\theta_{x})\mathrm{Re}(\langle x|\psi\rangle)=|\langle
x|\psi\rangle|$

and therefore

$\left(\sum_{x\in\\{0,1\\}^{n}}|x\rangle\langle x|\otimes|v_{x}\rangle\langle
v_{x}|\right)|\psi_{1}\rangle=\sum_{x}|x\rangle|v_{x}\rangle|\langle
x|\psi\rangle|=\sum_{x}|x\rangle(\langle x|\otimes
I)|\psi_{1}\rangle=|\psi_{1}\rangle$

From this we see that $\psi_{1}$ is an eigenvector of $H^{\prime\prime}$ with
eigenvalue $\lambda_{1}(H)+1$. Since $H^{\prime\prime}\geq H^{\prime}$, all
other eigenvalues of $H^{\prime\prime}$ are at least $\lambda_{2}(H)$.
Therefore $H^{\prime\prime}$ has unique ground state $\psi_{0}$ and its
spectral gap is at least as large as

$\min\\{1,\lambda_{2}(H)-\lambda_{1}(H)\\}=\min\\{1,\gamma\\}.$

∎

###### Proof of Theorem 1.

Our main technical tool is the so-called effective fixed-node Hamiltonian
proposed by Ceperley et al [8]. It can be viewed as a method of “curing” the
sign problem in Quantum Monte Carlo simulations. The method is applicable
whenever amplitudes of the ground state can be efficiently computed.

As discussed above, in light of Lemma 1 we shall assume without loss of
generality that $H$ has real matrix elements and a unique ground state with
real entries in the standard basis. Define sets

$S^{+}=\\{(x,y)\,:\,x\neq y\quad\mbox{and}\quad\langle\psi|x\rangle\langle
x|H|y\rangle\langle y|\psi\rangle>0\\},$ (10) $S^{-}=\\{(x,y)\,:\,x\neq
y\quad\mbox{and}\quad\langle\psi|x\rangle\langle x|H|y\rangle\langle
y|\psi\rangle\leq 0\\}.$ (11)

Here $x,y\in\\{0,1\\}^{n}$ are basis states. Define a fixed-node Hamiltonian
$F$ with matrix elements

$\langle x|F|y\rangle=\left\\{\begin{array}[]{rcl}0&\mbox{if}&(x,y)\in
S^{+},\\\ \langle x|H|y\rangle&\mbox{if}&(x,y)\in S^{-},\\\ &&\\\
\displaystyle\langle x|H|x\rangle+\sum_{(x,z)\in S^{+}}\langle
x|H|z\rangle\frac{\langle z|\psi\rangle}{\langle
x|\psi\rangle}&\mbox{if}&x=y.\\\ \end{array}\right.$ (12)

Note that $F$ is stoquastic (sign problem free) modulo a change of basis
$|x\rangle\to\mathrm{sign}(\langle x|\psi\rangle)|x\rangle$. The following
lemma is largely based on Ref. [8].

######  Lemma 2 (Fixed-node Hamiltonian).

The Hamiltonians $F$ and $H$ have the same unique ground state $\psi$ and the
same ground energy. The spectral gap of $F$ is at least as large as the one of
$H$.

###### Proof.

First we claim that

$F|\psi\rangle=H|\psi\rangle.$ (13)

Indeed,

$\langle x|F|\psi\rangle=\left(\langle x|H|x\rangle+\sum_{(x,z)\in
S^{+}}\langle x|H|z\rangle\frac{\langle z|\psi\rangle}{\langle
x|\psi\rangle}\right)\langle x|\psi\rangle+\sum_{(x,y)\in S^{-}}\langle
x|H|y\rangle\langle y|\psi\rangle=\langle x|H|\psi\rangle$

for any basis state $x$ which proves Eq. (13). Next we claim that

$\langle\phi|F|\phi\rangle\geq\langle\phi|H|\phi\rangle\quad\mbox{for any
state $\phi$}.$ (14)

Indeed,

$\langle\phi|F-H|\phi\rangle=\sum_{(x,y)\in S^{+}}\langle
x|H|y\rangle\left(\frac{\langle y|\psi\rangle|\langle
x|\phi\rangle|^{2}}{\langle x|\psi\rangle}-\langle\phi|x\rangle\langle
y|\phi\rangle\right)$

Using the definition of $S^{+}$ one gets

$\langle\phi|F-H|\phi\rangle=\sum_{(x,y)\in S^{+}}|\langle
x|H|y\rangle|\left(\frac{|\langle y|\psi\rangle|\cdot|\langle
x|\phi\rangle|^{2}}{|\langle x|\psi\rangle|}-s(x,y)\langle\phi|x\rangle\langle
y|\phi\rangle\right),$

where $s(x,y)=\mathrm{sign}(\langle x|H|y\rangle)$. This is equivalent to

$\langle\phi|F-H|\phi\rangle=\frac{1}{2}\sum_{(x,y)\in S^{+}}|\langle
x|H|y\rangle|\left|\sqrt{\frac{|\langle y|\psi\rangle|}{|\langle
x|\psi\rangle|}}\langle x|\phi\rangle-s(x,y)\sqrt{\frac{|\langle
x|\psi\rangle|}{|\langle y|\psi\rangle|}}\langle y|\phi\rangle\right|^{2}\geq
0.$

We claim that the smallest and second-smallest eigenvalues of $H$ and $F$
satisfy

$\lambda_{1}(F)=\lambda_{1}(H)\quad\mbox{and}\quad\lambda_{2}(F)\geq\lambda_{2}(H).$
(15)

Indeed, from Eq. (13) one infers
$\lambda_{1}(F)\leq\langle\psi|F|\psi\rangle=\langle\psi|H|\psi\rangle=\lambda_{1}(H)$.
Suppose $\phi$ is a ground state of $F$. From Eq. (14) one gets
$\lambda_{1}(F)=\langle\phi|F|\phi\rangle\geq\langle\phi|H|\phi\rangle\geq\lambda_{1}(H)$.
Thus $\lambda_{1}(H)=\lambda_{1}(F)$ and $\psi$ is a ground state of both $H$
and $F$. Let $\phi$ be an eigenvector of $F$ orthogonal to $\psi$ such that
$F|\phi\rangle=\lambda_{2}(F)|\phi\rangle$. Since $H$ has the unique ground
state $\psi$, one has
$\lambda_{2}(H)\leq\langle\phi|H|\phi\rangle\leq\langle\phi|F|\phi\rangle=\lambda_{2}(F)$.
Here the second inequality uses Eq. (14). Thus
$\lambda_{2}(F)\geq\lambda_{2}(H)$. ∎

The next step is to convert the fixed-node Hamiltonian $F$ to a Markov chain
with the state space $\\{0,1\\}^{n}$ such that the ground state distribution
$\pi$ is the unique steady state of the chain. One technical difficulty that
prevents us from applying the standard quantum-to-classical mapping commonly
used in Quantum Monte Carlo simulations (see e.g. Section 8 of Ref. [6]) is
that the norm of $F$ may be unbounded. Indeed, a diagonal matrix element
$\langle x|F|x\rangle$ depends on quantities $\langle x|H|z\rangle\langle
z|\psi\rangle/\langle x|\psi\rangle$ which can be arbitrarily large even if
the original Hamiltonian $H$ has a bounded norm. Instead, we shall convert $F$
to a generator matrix describing a continuous-time Markov chain. We will see
that the latter can be simulated efficiently using the well-known Gillespie’s
algorithm [12]. In the proof of Theorem 2 below we establish that the average
number of iterations in Gillespie’s algorithm can be bounded by a quantity
that depends only on the spectral gap of $F$ and a certain “off-diagonal norm”
of $F$ which is at most $n^{k}\|H\|$ even if diagonal matrix elements of $F$
are unbounded.

Define a continuous-time Markov chain with the state space $\\{0,1\\}^{n}$ as
a family of stochastic matrices of the form $e^{Gt}$, where we use the matrix
exponential, $t\geq 0$ is the evolution time, and $G$ is a generator matrix of
size $2^{n}\times 2^{n}$. By definition, the probability that the chain
evolved for time $t$ makes a transition from a state $x$ to a state $y$ is
given by $\langle y|e^{Gt}|x\rangle$. Equivalently, the probability of a
transition from $x$ to $y\in\\{0,1\\}^{n}\setminus\\{x\\}$ between time $t$
and $t+dt$ in the limit $dt\to 0$ is given by $\langle y|G|x\rangle dt$. A
valid generator matrix $G$ must have non-negative off-diagonal matrix elements
and each column of $G$ must sum to zero. This ensures that $e^{Gt}$ is a
stochastic matrix for any $t\geq 0$. A probability distribution $\eta$ is a
steady state of the chain iff

$\sum_{x}\langle y|G|x\rangle\eta(x)=0$ (16)

for any $y$. Equivalently, $\sum_{x}\langle y|e^{Gt}|x\rangle\eta(x)=\eta(y)$
for all $y$ and all $t\geq 0$.

Let $F$ be the fixed-node Hamiltonian constructed above. Define a generator
matrix $G$ such that

$\langle x|G|y\rangle=\lambda_{1}(F)\delta_{x,y}-\langle
x|F|y\rangle\frac{\langle x|\psi\rangle}{\langle y|\psi\rangle}$ (17)

for all states $x,y$. Using Eq. (12) we see that $G$ can equivalently be
expressed as in Eq. (5). Note that $\lambda_{1}(F)$ can be efficiently
computed by making $O(n^{k})$ calls to the amplitude computation subroutine.
Indeed, Lemma 2 implies that $\lambda_{1}(F)=\lambda_{1}(H)$. Furthermore,

$\lambda_{1}(H)=\frac{\langle x|H|\psi\rangle}{\langle
x|\psi\rangle}=\sum_{y}\langle x|H|y\rangle\cdot\frac{\langle
y|\psi\rangle}{\langle x|\psi\rangle}$

for any basis state $x$. It remains to note that each row of $H$ has at most
$O(n^{k})$ non-zeros.

The following lemma implies that $G$ is a valid generator matrix for a
continuous-time Markov chain with the unique steady state $\pi$.

######  Lemma 3 (Generator matrix).

The generator matrix $G$ has real non-positive eigenvalues. Its largest and
second-largest eigenvalues are $0$ and $-\gamma_{F}$ respectively, where
$\gamma_{F}$ is the spectral gap of the fixed-node Hamiltonian $F$. The matrix
exponential $e^{Gt}$ is a stochastic matrix for any $t\geq 0$. The
distribution $\pi$ is the unique steady state of $e^{Gt}$.

###### Proof.

Let $D$ be a diagonal matrix such that $\langle x|D|x\rangle=\langle
x|\psi\rangle$. By definition, $G=D(\lambda_{1}(F)I-F)D^{-1}$. Thus
eigenvalues of $G$ coincide with eigenvalues of $\lambda_{1}(F)I-F$. This
proves the first and the second claims of the lemma.

Let us check that $e^{Gt}$ is a stochastic matrix. It follows directly from
the definitions that $G$ has non-negative off-diagonal elements. Thus all
matrix elements of $e^{Gt}$ are non-negative. Any column of $G$ sums to zero
since

$\sum_{x}\langle
x|G|y\rangle=\lambda_{1}(F)-\frac{\langle\psi|F|y\rangle}{\langle
y|\psi\rangle}=\lambda_{1}(F)-\frac{\lambda_{1}(F)\langle\psi|y\rangle}{\langle\psi|y\rangle}=0.$

Here we noted that $\langle\psi|F=\lambda_{1}(F)\langle\psi|$ due to Lemma 2.
This implies that any column of $e^{Gt}$ sums to one, that is, $e^{Gt}$ is a
stochastic matrix. Finally,

$\displaystyle\sum_{x}\langle y|G|x\rangle\pi(x)$
$\displaystyle=\sum_{x}\langle
y|D(\lambda_{1}(F)I-F)D^{-1}|x\rangle\pi(x)=\langle y|\psi\rangle\cdot\langle
y|\lambda_{1}(F)I-F|\psi\rangle=0$

for any state $y$. Thus $\pi$ is a steady state. This is a unique steady state
since the zero eigenvalue of $G$ is non-degenerate. ∎

To complete the proof of Theorem 1 we now establish the bound Eq. (7).

Our proof strategy closely follows Ref. [13], see Proposition 3 thereof. First
we claim the Markov chain $e^{Gt}$ obeys the detailed balance condition

$\pi(z)\langle y|e^{Gt}|z\rangle=\pi(y)\langle z|e^{Gt}|y\rangle$ (18)

for all states $y,z$. Indeed, let $D$ be a diagonal matrix such that $\langle
z|D|z\rangle=\langle z|\psi\rangle$ for all $z$. We have

$G=DMD^{-1}\quad\mbox{where}\quad M:=\lambda_{1}(F)I-F.$

Note that $M$ is a symmetric matrix since the fixed-node Hamiltonian $F$ is
symmetric. The identity $e^{Gt}=De^{Mt}D^{-1}$ gives

$\pi(z)\langle y|e^{Gt}|z\rangle=\langle y|\psi\rangle\langle
z|\psi\rangle\langle y|e^{Mt}|z\rangle.$

Clearly this expression is symmetric under the exchange of $y$ and $z$, which
proves Eq. (18). Using Cauchy-Schwartz one gets

$\|\pi_{t}-\pi\|_{1}^{2}=\left(\sum_{y}\sqrt{\frac{\pi(y)}{\pi(y)}}|\pi_{t}(y)-\pi(y)|\right)^{2}\leq\sum_{y}\frac{1}{\pi(y)}|\pi_{t}(y)-\pi(y)|^{2}=-1+\sum_{y}\frac{(\pi_{t}(y))^{2}}{\pi(y)}.$
(19)

The detailed balance condition Eq. (18) gives

$\frac{(\pi_{t}(y))^{2}}{\pi(y)}=\frac{\langle
y|e^{Gt}|x\rangle^{2}}{\pi(y)}=\frac{\langle x|e^{Gt}|y\rangle\langle
y|e^{Gt}|x\rangle}{\pi(x)}$ (20)

for all $y$. Combining Eqs. (19,20) gives

$\|\pi_{t}-\pi\|_{1}^{2}\leq-1+\frac{\langle
x|e^{2Gt}|x\rangle}{\pi(x)}=-1+\frac{\langle x|e^{2Mt}|x\rangle}{\pi(x)}.$
(21)

Here we used an identity $e^{2Gt}=De^{2Mt}D^{-1}$. Consider an eigenvalue
decomposition

$M=\sum_{i=1}^{2^{n}}\lambda_{i}(M)|\phi_{i}\rangle\langle\phi_{i}|$

such that $\lambda_{i}(M)$ is the $i$-th smallest eigenvalue of $M$ and
$\phi_{i}$ is the corresponding eigenvector. We have
$\langle\phi_{i}|\phi_{j}\rangle=\delta_{i,j}$. Since the matrices $M$ and $G$
are related by a similarity transformation, Lemma 3 implies that
$\lambda_{2^{n}}(M)=0$ is the largest eigenvalue of $M$ and
$\lambda_{i}(M)\leq-\gamma_{F}$ for all $i<2^{n}$. Furthermore, Lemma 2
implies $M|\psi\rangle=0$, that is, the only zero eigenvector of $M$ is
$\phi_{2^{n}}=\psi$. Thus we can write

$e^{2Mt}=|\psi\rangle\langle\psi|+R,\quad\mbox{where}\quad
R:=\sum_{i=1}^{2^{n}-1}e^{2\lambda_{i}(M)t}|\phi_{i}\rangle\langle\phi_{i}|$

obeys $\|R\|\leq e^{-2\gamma_{F}t}$. Substituting
$e^{2Mt}=|\psi\rangle\langle\psi|+R$ into Eq. (21) gives

$\|\pi_{t}-\pi\|_{1}^{2}\leq-1+1+\frac{\langle
x|R|x\rangle}{\pi(x)}\leq\frac{\|R\|}{\pi(x)}\leq\frac{e^{-2\gamma_{F}t}}{\pi(x)}.$
(22)

It remains to note that $\gamma_{F}\geq\gamma$, see Lemma 2. ∎

###  1.2 Sampling algorithm

Let $x_{in}\in\\{0,1\\}^{n}$ be the starting state. We would like to sample
$x$ from a distribution

$\pi_{t}(x):=\langle x|e^{Gt}|x_{in}\rangle.$

To this end we use Gillespie’s algorithm [12]. The algorithm takes as input
the starting state $x_{in}$, evolution time $t\geq 0$, and returns a sample
$x$ from the distribution $\pi_{t}(x)$.

Algorithm 1 Gillespie’s algorithm

1:$x\leftarrow x_{in}$

2:$\tau\leftarrow 0$

3:$\xi(0)\leftarrow x_{in}$

4:Sample $u\in[0,1]$ from the uniform distribution

5:$\Delta\tau\leftarrow\frac{\log{(1/u)}}{|\langle x|G|x\rangle|}$

6:Set $\xi(s)=x$ for all $s\in(\tau,\tau+\Delta\tau]$

7:$\tau\leftarrow\tau+\Delta\tau$

8:if $\tau\geq t$ then

9: return $x$

10:end if

11:Sample $y\in\\{0,1\\}^{n}\setminus\\{x\\}$ from the probability
distribution $\frac{\langle y|G|x\rangle}{|\langle x|G|x\rangle|}$

12:$x\leftarrow y$

13:go to line 4

Line 6 can be safely ignored as far as the implementation is concerned.
However it helps us to prove certain properties of the algorithm. Namely, each
run of the algorithm generates a continuous-time random walk on the set of bit
strings described by the piecewise constant function $\xi(s)$ with
$s\in[0,t]$.

######  Fact 1 (Output of Gillespie’s algorithm).

For each $s\in[0,t]$, the random variable $\xi(s)\in\\{0,1\\}^{n}$ generated
by Gillespie’s algorithm is distributed according to $\pi_{s}(x)$.

For completeness, we include a proof sketch for the correctness of Gillespie’s
algorithm.

###### Proof sketch.

The following derivations are based on Ref. [14]. Let
$P_{x,y}(\tau):=\Pr(\xi(\tau)=y|\xi(0)=x)=:\langle y|P(\tau)|x\rangle$ denote
the probability of being in state $y$ at time $\tau$ given that the initial
state is $\xi(0)=x$. We claim that $P(\tau)$ satisfies the differential
equation $P^{\prime}(\tau)=P(\tau)G$. Let $h>0$ where we think of $h$ as a
small quantity tending to zero. Let $m(x,h)$ be a random variable that is
equal to the number of transitions that occur in the time interval $[0,h]$
starting with initial state $x$. Since the time $\Delta\tau$ at which the
first transition occurs is exponentially distributed with rate $|\langle
x|G|x\rangle|$,

$\Pr(m(x,h)=0)=e^{-|\langle x|G|x\rangle|h}=1-|\langle
x|G|x\rangle|h+O(h^{2}).$ (23)

We can also show that

$\Pr(m(x,h)\geq 2)=\sum_{y\neq x}\Pr\left(\Delta\tau_{x}+\Delta\tau_{y}\leq
h\right)\frac{\langle y|G|x\rangle}{|\langle x|G|x\rangle|}=O(h^{2}).$

where $\Delta\tau_{x},\Delta\tau_{y}$ are exponential random variables with
rates $|\langle x|G|x\rangle|,|\langle y|G|y\rangle|$ respectively. Thus,

$\Pr(m(x,h)=1)=|\langle x|G|x\rangle|h+O(h^{2}).$ (24)

Eqs. (23,24) imply that for every $x\in\\{0,1\\}^{n}$ and every $y\neq x$,

$P_{x,x}(h)=1-|\langle x|G|x\rangle|h+O(h^{2})\qquad\text{and}\qquad
P_{x,y}(h)=\langle y|G|x\rangle h+O(h^{2}).$

For every $\tau\geq 0$ and $x,y\in\\{0,1\\}^{n}$, we arrive at

$\displaystyle P_{x,y}(\tau+h)$
$\displaystyle=P_{x,y}(\tau)P_{x,x}(h)+\sum_{z\neq x}P_{z,y}(\tau)P_{x,z}(h)$
$\displaystyle=P_{x,y}(\tau)-|\langle x|G|x\rangle|P_{x,y}(\tau)h+\sum_{z\neq
x}P_{z,y}(\tau)\langle z|G|x\rangle h+O(h^{2}).$

Hence, we conclude that

$\langle y|P^{\prime}(\tau)|x\rangle=\lim_{h\rightarrow
0}\frac{P_{x,y}(\tau+h)-P_{x,y}(\tau)}{h}=\sum_{z\in\\{0,1\\}^{n}}\langle
y|P(\tau)|z\rangle\langle z|G|x\rangle=\langle y|P(\tau)G|x\rangle.$

Solving this differential equation known as Kolmogorov’s forward equation, we
obtain the solution $P(s)=e^{Gs}$. This establishes the desired equivalence
$\Pr(\xi(s)=x|\xi(0)=x_{in})=\langle x|e^{Gs}|x_{in}\rangle=\pi_{s}(x)$. ∎

Let $m(x_{in},t)$ be the number of flips performed by the algorithm, that is,
the number of times the function $\xi(s)$ changes its value. One can easily
check that the total number of calls to the amplitude computation subroutine
made at lines 5 and 11 is at most $O(n^{k})\cdot m(x_{in},t)$. Thus the
algorithm is efficient as long as the number of flips is not too large. Our
key observation is the following.

######  Lemma 4 (Average number of flips).

Let $F_{\mathsf{od}}$ be the off-diagonal part of the fixed-node Hamiltonian
obtained from $F$ by setting to zero all diagonal matrix elements. Then

$\sum_{x_{in}}\pi(x_{in})m(x_{in},t)=-t\langle\psi|F_{\mathsf{od}}|\psi\rangle.$
(25)

We note that $\langle\psi|F_{\mathsf{od}}|\psi\rangle\leq 0$ by definition of
$F$, see Eq. (12).

Eq. (25) is closely related to a known estimator for the off-diagonal part of
the Hamiltonian in many-body simulations based on the continuous-time
worldline quantum Monte Carlo method [15] (see, e.g., Eq. (19) of [16]).

###### Proof.

Let $\rho(\xi)$ be the probability distributions over paths
$\xi\,:\,[0,t]\to\\{0,1\\}^{n}$ generated by the algorithm with the starting
state $x_{in}$ sampled from $\pi(x_{in})$. Since $\pi$ is the steady
distribution of the considered Markov chain, one infers that Gillespie’s
algorithm terminated at any intermediate time $s$ samples the steady
distribution, that is,

$\mathrm{Pr}_{\xi\sim\rho}[\xi(s)=x]=\pi(x)$ (26)

for any fixed $s\in[0,t]$. Now let us consider the probability $P_{x,s}(\xi)$
for the path $\xi$ with a fixed value $\xi(s)=x$ to contain a flip between
time $s$ and $s+ds$ in the limit $ds\to 0$. This probability is

$\displaystyle P_{x,s}(\xi)ds=\sum_{y\neq x}\langle y|G|x\rangle ds.$

Thus the average number of flips between time $s$ and $s+ds$ is given by

$P_{s}ds:=\mathbb{E}_{\xi\sim\rho}P_{\xi(s),s}ds=\sum_{y\neq x}\langle
y|G|x\rangle\pi(x)ds=-\langle\psi|F_{\mathsf{od}}|\psi\rangle ds.$ (27)

Here we used Eq. (26) and the definitions of $F$ and $G$, see Eqs. (12,17).
Integrating Eq. (27) over $s\in[0,t]$ completes the proof. ∎

By definition of $F$, see Eq. (12), one has

$-\langle\psi|F_{\mathsf{od}}|\psi\rangle=\sum_{(x,y)\in
S^{-}}|\langle\psi|x\rangle\langle x|H|y\rangle\langle
y|\psi\rangle|\leq\sum_{x,y}|\langle\psi|x\rangle\langle x|H|y\rangle\langle
y|\psi\rangle|\leq n^{k}\|H\|.$ (28)

Here we noted that each row of $H$ contains at most $n^{k}$ non-zero matrix
elements and each matrix element has magnitude at most $\|H\|$. Thus the
average number of flips performed by Gillespie’s algorithm is at most
$tn^{k}\|H\|$, assuming that the starting state $x_{in}$ is sampled from the
steady distribution $\pi(x_{in})$.

###### Proof of Theorem 2.

We shall use a truncated version of Gillespie’s algorithm which terminates
whenever the condition $\tau\geq t$ at line 8 is satisfied or the number of
flips exceeds a cutoff value $4\epsilon^{-1}tn^{k}\|H\|$. In the latter case
the truncated algorithm declares an error. By Markov’s inequality, there
exists at least one starting state $x^{\star}_{in}$ such that the algorithm
errs with the probability at most $\epsilon/4$. The truncated Gillespie’s
algorithm makes at most $O(t\epsilon^{-1}n^{2k}\|H\|)$ calls to the amplitude
computation subroutine. Let us write $\tilde{\pi}_{t}(x)$ for the probability
distribution sampled by the truncated algorithm, with some starting state
$x_{in}$, and conditioned on the event that no error occurs. Let $E$ be the
event $m(x_{in},t)\geq 4\epsilon^{-1}tn^{k}\|H\|$ and let $P(E,x_{in},t)$ be
the probability of this event, i.e., the probability that the truncated
Gillespie algorithm outputs an error. Then

$\pi_{t}(x)=\tilde{\pi}_{t}(x)(1-P(E,x_{in},t))+\eta_{t}(x)P(E,x_{in},t)$

where $\eta_{t}(x)$ is the probability that Gillespie’s algorithm outputs $x$
conditioned on event $E$. From the above we get

$\displaystyle\|\tilde{\pi}_{t}-\pi_{t}\|_{1}\leq\sum_{x}\left|\eta_{t}(x)-\tilde{\pi}_{t}(x)\right|\cdot
P(E,x_{in},t)\leq 2P(E,x_{in},t).$

For the special starting state $x^{\star}_{in}$ we see that the distribution
sampled by the truncated Gillespie algorithm is $\epsilon/2$-close to the
distribution $\pi_{t}$ sampled by the Gillespie algorithm without truncation.

Finally, we use Theorem 1 to choose the evolution time $t$ in Gillespie’s
algorithm large enough such that the distribution $\pi_{t}(x)$ is
$\epsilon/2$-close to the steady distribution $\pi(x)$ in total variation
distance. In particular, we conclude that the distributions
$\pi_{t}(x)=\langle x|e^{Gt}|x_{in}\rangle$ and $\pi(x)$ are
$\epsilon/2$-close in the total variation distance for all $t\geq
t_{\epsilon/2}$, where

$t_{\epsilon}=\frac{1}{\gamma}\log{\left(\frac{1}{\epsilon\sqrt{\pi(x_{in})}}\right)}\leq\frac{1}{\gamma}\log{\left(\frac{1}{\epsilon\sqrt{\pi^{*}}}\right)}$

This concludes the proof of Theorem 2. ∎

##  2 Application to the Haldane-Shastry model

In this section we further investigate the properties of the continuous-time
Markov Chain (CTMC) based on the fixed node Hamiltonian. Here we focus on a
specific example of a quantum spin system where the amplitudes of the ground
state can be computed exactly and efficiently. We implement the CTMC in
software and use it to approximately sample from the resulting probability
distribution.

The example we consider is a system of $L\geq 2$ qubits interacting according
to the Haldane-Shastry Hamiltonian

$H=\sum_{1\leq i<j\leq
L}\frac{X_{i}X_{j}+Y_{i}Y_{j}+Z_{i}Z_{j}}{4\left(\frac{L}{\pi}\sin\frac{\pi(i-j)}{L}\right)^{2}}$
(29)

as defined in, e.g., Refs. [17] and [18]. This Hamiltonian describes a spin
chain with periodic boundary conditions (i.e., a ring) and long-range two-
qubit interactions that depend on the shortest distance between the qubits on
the ring. The Hamiltonian Eq. (29) has a unique ground state $\psi$ with the
following properties [17, 18]. For every $x\in\\{0,1\\}^{n}$ such that
$|x|=\frac{L}{2}$,

$\langle x|\psi\rangle\propto\prod_{k=1}^{L}(-1)^{(k-1)x_{k}}\prod_{1\leq
i<j\leq L}\left(\sin\frac{\pi(i-j)}{L}\right)^{2x_{i}x_{j}}$ (30)

where $\prod_{k=1}^{L}(-1)^{(k-1)x_{k}}$ determines the relative sign, and for
every $x\in\\{0,1\\}^{n}$ such that $|x|\neq\frac{L}{2}$, $\langle
x|\psi\rangle=0$. It is also known [17] that the spectral gap of $H$ scales
inversely with system size $L$, and in particular satisfies

$\gamma=\frac{2\pi c}{L}$ (31)

for some constant $c$.

We chose this example because the entries of the ground state Eq. (30) are
efficiently computable, and because it is not immediately clear how to sample
from the probability distribution $\pi(x)=|\langle x|\psi\rangle|^{2}$ using
other, more standard, techniques. For example, in Appendix A we show that the
state $\psi$ is not a fermionic Gaussian state and therefore cannot be
directly sampled using a naive application of free-fermion based methods. One
can easily check that Haldane-Shastry Hamiltonian $H$ has a sign problem (that
is, $H$ is not stoquastic). Indeed, each two-qubit term in $H$ has a positive
matrix element between basis states $10$ and $01$. Moreover, in Appendix B we
use techniques of Ref. [19] to show that the sign problem in $H$ cannot be
“cured” by a local change of basis, that is, a Hamiltonian $U^{\dagger}HU$ has
a sign problem for any product unitary $U=U_{1}\otimes
U_{2}\otimes\cdots\otimes U_{L}$. This suggests that standard Quantum Monte
Carlo methods applicable to sign problem free Hamiltonians cannot be directly
used to sample the ground state of $H$.

The first step in our method for sampling from the probability distribution
$\pi(x)$ is to construct the associated fixed-node Hamiltonian defined by Eq.
(12). Note that all matrix elements of $H$ and all entries of $|\psi\rangle$
are real numbers, so in this example there is no need to preprocess the
Hamiltonian using Lemma 1. The associated CTMC is generated by the matrix $G$
from Eq. (17). The mixing time of this CTMC scales inversely with the spectral
gap $\gamma_{F}$ of $F$ (see Eq. (22)). We have shown that
$\gamma_{F}\geq\gamma$ holds in general, which, combined with Eq. (31) implies
that the mixing time is upper bounded as $O(L)$. Figure 1 shows the spectral
gaps $\gamma$ and $\gamma_{F}$ for the Haldane-Shastry Hamiltonian and small
values of $L$, computed using exact numerical diagonalization. We find
empirically that $\gamma_{F}$ has a milder scaling with system size $L$ than
$\gamma$ and we therefore expect the CTMC to mix more rapidly than the
rigorous bounds suggest.

(a)

(b)

Figure 1: (a) The inverse of the spectral gap of the fixed-node Hamiltonian
$F$ increases more slowly than that of the Haldane-Shastry Hamiltonian for
small values of the number of qubits $L$ (computed using exact numerical
diagonalization). (b) A linear fit for $L\geq 4$ on a log-log plot suggests a
scaling of $\gamma_{F}\sim L^{-0.7\ldots}$.

We implement Gillespie’s algorithm using the fixed-node construction to
approximately sample from the distribution $\pi(x)=|\langle
x|\psi\rangle|^{2}$, $x\in\\{0,1\\}^{L}$. The starting state is chosen
randomly from $\text{Supp}(|\psi\rangle)=\\{x\in\\{0,1\\}^{L}:\langle
x|\psi\rangle\neq 0\\}$. By our rigorous mixing time upper bound and
numerically determined spectral gaps, the Markov chain is expected to converge
rapidly. To assess convergence numerically in practice, we use the vanilla
$\hat{R}$ statistic discussed in [20] to select an appropriate “burn-in”
period for our chain.

After allowing the CTMC to (approximately) converge, we are then able to
estimate physical quantities in the ground state $\psi$. These can be compared
with exact formulas that are available for the Haldane-Shastry ground state.
For example, it is known from [17] that for every $1\leq i<j\leq L$, the two-
point $ZZ$ correlation admits the formula

$\langle\psi|Z_{i}Z_{j}|\psi\rangle=\left(\sum_{k=1}^{L/2}\frac{\sin\frac{(2k-1)\pi(i-j)}{L}}{2k-1}\right)\left(\frac{(-1)^{i-j}}{2L\sin\frac{\pi(i-j)}{L}}\right).$
(32)

For every $d\in\\{1,\ldots,L-1\\}$, we consider a two-point correlator which
is averaged over all pairs of spins at a distance $d$ on the ring. In
particular, we let $M_{d}\equiv L^{-1}\sum_{i=1}^{L}Z_{i}Z_{\text{Mod}(i+d)}$
where $\text{Mod}(i+d)=((i+d-1)\text{ mod }L)+1$. Then the ground state
satisfies

$\langle\psi|M_{d}|\psi\rangle=\langle\psi|Z_{1}Z_{1+d}|\psi\rangle$ (33)

and the exact value is therefore given by Eq. (32). We have tested our
implementation of the CTMC by computing the expected value of $M_{d}$ for
$d=1,5,10$. Since these observables are diagonal in the computational basis,
they can be straightforwardly estimated from the output
$\\{\xi(\tau):0\leq\tau\leq T\\}$ of Gillespie’s algorithm. In particular, let
$f:\\{0,1\\}^{L}\rightarrow\mathbb{R}$ be such that $f(z)=\langle
z|M_{d}|z\rangle$. Letting $\tau_{0}\leq T$ be an initial “burn-in” time used
to equilibrate the CTMC, we compute an approximation to
$\mathop{\mathbb{missing}}{E}_{x\sim\pi}[f(x)]=\langle\psi|M_{d}|\psi\rangle$
using an estimator

$\frac{1}{T-\tau_{0}}\int_{\tau=\tau_{0}}^{T}f(\xi(\tau))d\tau.$ (34)

(note that this integral can be equivalently be expressed as a finite sum
since $\xi$ is piecewise constant). The data from our CTMC is shown in Fig. 2
and compared with the exact value computed using Eq. (32). See Appendix C for
a description of the analysis used to estimate the error bars in this plot.

(a)

(b)

(c)

Figure 2:  Two-point correlation functions
$\langle\psi|Z_{i}Z_{i+d}|\psi\rangle$ with $d=1,5,10$ estimated using the
continuous-time Markov Chain and compared with the exact formula Eq. (32).
Here the CTMC was run for a total time $T=10^{6}$ and an initial “burn-in”
time $\tau_{0}=100$ was used for equilibration. The estimator $\hat{\mu}$ is
computed according to Eq. (34) and the estimate of standard deviation
$\hat{\sigma}$ is constructed as described in Appendix C.

In order to compare our CTMC with existing methods, we implemented a simple
Metropolis-Hastings Markov Chain to sample from the ground state $\psi$. Fig.
3 shows the nearest-neighbor two-point correlator $M_{1}$ estimated using the
MH method. To use a Metropolis-Hastings Markov chain to sample from the ground
state distribution $\pi$, we define the proposal distribution by viewing $H$
as the unweighted adjacency matrix of a transition graph with the state space
$\\{x\in\\{0,1\\}^{L}:|x|=\frac{L}{2}\\}$. When the chain is in state $x$, it
will propose to move to a neighbour $y$ of $x$ chosen uniformly at random.
More explicitly, for every $x,y\in\\{0,1\\}^{L}$ such that
$|x|=|y|=\frac{L}{2}$ and $x$ and $y$ differ in exactly two bits, the chain
proposes to transition from $x$ to $y$ with probability
$Q(y|x)=1/(L/2)^{2}=\frac{4}{L^{2}}$. The acceptance probabilities
$A(y|x)=\min\\{1,\frac{\pi(y)Q(x|y)}{\pi(x)Q(y|x)}\\}$ are defined so that $Q$
and $A$ together satisfy the detailed balance condition. The fixed-node MH
chain is defined by selecting a different proposal distribution, namely the
one induced by viewing $F$ as an unweighted adjacency matrix. We note that
both MH chains are irreducible and aperiodic 999Clearly the MH chain based on
$H$ is irreducible. The one based on $F$ is also irreducible –otherwise $F$
would have a non-unique ground state, contradicting Lemma 2. Both MH chains
are aperiodic when $L\geq 4$. This can be seen to follow from the fact that
they are irreducible and that there is at least one acceptance probability
$A(y|x)$ which is less than one. The latter implies $\sum_{y}A(y|x)Q(y|x)<1$
and thus the self-loop probability for the state $x$ is strictly positive. and
therefore the distribution obtained by running the chain with any initial
state converges to the limiting distribution $\pi$ as the number of steps
grows.

(a)

(b)

Figure 3:  Nearest-neighbor two-point correlator computed using a Metropolis-
Hastings Markov Chain with $5\cdot 10^{6}$ steps of the Markov Chain and
$10^{6}$ used for equilibration. In (a) the Markov Chain proposal distribution
was determined by $H$, and in (b) the proposal distribution was determined by
$F$.

To enable a comparison between the Metropolis-Hastings Markov Chain and the
CTMC we have tried to assess the computational cost of generating independent
samples using each method. For the MH chain, this is determined by the
autocorrelation time as measured in the number of steps of the (discrete-time)
Markov Chain. For the CTMC, we can also compute the autocorrelation time but
the computational cost of running the chain for a given interval of time
(using Gillespie’s algorithm) is determined by the number of transitions
during the interval. For this reason we choose to normalize the
autocorrelation time of the CTMC so that, roughly speaking, it is measured in
units of transitions rather than time, see Appendix C for details.

(a)

(b)

Figure 4: (a) Autocorrelation time of observables estimated using CTMC,
measured in units of the number of transitions. (b) Autocorrelation time of
observables estimated using Metropolis-Hastings with proposal distributions
generated by either the Haldane-Shastry Hamiltonian $H$ or the associated
fixed-node Hamiltonian $F$.

We note that the normalized autocorrelation times in Fig. 4 (a) are
significantly lower than the ones for the Metropolis-Hastings algorithm
reported in Fig. 4 (b). Intriguingly, we also observe in Fig. 4 (b) that the
Metropolis-Hastings algorithm using the fixed-node Hamiltonian $F$ has shorter
autocorrelation times for each number of qubits and each of the three
observables. In other words, using the fixed-node Hamiltonian to generate the
proposal distribution appears to have a marked advantage over the naive
strategy based on the Hamiltonian itself. This is despite the fact that the
nonzero matrix elements of $F$ are a subset of those of $H$. It appears that
the information about the sign structure of $\psi$ that determines which
entries of $F$ are set to zero may help the MH chain converge more quickly.

To investigate this further, we looked at how the autocorrelation times behave
when we use a corrupted ground state distribution $\tilde{\pi}$ as opposed to
the true ground state distribution $\pi$. For some noise strength $\kappa>0$,
we construct a perturbed distribution $\tilde{\pi}$ by defining

$\langle x|\tilde{\psi}\rangle=\langle
x|\psi\rangle+N\left(0,\frac{\kappa}{2^{L}}\right)$

and setting $\tilde{\pi}(x)=|\langle x|\tilde{\psi}\rangle|^{2}$ for every
$x\in\\{0,1\\}^{L}$ where $N(\mu,\sigma)$ is a Gaussian random variable with
mean $\mu$ and standard deviation $\sigma$. $\tilde{\pi}$ is normalized so
that $\sum_{x\in\\{0,1\\}^{L}}\tilde{\pi}(x)=1$. Fig. 5 shows that as the
total variance distance between $\pi$ and $\tilde{\pi}$ increases, the
differences between the autocorrelation times based on $H$ (Haldane-Shastry
Hamiltonian) and $F$ (defined by Eqs. (10,11,12) with $\psi$ replaced by
$\tilde{\psi}$) disappear as expected.

Figure 5: Autocorrelation times of Metropolis-Hastings chains defined using
corrupted ground state distributions. Each datapoint is computed using a
length 1000000 chain with a 1000-sample burn in period, all for $L=20$. When
the noise strength $\kappa$ is small, applying the fixed-node transformation
reduces autocorrelation times for each of the chosen observables. As the error
strength increases (as the corrupted distributions become dominated by noise),
the differences between the autocorrelation times of F and H vanish.

##  3 Acknowledgments

SB thanks Vojtech Havlicek for helpful discussions. DG and YL acknowledge the
support of the Natural Sciences and Engineering Research Council of Canada
through grant number RGPIN-2019-04198. DG also acknowledges the support of the
Canadian Institute for Advanced Research, and IBM Research. Research at
Perimeter Institute is supported in part by the Government of Canada through
the Department of Innovation, Science and Economic Development Canada and by
the Province of Ontario through the Ministry of Colleges and Universities.

##  Appendix A Haldane-Shastry state: lack of free fermion representation

Given a bit string $x\in\\{0,1\\}^{n}$, let $|x|$ be the number of 1s in $x$
(the Hamming weight). Suppose $n$ is an even integer. Consider an $n$-qubit
state $|\psi_{n}\rangle$ with amplitudes

$\langle x|\psi_{n}\rangle=e^{i\theta_{x}}\prod_{1\leq a<b\leq
n}\left[\sin{\left(\frac{\pi(a-b)}{n}\right)}\right]^{2x_{a}x_{b}}$ (35)

if $|x|=n/2$ and $\langle x|\psi_{n}\rangle=0$ otherwise. Here
$\theta_{x}\in\mathbb{R}$ are arbitrary angles. For example,

$|\psi_{4}\rangle=\frac{e^{i\theta_{1100}}}{2}|1100\rangle+\frac{e^{i\theta_{0011}}}{2}|0011\rangle+\frac{e^{i\theta_{1001}}}{2}|1001\rangle+\frac{e^{i\theta_{0110}}}{2}|0110\rangle+e^{i\theta_{1010}}|1010\rangle+e^{i\theta_{0101}}|0101\rangle.$
(36)

The state defined in Eq. (35) coincides with the Haldane-Shastry state for a
suitable choice of the phase factors $e^{i\theta_{x}}=\pm 1$.

Define Majorana fermion operators $\gamma_{1},\gamma_{2},\ldots,\gamma_{2n}$
such that $\gamma_{1}=X_{1}$, $\gamma_{2}=Y_{1}$,

$\gamma_{2p-1}=Z_{1}\cdots
Z_{p-1}X_{p}\quad\mbox{and}\quad\gamma_{2p}=Z_{1}\cdots Z_{p-1}Y_{p}$ (37)

for $p=2,\ldots,n$. Recall that $n$-qubit state $|\phi\rangle$ is called a
free fermion state [21] if it obeys fermionic Wick’s theorem, that is, for any
$k$-tuple of Majorana operators $1\leq p_{1}<p_{2}<\ldots<p_{k}\leq 2n$ one
has

$\langle\phi|\gamma_{p_{1}}\gamma_{p_{2}}\cdots\gamma_{p_{k}}|\phi\rangle=\left\\{\begin{array}[]{rcl}0&\mbox{if}&\mbox{$k$
is odd}\\\ {\cal
A}\left(\langle\phi|\gamma_{p_{1}}\gamma_{p_{2}}|\phi\rangle\cdot\langle\phi|\gamma_{p_{3}}\gamma_{p_{4}}|\phi\rangle\cdots\langle\phi|\gamma_{p_{k-1}}\gamma_{p_{k}}|\phi\rangle\right)&\mbox{if}&\mbox{$k$
is even}\\\ \end{array}\right.$ (38)

Here ${\cal A}$ denotes anti-symmetrization over all permutations of indices
$p_{1},p_{2},\ldots,p_{k}$. For example,

$\langle\phi|\gamma_{1}\gamma_{2}\gamma_{3}\gamma_{4}|\phi\rangle=\langle\phi|\gamma_{1}\gamma_{2}|\phi\rangle\cdot\langle\phi|\gamma_{3}\gamma_{4}|\phi\rangle-\langle\phi|\gamma_{1}\gamma_{3}|\phi\rangle\cdot\langle\phi|\gamma_{2}\gamma_{4}|\phi\rangle+\langle\phi|\gamma_{1}\gamma_{4}|\phi\rangle\cdot\langle\phi|\gamma_{2}\gamma_{3}|\phi\rangle.$

The task of sampling the probability distribution $|\langle
x|\phi\rangle|^{2}$ with a free fermion state $|\phi\rangle$ admits an
efficient classical algorithm with the runtime $O(n^{3})$, see for instance
[21] or Appendix B of [22]. However, these free fermion simulation methods are
not directly applicable to simulating measurement of the Haldane-Shastry
state, as follows from the following lemma.

######  Lemma 5.

Let $|\psi_{n}\rangle$ be the state defined in Eq. (35) with an arbitrary
choice of the angles $\theta_{x}$. Suppose $n\geq 4$ is an integer multiple of
four. Then $|\psi_{n}\rangle$ is not proportional to a free fermion state.

###### Proof.

We shall need two well-known facts.

######  Fact 2 ([23]).

Suppose a four-qubit state $|\phi\rangle$ has support only on even-weight
basis vectors. Then $|\phi\rangle$ is a free fermion state if and only if

$-\langle 0000|\phi\rangle\cdot\langle 1111|\phi\rangle+\langle
1100|\phi\rangle\cdot\langle 0011|\phi\rangle-\langle
1010|\phi\rangle\cdot\langle 0101|\phi\rangle+\langle
1001|\phi\rangle\cdot\langle 0110|\phi\rangle=0.$ (39)

######  Fact 3 ([21]).

Suppose $|\phi\rangle$ is an $n$-qubit free fermion state and $D$ is a tensor
product of diagonal single-qubit operators. Then $D|\phi\rangle$ is
proportional to a free fermion state.

Assume that $|\psi_{4}\rangle\sim|\phi\rangle$ for some free fermion state
$|\phi\rangle$. Then Eq. (39) gives

$|\langle 1010|\psi_{4}\rangle\cdot\langle 0101|\psi_{4}\rangle|\leq|\langle
1100|\psi_{4}\rangle\cdot\langle 0011|\psi_{4}\rangle|+|\langle
1001|\psi_{4}\rangle\cdot\langle 0110|\psi_{4}\rangle|$ (40)

which contradicts to Eq. (36). This proves the lemma for $n=4$. Consider an
integer $n=4m\geq 4$. Partition the set of $n$ qubits as $[n]=AB$. Given a bit
string $x\in\\{0,1\\}^{n}$, let $x_{A}$ and $x_{B}$ be the restrictions of $x$
onto $A$ and $B$. Fix an arbitrary bit string $x_{B}$ such that
$|x_{B}|=n/2-2$ and ${}_{B}\langle x_{B}|\psi_{n}\rangle\neq 0$. By definition
of $|\psi_{n}\rangle$, one has

$\langle x|\psi_{n}\rangle=e^{i\theta_{x}}\prod_{j\in
A}f_{j}(x_{j})\prod_{\begin{subarray}{c}1\leq q<r\leq n\\\ q,r\in A\\\
\end{subarray}}\left[\sin{\left(\frac{\pi(q-r)}{n}\right)}\right]^{2x_{q}x_{r}}$
(41)

for all $x_{A}$ with $|x_{A}|=2$ and $\langle x|\psi_{n}\rangle=0$ if
$|x_{A}|\neq 2$. Here $f_{j}(x_{j})$ are some functions that may depend on
$x_{B}$. Choose

$A=\\{1,m+1,2m+1,3m+1\\}.$

Using the identity $n=4m$ and Eq. (41) one gets

$\langle x|\psi_{n}\rangle=\langle x_{A}|D_{1}D_{2}D_{3}D_{4}|\psi_{4}\rangle$
(42)

for all $x_{A}$ where $|\psi_{4}\rangle$ is defined in Eq. (36) and $D_{j}$
are diagonal single-qubit operators that may depend of $x_{B}$. Assume that
$|\psi_{n}\rangle$ is free. Then Fact 3 implies that a four-qubit state
${}_{B}\langle x_{B}|\psi_{n}\rangle$ is free for any $x_{B}$. Using Fact 3
again and Eq. (42) one infers that

$|\psi_{4}\rangle=D_{1}^{-1}D_{2}^{-1}D_{3}^{-1}D_{4}^{-1}{}_{B}\langle
x_{B}|\psi_{n}\rangle$

is also free. However this is a contradiction since we have already proved
that $|\psi_{4}\rangle$ is not free. ∎

##  Appendix B Sign problem in Haldane-Shastry Hamiltonian

Recall that an $n$-qubit Hamiltonian $H$ is called stoquastic (sign problem
free) if $H$ has real matrix elements in the standard basis and $\langle
x|H|y\rangle\leq 0$ for all $x\neq y$. Consider a Hamiltonian

$H=\sum_{1\leq i<j\leq n}J_{i,j}(X_{i}X_{j}+Y_{i}Y_{j}+Z_{i}Z_{j}),$ (43)

where $J_{i,j}>0$ are arbitrary coefficients. This includes Haldane-Shastry
Hamiltonian Eq. (29) as a special case. Choose any basis vectors
$x,y\in\\{0,1\\}^{n}$ that differ only on two qubits $i,j$ such that
$x_{i}x_{j}=10$, $y_{i}y_{j}=01$, and $x_{\ell}=y_{\ell}$ for all
$\ell\notin\\{i,j\\}$. A simple calculation gives $\langle
x|H|y\rangle=2J_{i,j}>0$, that is, $H$ is not stoquastic. Suppose there exist
single-qubit unitary operators $U_{1},U_{2},\ldots,U_{n}$ such that

$H_{U}:=(U_{1}\otimes U_{2}\otimes\cdots\otimes U_{n})^{\dagger}H(U_{1}\otimes
U_{2}\otimes\cdots\otimes U_{n})$

is stoquastic. The question of whether Hamiltonians of the form Eq. (43) can
be made stoquastic by a local change of basis has been studied by Klassen and
Terhal [19]. Lemma 22 of Ref. [19] implies that without loss of generality the
unitaries $U_{j}$ can be chosen from a finite subgroup of the unitary group
(known as the Clifford group). Hence the number of distinct unitaries $U_{j}$
is upper bounded by a constant independent of $n$. Thus for a sufficiently
large number of qubits $n$, there will be at least one pair of qubits $i<j$
such that $U_{i}=U_{j}$. Note that $U_{i}\otimes U_{j}$ commutes with
$X\otimes X+Y\otimes Y+Z\otimes Z$ if $U_{i}=U_{j}$. Thus we can write

$H_{U}=J_{i,j}(X_{i}X_{j}+Y_{i}Y_{j}+Z_{i}Z_{j})+H_{\mathrm{else}},$

where $H_{\mathrm{else}}$ is a sum of operators that act non-trivially on at
most one of the qubits $i,j$. Now the same calculation as above shows that
$\langle x|H_{U}|y\rangle=2J_{i,j}>0$, that is, $H_{U}$ is not stoquastic.

##  Appendix C Details of numerical implementation

In this appendix we describe how error bars are computed in Fig. 2. We also
describe the definition of the normalized autocorrelation time reported in
Fig. 4.

Let $\\{\xi(\tau):\tau\geq 0\\}$ denote the stochastic process induced by
running Gillespie’s algorithm on the HS fixed-node Hamiltonian where the
initial state $\xi(0)\sim\pi$, the true ground state distribution. The
following derivation is based on Ref. [24]. Let
$f:\\{0,1\\}^{L}\rightarrow\mathbb{R}$ be a function and suppose our goal is
to estimate the mean of $f$ with respect to the steady distribution $\pi$. Let
$T>0$ be large and $h>0$ be small. We can estimate
$\mu=\mathop{\mathbb{missing}}{E}_{x\sim\pi}[f(x)]$ using the estimator

$\hat{\mu}=\frac{1}{\lfloor\frac{T}{h}\rfloor}\sum_{j=1}^{\lfloor\frac{T}{h}\rfloor}f(\xi(jh)).$

Note that in the limit $h\rightarrow 0$ this can be represented as an integral
(cf. Eq. (34)). For the purposes of estimating error bars it will be
convenient to use the discretized representation however. Since
$\xi(0)\sim\pi$, $\xi(\tau)\sim\pi$ for every $\tau\geq 0$. Let
$\sigma^{2}_{f}=\mathop{\text{Var}}(f(\xi(0)))=\mathop{\text{Var}}(f(\xi(\tau)))$
for every $\tau\geq 0$. Recall that the Pearson correlation coefficient of two
random variables $A$ and $B$ is defined by

$\text{cor}(A,B)=\frac{\text{cov}(A,B)}{\sigma_{A}\sigma_{B}}$

where $\sigma_{A}^{2}$ and $\sigma_{B}^{2}$ are the variances of $A$ and $B$
respectively. Then, the variance of the estimator $\hat{\mu}$ is

$\displaystyle\hat{\sigma}^{2}=\mathop{\text{Var}}(\hat{\mu})$
$\displaystyle=\frac{1}{\lfloor\frac{T}{h}\rfloor^{2}}\mathop{\text{Var}}\left(\sum_{j=1}^{\lfloor\frac{T}{h}\rfloor}f(\xi(jh))\right)$
$\displaystyle=\frac{1}{\lfloor\frac{T}{h}\rfloor^{2}}\sum_{i=1}^{\lfloor\frac{T}{h}\rfloor}\sum_{j=1}^{\lfloor\frac{T}{h}\rfloor}\text{cov}(f(\xi(ih)),f(\xi(jh)))$
$\displaystyle=\frac{\sigma^{2}_{f}}{\lfloor\frac{T}{h}\rfloor^{2}}\sum_{i=1}^{\lfloor\frac{T}{h}\rfloor}\sum_{j=1}^{\lfloor\frac{T}{h}\rfloor}\text{cor}(f(\xi(ih)),f(\xi(jh))).$

We assume for every $i\in\\{1,\ldots,\lfloor\frac{T}{h}\rfloor\\}$, each
correlation term in
$\tau_{f}=\sum_{j=1}^{\lfloor\frac{T}{h}\rfloor}\text{cor}(f(\xi(ih)),f(\xi(jh)))$
only depends on $|i-j|$ and is otherwise independent of $i$. Thus,

$\hat{\sigma}^{2}=\frac{\sigma^{2}_{f}\tau_{f}}{\lfloor\frac{T}{h}\rfloor}.$

$\tau_{f}$ is the integrated autocorrelation time w.r.t $f$. We use the emcee
library [25] to obtain an estimate $\hat{\tau}_{f}$ of $\tau_{f}$. We further
estimate $\hat{\sigma}^{2}$ using $\hat{\tau}_{f}$ and the sample variance of
$\\{f(\xi(1h)),f(\xi(2h)),\ldots,f(\xi(\lfloor\frac{T}{h}\rfloor h))\\}$.

Notice that $\hat{\tau}_{f}h$ estimates the autocorrelation time of the CTMC.
Let $r$ denote the total number of transitions it took for the CTMC to reach
time $T$. Then $\frac{r}{T}$ gives the average number of transitions needed
for the CTMC to advance time by $1$ unit. Thus, we infer that it takes an
average of $\hat{\tau}_{f}h\frac{r}{T}$ transitions for the CTMC to advance
time by $\hat{\tau}_{f}h$ units. In Figure 4 we plot the normalized
autocorrelation time $\hat{\tau}_{f}h\frac{r}{T}$ against the number of
qubits.

## References

  * [1] Francisco Barahona.  On the computational complexity of Ising spin glass models.  Journal of Physics A: Mathematical and General, 15(10):3241, 1982\.
  * [2] Philip M Long and Rocco A Servedio.  Restricted Boltzmann machines are hard to approximately evaluate or simulate.  Proceedings of the 27th International Conference on International Conference on Machine Learning. ICML’10. Haifa, Israel, page 703–710, 2010.
  * [3] W. K. Hastings.  Monte Carlo sampling methods using Markov chains and their applications.  Biometrika, 57(1):97–109, April 1970.
  * [4] David A Levin and Yuval Peres.  Markov chains and mixing times, volume 107.  American Mathematical Soc., 2017.
  * [5] Sergey Bravyi, David Gosset, and Yinchen Liu.  How to simulate quantum measurement without computing marginals.  Physical Review Letters, 128(22):220503, 2022.
  * [6] Dorit Aharonov and Amnon Ta-Shma.  Adiabatic quantum state generation and statistical zero knowledge.  In Proceedings of the thirty-fifth annual ACM symposium on Theory of computing, pages 20–29, 2003.
  * [7] Sergey Bravyi and Barbara Terhal.  Complexity of stoquastic frustration-free hamiltonians.  SIAM Journal on Computing, 39(4):1462–1485, 2010.
  * [8] DFB Ten Haaf, HJM Van Bemmel, JMJ Van Leeuwen, W Van Saarloos, and DM Ceperley.  Proof for an upper bound in fixed-node Monte Carlo for lattice fermions.  Physical Review B, 51(19):13039, 1995.
  * [9] WMC Foulkes, Lubos Mitas, RJ Needs, and Guna Rajagopal.  Quantum monte carlo simulations of solids.  Reviews of Modern Physics, 73(1):33, 2001.
  * [10] Federico Becca and Sandro Sorella.  Quantum Monte Carlo Approaches for Correlated Systems.  Cambridge University Press, 2017.
  * [11] Vojtech Havlicek.  Amplitude ratios and neural network quantum states.  arXiv preprint arXiv:2201.09128, 2022.
  * [12] Daniel T Gillespie.  Exact stochastic simulation of coupled chemical reactions.  The journal of physical chemistry, 81(25):2340–2361, 1977.
  * [13] Persi Diaconis and Daniel Stroock.  Geometric bounds for eigenvalues of Markov chains.  The Annals of Applied Probability, pages 36–61, 1991.
  * [14] Glen Takahara.  STAT 455 Stochastic Process Lecture Notes.  2017\.
  * [15] NV Prokof’Ev, BV Svistunov, and IS Tupitsyn.  Exact, complete, and universal continuous-time worldline monte carlo approach to the statistics of discrete quantum systems.  Journal of Experimental and Theoretical Physics, 87(2):310–321, 1998.
  * [16] Edward Farhi, Jeffrey Goldston, David Gosset, Sam Gutmann, Harvey B. Meyer, and Peter Shor.  Quantum adiabatic algorithms, small gaps, and different paths.  Quantum Info. Comput., 11(3):181–214, mar 2011.
  * [17] Jean-Marie Stephan and Frank Pollmann.  Full counting statistics in the haldane-shastry chain.  Physical Review B, 95(3):035119, 2017.
  * [18] Shriya Pai, NS Srivatsa, and Anne EB Nielsen.  Disordered haldane-shastry model.  Physical Review B, 102(3):035117, 2020.
  * [19] Joel Klassen and Barbara M Terhal.  Two-local qubit hamiltonians: when are they stoquastic?  Quantum, 3:139, 2019.
  * [20] Aki Vehtari, Andrew Gelman, Daniel Simpson, Bob Carpenter, and Paul-Christian Bürkner.  Rank-normalization, folding, and localization: An improved r for assessing convergence of mcmc (with discussion).  Bayesian analysis, 16(2):667–718, 2021.
  * [21] Barbara M Terhal and David P DiVincenzo.  Classical simulation of noninteracting-fermion quantum circuits.  Physical Review A, 65(3):032325, 2002.
  * [22] Sergey Bravyi, Matthias Englbrecht, Robert König, and Nolan Peard.  Correcting coherent errors with surface codes.  npj Quantum Information, 4(1):1–6, 2018.
  * [23] Sergey Bravyi.  Contraction of matchgate tensor networks on non-planar graphs.  Contemp. Math, 482:179–211, 2009.
  * [24] Tom Kennedy.  Monte Carlo Methods - a special topics course.  2016\.
  * [25] Daniel Foreman-Mackey, David W Hogg, Dustin Lang, and Jonathan Goodman.  emcee: the mcmc hammer.  Publications of the Astronomical Society of the Pacific, 125(925):306, 2013.