Quantum-Inspired Tempering for Ground State Approximation using Artificial
Neural Networks

T. Albash1,2,3*, C. Smith1,3, Q. Campbell3 A. D. Baczewski2,3,4

1 Department of Electrical and Computer Engineering, University of New Mexico,
Albuquerque, NM 87131, USA

2 Department of Physics and Astronomy, University of New Mexico, Albuquerque,
NM 87131, USA

3 Center for Quantum Information and Control (CQuIC), University of New
Mexico, Albuquerque, NM 87131, USA

4 Quantum Algorithms and Applications Collaboratory, Sandia National
Laboratories, Albuquerque NM 87185, USA

*<EMAIL_ADDRESS>
## Abstract

A large body of work has demonstrated that parameterized artificial neural
networks (ANNs) can efficiently describe ground states of numerous interesting
quantum many-body Hamiltonians. However, the standard variational algorithms
used to update or train the ANN parameters can get trapped in local minima,
especially for frustrated systems and even if the representation is
sufficiently expressive. We propose a parallel tempering method that
facilitates escape from such local minima. This methods involves training
multiple ANNs independently, with each simulation governed by a Hamiltonian
with a different “driver” strength, in analogy to quantum parallel tempering,
and it incorporates an update step into the training that allows for the
exchange of neighboring ANN configurations. We study instances from two
classes of Hamiltonians to demonstrate the utility of our approach using
Restricted Boltzmann Machines as our parameterized ANN. The first instance is
based on a permutation-invariant Hamiltonian whose landscape stymies the
standard training algorithm by drawing it increasingly to a false local
minimum. The second instance is four hydrogen atoms arranged in a rectangle,
which is an instance of the second quantized electronic structure Hamiltonian
discretized using Gaussian basis functions. We study this problem in a minimal
basis set, which exhibits false minima that can trap the standard variational
algorithm despite the problem’s small size. We show that augmenting the
training with quantum parallel tempering becomes useful to finding good
approximations to the ground states of these problem instances.

###### Contents

  1. 1 Introduction
  2. 2 Overview of the Methodology
  3. 3 Results
    1. 3.1 Precipice Problem
    2. 3.2 H4 Rectangle
      1. 3.2.1 Fixed Hamming Weight
      2. 3.2.2 No Fixed Hamming Weight
  4. 4 Conclusions
  5. 5 Code
  6. A Review of Stochastic Reconfiguration
  7. B Adaptive Learning Rate Scheme
  8. C Calculating the expectation value of operators with respect to an RBM ansatz
  9. D Alternative Parallel Tempering Schemes
    1. D.1 Standard deviation of $H^{2}$ as a measure of temperature
    2. D.2 Almost constant swap probability
  10. E Short Note On the Bootstrap Procedure
  11. F The H4 rectangle

##  1 Introduction

Approximating the ground state of a many-body quantum Hamiltonian is a
computational task at the heart of many problems in the physical sciences.
Tackling this problem requires both an efficient representation of the wave
function and an algorithm for optimizing the wave function to minimize its
energy whose individual update steps are efficient. A variety of methods have
been developed, each with its advantages and limitations. Quantum Monte Carlo
(QMC) [1, 2] methods are generically useful in the absence of the sign problem
[3], although there are QMC approaches that try to mitigate the bottlenecks
associated with the sign problem [4, 5, 6]. Methods based on Matrix Product
States (MPS) [7, 8, 9, 10, 11, 12] have been extremely successful at
addressing systems in one dimension and methods based on Projected Entangled
Pair States (PEPS) [13, 14] have allowed the study of certain systems in two
dimensions. More recently, artificial neural networks (ANNs) have been used as
ansätze for quantum ground states [15], with the advantages of being agnostic
to the sign problem and able to take as input Hamiltonians with arbitrary
connectivity. While there have been important demonstrations of the utility of
this approach [16, 17, 18], it still suffers when studying frustrated systems
by getting stuck in the many local minima of a rugged optimization landscape
[19].

Parallel tempering (PT) [20, 21, 22] is a Monte Carlo method that has been
extensively used to help explore the rugged energy landscapes of spin glasses
[23]. In the typical PT approach, multiple “replicas” of the system evolve
independently at different temperatures and hence in different free energy
landscapes. The exchange of configurations between higher and lower
temperature replicas facilitates escape from local minima. There is nothing
special about using temperature to alter the landscape, and other physically
motivated choices can be made. In quantum parallel tempering (QPT) [24, 25],
replicas are evolved according to different Hamiltonians, where the strength
of a driver Hamiltonian, e.g. a uniform transverse field, is used to control
the landscape the replicas are evolving in instead of the temperature.

In this paper we combine the methods of QPT and the training of ANNs to
approximate ground states of quantum spin systems. For simplicity of
presentation, we focus on a specific choice of ANN, the restricted Boltzmann
machine (RBM), and we demonstrate how a QPT method implemented in addition to
the standard training can help overcome barriers in energy landscape. We do
this using two examples. The first is based on a permutation-invariant
Hamiltonian that is designed to trap the standard training approach in a false
minimum. In this case, the QPT provides a clear advantage over simply
repeating the standard algorithm starting from independent initial points. The
second is based on a molecular system of four hydrogen atoms arranged in a
rectangle with an angle $\theta$ between adjacent hydrogen atoms that we can
vary. The fermionic second quantized Hamiltonian is mapped to a system of
qubits using the Jordan-Wigner transformation [26], and we apply the ANN
training methods to the qubit Hamiltonian [27]. We find that the hardness of
approximating the ground state of this system depends on different aspects,
such as the choice of angles and whether we constrain the training to a fixed
magnetization sector corresponding to a fixed particle number, and we find
that QPT can provide a computational advantage in the cases where the standard
approach struggles to find a good approximation of the ground state.

The paper is organized as follows. In Section 2, we review the standard
training approach of the RBM ansatz for quantum wave functions in addition to
how we introduce QPT to it. In Section 3, we introduce our two problem classes
and compare the performance of the standard training approach to our QPT
augmented approach. In Section 4, we conclude with open questions and future
directions of our work.

##  2 Overview of the Methodology

We consider an unnormalized parameterized wave function expressed in the
computational basis
$|\psi(\alpha)\rangle=\sum_{x\in\left\\{0,1\right\\}^{n}}\psi_{x}(\alpha)|x\rangle$,
where $|x\rangle\equiv|x_{n}\rangle\otimes\dots\otimes|x_{1}\rangle$ denotes
the $n$-qubit computational basis state with
$\sigma_{i}^{z}|x_{i}\rangle=(1-2x_{i})|x_{i}\rangle=v_{i}|x_{i}\rangle$ and
$\alpha$ denotes the set of $K$ complex variational parameters. Here we
restrict our attention to the case where $\psi_{x}(\alpha)$ is the weight
associated with a RBM with parameters
$\alpha=\left\\{a_{i},b_{\mu},W_{i\mu}\right\\}_{i=1,\mu=1}^{n,m}$ and with
$x$ as the value of its input layer,

$\displaystyle\psi_{x}(\alpha)$ $\displaystyle=$ $\displaystyle
e^{\sum_{i=1}^{n}a_{i}v_{i}}\prod_{\mu=1}^{m}\cosh\left(b_{\mu}+\sum_{i=1}^{n}W_{i\mu}v_{i}\right),$
(1)

where $v_{i}\in\left\\{-1,1\right\\}$. As long as $m$ scales polynomially with
$n$, the RBM representation provides an efficient representation of the
quantum state $|\psi(\alpha)\rangle$.

Given a many-body Hamiltonian $H$ defined on $n$ qubits, our aim is to find an
approximate RBM representation to the ground state of $H$. Therefore, our
optimization problem is to identify a set of parameters $\alpha_{\ast}$ such
that

$\alpha_{\ast}\in\operatorname*{argmin}_{\alpha}\
\frac{\langle\psi(\alpha)|H|\psi(\alpha)\rangle}{\langle\psi(\alpha)|\psi(\alpha)\rangle}\
.$ (2)

The cost function on the right hand side defines a multi-dimensional parameter
landscape, and if the number of variational parameters is sufficiently large,
the global minimum of this landscape gives an excellent approximation of the
ground state.

This optimization problem is solved using Stochastic Reconfiguration (SR) [28,
29], which updates the parameters as

$\displaystyle\alpha^{\prime}_{k}$ $\displaystyle=$
$\displaystyle\alpha_{k}+\delta\alpha_{k}=\alpha_{k}-\gamma\sum_{k^{\prime}}\left(S(\alpha)\right)^{-1}_{kk^{\prime}}F_{k^{\prime}}(\alpha)\
,$ (3)

where $S_{kk^{\prime}}(\alpha)$ are the matrix elements of the covariance
matrix and $F_{k}(\alpha)$ are the elements of the force vector defined as,

$\displaystyle S_{kk^{\prime}}(\alpha)$ $\displaystyle=\langle
O_{k}(\alpha)^{\dagger}O_{k^{\prime}}(\alpha)\rangle-\langle
O_{k}(\alpha)^{\dagger}\rangle\langle O_{k^{\prime}}(\alpha)\rangle$ (4a)
$\displaystyle F_{k}(\alpha)$ $\displaystyle=\langle
O_{k}(\alpha)^{\dagger}H\rangle-\langle O_{k}(\alpha)^{\dagger}\rangle\langle
H\rangle\ ,$ (4b)

where the expectation values are taken with respect to the state
$|\psi(\alpha)\rangle$ and where $O_{k}(\alpha)$ is the operator associated
with the gradient of the state $|\psi(\alpha)\rangle$ with respect to the
variational parameter $\alpha_{k}$,

$O_{k}(\alpha)=\sum_{x=\left\\{0,1\right\\}^{n}}\frac{1}{\psi_{x}(\alpha)}\frac{\partial}{\partial\alpha_{k}}\psi_{x}(\alpha)|x\rangle\langle
x|\ .$ (5)

Since $S$ may not be invertible, $S^{-1}$ is strictly speaking the Moore-
Penrose pseudoinverse [30, 31]. We review the derivations of these formulas in
Appendix A.

The parameter $\gamma$ sets the learning rate, and we choose it adaptively
using a second-order Runge-Kutta (RK) method [32, 17, 19]. We review this
adaptive procedure in Appendix B. While the RK method incurs a computational
overhead because it requires multiple estimates of the covariance matrix and
force vector, the combination of RK with SR generally results in more accurate
results with fewer updates in total compared to other optimization schemes
[19] such as ADAM [33] with SR or stochastic gradient descent, and it
eliminates the need to specify a learning rate. In this regard, the RK with SR
optimization approach to train the RBM sets a competitive baseline for us to
improve upon.

Calculating the covariance matrix $S$ and force vector $F$ requires
expectation values of the observables associated with the set
$\left\\{O_{k}\right\\}$ and $H$. These can be estimated using Monte Carlo
importance sampling [15], but in this work we will restrict ourselves to
evaluating these expectation values exactly. The reason for using exact
sampling is to eliminate the role of finite sampling in the training. We give
details of how to calculate the necessary expectation values in Appendix C.

With $S$ and $F$ in hand, the final step to implement Eq. (3) is to calculate
the pseudoinverse of $S$. In what follows we calculate the inverse using an
explicit regularization where we replace the diagonal elements of $S$ by [15]

$S_{kk}\to S_{kk}^{\mathrm{reg}}=S_{kk}\left(1+\lambda(p)\right)\ ,$ (6)

where $p$ denotes the update step and
$\lambda(p)=\max\left(\lambda_{0}b^{p},\lambda_{\mathrm{min}}\right)$ with
$\lambda_{0}=100,b=0.9,\lambda_{\mathrm{min}}=10^{-4}$. We use the MINRES-QLP
[34] with a maximum of $10^{3}$ iterations and a stopping tolerance of
$10^{-6}$ to solve the linear system of equations. The choice of these
parameters can affect the performance of SR to find a description of the
ground state, and we have not tried to optimize these parameters for our
problem; our aim is to compare the _relative_ performance of the overall
algorithm with and without QPT using the same set of parameters otherwise.

The SR algorithm can be shown [35] to be a generalization of the Riemannian
natural gradient [36], so like other gradient descent algorithms it is
susceptible to being trapped in local minima of the energy landscape. Finite
sampling and the resulting fluctuations in $S$ and $F$ can be a source of
noise that can help the training escape local minima [37], but it can also
hinder an accurate determination of the gradient and hence prevent the
training from finding high-quality solutions. Motivated by the success of PT
[20, 21, 22] in the study of equilibration in spin systems, we propose the
following optimization algorithm that is inspired by QPT [24, 25]. We perform
the usual training on $N$ RBMs, each called a replica in the parlance of PT.
Each training is done with a different training Hamiltonian, where the
Hamiltonian for the $r$th replica is given as

$H(\Gamma_{r})=H_{\mathrm{T}}+\Gamma_{r}H_{\mathrm{M}},$ (7)

where $H_{\mathrm{T}}$ is the target Hamiltonian whose ground state we are
actually after and $H_{\mathrm{M}}$ is a suitably chosen mixer Hamiltonian
with a trivial ground state, preferably the one with a uniform superposition
over the relevant basis states. This choice of training Hamiltonian is not
unique, and we give an alternative formulation in Appendix D. For example, in
the absence of any constraints, the mixer can be taken to be the uniform
transverse field Hamiltonian,
$H_{\mathrm{M}}=-H_{0}\sum_{i=1}^{n}\frac{1}{2}\left(\mathds{1}-\sigma_{i}^{x}\right)$,
where $H_{0}$ is a characteristic energy scale of $H_{\mathrm{T}}$. In this
case, for replicas with $\Gamma_{r}\gg 1$ the training Hamiltonian is
dominated by the transverse field, and the ground state is the uniform
superposition state. This is analogous to the high-temperature limit in
“standard” PT, where the different replicas operate at different temperatures
but with the same Hamiltonian. For replicas with $\Gamma_{r}$ close to 0, the
training Hamiltonian is dominated by the target Hamiltonian, and the ground
state is close to the target ground state. This is analogous to the low-
temperature limit in standard parallel tempering. The parameter $\Gamma_{r}$
thus allows us to interpolate between the target parameter landscape at
$\Gamma_{r}=0$ and a much simpler parameter landscape at $\Gamma_{r}\gg 1$. We
can therefore think of $\Gamma_{r}$ as controlling the effective temperature
of the simulation. We note that using Eq. (7) as the training Hamiltonian and
using $\Gamma_{r}$ as the temperature is similar to schemes that use the free
energy to train the ANN [38, 39, 18], where instead of using the entropy we
use $\langle H_{\mathrm{T}}\rangle$. We leave it to future work to directly
compare these two approaches.

We choose our distribution of $\Gamma_{r}$ to be a cubic function of the
replica indices,

$\Gamma_{r}=10\left(\frac{r}{N-1}\right)^{3}\ ,\quad r=0,\dots,N-1\ .$ (8)

For our simulations, a $\Gamma_{N-1}=10$ value was sufficient for the
$H_{\mathrm{M}}$ to dominate over $H_{\mathrm{T}}$. We emphasize that the
choice of the distribution of $\Gamma_{r}$ values can strongly influence the
performance of the QPT algorithm; for example, choosing a linear function for
our simulations resulted in very poor performance.

Our algorithm proceeds as follows. The $N$ replicas (indexed by
$r=0,\dots,N-1$) are initialized randomly with both the real and imaginary
components of the parameters $\alpha_{r}$ chosen from a Gaussian distribution
with mean 0 and standard deviation $10^{-2}$. We perform 10 updates of SR for
all replicas, configuration swap between replicas $r$ and $r+1$ for $r$ even,
10 updates of SR for all replicas, and finally a configuration swap between
replicas $r$ and $r+1$ for $r$ odd. The choice of the number of updates of SR
between configuration swaps was not optimized, and other choices may improve
performance. This process is repeated until a total of $u$ updates are
performed. The deterministic even-odd pattern for configuration swaps is a
typical update pattern in PT that helps improve the diffusion of replicas [40,
41]. For every update, we calculate the expectation value of the target
Hamiltonian for each replica, which we denote by

$E_{r}\equiv\langle\psi(\alpha_{r})|H_{\mathrm{T}}|\psi(\alpha_{r})\rangle\ .$
(9)

We denote the running average of an expectation value of an observable $O$
over 10 updates of SR by $\overline{\langle O\rangle}$, which we will use to
calculate the probabilities for configuration swaps.

The final issue we must address about the QPT algorithm is how to perform the
swap updates between the $N$ replicas. We begin by reminding ourselves how the
swaps are performed in the case where each of the $N$ replicas is evolving
according to a Markov Chain Monte Carlo simulation that is ergodic and has a
unique equilibrium distribution. Assuming a finite state space, the $i$-th
configuration of replica $r$ has an equilibrium probability given by
$\Pi(\Gamma_{r})=W_{i}(\Gamma_{r})/\sum_{j}W_{j}(\Gamma_{r})$, where
$W_{i}(\Gamma_{r})$ is the unnormalized weight of the $i$-th configuration.
The swap updates between replicas then occurs with a probability $p$ chosen to
satisfy the detailed balance condition; specifically, if replica $r$ is in
configuration $i$ and replica $r+1$ is in configuration $i^{\prime}$, then the
probability of a swap between the two is taken to be:

$p_{r\leftrightarrow
r+1}=\min\left(1,\frac{W_{i}(\Gamma_{r+1})W_{i^{\prime}}(\Gamma_{r})}{W_{i}(\Gamma_{r})W_{i^{\prime}}(\Gamma_{r+1})}\right)\
.$ (10)

While we may think of the SR algorithm as implementing imaginary-time
evolution in a fixed subspace (see Appendix A), we are not aware of a way to
implement this evolution in terms of an ergodic MCMC simulation that
eventually converges to an equilibrium distribution. Therefore we are not able
to use the form in Eq. (10), and instead we will use physical intuition to
pick a swap update probability. This choice is not unique, and we expect
different choices to perform differently.

With these concerns in mind, we propose the following update. For
$\Gamma_{r}=0(r=0)$, we expect the ANN to reach a good approximation of an
eigenstate of $H_{\mathrm{T}}$, assuming the number ANN parameters is
sufficiently large and the ANN is sufficiently expressive. Similarly, for
replicas with $\Gamma_{r}\gg 1$ we expect the ANN to be a good approximation
of an eigenstate of $H_{\mathrm{M}}$. Therefore, as we have already mentioned,
we can think of $\Gamma_{r}$ as being a measure of a “noise” introduced, and
we can think of $\Gamma_{r}$ as being proportional to the temperature of the
replica. We thus propose the following probability rule for swap updates
between neighboring pairs $r$ and $r+1$:

$p_{r\leftrightarrow
r+1}=\min\left(1,\exp\left[\left(\overline{E_{r}}-\overline{E_{r+1}}\right)\left(\frac{1}{H_{0}\Gamma_{r}}-\frac{1}{H_{0}\Gamma_{r+1}}\right)\right]\right)\
.$ (11)

This choice for the probability rule is a heuristic choice motivated by the
Metropolis update [42] for the standard PT algorithm [20, 21, 22], but in no
way do we claim that it satisfies detailed balance [43]. The use of the target
energy comes from the desire to move configurations with lower target energies
to smaller $r$. For our temperature, we use $\Gamma_{r}$ multiplied by some
relevant energy scale $H_{0}$ as motivated by our discussion above. We use
$\overline{E_{r}}$ as opposed to simply ${E_{r}}$ in order to average out any
fluctuations in the estimation of the target energy during the ANN training.

In what follows, we define success as finding a configuration below some
relative error $\epsilon$ of the energy given by

$\epsilon=\left\lvert\frac{\langle E\rangle-
E_{\mathrm{GS}}}{E_{\mathrm{GS}}}\right\rvert\ ,$ (12)

where $E_{\mathrm{GS}}$ is the true ground state energy. To determine whether
QPT can provide practical advantages over its standard counterpart, we use the
time-to-epsilon (TT$\epsilon$) metric [44, 45], which is the time required to
find a solution below a target relative error $\epsilon$ at least once with
probability 0.99. The TT$\epsilon$ is given by

$\mathrm{TT}\epsilon=\frac{N}{N_{\mathrm{max}}}u\frac{\ln(1-0.99)}{\ln(1-p_{\mathrm{S}}(u))}\
,$ (13)

where $p_{\mathrm{S}}(u)$ is the probability of reaching the target $\epsilon$
within $u$ steps of the algorithms, such that the factor
$\frac{\ln(1-0.99)}{\ln(1-p_{\mathrm{S}}(u))}$ counts the number of
repetitions of the algorithm needed to ensure the target $\epsilon$ is reached
with probability $0.99$ at least once. $N_{\max}$ is the largest number of QPT
replicas we will use in our simulations, so that the factor $N/N_{\max}$ takes
into account the reduction in computational cost possible by parallelization
by exchanging replicas $N$ with independent runs of the algorithm. For every
pair of values $(n,N)$, we identify the value of $u$ that minimizes the
TT$\epsilon$, and we choose this as the cost of running the algorithm. Using
this approach we can determine unambiguously when it becomes more advantageous
to use QPT versus independent runs of the algorithm. We emphasize that this
metric does not account for the extra overhead associated with the
communication needed between replicas to implement the swap update, but we
expect this to be minimal compared to the cost of finding a high quality
solution for hard problems.

##  3 Results

###  3.1 Precipice Problem

In order to illustrate the computational viability of the QPT approach, we
study a problem class of qubit Hamiltonians that are invariant under any
permutation of the qubits. Specifically, we take our target Hamiltonian to be
of the form

$H_{\mathrm{T}}=\frac{1}{2}(1-s)\sum_{i=1}^{n}\left(\mathds{1}-\sigma_{i}^{x}\right)+s\sum_{x\in\left\\{0,1\right\\}^{n}}f(x)|x\rangle\langle
x|\ ,$ (14)

where $n$ is the total number of qubits and $f(x)$ is a function that only
depends on the Hamming weight of the bit-string $x$. We can therefore label
the energy eigenstates of $H_{\mathrm{T}}$ in terms of their symmetry
properties under qubit permutation. We take the mixer Hamiltonian to be given
by the uniform transverse field Hamiltonian,

$H_{\mathrm{M}}=-\frac{H_{0}}{2}\sum_{i=1}^{n}\left(\mathds{1}-\sigma_{i}^{x}\right),$
(15)

which is also qubit permutation invariant and we set $H_{0}=1$. Because the
Hamiltonian is stoquastic [46, 47], the ground state of this Hamiltonian can
be expressed entirely with non-negative amplitudes, so it must be in the
completely symmetric subspace, which has dimension $n+1$. We can enforce this
symmetry in our RBM ansatz

$\displaystyle\psi_{x}(\alpha)$ $\displaystyle=$ $\displaystyle
e^{a\sum_{i=1}^{n}v_{i}}\prod_{\mu=1}^{m}\cosh\left(b_{\mu}+W_{\mu}\sum_{i=1}^{n}v_{i}\right),$
(16)

where $v_{i}\equiv(1-2x_{i})\in\left\\{-1,1\right\\}$. This ansatz only
depends on the total magnetization of the spin configuration $v$, or
equivalently on the Hamming weight of the bit-string $x$. This will enforce
that all equal Hamming weight bit-strings will have the same amplitude. This
representation has only $K=1+2m$ parameters. In principle these parameters
need only be purely real to capture the ground state, but we will allow the
training to treat them as complex parameters.

In what follows, we will focus on a function $f(x)$ first introduced in Ref.
[48]

$f(x)=\left\\{\begin{array}[]{ll}-1&\text{if }x=1\cdots 1\\\
w(x)&\text{otherwise}\ ,\end{array}\right.$ (17)

where $w(x)$ denotes the Hamming weight of the bit-string $x$. We call this
problem the “Precipice” problem because of the form of $f(x)$, since it
increases linearly with Hamming weight until reaching Hamming weight $n$ where
it plummets to a value of $-1$.

For concreteness, we pick $s=4/5$ for our target Hamiltonian in Eq. (14). For
this value of $s$, $H_{T}$ is dominated by its diagonal component,
corresponding to the second term in Eq. (14). The ground state of the
corresponding target Hamiltonian $H_{\mathrm{T}}$ has most but not all of its
weight on the all-one bit-string, which is the bit-string that minimizes
$f(x)$. On the other hand, the first excited state of $H_{\mathrm{T}}$ has
most of its weight on the false minimum of $f(x)$ given by the all-zero bit-
string. The Hamiltonian $H_{\mathrm{T}}$ has the feature that annealing from
$\Gamma\gg 1$ to $\Gamma=0$ (Eq. (7)) encounters an exponentially closing
ground state energy gap as a function of the system size. One can understand
this intuitively as arising from the system requiring to tunnel between the
two configurations above (the state with most of its weight on the all-one bit
string and the state with most of its weight on the all-zero bit string) that
requires flipping $n$ spins [49].

This feature of the energy spectrum stymies the SR algorithm. Whether it will
find the the correct global minimum depends entirely on where the RBM
parameters are initialized, because once it reaches the false minimum
associated with the first excited state, the SR updates do not provide a means
to escape this local minimum and reach the global minimum. This is where our
QPT implementation helps, as it provides a mechanism for RBM configurations
with weight on the all-one bit-string to be reached and guide the training to
the true ground state. Since our aim is to distinguish between the two minima
that the algorithm can get stuck in, the choice of $\epsilon$ is not crucial
as long as it is sufficiently small. We therefore fix it to
$\epsilon=10^{-1}$, which is enough to distinguish between the two lowest
energy eigenstates.

(a)  (b)

Figure 1: Relative frequency of finding the ground state for the Precipice
problem to within a relative error of $\epsilon=0.1$ for a random
initialization of a trial RBM ansatz without using QPT as a function of (a)
number of updates and (b) system size $n$. We used $10^{3}$ independent trials
for $n=6,7,8$, $4\times 10^{3}$ for $n=9$ and $4\times 10^{4}$ for $n=10$. The
error bars correspond to the 95% confidence interval calculated using a
bootstrap with $10^{3}$ resamplings of the data. We describe this procedure in
Appendix E. Inset: Median number of updates to reach a relative error of
$\epsilon=0.1$ when the algorithm succeeds.

We begin by analyzing the behavior of the standard algorithm without QPT. We
show in Fig. 1 how the probability of finding the ground state within a
relative error of $\epsilon$ behaves as a function of system size. We see that
the probability of finding the ground state saturates to a value less than 1
as a function of the number of SR updates, and the saturation value decreases
faster than exponentially with system size. This behavior indicates that once
the system reaches the false minimum it is trapped there without the
possibility of escape, and the probability of initializing at a point where
the SR updates guide you to the correct minimum decays faster than
exponentially.

To show how QPT boosts the probability of finding the ground state, we show in
Fig. 2 how with increasing number of updates in QPT we can increase the
probability of finding the ground state. In order to compare these results to
the standard approach without QPT, we use the results at $n=10$ as a
benchmark. In the standard approach, the probability of finding the ground
state is approximately $2.5\times 10^{-3}$ with a median number of steps of
$440$ (Fig. 1). Therefore, to guarantee that we find the ground state at least
once with probability $0.99$, we would need to perform $1840$ independent
trials. With QPT, we can achieve almost a probability of $1$ of finding the
ground state using less than 500 steps and 10 replicas with a single
simulation. Therefore, we already see a significant advantage for using QPT in
this model.

Figure 2: Relative frequency of finding the ground state for the Precipice
problem to within $\epsilon=10^{-1}$ with increasing number of updates in QPT
for 50 independent simulations. Here we use $N=10$ replicas for all QPT
simulations.

To illustrate how the QPT algorithm is helping, we show in Fig. 3 an example
of a random walk for one of our simulation trials. In Fig. 3a, we see that
both replicas at $r=0$ and $r=1$ very quickly converge to the first excited
state and only after approximately $1500$ network updates does the $r=0$
replica find the ground state. As we can see in Fig. 3b, the transition to the
ground state is precipitated by a swap between the $r=0$ and $r=1$
configurations. While the configuration at $r=0$ remains unchanged for a large
fraction of the earlier updates, the configuration at $r=1$ is swapped
multiple times with other configurations, allowing it to explore the landscape
of possible solutions.

(a)

(b)

Figure 3: An example of a random walk for a single simulation of the Precipice
problem with $n=12$ and $N=10$. In (a), we show the target energy (Eq. (9)) of
the configuration at $r=0$ and $r=1$ during the simulation. The dotted lines
correspond to the ground state energy and first excited state energy for the
target Hamiltonian. In (b), we label each configuration by its original
replica index at initialization, and we show which index is at position $r=0$
and $r=1$ during the simulation. This allows us to monitor how configurations
walk on the space of $\Gamma_{r}$.

Finally, we show our results for the $\mathrm{TT}\epsilon$ in Fig. 4 with
$N_{\max}=40$. We see that at small system sizes ($<8$), using the standard
approach ($N=1$) is more advantageous than incurring the added computational
cost of QPT. Here the probability of reaching the ground state is sufficiently
high (even if it is not close to 1 as seen in Fig. 1) that running multiple
independent runs of the algorithm is the optimal strategy. However beyond
system sizes of $>8$, we see a significant advantage to using QPT. In this
case, the optimal strategy is to use QPT with more updates to achieve higher
probabilities of reaching the ground state.

We note that for intermediate values of $n$ we see similar performance for
$N=10$ and $N=20$, indicating that increasing the number of replicas does not
improve the success probability enough to outweigh the cost of simulating
additional replicas. However, we also see that at larger system sizes it does
become more advantageous to use more replicas. For example, at $n=17$ our
simulations with $N=10$ replicas fail to find the ground state, and at $n=18$
our simulations with $N=20$ replicas fail to find the ground state. Because
increasing $N$ changes the distribution of $\Gamma_{r}$ values and includes
more replicas at smaller $\Gamma_{r}$ values, the need to increase $N$ to find
the ground state at larger system sizes indicates the importance of optimizing
this distribution to get the best results.

Figure 4: Optimum time to epsilon to reach a relative error of
$\epsilon=10^{-1}$ for the Precipice problem as a function of system size $n$
for the Precipice problem with different number of QPT replicas $N$ with
$N_{\max}=40$. The error bars correspond to the 95% confidence interval
calculated using a bootstrap with $10^{3}$ resamplings of the data.

###  3.2 H4 Rectangle

The second system we study is that of four hydrogen atoms arranged in a
rectangle with an angle $\theta$ between adjacent atoms. This is an example of
a molecular system that is small enough to be exactly solved that includes a
type of configurational degeneracy that gives some approximate quantum
chemistry algorithms trouble [50, 51]. We consider a particular encoding of
the second quantized Hamiltonian into 8 qubits; further details of this system
and how we construct the qubit Hamiltonian are given in Appendix F.

The ground state of the 8-qubit Hamiltonian has non-zero amplitudes only on
computational basis states with Hamming Weight four or equivalently with four
fermions or zero $z$-magnetization. This observation motivates two different
RBM ansätze: (1) we restrict the input of the RBM to Hamming-weight-four
computational basis states and assume that all other inputs give zero
amplitude, and (2) we allow all inputs to the RBM and the training of the RBM
must find the Hamming weight-four subspace. While the first case is a more
efficient representation for the H4 Rectangle, we still consider both ansätze
since they could be relevant in different contexts for fermionic simulation
beyond the H4 Rectangle: the first applies to cases where particle number is
conserved, whereas the second applies to cases where the chemical potential is
fixed.

In what follows, we fix our desired error to $\epsilon=10^{-3}$. The reason
for this smaller choice of relative error, compared to the Precipice example,
is that this problem requires a relative error less than $6\times 10^{-3}$ for
$\theta=90^{\circ}$ to begin distinguishing the two lowest energy eigenstates.
So we choose an $\epsilon$ that is sufficiently small to ensure that we are
finding a high quality approximation of the ground state. For the simulations,
we use a up to $10^{4}$ updates.

####  3.2.1 Fixed Hamming Weight

We first consider the case where the RBM is only sampled with Hamming-weight-
four states. In this case, we choose a mixer Hamiltonian given by the
ferromagnetic Heisenberg model with all-to-all connectivity,

$H_{\mathrm{M}}=\frac{H_{0}}{2n}\sum_{i=1}^{n}\sum_{j=i+1}^{n}\left(\mathds{1}-\sigma_{i}^{x}\sigma_{j}^{x}-\sigma_{i}^{y}\sigma_{j}^{y}-\sigma_{i}^{z}\sigma_{j}^{z}\right)\
,$ (18)

for which the (degenerate) ground state is the uniform superposition of fixed
Hamming weight states with energy 0 and $H_{0}$ is the absolute magnitude of
the largest coupling strength in the target Hamiltonian. This choice is made
to ensure that the energy scales of $H_{\mathrm{M}}$ and $H_{\mathrm{T}}$ are
comparable.

For the simulations, we find that using $m=8$ hidden units is sufficient to
represent the ground state with the required accuracy, so we use this value
for all our training simulations with a fixed Hamming weight exact sampling.
We show in Fig. 5 the results for the optimal TT$\epsilon$ for three different
$\theta$ values, where we find rich behavior as a function of $\theta$. We
find that for a sufficient number of replicas, the QPT approach demonstrates
an advantage over the standard approach for angles $\theta=80^{\circ}$ and
$85^{\circ}$. For $\theta=90^{\circ}$, the ground state is found readily by
the standard approach that the additional cost of running $N$ replicas does
not introduce any benefit. The $\theta=90^{\circ}$ instance was found to be
pathological for the coupled cluster doubles (CCD) approximation but not its
variational counterpart in Ref. [51]. It is therefore unsurprising that a
variational RBM ansatz that explores a similar Hilbert space captures the
ground state well. Because the root of the CCD pathology is ultimately
degeneracy of two of the dominant configurations, it seems like this
degeneracy might provide an advantage for the RBM ansatz training that is
worthy of further investigation.

We also note that this particular radius provides a uniquely easy energy
landscape for the single-replica approach at $90^{\circ}$, and at different
radii the $90^{\circ}$ case is more difficult. In fact, among the instances
considered, the instance that we are highlighting was found to be the one for
which achieving an advantage with QPT was most difficult. At
$\theta=80^{\circ}$ and $85^{\circ}$, we observe that $N=5$, $N=10$ and $N=20$
replicas exhibit a definitive advantage over the single-replica approach.

Figure 5: Optimum time to epsilon to reach a relative error of
$\epsilon=10^{-3}$ for the H4 Rectangle at different angles with fixed Hamming
weight exact sampling. We use $N_{\max}=20$. The error bars correspond to the
95% confidence interval calculated using a bootstrap with $10^{3}$ resamplings
of the data. The lines are to guide the eye.

####  3.2.2 No Fixed Hamming Weight

We now consider the case where the RBM is sampled with all possible
computational basis states. In this case, we choose a mixer Hamiltonian to be
the uniform transverse field:

$H_{\mathrm{M}}=\frac{H_{0}}{2}\sum_{i=1}^{n}\left(\mathds{1}-\sigma_{i}^{x}\right)\
,$ (19)

where $H_{0}$ is again the absolute magnitude of the largest coupling strength
in the target Hamiltonian.

For the simulations, we find that using $m=48$ hidden units is sufficient to
represent the ground state with the required accuracy but also allow the
simulations with $N=1$ to find the ground state with sufficient frequency, so
we use this value for these simulations111We can actually use fewer hidden
units (as low as $m=16$) and still have an approximation of the ground state
to the desired accuracy, but in this case we only find a good approximation to
the ground state using our QPT method when using our simulation parameters..
We show in Fig. 6 the results for the optimal TT$\epsilon$ for the same three
values of $\theta\in\left\\{80^{\circ},85^{\circ},90^{\circ}\right\\}$. While
generally the problem of approximating the ground state is harder in this case
compared to the fixed Hamming weight case, the advantage for QPT training with
a sufficient number of replicas for all three angles is clear, indicating that
the energy landscape is more efficiently searched using the QPT training in
this case. This is a promising indication of the advantages of the QPT
approach for more generic problems that may not exhibit any symmetries that
can be used to restrict the wave function ansätze.

Figure 6: Optimum time to epsilon to reach a relative error of
$\epsilon=10^{-3}$ for the H4 Rectangle at different angles with full exact
sampling. We use $N_{\max}=20$. The error bars correspond to the 95%
confidence interval calculated using a bootstrap with $10^{3}$ resamplings of
the data. The lines are to guide the eye.

##  4 Conclusions

We have shown that using a QPT method as part of the training of ANNs to
approximate ground states of quantum Hamiltonians can help overcome barriers
in the energy landscape. To illustrate this, we simplify our analysis by using
exact sampling in order to eliminate the role of fluctuations introduced by
finite sampling. We demonstrated the utility of the QPT approach using two
Hamiltonians with very different profiles. For the permutation invariant
problem, the QPT approach demonstrated a clear advantage over repetitions of
the standard algorithm even for small system sizes. A possible criticism of
this problem is that it is described by a stoquastic Hamiltonian [46, 47].
Stoquastic Hamiltonians have ground state that can be described in terms of
non-negative real amplitudes, and they admit efficient quantum-to-classical
mappings [52]. While these properties do not necessarily mean there is an
efficient algorithm to find the ground state, the existence of an efficient
classical representation is often taken to mean that the ground states of
stoquastic Hamiltonians are not as difficult to find as those of non-
stoquastic Hamiltonians. We therefore consider a second problem based on four
hydrogen atoms arranged in a rectangle with a variable angle between adjacent
atoms, and we observe a clear advantage for the QPT approach at some angles
where the standard approach struggled to find the ground state. These
promising initial results indicate that using the QPT approach to train ANNs
may be useful for a broad range of Hamiltonians with difficult to explore
parameter landscapes.

We focused on using an RBM as our ANN ansatz, but the implementation of QPT is
independent of the choice of ANN architecture and can be applied equally well
to convolutional neural networks (CNNs) [16], to two disconnected real ANNs to
represent the magnitude and phase of wave function amplitudes [19], or to
group convolutional neural networks (GCNNs) [18], which have also been shown
to be efficient at representing quantum states. This makes our proposed
approach quite versatile, with only the added complexity of training multiple
replicas of the ANNs and implementing configuration swaps between them.

We emphasize that we made several choices in this work for simplicity that
might have been suboptimal for the performance of the QPT approach. For
example, we chose a cubic distribution of the parameter $\Gamma$ for our
replicas (Eq. (8)), and a suboptimal distribution of $\Gamma$ values around a
bottleneck for the simulation does hinder the random walk of replicas [53]. We
conjecture that optimizing this distribution will further increase the
advantage of the QPT approach in our problems. We leave it to future work to
develop a methodology to identify more optimal distributions.

Furthermore, our choice of swap probability (Eq. (11)) was based on an
intuitive argument, and we provided examples of other possible choices in
Appendix D. It would be ideal to identify a choice that can satisfy some
appropriate notion of detailed balance. We believe that this is the most
important open theoretical problem with our proposal. Such a formalism would
also allow us to adopt some of the strategies and methods using in conjunction
with QPT to characterize the hardness of the landscape, such as using the
autocorrelation time of the replica random walk as a proxy for the mixing time
[23].

##  5 Code

Code for this project can be found at https://github.com/talbash/QPTforRBM.

## Acknowledgements

We would like to thank the UNM Center for Advanced Research Computing,
supported in part by the National Science Foundation, for providing the high
performance computing resources used in this work.

##### Funding information

This material is based upon work supported by the National Science Foundation
under Grant No. 2037755. This material is based upon work supported by the Air
Force Office of Scientific Research under award number FA9550-22-1-0498. Any
opinions, findings, and conclusions or recommendations expressed in this
material are those of the author(s) and do not necessarily reflect the views
of the United States Air Force. CS and QC were supported by Sandia National
Laboratories’ Laboratory Directed Research and Development program (Project
222396). Sandia National Laboratories is a multi-mission laboratory managed
and operated by National Technology and Engineering Solutions of Sandia, LLC,
a wholly owned subsidiary of Honeywell International, Inc., for DOE’s National
Nuclear Security Administration under contract DE-NA0003525. This paper
describes objective technical results and analysis. Any subjective views or
opinions that might be expressed in the paper do not necessarily represent the
views of the U.S. Department of Energy or the United States Government.

##  Appendix A Review of Stochastic Reconfiguration

We consider a (possibly unnormalized) wave function $|\psi(\alpha)\rangle$
parameterized by $K$ complex parameters denoted by
$\alpha=\left(\alpha_{1},\dots,\alpha_{K}\right)^{T}$. We are interested in
finding the parameter values ${\alpha}$ that minimize

$E(\alpha)\equiv\frac{\langle\psi(\alpha)|H|\psi(\alpha)\rangle}{\langle\psi(\alpha)|\psi(\alpha)\rangle}=\langle
H\rangle\ .$ (20)

The true minimum is given by the ground state energy $E_{\mathrm{GS}}$.
Because the parameterization of the wave function $|\psi(\alpha)\rangle$ is
unlikely to be general enough to capture the ground state _exactly_ , we can
only hope for $|\psi(\alpha)\rangle$ (when normalized) to provide an
approximation to the ground state.

Assume the wave function has an expansion in a basis of the form:
$|\psi(\alpha)\rangle=\sum_{x}\psi_{x}(\alpha)|x\rangle$. For simplicity, we
will assume that $\psi_{x}(\alpha)$ only depends on $\alpha_{k}$ and not on
both $(\alpha_{k},\alpha_{k}^{\ast})$. Let us consider a small variation of
the parameters $\alpha_{k}^{\prime}=\alpha_{k}+\delta\alpha_{k}$ Let us expand
the wave function:

$\displaystyle|\psi(\alpha^{\prime})\rangle$ $\displaystyle=$
$\displaystyle|\psi(\alpha)\rangle$ (21)
$\displaystyle+\sum_{k}\delta\alpha_{k}\sum_{x}\left(\frac{1}{\psi_{x}(\alpha)}\frac{\partial}{\partial\alpha_{k}}\psi_{x}(\alpha)\right)\psi_{x}(\alpha)|x\rangle$
$\displaystyle+\frac{1}{2}\sum_{k,k^{\prime}}\delta\alpha_{k}\delta\alpha_{k^{\prime}}\left(\frac{1}{\psi_{x}(\alpha)}\frac{\partial}{\partial\alpha_{k}}\frac{\partial}{\partial\alpha_{k^{\prime}}}\psi_{x}(\alpha)\right)\psi_{x}(\alpha)|x\rangle$
$\displaystyle+\dots..$

It is convenient to define the operators:

$\displaystyle O_{k}(\alpha)$
$\displaystyle=\sum_{x}\frac{1}{\psi_{x}(\alpha)}\frac{\partial}{\partial\alpha_{k}}\psi_{x}(\alpha)|x\rangle\langle
x|\ ,$ (22a) $\displaystyle M_{kk^{\prime}}(\alpha)$
$\displaystyle=\sum_{x}\frac{1}{\psi_{x}(\alpha)}\frac{\partial}{\partial\alpha_{k}}\frac{\partial}{\partial\alpha_{k^{\prime}}}\psi_{x}(\alpha)|x\rangle\langle
x|\ .$ (22b)

(We will suppress their $\alpha$ dependence for notational brevity) such that
we can write

$\displaystyle|\psi(\alpha^{\prime})\rangle$ $\displaystyle=$
$\displaystyle|\psi(\alpha)\rangle+\sum_{k}\delta\alpha_{k}O_{k}|\psi(\alpha)\rangle$
(23)
$\displaystyle+\frac{1}{2}\sum_{k,k^{\prime}}\delta\alpha_{k}\delta\alpha_{k^{\prime}}M_{kk^{\prime}}|\psi(\alpha)\rangle+\dots$
$\displaystyle=$ $\displaystyle\left(1+\sum_{k}\delta\alpha_{k}\langle
O_{k}\rangle+\dots\right)|\psi(\alpha)\rangle$
$\displaystyle+\sum_{k}\delta\alpha_{k}\left(O_{k}-\langle
O_{k}\rangle\right)|\psi(\alpha)\rangle$
$\displaystyle+\frac{1}{2}\sum_{k,k^{\prime}}\delta\alpha_{k}\delta\alpha_{k^{\prime}}\left(M_{kk^{\prime}}-\langle
M_{kk^{\prime}}\rangle\right)|\psi(\alpha)\rangle$ $\displaystyle+\dots\ ,$

where

$\displaystyle\langle O_{k}(\alpha)\rangle$
$\displaystyle=\frac{\langle\psi(\alpha)|O_{k}(\alpha)|\psi(\alpha)\rangle}{\langle\psi(\alpha)|\psi(\alpha)\rangle}\
,$ (24a) $\displaystyle\langle M_{kk^{\prime}}(\alpha)\rangle$
$\displaystyle=\frac{\langle\psi(\alpha)|M_{kk^{\prime}}(\alpha)|\psi(\alpha)\rangle}{\langle\psi(\alpha)|\psi(\alpha)\rangle}\
.$ (24b)

We will denote the coefficient of the $|\psi(\alpha)\rangle$ term in Eq. (23)
by $\delta\alpha_{0}$, such that
$\langle\psi(\alpha)|\psi(\alpha^{\prime})\rangle=\delta\alpha_{0}\langle\psi(\alpha)|\psi(\alpha)\rangle$
and

$\displaystyle\langle\psi(\alpha^{\prime})|\psi(\alpha^{\prime})\rangle$
$\displaystyle=$
$\displaystyle|\delta\alpha_{0}|^{2}\langle\psi(\alpha)|\psi(\alpha)\rangle$
(25)
$\displaystyle+\sum_{k,k^{\prime}}\delta\alpha_{k}^{\ast}\delta\alpha_{k^{\prime}}\langle\psi(\alpha)|\left(O_{k}^{\dagger}-\langle
O_{k}^{\dagger}\rangle\right)\left(O_{k^{\prime}}-\langle
O_{k^{\prime}}\rangle\right)|\psi(\alpha)\rangle$ $\displaystyle+\dots$
$\displaystyle=|\delta\alpha_{0}|^{2}\left(1+\frac{1}{|\delta\alpha_{0}|^{2}}\sum_{k,k^{\prime}}\delta\alpha_{k}^{\ast}\delta\alpha_{k^{\prime}}S_{k,k^{\prime}}+\dots\right)$
$\displaystyle\times\langle\psi(\alpha)|\psi(\alpha)\rangle\ ,$

where we have defined the covariance matrix

$\displaystyle S_{kk^{\prime}}$ $\displaystyle=$
$\displaystyle\frac{\langle\psi(\alpha)|\left(O_{k}^{\dagger}-\langle
O_{k}^{\dagger}\rangle\right)\left(O_{k^{\prime}}-\langle
O_{k^{\prime}}\rangle\right)|\psi(\alpha)\rangle}{\langle\psi(\alpha)|\psi(\alpha)\rangle}$
(26) $\displaystyle=$ $\displaystyle\langle
O_{k}^{\dagger}O_{k^{\prime}}\rangle-\langle O_{k}^{\dagger}\rangle\langle
O_{k^{\prime}}\rangle\ .$

Here we have suppressed the dependence of $S$ on $\alpha$. We see from Eq.
(25) that renormalizing $|\psi(\alpha^{\prime})\rangle$ by $\delta\alpha_{0}$
amounts to rescaling $\delta\alpha_{k}$ by $\delta\alpha_{0}$. We will exploit
this freedom later. Finally, for completeness, we note that we can expand the
energy as:

$\displaystyle E(\alpha^{\prime})$ $\displaystyle=$ $\displaystyle E(\alpha)$
(27)
$\displaystyle+\frac{1}{\langle\psi(\alpha)|\psi(\alpha)\rangle}\sum_{k}\left(\delta\alpha_{k}^{\ast}\langle\psi(\alpha)|\left(O_{k}^{\dagger}-\langle
O_{k}^{\dagger}\rangle\right)H|\psi(\alpha)\rangle\right.$
$\displaystyle\left.+\delta\alpha_{k}\langle\psi(\alpha)|H\left(O_{k}-\langle
O_{k}\rangle\right)|\psi(\alpha)\rangle\right)+\dots\ ,$ $\displaystyle=$
$\displaystyle
E(\alpha)+\sum_{k}\left(\delta\alpha_{k}^{\ast}F_{k}+\delta\alpha_{k}F_{k}^{\ast}\right)+\dots\
,$

where we have defined

$\displaystyle F_{k}$ $\displaystyle=$
$\displaystyle\frac{\langle\psi(\alpha)|\left(O_{k}^{\dagger}-\langle
O_{k}^{\dagger}\rangle\right)\left(H-E(\alpha)\right)|\psi(\alpha)\rangle}{\langle\psi(\alpha)|\psi(\alpha)\rangle}$
(28) $\displaystyle=$ $\displaystyle\langle O_{k}^{\dagger}H\rangle-\langle
O_{k}^{\dagger}\rangle\langle H\rangle\ .$

Here we have suppressed the dependence of $F$ on $\alpha$.

For the SR algorithm we further demand that:

$|\psi(\alpha^{\prime})\rangle=P\left(\Lambda-H\right)|\psi(\alpha)\rangle\ ,$
(29)

where $\Lambda$ is taken to be a sufficiently large (larger than the largest
eigenvalue of $H)$ positive number and $P$ is the projector onto the subspace
of the parameterization. By acting with $\langle\psi(\alpha)|$ and
$\langle\psi(\alpha)|\left(O_{k}^{\dagger}-\langle
O_{k}^{\dagger}\rangle\right)$, and by restricting to linear order in
$\delta\alpha$, we get $K+1$ linear equations:

$\displaystyle\delta\alpha_{0}=\Lambda-E(\alpha)$ (30a)
$\displaystyle\sum_{k^{\prime}}\delta\alpha_{k^{\prime}}\langle\psi(\alpha)|\left(O_{k}^{\dagger}-\langle
O_{k}^{\dagger}\rangle\right)\left(O_{k^{\prime}}-\langle
O_{k^{\prime}}\rangle\right)|\psi(\alpha)\rangle$ $\displaystyle=$
$\displaystyle-\langle\psi(\alpha)|\left(O_{k}^{\dagger}-\langle
O_{k}^{\dagger}\rangle\right)H|\psi(\alpha)\rangle$
$\displaystyle=-\langle\psi(\alpha)|\left(O_{k}^{\dagger}-\langle
O_{k}^{\dagger}\rangle\right)\left(H-E(\alpha)\right)|\psi(\alpha)\rangle\ ,$
(30b)

we get the linear equation:

$-F_{k}=\sum_{k^{\prime}}S_{kk^{\prime}}\delta\alpha_{k^{\prime}}\Rightarrow\delta\alpha_{k}=-\sum_{k^{\prime}}S_{kk^{\prime}}^{-1}F_{k^{\prime}}\
.$ (31)

Since $S$ may not be invertible, $S^{-1}$ is strictly speaking the Moore-
Penrose pseudoinverse. Therefore, the update rule for the parameters $\alpha$
is given by:

$\alpha^{\prime}_{k}=\alpha_{k}+\delta\alpha_{k}=\alpha_{k}-\sum_{k^{\prime}}S_{kk^{\prime}}^{-1}F_{k^{\prime}}\
.$ (32)

Furthermore, if we renormalize the state $\psi(\alpha^{\prime})$ by
$\delta\alpha_{0}$, we effectively are rescaling $\delta\alpha_{k}$ by
$\delta\alpha_{0}$. Therefore, our update rule after this renormalization is
given by:

$\alpha^{\prime}_{k}=\alpha_{k}+\frac{\delta\alpha_{k}}{\delta\alpha_{0}}=\alpha_{k}-\frac{1}{\delta\alpha_{0}}\sum_{k^{\prime}}S_{kk^{\prime}}^{-1}F_{k^{\prime}}\
,$ (33)

and we have (from Eq. (25)):

$\displaystyle\frac{\langle\psi(\alpha^{\prime})|\psi(\alpha^{\prime})\rangle-\langle\psi(\alpha)|\psi(\alpha)\rangle}{\langle\psi(\alpha)|\psi(\alpha)\rangle}$
$\displaystyle=$ (34)
$\displaystyle\frac{1}{|\delta\alpha_{0}|^{2}}\sum_{k,k^{\prime}}\delta\alpha_{k}^{\ast}s_{k,k^{\prime}}\delta\alpha_{k^{\prime}}+\dots\
.$

By repeating this iterative scheme, convergence is reached when the ratio
$\delta\alpha_{k}/\delta\alpha_{0}\to 0$: we have from Eq. (27):

$\displaystyle E(\alpha^{\prime})-E(\alpha)=$
$\displaystyle-\frac{1}{|\delta\alpha_{0}|^{2}}\sum_{k,k^{\prime}}\left(F_{k}^{\ast}S_{kk^{\prime}}^{-1}F_{k^{\prime}}+\left(F_{k}^{\ast}S_{kk^{\prime}}^{-1}F_{k^{\prime}}\right)^{\ast}\right)\
.$ (35)

Because the covariance matrix is a positive semi-definite matrix, the term in
the sum is positive, which says that on every iteration the energy decreases,
and in the limit of $\delta\alpha_{k}/\delta\alpha_{0}\to 0$, the energy
ceases to change.

As we noted earlier, there is some freedom in picking the value of $\Lambda$.
Since $\Lambda$ controls the value of $\delta\alpha_{k}$, which in turn
controls both the parameter update rate (Eq. (33)) and the wave function
change rate (Eq. (34)), in practice we treat
$\delta\alpha_{0}^{-1}\equiv\gamma$ as a free parameter that is chosen to get
a stable convergence. The value of $\gamma$ is the learning rate. As we
describe in Appendix B, we will use an adaptive scheme to update this learning
rate.

In the case of the RBM ansatz, the operators $O_{k}$ can be calculated
analytically. The relevant quantities are given by:

$\displaystyle\frac{1}{\psi_{x}(\alpha)}\frac{\partial}{\partial
a_{i}}\psi_{x}(\alpha)$ $\displaystyle=(1-2x_{i})\ ,$ (36a)
$\displaystyle\frac{1}{\psi_{x}(\alpha)}\frac{\partial}{\partial
b_{\mu}}\psi_{x}(\alpha)$ $\displaystyle=\tanh\left(\theta(x)\right)\ ,$ (36b)
$\displaystyle\frac{1}{\psi_{x}(\alpha)}\frac{\partial}{\partial
W_{i\mu}}\psi_{x}(\alpha)$
$\displaystyle=(1-2x_{i})\tanh\left(\theta(x)\right)\ .$ (36c)

##  Appendix B Adaptive Learning Rate Scheme

We begin by interpreting Eq. (3) as being a discretization of a dynamical
equation for $\alpha$ with $\gamma$ being the step size; therefore, let us
recast it as:

$\frac{d}{d\gamma}\alpha_{k}(\gamma)=-\sum_{k^{\prime}}(S(\alpha(\gamma))^{-1})_{kk^{\prime}}F_{k^{\prime}}(\alpha(\gamma))\equiv
f_{k}(\alpha(\gamma))\ ,$ (37)

where we have introduced the dependence of $\alpha$ on $\gamma$. In this form,
we can use a Heun second order consistent integrator to adaptively update
$\gamma$. We calculate two quantities:

$\displaystyle\vec{\kappa}_{1}$ $\displaystyle=\vec{f}(\vec{\alpha}(\gamma))\
,$ (38a) $\displaystyle\vec{\kappa}_{2}$
$\displaystyle=\vec{f}(\vec{\alpha}(\gamma)+\Delta\gamma\vec{\kappa}_{1})\ ,$
(38b)

and then update using:

$\vec{\alpha}(\gamma+\Delta\gamma)=\vec{\alpha}(\gamma)+\frac{\Delta\gamma}{2}\left(\vec{\kappa}_{1}+\vec{\kappa}_{2}\right)$
(39)

In order to dynamically adjust the step $\Delta\gamma$, we calculate the ANN
parameters at the same time $\gamma+\Delta\gamma$ but now using two
$\Delta\gamma/2$ steps. Therefore, we define:

$\displaystyle\vec{\kappa}_{3}$
$\displaystyle=\vec{f}(\vec{\alpha}(\gamma)+\frac{\Delta\gamma}{2}\vec{\kappa}_{1})\
,$ (40a) $\displaystyle\vec{\kappa}_{4}$
$\displaystyle=\vec{f}(\vec{\alpha}(\gamma)+\frac{\Delta\gamma}{4}\left(\vec{\kappa}_{1}+\vec{\kappa_{3}}\right))\
,$ (40b) $\displaystyle\vec{\kappa}_{5}$
$\displaystyle=\vec{f}(\vec{\alpha}(\gamma)+\frac{\Delta\gamma}{4}\left(\vec{\kappa}_{1}+\vec{\kappa_{3}}\right)+\frac{\Delta\gamma}{2}\vec{\kappa}_{4})\
,$ (40c)

such that

$\vec{\alpha}^{\prime}(\gamma+\Delta\gamma)=\vec{\alpha}(\gamma)+\frac{\Delta\gamma}{4}\left(\vec{\kappa}_{1}+\vec{\kappa_{3}}\right)+\frac{\Delta\gamma}{4}\left(\vec{\kappa}_{4}+\vec{\kappa}_{5}\right)\
.$ (41)

We define
$\Delta\alpha=\lVert\vec{\alpha}^{\prime}(\gamma+\Delta\gamma)-\vec{\alpha}(\gamma+\Delta\gamma)\rVert_{S(\alpha(\gamma))}$,
where the norm here is defined as:

$\lVert\vec{x}\rVert_{S(\alpha(\gamma))}=\frac{1}{M}\sqrt{\sum_{k,k^{\prime}=1}^{M}x_{k}^{\ast}S_{kk^{\prime}}(\alpha(\gamma))x_{k^{\prime}}}\
.$ (42)

The adjusted time step is then taken to be [17]:

$\Delta\gamma^{\prime}=\Delta\gamma\left(\frac{6\epsilon}{\Delta{\alpha}}\right)^{1/3}\
,$ (43)

for a fixed tolerance $\epsilon$. In our simulations, we fix
$\epsilon=10^{-6}$. We note that in a single update step in this framework
requires performing 5 different calculations of $S$ and $F$, one for each of
the different $\kappa$ values.

##  Appendix C Calculating the expectation value of operators with respect to
an RBM ansatz

We begin by first considering a diagonal operator $D$, meaning it is diagonal
in the computational basis. In this case, the expectation value of the
operator can be straightforwardly calculated using the RBM ansatz in Eq. (1)
(we will suppress the $\alpha$ dependence):

$\displaystyle\langle D\rangle$ $\displaystyle=$
$\displaystyle\frac{1}{\mathcal{Z}}\sum_{x}|\psi_{x}|^{2}\langle x|D|x\rangle\
,$ (44)

where $\mathcal{Z}=\langle\psi|\psi\rangle=\sum_{x}|\psi_{x}|^{2}$. Note that
$|\psi\rangle/\sqrt{\mathcal{Z}}$ defines a normalized state.

Let us now consider the expectation value of an arbitrary Pauli operator
acting on $n$ qubits, which we denote by $P$. We can express the expectation
value as:

$\displaystyle\langle P\rangle$ $\displaystyle=$
$\displaystyle\frac{1}{\mathcal{Z}}\sum_{x,x^{\prime}}\psi_{x}^{\ast}\psi_{x^{\prime}}\langle
x|P|x^{\prime}\rangle$ (45) $\displaystyle=$
$\displaystyle\frac{1}{\mathcal{Z}}\sum_{x}|\psi_{x}|^{2}\left(\sum_{x^{\prime}}\frac{\psi_{x^{\prime}}}{\psi_{x}}\langle
x|P|x^{\prime}\rangle\right)\ ,$

where the term in parenthesis defines the local Pauli operator $P(x)$:

$P(x)=\sum_{x^{\prime}}\frac{\psi_{x^{\prime}}}{\psi_{x}}\langle
x|P|x^{\prime}\rangle\ .$ (46)

For a given $P$, let $\mathcal{I}_{\alpha}$, $\alpha=x,y,z$, denote the qubit
indices where a single-qubit Pauli operator $\sigma^{\alpha}$ acts. For
example, consider the Pauli operator acting on 3 qubits
$P=\sigma_{2}^{x}\otimes\sigma_{1}^{y}\otimes\sigma^{z}_{0}$; in this case we
would have $\mathcal{I}_{x}=\left\\{2\right\\}$,
$\mathcal{I}_{y}=\left\\{1\right\\}$, $\mathcal{I}_{z}=\left\\{0\right\\}$.
For the RBM ansatz in Eq. (1), we can calculate the local Pauli operator
straightforwardly:

$\displaystyle P(x)$ $\displaystyle=$
$\displaystyle\left[\prod_{j\in\mathcal{I}_{z}}\left(-1\right)^{x_{j}}\right]\left[i^{\lvert\mathcal{I}_{y}\rvert}\prod_{j\in\mathcal{I}_{y}}\left(-1\right)^{1-x_{j}}\right]$
$\displaystyle\times\exp\left(-2\sum_{j\in\mathcal{I}_{x}\cup\mathcal{I}_{y}}a_{j}(1-2x_{j})\right)$
$\displaystyle\times\prod_{\mu=1}^{m}\frac{\cosh\left(\theta_{\mu}(x)-2\sum_{j\in\mathcal{I}_{x}\cup\mathcal{I}_{y}}W_{j\mu}(1-2x_{j})\right)}{\cosh\left(\theta_{\mu}(x)\right)}\
,$

where $\theta_{\mu}(x)=b_{\mu}+\sum_{j=1}^{n}W_{j\mu}(1-2x_{j})$.

The above prescription allows us to calculate the expectation values in the
covariance matrix and force vector in Eq. (4). For example, let us consider
the term $\langle O_{k}^{\dagger}H\rangle$. We can write:

$\displaystyle\langle O_{k}^{\dagger}H\rangle$ $\displaystyle=$
$\displaystyle\frac{1}{\mathcal{Z}}\sum_{x,x^{\prime}}\psi_{x}^{\ast}\psi_{x^{\prime}}\langle
x|O_{k}^{\dagger}H|x^{\prime}\rangle$ (48) $\displaystyle=$
$\displaystyle\frac{1}{\mathcal{Z}}\sum_{x,x^{\prime},x^{\prime\prime}}\psi_{x}^{\ast}\psi_{x^{\prime}}\langle
x|O_{k}^{\dagger}|x^{\prime\prime}\rangle\langle
x^{\prime\prime}|H|x^{\prime}\rangle$ $\displaystyle=$
$\displaystyle\frac{1}{\mathcal{Z}}\sum_{x}|\psi_{x}|^{2}\left(\sum_{x^{\prime\prime}}\frac{\psi_{x^{\prime\prime}}}{\psi_{x}}\langle
x|O_{k}^{\dagger}|x^{\prime\prime}\rangle\right.$
$\displaystyle\times\left.\sum_{x^{\prime}}\frac{\psi_{x^{\prime}}}{\psi_{x^{\prime\prime}}}\langle
x^{\prime\prime}|H|x^{\prime}\rangle\right)$ $\displaystyle=$
$\displaystyle\frac{1}{\mathcal{Z}}\sum_{x}|\psi_{x}|^{2}\langle
x|O_{k}^{\dagger}|x\rangle H(x)\ ,$

where we have used that the operators $O_{k}$ are diagonal and where we
defined a local Hamiltonian:

$H(x)=\sum_{x^{\prime}}\frac{\psi_{x^{\prime}}}{\psi_{x}}\langle
x|H|x^{\prime}\rangle\ .$ (49)

Since any Hamiltonian can be expressed in terms of diagonal operators and
arbitrary Pauli operators, we can calculate the local Hamiltonian using our
prescription above for diagonal and local Pauli operators.

In the case of the completely symmetric RBM ansatz, the expectation values
take a simple form. For example, we can express the contribution of the
transverse field as:

$\displaystyle\sum_{i=1}^{n}\left\langle\sigma_{i}^{x}\right\rangle$
$\displaystyle=$ $\displaystyle\sum_{w=0}^{n}{n\choose
w}|\psi_{w}(\alpha)|^{2}$ (50)
$\displaystyle\times\left(n\delta_{w0}\frac{\psi_{w+1}(\alpha)}{\psi_{w}(\alpha)}+n\delta_{wn}\frac{\psi_{w-1}(\alpha)}{\psi_{w}(\alpha)}\right.$
$\displaystyle\left.+\delta_{0<w<n}\left((n-w)\frac{\psi_{w+1}(\alpha)}{\psi_{w}(\alpha)}+w\frac{\psi_{w-1}(\alpha)}{\psi_{w}(\alpha)}\right)\right)\
,$

where we have defined $\psi_{w}(\alpha)$ as

$\displaystyle\psi_{w}(\alpha)$ $\displaystyle=$ $\displaystyle
e^{a(n-2w)}\prod_{\mu=1}^{m}\cosh\left(b_{\mu}+W_{\mu}(n-2w)\right)\ .$ (51)

We can also calculate:

$\displaystyle\sum_{i\neq j}^{n}\langle\sigma_{i}^{x}\sigma_{j}^{x}\rangle$
$\displaystyle=$ $\displaystyle\sum_{w=0}^{n}{n\choose
w}|\psi_{w}(\alpha)|^{2}$ (52)
$\displaystyle\times\left[\delta_{w0}n(n-1)\frac{\psi_{w+2}}{\psi_{w}}+\delta_{wn}n(n-1)\frac{\psi_{w-2}}{\psi_{w}}\right.$
$\displaystyle\left.+\delta_{w1}\left(2(n-1)+(n-1)(n-2)\frac{\psi_{w+2}}{\psi_{w}}\right)\right.$
$\displaystyle\left.+\delta_{w,n-1}\left(2(n-1)+(n-1)(n-2)\frac{\psi_{w-2}}{\psi_{w}}\right)\right.$
$\displaystyle\left.+\delta_{1<w<n-1}\left(2w(n-w)+(n-w)(n-w-1)\frac{\psi_{w+2}}{\psi_{w}}\right.\right.$
$\displaystyle\left.\left.+w(w-1)\frac{\psi_{w-2}}{\psi_{w}}\right)\right]\ .$

##  Appendix D Alternative Parallel Tempering Schemes

###  D.1 Standard deviation of $H^{2}$ as a measure of temperature

As we indicated in the main text, there is a lot of freedom in choosing the
training Hamiltonian and the parallel tempering update scheme. Here we present
an alternative approach. We can choose the Hamiltonian for the $r$th replica
to be given by

$H(\Gamma_{r})=(1-\Gamma_{r})H_{\mathrm{M}}+\Gamma_{r}H_{\mathrm{T}}\ ,$ (53)

where $H_{\mathrm{T}}$ is the target Hamiltonian whose ground state we are
actually after and $H_{\mathrm{M}}$ is a suitably chosen mixer Hamiltonian
with a trivial ground state. For example, in the absence of any constraints,
the mixer can be taken to be the uniform transverse field Hamiltonian,
$H_{\mathrm{M}}=-H_{0}\sum_{i=1}^{n}\sigma_{i}^{x}$, where $H_{0}$ is a
characteristic energy scale of $H_{\mathrm{T}}$. In this case, for replicas
with $\Gamma_{r}$ close to 0 the training Hamiltonian is dominated by the
transverse field, and the ground state is the uniform superposition state.
This is analogous to the high-temperature limit in “standard” PT, where the
different replicas operate at different temperatures but with the same
Hamiltonian. For replicas with $\Gamma_{r}$ close to 1, the training
Hamiltonian is dominated by the target Hamiltonian, and the ground state is
close to the target ground state. This is analogous to the low-temperature
limit in standard parallel tempering. The parameter $\Gamma_{r}$ thus allows
us to interpolate between the target parameter landscape at $\Gamma_{r}=1$ and
a much simpler parameter landscape at $\Gamma_{r}=0$. For simplicity we choose
our distribution of $\Gamma_{r}$ to be linearly spaced,

$\Gamma_{r}=1-\frac{r}{N-1}\ ,\quad r=0,\dots,N-1\ .$ (54)

We then propose the following PT update. Since $\Gamma_{r}$ controls the
relative energy scales of $H_{\mathrm{T}}$ and $H_{\mathrm{M}}$, it is not a
suitable candidate for the temperature as it was in the main text, so we
propose a different choice. We expect the standard deviation $\sigma_{r}$ to
be 0 for $\Gamma_{r}=1~{}(r=0)$ because we expect the ANN to reach a good
approximation of an eigenstate of $H_{\mathrm{T}}$, assuming the number ANN
parameters is sufficiently large. Similarly, for larger $r$ replicas this
quantity should be non-zero since we expect the ANN to be a good approximation
of an eigenstate of $H(\Gamma_{r})$. Therefore, we can think of $\sigma_{r}$
as being a measure of the broadening of the eigenstates of $H_{\mathrm{T}}$
due to the “noise” introduced by $\Gamma_{r}<1$. In this sense, we can think
of $\sigma_{r}$ as being a kind of temperature. We thus propose the following
probability rule for swap updates between neighboring pairs $r$ and $r+1$:

$p_{r\leftrightarrow r+1}=\left\\{\begin{array}[]{lr}1&\text{if }r=0\text{ and
}\overline{E_{r+1}}<\overline{E_{r}}\\\ 0&\text{if }r=0\text{ and
}\overline{E_{r+1}}>\overline{E_{r}}\\\
\min\left(1,\exp\left[\left(\overline{E_{r}}-\overline{E_{r+1}}\right)\left(\frac{2}{\overline{\sigma_{r}}}-\frac{2}{\overline{\sigma_{r+1}}}\right)\right]\right)&\text{otherwise}\
.\end{array}\right.$ (55)

While it is not strictly necessary, the $r=0$ replica only swaps its
configuration with its neighbor if its neighbor has a better approximation of
the ground state. This is equivalent to what happens in the main text at zero
temperature. The factor of 2 in $2/\overline{\sigma_{r}}$ is a choice we make
to ensure that enough replica swaps happen in our simulations. This is a
hyper-parameter that could be optimized.

We show the QPT simulation results using this rule for the Precipice problem
in Fig. 7 and for the H4 Rectangle in Fig. 8. The results are comparable to
our results from the main text (Figs. 4, 5, and 6), although the performance
using the rule in the main text is better. We also emphasize that our
simulations with this rule did not succeed at finding a good approximation of
the ground state with $N=5$ replicas, while the approach in the main text
does.

Figure 7: Optimum time to epsilon to reach a relative error of
$\epsilon=10^{-1}$ for the Precipice problem as a function of system size $n$
using the QPT swap probability rule in Eq. (55) with different number of QPT
replicas $N$ with $N_{\max}=40$. The error bars correspond to the 95%
confidence interval calculated using a bootstrap with $10^{3}$ resamplings of
the data.

(a)  (b)

Figure 8: Optimum time to epsilon to reach a relative error of
$\epsilon=10^{-3}$ for the H4 Rectangle at different angles with (a) fixed
Hamming weight exact sampling and (b) full exact sampling using the QPT swap
probability rule in Eq. (55). We use $N_{\max}=20$. The error bars correspond
to the 95% confidence interval calculated using a bootstrap with $10^{3}$
resamplings of the data. The lines are to guide the eye.

###  D.2 Almost constant swap probability

In what follows we use the same the same definition of the training
Hamiltonian as in the main text (Eq. (7)), as well as the same distribution of
$\Gamma_{r}$ values (Eq. (8)). We consider a swap rule that abandons the idea
of a temperature altogether; instead we consider the following swap
probability rule:

$p_{r\leftrightarrow r+1}=\left\\{\begin{array}[]{lr}1&\text{if }r=0\text{ and
}\overline{E_{r+1}}<\overline{E_{r}}\\\ 0&\text{if }r=0\text{ and
}\overline{E_{r+1}}>\overline{E_{r}}\\\
1/4&\text{otherwise}\end{array}\right.$ (56)

The motivation for this choice is as follows. The $r=0$ replica only swaps its
configuration with its neighbor if its neighbor has a better approximation of
the ground state. This is equivalent to what happens in the main text at zero
temperature. For the remaining replicas, we swap the replicas with probability
$1/4$ irrespective of the energies of the configurations. We show the QPT
simulation results using this rule for the Precipice problem in Fig. 9 and for
the H4 Rectangle in Fig. 10. The results are comparable to our results from
the main text (Figs. 4, 5, and 6), although the performance using the rule in
the main text is better in all cases, although there are some qualitative
differences. For the Precipice problem, we see a clear separation between the
different $N$ values for intermediate $n$ values, which suggests that the
constant swap probability with a smaller separation between replicas is
hindering performance compared to the probability rule in the main text. For
the H4 Rectangle, the larger error bars suggest less consistency in the
performance using this scheme.

Figure 9: Optimum time to epsilon to reach a relative error of
$\epsilon=10^{-1}$ for the Precipice problem as a function of system size $n$
using the QPT swap probability rule in Eq. (56) with different number of QPT
replicas $N$ with $N_{\max}=40$. The error bars correspond to the 95%
confidence interval calculated using a bootstrap with $10^{3}$ resamplings of
the data.

(a)  (b)

Figure 10: Optimum time to epsilon to reach a relative error of
$\epsilon=10^{-3}$ for the H4 Rectangle at different angles with (a) fixed
Hamming weight exact sampling and (b) full exact sampling using the QPT swap
probability rule in Eq. (56). We use $N_{\max}=20$. The error bars correspond
to the 95% confidence interval calculated using a bootstrap with $10^{3}$
resamplings of the data. The lines are to guide the eye.

##  Appendix E Short Note On the Bootstrap Procedure

Let us assume we have repeated our stochastic simulations $n$ times and gotten
the set of outcomes $\left\\{x_{1},x_{2},\dots,x_{n}\right\\}$. The outcome of
the $i$-th experiment can be treated as a random variable $\mathbf{X_{i}}$,
with the $i$-th experimental outcome given $\mathbf{X_{i}}=x_{i}$. We will
assume that all the random variables are independent and identically
distributed with mean $\mu_{x}$ and standard deviation $\sigma_{x}$. Our goal
is to report the mean and uncertainty of some property of the random
variables. For example, if we are interested in estimating $\mu_{x}$, we may
use the the sample mean $\mathbf{\bar{X}}$ as our function:

$\mathbf{\bar{X}}=\frac{1}{n}\sum_{i=1}^{n}\mathbf{X_{i}}$ (57)

From the Central Limit Theorem, we expect that
$\mathbf{\bar{X}}\sim\mathcal{N}(\mu_{x},\sigma_{x}^{2}/n)$ for sufficiently
large $n$. So for sufficiently large $n$, we can report our 95% confidence
estimate of $\mu_{x}$ as $\bar{x}\pm 2\frac{\sigma_{x}}{\sqrt{n}}$ where
$\bar{x}=\frac{1}{n}\sum_{i=1}^{n}x_{i}$. We can estimate $\sigma_{x}^{2}$ by
using the sample variance for example.

However, if instead of estimating $\mu_{x}$, we were interested in the median
of $\mathbf{X_{i}}$, then it is not obvious how to use the above method to
give confidence estimates. The bootstrap method allows us to mitigate this
problem in a simple way.

Bootstrapping requires resampling (with replacement) the data
$(x_{1},\dots,x_{n})$ uniformly to acquire $n_{b}$ bootstrap resamples. Each
bootstrap resample $S_{i}$ picks $n$ of the data outcomes randomly with
replacement, i.e. $S_{i}=\left\\{x_{i_{1}},x_{i_{2}},\dots
x_{i_{n}}\right\\}$. So a bootstrap resample effectively amounts to randomly
picking indices from $1$ to $n$ with replacement. For each $S_{i}$, we
calculate the property we are interested in, like the median, giving a value
$F_{i}$. We can treat the set $\left\\{F_{1},F_{2},\dots,F_{n_{b}}\right\\}$
as $n_{b}$ samples of the random variables $\mathbf{F_{i}}$. We now calculate
the sample mean of $\mathbf{F}$:

$\mathbf{\bar{F}}=\frac{1}{n_{b}}\sum_{i=1}^{n_{b}}\mathbf{F_{i}}$ (58)

Now we can use the Central Limit Theorem (for sufficiently large $n_{b}$) that
$\mathbf{\bar{F}}\sim\mathcal{N}$. We can now report our estimate with $95\%$
confidence interval as:

$\bar{F}\pm 2\sigma_{F}$ (59)

where $\bar{F}$ is the sample mean:

$\bar{F}=\frac{1}{n_{b}}\sum_{i=1}^{n_{b}}{F_{i}}$ (60)

and $\sigma_{F}$ is the standard deviation:

$\sigma_{F}^{2}=\frac{1}{n_{b}-1}\sum_{i=1}^{n_{b}}(F_{i}-\bar{F})^{2}$ (61)

##  Appendix F The H4 rectangle

Figure 11: The arrangement of hydrogen atoms in the H4 rectangle. The
configuration is controlled by the radius $R$ and the angle $\theta$.

The H4 rectangle, illustrated in Fig. 11, consists of four hydrogen atoms
placed in a rectangle, each with a distance to the center of $R$. The
separation between these atoms can be further controlled by the angle between
adjacent hydrogen atoms $\theta$. When $\theta$ is far away from 90∘, the
system resembles two disconnected H2 molecules and the energy can be
approximately calculated as twice the energy of an equivalent H2 molecule. As
the angle approaches 90∘, however, the distance between the different hydrogen
atoms becomes equivalent and the system becomes degenerate, leading to a
frustrated system with a peak in energy at 90∘. By moving from a rectangular
system to a square one, the system can be brought arbitrarily close to a
strongly correlated, two configuration system. This configuration has been
shown to lead to difficulty for many atomic modeling algorithms, notably
including Coupled Cluster methods, which incorrectly predict a cusp and
minimum at 90∘ for certain radii [51]. Throughout this work, following Pfau et
al. and others [54, 55], we have used a radius of 1.738 Å, which places the H4
rectangle in a regime where the atoms are nearing full dissociation. Coupled
Cluster methods have been shown to fail at this radius [51], and single
replica RBM runs often struggle to find the ground state energy as seen in the
main text, providing a useful benchmark for our parallel tempering scheme.

To create inputs for our calculations, we consider the Fermionic problem in
the second quantitized formalism

$H=\sum_{i,j}t_{ij}c_{i}^{\dagger}c_{j}+\sum_{i,j,k,m}u_{ijkm}c_{i}^{\dagger}c_{k}^{\dagger}c_{m}c_{j}\
,$ (62)

where we define the fermionic annihilation and creation operators with the
anticommutation relation $\\{c_{i}^{\dagger},c_{j}\\}=\delta_{i,j}$. We define
$t_{ij}$ and $u_{ijkm}$ as one- and two-body integrals, respectively. We then
map this fermionic Hamiltonian onto a spin Hamiltonian with the form

$H=\sum_{j=1}^{r}h_{j}\sigma_{j}\ ,$ (63)

where $h_{j}$ are coefficients and $\sigma_{j}$ is an N-fold tensor product of
single-qubit Pauli operators $I$, $\sigma^{x}$, $\sigma^{y}$, and
$\sigma^{z}$. We use the canonical Jordan-Wigner transformation [26] to map
Hamiltonians of the form of Eq. 62 to Eq. 63. We use a minimal STO-3G basis
set [56] and the OpenFermion software package [57] to generate and map
Hamiltonians for the H4 problem.

For $\theta=80^{\circ}$, the ground state energy of the 8 qubit Hamiltonian is
$-1.88016686$; for $\theta=85^{\circ}$, the ground state energy is
$-1.87584118$; for $\theta=90^{\circ}$, the ground state energy is
$-1.87420093$.

## References

  * [1] D. Landau and K. Binder,  _A Guide to Monte Carlo Simulations in Statistical Physics_ ,  Cambridge University Press, USA,  ISBN 0521842387 (2005).
  * [2] F. Becca and S. Sorella,  _Quantum Monte Carlo Approaches for Correlated Systems_ ,  Cambridge University Press,  10.1017/9781316417041 (2017).
  * [3] E. Y. Loh, J. E. Gubernatis, R. T. Scalettar, S. R. White, D. J. Scalapino and R. L. Sugar,  _Sign problem in the numerical simulation of many-electron systems_ ,  Phys. Rev. B 41, 9301 (1990),  10.1103/PhysRevB.41.9301.
  * [4] P. J. Reynolds, D. M. Ceperley, B. J. Alder and W. A. Lester,  _Fixed‐node quantum monte carlo for molecules_ ,  The Journal of Chemical Physics 77(11), 5593 (1982),  10.1063/1.443766.
  * [5] G. Ortiz, D. M. Ceperley and R. M. Martin,  _New stochastic method for systems with broken time-reversal symmetry: 2d fermions in a magnetic field_ ,  Phys. Rev. Lett. 71, 2777 (1993),  10.1103/PhysRevLett.71.2777.
  * [6] S. Zhang, J. Carlson and J. E. Gubernatis,  _Constrained path quantum monte carlo method for fermion ground states_ ,  Phys. Rev. Lett. 74, 3652 (1995),  10.1103/PhysRevLett.74.3652.
  * [7] I. Affleck, T. Kennedy, E. H. Lieb and H. Tasaki,  _Valence bond ground states in isotropic quantum antiferromagnets_ ,  Communications in Mathematical Physics 115(3), 477 (1988),  10.1007/BF01218021.
  * [8] S. R. White,  _Density matrix formulation for quantum renormalization groups_ ,  Phys. Rev. Lett. 69, 2863 (1992),  10.1103/PhysRevLett.69.2863.
  * [9] S. R. White,  _Density-matrix algorithms for quantum renormalization groups_ ,  Phys. Rev. B 48, 10345 (1993),  10.1103/PhysRevB.48.10345.
  * [10] M. Fannes, B. Nachtergaele and R. F. Werner,  _Finitely correlated states on quantum spin chains_ ,  Communications in Mathematical Physics 144(3), 443 (1992),  10.1007/BF02099178.
  * [11] U. Schollwöck,  _The density-matrix renormalization group_ ,  Rev. Mod. Phys. 77, 259 (2005),  10.1103/RevModPhys.77.259.
  * [12] U. Schollwöck,  _The density-matrix renormalization group in the age of matrix product states_ ,  Annals of Physics 326(1), 96 (2011),  https://doi.org/10.1016/j.aop.2010.09.012,  January 2011 Special Issue.
  * [13] F. Verstraete and J. I. Cirac,  _Renormalization algorithms for quantum-many body systems in two and higher dimensions_ ,  arXiv e-prints cond-mat/0407066 (2004).
  * [14] R. Orús,  _A practical introduction to tensor networks: Matrix product states and projected entangled pair states_ ,  Annals of Physics 349, 117 (2014),  https://doi.org/10.1016/j.aop.2014.06.013,  1306.2164.
  * [15] G. Carleo and M. Troyer,  _Solving the quantum many-body problem with artificial neural networks_ ,  Science 355(6325), 602 (2017),  10.1126/science.aag2302.
  * [16] K. Choo, T. Neupert and G. Carleo,  _Two-dimensional frustrated ${J}_{1}\text{$-$}{J}_{2}$ model studied with neural network quantum states_,  Phys. Rev. B 100, 125124 (2019),  10.1103/PhysRevB.100.125124.
  * [17] M. Schmitt and M. Heyl,  _Quantum many-body dynamics in two dimensions with artificial neural networks_ ,  Phys. Rev. Lett. 125, 100503 (2020),  10.1103/PhysRevLett.125.100503.
  * [18] C. Roth, A. Szabó and A. MacDonald,  _High-accuracy variational Monte Carlo for frustrated magnets with deep neural networks_ ,  arXiv e-prints arXiv:2211.07749 (2022),  2211.07749.
  * [19] M. Bukov, M. Schmitt and M. Dupont,  _Learning the ground state of a non-stoquastic quantum Hamiltonian in a rugged neural network landscape_ ,  SciPost Phys. 10, 147 (2021),  10.21468/SciPostPhys.10.6.147.
  * [20] R. H. Swendsen and J.-S. Wang,  _Replica monte carlo simulation of spin-glasses_ ,  Phys. Rev. Lett. 57, 2607 (1986),  https://doi.org/10.1103/PhysRevLett.57.2607.
  * [21] C. J. Geyer,  _Parallel tempering_ ,  In E. M. Keramidas and S. M. Kaufman, eds., _Computing Science and Statistics Proceedings of the 23rd Symposium on the Interface_ , p. 156. American Statistical Association, New York (1991).
  * [22] K. Hukushima and K. Nemoto,  _Exchange monte carlo method and application to spin glass simulations_ ,  Journal of the Physical Society of Japan 65(6), 1604 (1996),  https://doi.org/10.1143/JPSJ.65.1604.
  * [23] R. A. Baños, A. Cruz, L. A. Fernandez, J. M. Gil-Narvion, A. Gordillo-Guerrero, M. Guidetti, A. Maiorano, F. Mantovani, E. Marinari, V. Martin-Mayor, J. Monforte-Garcia, A. M. Sudupe _et al._ ,  _Nature of the spin-glass phase at experimental length scales_ ,  Journal of Statistical Mechanics: Theory and Experiment 2010(06), P06026 (2010),  10.1088/1742-5468/2010/06/p06026.
  * [24] E. Farhi, D. Gosset, I. Hen, A. W. Sandvik, P. Shor, A. P. Young and F. Zamponi,  _Performance of the quantum adiabatic algorithm on random instances of two optimization problems on regular hypergraphs_ ,  Phys. Rev. A 86, 052334 (2012),  10.1103/PhysRevA.86.052334.
  * [25] T. Albash, G. Wagenbreth and I. Hen,  _Off-diagonal expansion quantum monte carlo_ ,  Phys. Rev. E 96, 063309 (2017),  10.1103/PhysRevE.96.063309.
  * [26] P. Jordan and E. P. Wigner,  _Über das paulische äquivalenzverbot_ ,  In _The Collected Works of Eugene Paul Wigner_ , pp. 109–129. Springer (1993).
  * [27] K. Choo, A. Mezzacapo and G. Carleo,  _Fermionic neural-network states for ab-initio electronic structure_ ,  Nature Communications 11(1), 2368 (2020),  10.1038/s41467-020-15724-9.
  * [28] S. Sorella,  _Green function monte carlo with stochastic reconfiguration_ ,  Phys. Rev. Lett. 80, 4558 (1998),  10.1103/PhysRevLett.80.4558.
  * [29] S. Sorella,  _Generalized lanczos algorithm for variational quantum monte carlo_ ,  Phys. Rev. B 64, 024512 (2001),  10.1103/PhysRevB.64.024512.
  * [30] E. H. Moore,  _On the reciprocal of the general algebraic matrix_ ,  In _Bulletin of the American Mathematical Society_ , vol. 26, pp. 394–395 (1920).
  * [31] R. Penrose,  _A generalized inverse for matrices_ ,  Mathematical Proceedings of the Cambridge Philosophical Society 51(3), 406 (1955),  10.1017/S0305004100030401.
  * [32] G. Fabiani and J. H. Mentink,  _Investigating ultrafast quantum magnetism with machine learning_ ,  SciPost Phys. 7, 4 (2019),  10.21468/SciPostPhys.7.1.004.
  * [33] D. P. Kingma and J. Ba,  _Adam: A Method for Stochastic Optimization_ ,  arXiv e-prints arXiv:1412.6980 (2014),  1412.6980.
  * [34] S.-C. T. Choi and M. A. Saunders,  _Algorithm 937: Minres-qlp for symmetric and hermitian linear equations and least-squares problems_ ,  ACM Trans. Math. Softw. 40(2) (2014),  10.1145/2527267.
  * [35] C.-Y. Park and M. J. Kastoryano,  _Geometry of learning neural quantum states_ ,  Phys. Rev. Research 2, 023232 (2020),  10.1103/PhysRevResearch.2.023232.
  * [36] S.-i. Amari,  _Natural Gradient Works Efficiently in Learning_ ,  Neural Computation 10(2), 251 (1998),  10.1162/089976698300017746.
  * [37] M. Casula, C. Attaccalite and S. Sorella,  _Correlated geminal wave function for molecules: An efficient resonating valence bond approach_ ,  The Journal of Chemical Physics 121(15), 7110 (2004),  10.1063/1.1794632.
  * [38] C. Roth,  _Iterative Retraining of Quantum Spin Models Using Recurrent Neural Networks_ ,  arXiv e-prints arXiv:2003.06228 (2020),  2003.06228.
  * [39] M. Hibat-Allah, R. G. Melko and J. Carrasquilla,  _Supplementing Recurrent Neural Network Wave Functions with Symmetry and Annealing to Improve Accuracy_ ,  arXiv e-prints arXiv:2207.14314 (2022),  2207.14314.
  * [40] T. Okabe, M. Kawata, Y. Okamoto and M. Mikami,  _Replica-exchange monte carlo method for the isobaric–isothermal ensemble_ ,  Chemical Physics Letters 335(5), 435 (2001),  https://doi.org/10.1016/S0009-2614(01)00055-0.
  * [41] S. Syed, A. Bouchard-Côté, G. Deligiannidis and A. Doucet,  _Non-reversible parallel tempering: a scalable highly parallel mcmc scheme_ ,  Journal of Royal Statistical Society, Series B 84, 321 (2021),  10.1111/rssb.12464.
  * [42] N. Metropolis, A. W. Rosenbluth, M. N. Rosenbluth, A. H. Teller and E. Teller,  _Equation of state calculations by fast computing machines_ ,  The Journal of Chemical Physics 21(6), 1087 (1953),  10.1063/1.1699114.
  * [43] W. K. Hastings,  _Monte Carlo sampling methods using Markov chains and their applications_ ,  Biometrika 57(1), 97 (1970),  10.1093/biomet/57.1.97.
  * [44] T. F. Rønnow, Z. Wang, J. Job, S. Boixo, S. V. Isakov, D. Wecker, J. M. Martinis, D. A. Lidar and M. Troyer,  _Defining and detecting quantum speedup_ ,  Science 345(6195), 420 (2014),  10.1126/science.1252319.
  * [45] J. King, S. Yarkoni, M. M. Nevisi, J. P. Hilton and C. C. McGeoch,  _Benchmarking a quantum annealing processor with the time-to-target metric_ ,  arXiv e-prints arXiv:1508.05087 (2015),  1508.05087.
  * [46] S. Bravyi, D. P. DiVincenzo, R. I. Oliveira and B. M. Terhal,  _The complexity of stoquastic local hamiltonian problems_ ,  Quant. Inf. Comp. 8(5), 0361 (2008).
  * [47] S. Bravyi and B. Terhal,  _Complexity of stoquastic frustration-free hamiltonians_ ,  SIAM Journal on Computing 39(4), 1462 (2009),  10.1137/08072689X.
  * [48] W. van Dam, M. Mosca and U. Vazirani,  _How powerful is adiabatic quantum computation?_ ,  In _Proceedings 42nd IEEE Symposium on Foundations of Computer Science_ , pp. 279–287,  10.1109/SFCS.2001.959902 (2001).
  * [49] S. Muthukrishnan, T. Albash and D. A. Lidar,  _Tunneling and speedup in quantum optimization for permutation-symmetric problems_ ,  Phys. Rev. X 6, 031010 (2016),  10.1103/PhysRevX.6.031010.
  * [50] J. Paldus, P. Piecuch, L. Pylypow and B. Jeziorski,  _Application of hilbert-space coupled-cluster theory to simple (h 2) 2 model systems: Planar models_ ,  Physical Review A 47(4), 2738 (1993),  10.1103/PhysRevA.47.2738.
  * [51] T. Van Voorhis and M. Head-Gordon,  _Benchmark variational coupled cluster doubles results_ ,  The Journal of Chemical Physics 113(20), 8873 (2000),  10.1063/1.1319643.
  * [52] M. Suzuki,  _Relationship between d-dimensional quantal spin systems and (d+ 1)-dimensional Ising systems: Equivalence, critical exponents and systematic approximants of the partition function and spin correlations_ ,  Progress of theoretical physics 56(5), 1454 (1976),  10.1143/PTP.56.1454.
  * [53] H. G. Katzgraber, S. Trebst, D. A. Huse and M. Troyer,  _Feedback-optimized parallel tempering monte carlo_ ,  Journal of Statistical Mechanics: Theory and Experiment 2006(03), P03018 (2006),  10.1088/1742-5468/2006/03/p03018.
  * [54] D. Pfau, J. S. Spencer, A. G. D. G. Matthews and W. M. C. Foulkes,  _Ab initio solution of the many-electron schrödinger equation with deep neural networks_ ,  Phys. Rev. Research 2, 033429 (2020),  10.1103/PhysRevResearch.2.033429.
  * [55] H. G. A. Burton and A. J. W. Thom,  _Holomorphic hartree–fock theory: An inherently multireference approach_ ,  Journal of Chemical Theory and Computation 12(1), 167 (2016),  10.1021/acs.jctc.5b01005.
  * [56] W. J. Hehre, R. F. Stewart and J. A. Pople,  _Self‐consistent molecular‐orbital methods. i. use of gaussian expansions of slater‐type atomic orbitals_ ,  The Journal of Chemical Physics 51(6), 2657 (1969),  10.1063/1.1672392.
  * [57] J. R. McClean, N. C. Rubin, K. J. Sung, I. D. Kivlichan, X. Bonet-Monroig, Y. Cao, C. Dai, E. S. Fried, C. Gidney, B. Gimby, P. Gokhale, T. Häner _et al._ ,  _OpenFermion: the electronic structure package for quantum computers_ ,  Quantum Science and Technology 5(3), 034014 (2020),  10.1088/2058-9565/ab8ebc.