# Purifying Deep Boltzmann Machines for Thermal Quantum States

Yusuke Nomura1<EMAIL_ADDRESS>1RIKEN Center for Emergent Matter
Science, 2-1 Hirosawa, Wako, Saitama 351-0198, Japan  Nobuyuki Yoshioka2,3
<EMAIL_ADDRESS>Franco Nori3,4,5 2Department of Applied Physics,
University of Tokyo, 7-3-1 Hongo, Bunkyo-ku, Tokyo 113-8656, Japan
3Theoretical Quantum Physics Laboratory, RIKEN Cluster for Pioneering Research
(CPR), Wako-shi, Saitama 351-0198, Japan  4RIKEN Center for Quantum Computing
(RQC), Wako-shi, Saitama 351-0198, Japan  5Physics Department, University of
Michigan, Ann Arbor, Michigan 48109-1040, USA

###### Abstract

We develop two cutting-edge approaches to construct deep neural networks
representing the purified finite-temperature states of quantum many-body
systems. Both methods commonly aim to represent the Gibbs state by a highly
expressive neural-network wave function, exemplifying the idea of
purification. The first method is an entirely deterministic approach to
generate deep Boltzmann machines representing the purified Gibbs state
exactly. This strongly assures the remarkable flexibility of the ansatz which
can fully exploit the quantum-to-classical mapping. The second method employs
stochastic sampling to optimize the network parameters such that the imaginary
time evolution is well approximated within the expressibility of neural
networks. Numerical demonstrations for transverse-field Ising models and
Heisenberg models show that our methods are powerful enough to investigate the
finite-temperature properties of strongly correlated quantum many-body
systems, even when the problematic effect of frustration is present.

Introduction.— The thermal behavior of quantum many-body systems is one of the
most fundamental problems in physics. Statistical mechanics states that the
density matrix describing a system in thermal equilibrium, governed by a
Hamiltonian $\mathcal{H}$ at an inverse temperature $\beta$, is given by the
Gibbs state $\rho=e^{-\beta\mathcal{H}}/{\rm Tr}[e^{-\beta\mathcal{H}}]$.
Computing and extracting physical properties from the Gibbs state is a
significant challenge to understand natural phenomena, which in reality all
occur at finite $\beta$.

One of the most celebrated numerical techniques in lattice systems is the
quantum Monte Carlo (QMC) method Suzuki (1976); Suzuki  _et al._ (1977);
Sandvik (1999); Gubernatis  _et al._ (2016), typically based on the path
integral formalism of the partition function. The QMC method yields
numerically exact solutions when the positive definiteness is assured;
otherwise, the infamous negative sign problem arises. Many physically
intriguing system falls into the latter category, and therefore various
efforts have been devoted to overcome this difficulty: tensor-network-based
algorithms Verstraete  _et al._ (2004); Zwolak and Vidal (2004); Feiguin and
White (2005); White (2009) mostly applied to one-dimensional (1D) systems,
dynamical mean-field theory Georges  _et al._ (1996) which becomes exact in
the infinite coordination-number limit, diagrammatic Monte Carlo methods
Prokof’ev and Svistunov (1998); Van Houcke  _et al._ (2010), to name a few
Motta  _et al._ (2020); Irikura and Saito (2020); Iitaka (2020). In another
notable approach Sugiura and Shimizu (2012, 2013); Takai  _et al._ (2016);
Iwaki  _et al._ (2020) using the thermal pure quantum (TPQ) states, one can
extract the ensemble property from a single pure state that represents the
thermal equilibrium. We point out, however, that it remains extremely
challenging to establish a methodology that is both reliable and scalable for
finite-temperature calculations in two-dimensional (2D) systems—the most
exotic and intriguing realm in quantum many-body problems.

Neural networks, initially developed for classical data processing in the
context of machine learning, offer a very strong methodology for quantum
physics Carleo and Troyer (2017); Torlai  _et al._ (2018); Wang and Sandvik
(2018); Carleo  _et al._ (2019a); Melko  _et al._ (2019); Melkani  _et al._
(2020); Ahmed  _et al._ (2020). As was firstly demonstrated by Carleo and
Troyer Carleo and Troyer (2017), neural networks applied as variational wave
functions, commonly dubbed as the neural quantum states, are capable of
simulating ground states Carleo and Troyer (2017); Nomura  _et al._ (2017);
Saito (2017); Cai and Liu (2018), excited states Choo  _et al._ (2018); Nomura
(2020); Yoshioka  _et al._ (2021), and even out-of-equilibrium property Carleo
and Troyer (2017); Schmitt and Heyl (2020); Xie and E (2021); Yoshioka and
Hamazaki (2019); Nagy and Savona (2019); Hartmann and Carleo (2019); Vicentini
_et al._ (2019); Yuan  _et al._ (2020) of strongly correlated systems up to
unprecedently large size. Among the tremendous variety of network structures,
Boltzmann machines with restricted connectivity, known to be universal
approximators for arbitrary real/complex amplitudes Le Roux and Bengio (2008);
Montufar and Ay (2011), are useful for statistical mechanics and quantum
information. The Boltzmann machines with the shallowest structure are already
powerful enough to compactly express complex quantum states with extensively-
growing quantum entanglement Deng  _et al._ (2017); Levine  _et al._ (2019).
Furthermore, deep Boltzmann machines (DBMs), i.e., the ones with multiple
hidden layers, are guaranteed to provide an efficient description for an even
wider range of quantum states Gao and Duan (2017). Strongly motivated by their
extremely high representability, the ground states in quantum many-body spin
systems have been successfully simulated by DBMs Carleo  _et al._ (2018).

In the present study, we provide two state-of-the-art methods to construct
DBMs that capture the finite-temperature behavior of quantum many-body
systems. Both methods share the strategy of employing DBMs to express the
purified Gibbs state. Namely, a mixed state under imaginary time evolution is
compactly encoded as a pure DBM wave function in the enlarged Hilbert space.
In the first method, we find a completely deterministic way to construct DBMs,
realizing the exact purified expression of finite temperature states. This
proves the remarkable flexibility and power of DBMs for investigating finite-
temperature many-body phenomena. In the second method, we provide a stochastic
way to simulate the imaginary time evolution which exploits the versatile
expressive power of DBMs as approximators. For demonstration, we apply these
methods to the 1D transverse-field Ising (TFI) model and the 2D
$J_{1}$–$J_{2}$ Heisenberg model on the square lattice to find surprisingly
high accuracy compared to numerically exact methods. We emphasize that only a
polynomial number of auxiliary spins suffices in both approaches, yielding a
huge computational advantage even under the presence of the problematic effect
of frustration.

Figure 1:  (a) Structure of the three-layer DBM used in the present study. The
visible layer (blue) corresponds to physical spins $\sigma_{i}$. Two hidden
layers are distinguished as hidden (green) and deep (red) layers; while the
hidden layer simply consists of hidden spins $h_{j}$, the deep layer is
composed of the deep spins $d_{k}$ and ancilla spins $\sigma_{i}^{\prime}$,
which are introduced to purify the Gibbs-state density matrix. Here, the
number of physical/ancilla, hidden, and deep spins are denoted as $N$, $M$,
$M^{\prime}$, respectively. (b) Two different ways of depicting the identical
DBM structure. A network with seemingly many hidden layers can always be
recast into that with only two hidden layers. (c) The DBM representing the
infinite-temperature state $\ket{\Psi(T=\infty)}$ for quantum spin-1/2
systems. (d) The DBM construction for the finite-temperature states of the 1D
TFI model [$|\Psi(T)\rangle$ in Eq. (3)]. The arrow denotes the growth of the
DBM structure along the imaginary-time $\tau$ propagation. Light blue, orange,
and purple bonds in (c) and (d) have couplings $\mathrm{i}\frac{\pi}{4}$,
$\frac{1}{2}{\rm arcosh}\left(\frac{1}{{\rm
tanh}(\Gamma\delta_{\tau})}\right)$, and $\frac{1}{2}{\rm
arcosh}(e^{2J\delta_{\tau}})$, respectively.

DBM for purification.— Throughout this work, we utilize the idea of
purification to represent the Gibbs state. Namely, the finite-temperature
density matrix $\rho$ of a target system $\mathcal{S}$ is encoded as a pure
state in an extended system $\mathcal{S}+\mathcal{A}$, such that $\rho$ is
obtained by tracing out the ancillary system $\mathcal{A}$ Nielsen and Chuang
(Cambridge University Press, (2000). For instance, the purification of an
infinite-temperature state can be generated by the superposition
$\sum_{x}\ket{x}_{\mathcal{S}}\ket{a_{x}}_{\mathcal{A}}$, where
$\\{\ket{x}_{\mathcal{S}}\\}$ is the complete orthonormal basis set of the
target system, while $\\{\ket{a_{x}}_{\mathcal{A}}\\}$ is an orthonormal but
not necessarily complete basis set of ancillary system.

For the sake of concreteness, let us consider a quantum many-body spin-1/2
system. We introduce $N_{\rm site}$ binary degrees of freedom
$\\{\sigma_{i}\\}$ so that $\ket{\sigma}=\ket{\sigma_{1},...,\sigma_{N_{\rm
site}}}$ spans the Hilbert space of the target system $\mathcal{S}$. Hereafter
we call them physical spins. As the ancillary system $\mathcal{A}$, we
introduce an identical number of ancilla spins $\\{\sigma_{i}^{\prime}\\}$ so
that an arbitrary mixed state can be purified in principle. While the
purification of a mixed state is not unique, here we exclusively take the
purified infinite-temperature state as
$|\Psi(T\\!=\\!\infty)\rangle=\bigotimes_{i=1}^{N_{\rm
site}}\left(|\uparrow\downarrow^{\prime}\rangle+|\downarrow\uparrow^{\prime}\rangle\right)_{i}$,
and perform the imaginary time evolution as
$\ket{\Psi(T)}=e^{-\beta\mathcal{H}/2}\otimes\mathbb{1}^{\prime}\ket{\Psi(T=\infty)}$,
with $\beta=1/T$ to simulate the finite-temperature Gibbs state. Note that the
infinite-temperature state $\rho_{\infty}=\mathbb{1}/2^{N_{\rm site}}$ and the
finite-temperature state $\rho_{T}=e^{-\beta\mathcal{H}}/{\rm
Tr}[e^{-\beta\mathcal{H}}]$ are reproduced by tracing out the ancilla spins.

Intriguingly, the purified Gibbs state at each temperature can be efficiently
expressed by the DBM. In particular, we use the DBM with two hidden layers
[see Fig. 1(a)] to represent a purified wave function whose amplitude
$\Psi(\sigma,\sigma^{\prime})=\braket{\sigma,\sigma^{\prime}}{\Psi}$ is
parametrized as

$\displaystyle\Psi(\sigma,\sigma^{\prime})$ $\displaystyle=$
$\displaystyle\sum_{h,d}\phi(\sigma,\sigma^{\prime};h,d),$ (1)
$\displaystyle\phi(\sigma,\sigma^{\prime};h,d)$ $\displaystyle=$
$\displaystyle\exp\biggl{[}\sum_{j}b_{j}h_{j}\\!+\\!\sum_{ji}\\!h_{j}(W_{ji}\sigma_{i}\\!+\\!W_{ji}^{\prime}\sigma^{\prime}_{i})\\!$
(2) $\displaystyle+\\!\sum_{jk}W_{jk}^{\prime}h_{j}d_{k}\biggr{]},$

where we have introduced hidden spins $\\{h_{j}\\}$ and deep spins
$\\{d_{k}\\}$, in addition to the physical spins $\\{\sigma_{i}\\}$
constituting the visible layer and ancilla spins $\\{\sigma_{i}^{\prime}\\}$
allocated in the second hidden layer (we define all spins to be either $+1$ or
$-1$). This structure is “universal” in terms of the DBM architecture; the
arbitrary multi-hidden layer structure can be rearranged to have only two
hidden layers as shown in Fig 1(b). The number of complex variational
parameters $\mathcal{W}=\\{b,W,W^{\prime}\\}$ 111The renormalization factor is
omitted in the expression. We ignore bias terms for visible and deep units
because they are irrelevant in the present case. are directly related to the
number of $h$ and $d$ spins, which controls the representability of the DBM
wave function. Note that the purification technique for neural networks has
been considered in the context of quantum tomography Torlai and Melko (2018)
and dissipative quantum physics Yoshioka and Hamazaki (2019); Nagy and Savona
(2019); Hartmann and Carleo (2019); Vicentini  _et al._ (2019).

Method (I): Analytic purification using DBM.— In Method (I), we analytically
construct DBMs that exactly reproduce the behaviours of Gibbs states. First,
to represent the infinite-temperature state, we introduce a DBM with $N_{\rm
site}$ hidden and ancilla spins [Fig. 1(c)]. By setting the parameters
$\mathcal{W}_{\infty}$ as $b_{j}=0$ and
$W_{ji}=W^{\prime}_{ji}=\mathrm{i}\frac{\pi}{4}\delta_{ji}$, we find that the
analytical expression of the DBM wave function becomes
$\Psi(\sigma,\sigma^{\prime})=\prod_{i}2\cosh\left[\mathrm{i}\frac{\pi}{4}(\sigma_{i}+\sigma^{\prime}_{i})\right]$.
This exactly reproduces the $|\Psi(T\\!=\\!\infty)\rangle$ described above.

Next, to express finite-temperature states $|\Psi(T)\rangle$ analytically, we
introduce the Suzuki-Trotter decomposition Suzuki (1976). Namely, given the
Hamiltonian $\mathcal{H}=\sum_{\nu}\mathcal{H}_{\nu}$, the purified Gibbs
state up to the Trotter error is expressed as

$\displaystyle|\Psi(T)\rangle:=\left(\left[\prod_{\nu}e^{-\delta_{\tau}\mathcal{H}_{\nu}}\right]^{N_{\tau}}\otimes\mathbb{1}^{\prime}\right)\ket{\Psi(T\\!=\\!\infty)},$
(3)

where $\delta_{\tau}=\beta/2N_{\tau}$ is the propagation time step.
Remarkably, thanks to the flexible representability of the DBM, we can find
analytical solutions for Eq. (3): Starting from the initial DBM state
$|\Psi(T\\!=\\!\infty)\rangle$ given by the $\mathcal{W}_{\infty}$ parameters,
each short-time propagation $e^{-\delta_{\tau}\mathcal{H}_{\nu}}$ can be
implemented exactly by modifying the parameters and structure of the DBM
222For spin-1/2 systems, imaginary-time propagators can be analytically
represented by the DBM. This has been discussed in the context of classical
many-body interacting systems as well yoshioka_transforming_2019., whose the
explicit expression depends on the form of the Hamiltonian.

As a concrete example, let us consider the TFI model on the $N_{\rm
site}$-spin chain under periodic boundary condition. The Hamiltonian is given
by ${\mathcal{H}}={\mathcal{H}_{1}}+{\mathcal{H}_{2}}$, with
${\mathcal{H}_{1}}=-J\sum_{i}\sigma^{z}_{i}\sigma^{z}_{i+1}$ and
${\mathcal{H}_{2}}=-\Gamma\sum_{i}\sigma^{x}_{i}$, where
$\sigma_{i}^{a}(a=x,y,z)$ denotes the Pauli matrix operating on the $i$-th
site. We take the Ising-type interaction $J$ as the energy unit ($J=1$), and
$\Gamma$ as the strength of the transverse magnetic field. Using Method (I),
we can analytically construct the finite-temperature state using the DBM.
Solutions to realize the propagation by updating the DBM parameters from
${\mathcal{W}}$ to $\bar{\mathcal{W}}$, i.e.,
$e^{-\delta_{\tau}\mathcal{H}_{\nu}}|\Psi_{\mathcal{W}}\rangle=C|\Psi_{\bar{\mathcal{W}}}\rangle$
($C$ : a constant) can be sketched as follows Carleo  _et al._ (2018): For
$\nu=1$ (interaction propagator), we add a single hidden spin for each
neighboring visible spins $\sigma_{i}$ and $\sigma_{i+1}$. For $\nu=2$
(transverse-field propagator), we add new hidden- and deep-spin layers between
the visible and neighboring hidden layers. All the DBM parameters are
determined analytically [see Supplementary Materials (SM) for detailed
implementation 333See Supplementary Materials for more details of the method
(URL to be added).]. As a consequence of propagarions, the DBM architecture
grows as in Fig. 1(d), and the number of hidden and deep spins scale as
$\mathcal{O}(N_{\tau}N_{\rm site})$,  which is also common among general local
Hamiltonians.

Now that the Gibbs states are represented, let us discuss how the physical
quantities are computed using the DBM framework in general.  Since the
expectation value of a physical observable
$\langle{\mathcal{O}}\rangle=\frac{\langle\Psi(T)|{\mathcal{O}}\otimes{\mathbb{1}}^{\prime}|\Psi(T)\rangle}{\langle\Psi(T)|\Psi(T)\rangle}$
becomes analytically intractable, we use the Monte Carlo (MC) method for its
numerical estimation. The sampling weight is based on the expression of the
normalization factor of the DBM state [Eq. (1)] given as

$\displaystyle\langle\Psi|\Psi\rangle=\sum_{\sigma,\sigma^{\prime}}\bigl{|}\Psi(\sigma,\sigma^{\prime})\bigr{|}^{2}$
$\displaystyle\ \ \
=\sum_{\sigma,\sigma^{\prime}}\sum_{h_{1},d_{1},h_{2},d_{2}}\phi^{\ast}(\sigma,\sigma^{\prime};h_{1},d_{1})\phi(\sigma,\sigma^{\prime};h_{2},d_{2}).$
(4)

Namely, we sample over the configurations of
$(\sigma,\sigma^{\prime},h_{1},h_{2},d_{1},d_{2})$ weighted by the product of
amplitudes as
$w(\sigma,\sigma^{\prime};h,d)\equiv\phi^{*}(\sigma,\sigma^{\prime};h_{1},d_{1})~{}\phi(\sigma,\sigma^{\prime};h_{2},d_{2})$.
Alternatively, it is possible to trace out the hidden spins $h$ analytically,
and use
$\sum_{h_{1},h_{2}}\phi^{*}(\sigma,\sigma^{\prime};h_{1},d_{1})\phi(\sigma,\sigma^{\prime};h_{2},d_{2})$
as the MC sampling weight over configurations of
$(\sigma,\sigma^{\prime},d_{1},d_{2})$ (or trace out $d$ spins and sample over
$h$ spins). See SM for more details of the method Note (3).

To verify the construction of the DBM and the proposed MC sampling strategy,
we apply the method to the 16-site TFI model. Figure 2 shows the DBM results
(symbols) for the temperature dependence of the (a) energy, (b) specific heat,
and (c) susceptibility. As expected, the DBM results follow the exact
temperature evolution (solid curves). This confirms the remarkable
representability of the DBM not only at zero temperature Carleo  _et al._
(2018) but also at finite temperatures, offering an intriguing quantum-to-
classical mapping.

Figure 2:  Finite-temperature calculations for the 1D TFI model on a 16-site
chain ($N_{\rm site}=16$) with periodic boundary condition: (a) Energy $E$,
(b) specific heat $C$, and (c) susceptibility
$\chi=\frac{1}{T}\sum_{i}\langle\sigma_{0}\sigma_{i}\rangle$. The symbols
denote the DBM results [Method (I) with $\delta_{\tau}=0.05$], which agree
well with the exact-diagonalization results (solid curves).

Method (II): Numerical purification using DBM.— When the weight of individual
spin configuration can be taken to be always positive, Method (I) is quite
useful and provides numerically exact finite-temperature results. However,
when the frustration exists in the spin Hamiltonian, for instance, we cannot
avoid the existence of negative weights as in other finite-temperature
calculations based on the QMC method. While it is possible to construct
finite-temperature states analytically, the estimation of physical quantities
becomes extremely difficult because of the negative sign problem.

To make the application to frustrated models possible, we propose an
alternative method which employs DBMs with only ancilla spins
$\sigma^{\prime}$ in the second hidden layer. By tracing out hidden spins $h$,
such a purified DBM wave function has a simple form:

$\displaystyle\Psi(\sigma,\sigma^{\prime})=\prod_{j}2\cosh\Bigl{[}b_{j}\\!+\\!\sum_{i}(W_{ji}\sigma_{i}\\!+\\!W^{\prime}_{ji}\sigma_{i}^{\prime})\Bigr{]}.$
(5)

Then, we can avoid negative signs by simply employing
$\bigl{|}\Psi(\sigma,\sigma^{\prime})\bigr{|}^{2}$ as the weight for sampling
over $\sigma$ and $\sigma^{\prime}$ spins. However, in this case, differently
from Method (I), the imaginary-time evolution cannot be followed analytically.
Instead, we need to update parameters numerically at each time so that the DBM
obeys the imaginary-time evolution starting from the infinite-temperature
state
$\Psi_{\infty}(\sigma,\sigma^{\prime})=\prod_{i}2\cosh\left[\mathrm{i}\frac{\pi}{4}(\sigma_{i}+\sigma^{\prime}_{i})\right]$
(recall that this infinite-temperature DBM does not require $d$ spins). For
this purpose, we employ the stochastic reconfiguration (SR) method Sorella
(1998, 2001); Note (3) 444The same optimization scheme has been developed
independently in artificial-intelligence community, called “natural gradient”
amari_1992; amari_1998.. The SR optimization is designed to minimize the
distance between quantum states following the exact and variational imaginary-
time evolution, as much as possible within the expressive power of the DBM
wave function in Eq. (5). The expressive power is systematically controlled by
the number of $h$ spins; It is ensured that any quantum states can be
represented exactly (universal approximation) by an infinitely wide network
structure Le Roux and Bengio (2008); Montufar and Ay (2011).

Figure 3:  Finite-temperature calculations for the 2D $J_{1}$-$J_{2}$
Heisenberg model ($J_{1}=1$) on the $6\times 6$ square lattice with periodic
boundary condition: (a) Energy $E$, (b) specific heat $C$, and (c) $z$
component of the spin structure factor $S^{z}({\bf q})=\frac{1}{N_{\rm
site}}\sum_{ij}e^{i{\bf q}\cdot({\bf R}_{i}-{\bf R}_{j})}\langle
S^{z}_{i}S^{z}_{j}\rangle$ at ${\bf q}=(\pi,\pi)$. The symbols denote Method
(II) results with $8N_{\rm site}$ hidden spins, which show a good agreement
with the numerically exact references (solid curves) obtained by Method (I)
with $\delta_{\tau}=0.005$ ($J_{2}=0$) and the TPQ method ($J_{2}=0.5$). The
TPQ calculations are performed with ${\mathcal{H}}\Phi$ Kawamura  _et al._
(2017). The shaded regions show the size of the error bars of the TPQ results
[in Method (I), the size of the error bars is small].

As a demonstration of Method (II), we turn to a highly challenging problem:
frustrated spin systems. As a representative, here, we consider the 2D
antiferromagnetic $J_{1}$–$J_{2}$ Heisenberg model on the $L\times L$ square
lattices with periodic boundary condition. The Hamiltonian reads
${\mathcal{H}}=J_{1}\sum_{\langle i,j\rangle}{\bf S}_{i}\cdot{\bf
S}_{j}+J_{2}\sum_{\langle\langle i,j\rangle\rangle}{\bf S}_{i}\cdot{\bf
S}_{j}$. Here, ${\bf S}_{i}$ is the spin-$1/2$ operator at site $i$, and
$J_{1}(=1)$ and $J_{2}$ are the nearest-neighbor and next-nearest-neighbor
couplings, respectively. When $J_{2}$ is finite, the spin configuration cannot
satisfy the energy gain by the $J_{1}$ and $J_{2}$ interactions simultaneously
(frustration). Around $J_{2}=0.5$, where the frustration is strong, an exotic
state of matter, quantum spin liquid without any symmetry breaking, might be
stabilized as the ground state Jiang  _et al._ (2012); Hu  _et al._ (2013);
Morita  _et al._ (2015); Wang and Sandvik (2018); Ferrari and Becca (2020);
Liu  _et al._ (2020); Nomura and Imada (2020). The model also attracts
attention because of its possible relevance to the physics of high-$T_{\rm c}$
cuprates Nori and Zimanyi (1991); Nori  _et al._ (1992, 1995). However,
because numerically exact QMC results are not available due to the sign
problem, the ground-state phase diagram is still under active debate.

For this challenging problem, the wave functions using neural networks have
started to be applied to the zero temperature calculation Liang  _et al._
(2018); Choo  _et al._ (2019); Ferrari  _et al._ (2019); Westerhout  _et al._
(2020); Szabó and Castelnovo (2020); Nomura and Imada (2020); Nomura (2021).
However, to detect a hallmark of the possible quantum spin liquid phase
experimentally, the finite-temperature behavior needs to be elucidated. Here,
we apply Method (II) to perform the finite-temperature calculations for
$J_{2}=0$ and 0.5. To check its accuracy, we compare the results with
numerically exact ones obtained by Method (I) for the non-frustrated case
($J_{2}=0$). For the frustrated case ($J_{2}=0.5$), Method (I) suffers from
the sign problem, but, up to a $6\times 6$ lattice, the TPQ results are
available, which are also numerically exact. Therefore, we perform
calculations using the $6\times 6$ lattice with the total magnetization
restricted to be zero $\bigl{(}\sum_{i}S^{z}_{i}=0\bigr{)}$. To further
improve the accuracy of the calculation, we utilize the translational and
point-group symmetry of the extended system Nomura (2021). See SM for the
practical details Note (3).

Figure 3 shows Method (II) results (symbols) for the temperature dependence of
the (a) energy, (b) specific heat, and (c) $z$ component of spin structure
factor $S^{z}(\pi,\pi)$, which quantifies the Néel-type antiferromagnetic
correlation. We can see that, by the frustration, the antiferromagnetic
correlation is largely suppressed, and the entropy release slows down. Method
(II) results accurately reproduce the exact imaginary-time evolution, showing
its reliability even in the frustrated regime.

In Method (II), by optimizing parameters numerically, we obtain a more compact
and dense network to represent finite-temperature states compared to the
analytically derived network in Method (I): the number of hidden units is
$N_{h}=8N_{\rm site}$ in this case, which is in contrast to
${\mathcal{O}}(N_{\tau}N_{\rm site})$ in Method (I). The compactness of the
network without $d$ spins results in the absence of the negative weights in
the MC sampling. Considering the observation that the $N_{h}$ scales
polynomially with respect to the system size, our results are strongly
encouraging to expect a computational advantage in even more dedicated
simulations for larger systems.

From the application to the Heisenberg model, we see that both the
implementation of symmetrization and increasing $N_{h}$ contributes to
enhancing the accuracy of calculations Note (3). It is an important future
task to check the reliability of Method (II) for other Hamiltonians, since the
convergence of the error in the imaginary-time evolution with respect to
$N_{h}$ should be model dependent.

Summary and outlook.— In summary, we have proposed two cutting-edge approaches
that utilize DBMs to simulate the finite-temperature properties of quantum
many-body systems. In the first approach, we provide a deterministic
construction of DBMs that exactly represents Gibbs states, which proves the
suitability and flexibility of neural networks for encoding thermal
properties. In the second approach, the DBM network parameters are optimized
stochastically so that the imaginary time evolution can be approximated
efficiently, even for one of the most challenging 2D problems, such as
frustrated systems.

Several future directions can be envisioned. It is an interesting question how
the neural-network quantum states perform under other schemes of finite-
temperature calculation such as TPQ methods. All variational ansatze are, by
construction, not powerful enough to express Haar random states, which are
taken as the initial states in TPQ calculations. Nonetheless, results by
tensor-network-based algorithms White (2009) imply that using the truncated
Hilbert space is sufficient in practical simulations. The major obstacle is
considered to be the entanglement growth along the time propagation, which we
expect to be simulated well by neural networks, based on previous works on
real-time evolution Carleo and Troyer (2017); Schmitt and Heyl (2020); Xie and
E (2021). Also, it is natural to explore the scalability of our methods in
larger and/or more complex systems, or ask whether other network structures
(e.g. deep feed-forward networks) are suited for the finite-temperature
calculations; the trade-off relationship between the representability and
trainability of shallow/deep neural networks remains to be open.

##  I Acknowledgements

Y.N. is supported by Grant-in-Aids for Scientific Research (JSPS KAKENHI)
(Grants No. 16H06345, No. 17K14336, No. 18H01158, and No. 20K14423) and MEXT
as “Program for Promoting Researches on the Supercomputer Fugaku” (Basic
Science for Emergence and Functionality in Quantum Matter —Innovative
Strongly-Correlated Electron Science by Integration of “Fugaku” and Frontier
Experiments—) (Project ID: hp200132). N.Y. is supported by the Japan Science
and Technology Agency (JST) (via the Q-LEAP program). F.N. is supported in
part by: NTT Research, Army Research Office (ARO) (Grant No.
W911NF-18-1-0358), Japan Science and Technology Agency (JST) (via the CREST
Grant No. JPMJCR1676), JSPS (via the KAKENHI Grant No. JP20H00134 and the
JSPS-RFBR Grant No. JPJSBP120194828), the Asian Office of Aerospace Research
and Development (AOARD) (via Grant No. FA2386-20-1-4069), and the Foundational
Questions Institute Fund (FQXi) via Grant No. FQXi-IAF19-06. Part of numerical
calculations were performed using NetKet Carleo  _et al._ (2019b). Some
calculations were performed using the computational resources of
supercomputers in RIKEN (Fugaku and HOKUSAI) and the Institute of Solid State
Physics at the University of Tokyo.

## References

  * Suzuki (1976) Masuo 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, 1454–1469 (1976).
  * Suzuki  _et al._ (1977) Masuo Suzuki, Seiji Miyashita, and Akira Kuroda, “Monte Carlo Simulation of Quantum Spin Systems. I,” Progress of Theoretical Physics 58, 1377–1387 (1977).
  * Sandvik (1999) Anders W. Sandvik, “Stochastic series expansion method with operator-loop update,” Phys. Rev. B 59, R14157–R14160 (1999).
  * Gubernatis  _et al._ (2016) James Gubernatis, Naoki Kawashima, and Philipp Werner,  _Quantum Monte Carlo Methods_ (Cambridge University Press, 2016).
  * Verstraete  _et al._ (2004) F. Verstraete, J. J. García-Ripoll, and J. I. Cirac, “Matrix product density operators: Simulation of finite-temperature and dissipative systems,” Phys. Rev. Lett. 93, 207204 (2004).
  * Zwolak and Vidal (2004) Michael Zwolak and Guifré Vidal, “Mixed-state dynamics in one-dimensional quantum lattice systems: A time-dependent superoperator renormalization algorithm,” Phys. Rev. Lett. 93, 207205 (2004).
  * Feiguin and White (2005) Adrian E. Feiguin and Steven R. White, “Finite-temperature density matrix renormalization using an enlarged Hilbert space,” Physical Review B 72, 220401 (2005).
  * White (2009) Steven R. White, “Minimally entangled typical quantum states at finite temperature,” Phys. Rev. Lett. 102, 190601 (2009).
  * Georges  _et al._ (1996) Antoine Georges, Gabriel Kotliar, Werner Krauth, and Marcelo J. Rozenberg, “Dynamical mean-field theory of strongly correlated fermion systems and the limit of infinite dimensions,” Rev. Mod. Phys. 68, 13–125 (1996).
  * Prokof’ev and Svistunov (1998) Nikolai V. Prokof’ev and Boris V. Svistunov, “Polaron Problem by Diagrammatic Quantum Monte Carlo,” Phys. Rev. Lett. 81, 2514â2517 (1998).
  * Van Houcke  _et al._ (2010) Kris Van Houcke, Evgeny Kozik, N. Prokof’ev, and B. Svistunov, “Diagrammatic Monte Carlo,” Physics Procedia 6, 95–105 (2010), computer Simulations Studies in Condensed Matter Physics XXI.
  * Motta  _et al._ (2020) Mario Motta, Chong Sun, Adrian T. K. Tan, Matthew J. O’Rourke, Erika Ye, Austin J. Minnich, Fernando G. S. L. Brandão, and Garnet Kin-Lic Chan, “Determining eigenstates and thermal states on a quantum computer using quantum imaginary time evolution,” Nature Physics 16, 205–210 (2020).
  * Irikura and Saito (2020) Naoki Irikura and Hiroki Saito, “Neural-network quantum states at finite temperature,” Phys. Rev. Research 2, 013284 (2020).
  * Iitaka (2020) Toshiaki Iitaka, “Random phase product sate for canonical ensemble,” (2020), arXiv:2006.14459 [cond-mat.str-el] .
  * Sugiura and Shimizu (2012) Sho Sugiura and Akira Shimizu, “Thermal pure quantum states at finite temperature,” Phys. Rev. Lett. 108, 240401 (2012).
  * Sugiura and Shimizu (2013) Sho Sugiura and Akira Shimizu, “Canonical thermal pure quantum state,” Phys. Rev. Lett. 111, 010401 (2013).
  * Takai  _et al._ (2016) Kensaku Takai, Kota Ido, Takahiro Misawa, Youhei Yamaji, and Masatoshi Imada, “Finite-Temperature Variational Monte Carlo Method for Strongly Correlated Electron Systems,” Journal of the Physical Society of Japan 85, 034601 (2016).
  * Iwaki  _et al._ (2020) Atsushi Iwaki, Akira Shimizu, and Chisa Hotta, “Thermal pure quantum matrix product states recovering a volume law entanglement,” (2020), arXiv:2005.06829 [cond-mat.str-el] .
  * Carleo and Troyer (2017) Giuseppe Carleo and Matthias Troyer, “Solving the quantum many-body problem with artificial neural networks,” Science 355, 602–606 (2017).
  * Torlai  _et al._ (2018) Giacomo Torlai, Guglielmo Mazzola, Juan Carrasquilla, Matthias Troyer, Roger Melko, and Giuseppe Carleo, “Neural-network quantum state tomography,” Nature Physics 14, 447–450 (2018).
  * Wang and Sandvik (2018) Ling Wang and Anders W. Sandvik, “Critical Level Crossings and Gapless Spin Liquid in the Square-Lattice Spin-$1/2$ $J_{1}$–$J_{2}$ Heisenberg Antiferromagnet,” Phys. Rev. Lett. 121, 107202 (2018).
  * Carleo  _et al._ (2019a) Giuseppe Carleo, Ignacio Cirac, Kyle Cranmer, Laurent Daudet, Maria Schuld, Naftali Tishby, Leslie Vogt-Maranto, and Lenka Zdeborová, “Machine learning and the physical sciences,” Rev. Mod. Phys. 91, 045002 (2019a).
  * Melko  _et al._ (2019) Roger G Melko, Giuseppe Carleo, Juan Carrasquilla, and J Ignacio Cirac, “Restricted Boltzmann machines in quantum physics,” Nature Physics 15, 887–892 (2019).
  * Melkani  _et al._ (2020) Abhijeet Melkani, Clemens Gneiting, and Franco Nori, “Eigenstate extraction with neural-network tomography,” Phys. Rev. A 102, 022412 (2020).
  * Ahmed  _et al._ (2020) Shahnawaz Ahmed, Carlos SÃ¡nchez MuÃ±oz, Franco Nori, and Anton Frisk Kockum, “Quantum state tomography with conditional generative adversarial networks,” (2020), arXiv:2008.03240 [quant-ph] .
  * Nomura  _et al._ (2017) Yusuke Nomura, Andrew S. Darmawan, Youhei Yamaji, and Masatoshi Imada, “Restricted Boltzmann machine learning for solving strongly correlated quantum systems,” Phys. Rev. B 96, 205152 (2017).
  * Saito (2017) Hiroki Saito, “Solving the Bose-Hubbard Model with Machine Learning,” J. Phys. Soc. Jpn. 86, 093001 (2017).
  * Cai and Liu (2018) Zi Cai and Jinguo Liu, “Approximating quantum many-body wave functions using artificial neural networks,” Phys. Rev. B 97, 035116 (2018).
  * Choo  _et al._ (2018) Kenny Choo, Giuseppe Carleo, Nicolas Regnault, and Titus Neupert, “Symmetries and many-body excitations with neural-network quantum states,” Phys. Rev. Lett. 121, 167204 (2018).
  * Nomura (2020) Yusuke Nomura, “Machine learning quantum states — extensions to fermion–boson coupled systems and excited-state calculations,” J. Phys. Soc. Jpn. 89, 054706 (2020).
  * Yoshioka  _et al._ (2021) Nobuyuki Yoshioka, Wataru Mizukami, and Franco Nori, “Solving quasiparticle band spectra of real solids using neural-network quantum states,” Communications Physics 4, 106 (2021).
  * Schmitt and Heyl (2020) Markus Schmitt and Markus Heyl, “Quantum many-body dynamics in two dimensions with artificial neural networks,” Phys. Rev. Lett. 125, 100503 (2020).
  * Xie and E (2021) Pinchen Xie and Weinan E, “Coarse-grained spectral projection: A deep learning assisted approach to quantum unitary dynamics,” Physical Review B 103, 024304 (2021).
  * Yoshioka and Hamazaki (2019) Nobuyuki Yoshioka and Ryusuke Hamazaki, “Constructing neural stationary states for open quantum many-body systems,” Phys. Rev. B 99, 214306 (2019).
  * Nagy and Savona (2019) Alexandra Nagy and Vincenzo Savona, “Variational Quantum Monte Carlo Method with a Neural-Network Ansatz for Open Quantum Systems,” Phys. Rev. Lett. 122, 250501 (2019).
  * Hartmann and Carleo (2019) Michael J. Hartmann and Giuseppe Carleo, “Neural-Network Approach to Dissipative Quantum Many-Body Dynamics,” Phys. Rev. Lett. 122, 250502 (2019).
  * Vicentini  _et al._ (2019) Filippo Vicentini, Alberto Biella, Nicolas Regnault, and Cristiano Ciuti, “Variational Neural-Network Ansatz for Steady States in Open Quantum Systems,” Phys. Rev. Lett. 122, 250503 (2019).
  * Yuan  _et al._ (2020) Dong Yuan, Heran Wang, Zhong Wang, and Dong-Ling Deng, “Solving the Liouvillian Gap with Artificial Neural Networks,” (2020), arXiv:2009.00019 [quant-ph] .
  * Le Roux and Bengio (2008) Nicolas Le Roux and Yoshua Bengio, “Representational Power of Restricted Boltzmann Machines and Deep Belief Networks,” Neural Computation 20, 1631–1649 (2008).
  * Montufar and Ay (2011) Guido Montufar and Nihat Ay, “Refinements of universal approximation results for deep belief networks and restricted boltzmann machines,” Neural Comput. 23, 1306–1319 (2011).
  * Deng  _et al._ (2017) Dong-Ling Deng, Xiaopeng Li, and S. Das Sarma, “Quantum entanglement in neural network states,” Phys. Rev. X 7, 021021 (2017).
  * Levine  _et al._ (2019) Yoav Levine, Or Sharir, Nadav Cohen, and Amnon Shashua, “Quantum entanglement in deep learning architectures,” Phys. Rev. Lett. 122, 065301 (2019).
  * Gao and Duan (2017) Xun Gao and Lu-Ming Duan, “Efficient representation of quantum many-body states with deep neural networks,” Nat. Commun. 8, 662 (2017).
  * Carleo  _et al._ (2018) Giuseppe Carleo, Yusuke Nomura, and Masatoshi Imada, “Constructing exact representations of quantum many-body systems with deep neural networks,” Nature Communications 9, 5322 (2018).
  * Nielsen and Chuang (Cambridge University Press, (2000) Michael A Nielsen and Isaac L Chuang, Quantum computation and quantum information (Cambridge University Press, (2000)).
  * Note (1) The renormalization factor is omitted in the expression. We ignore bias terms for visible and deep units because they are irrelevant in the present case.
  * Torlai and Melko (2018) Giacomo Torlai and Roger G. Melko, “Latent space purification via neural density operators,” Phys. Rev. Lett. 120, 240503 (2018).
  * Note (2) For spin-1/2 systems, imaginary-time propagators can be analytically represented by the DBM. This has been discussed in the context of classical many-body interacting systems as well yoshioka_transforming_2019.
  * Note (3) See Supplementary Materials for more details of the method (URL to be added).
  * Sorella (1998) Sandro Sorella, “Green function monte carlo with stochastic reconfiguration,” Phys. Rev. Lett. 80, 4558–4561 (1998).
  * Sorella (2001) Sandro Sorella, “Generalized Lanczos algorithm for variational quantum Monte Carlo,” Phys. Rev. B 64, 024512 (2001).
  * Note (4) The same optimization scheme has been developed independently in artificial-intelligence community, called “natural gradient” amari_1992; amari_1998.
  * Kawamura  _et al._ (2017) Mitsuaki Kawamura, Kazuyoshi Yoshimi, Takahiro Misawa, Youhei Yamaji, Synge Todo, and Naoki Kawashima, “Quantum lattice model solver $\mathcal{H}\Phi$,” Computer Physics Communications 217, 180–192 (2017).
  * Jiang  _et al._ (2012) Hong-Chen Jiang, Hong Yao, and Leon Balents, “Spin liquid ground state of the spin-$\frac{1}{2}$ square $J_{1}$–$J_{2}$ Heisenberg model,” Phys. Rev. B 86, 024424 (2012).
  * Hu  _et al._ (2013) Wen-Jun Hu, Federico Becca, Alberto Parola, and Sandro Sorella, “Direct evidence for a gapless ${Z}_{2}$ spin liquid by frustrating Néel antiferromagnetism,” Phys. Rev. B 88, 060402 (2013).
  * Morita  _et al._ (2015) Satoshi Morita, Ryui Kaneko, and Masatoshi Imada, “Quantum Spin Liquid in Spin 1/2 $J_{1}$–$J_{2}$ Heisenberg Model on Square Lattice: Many-Variable Variational Monte Carlo Study Combined with Quantum-Number Projections,” J. Phys. Soc. Jpn. 84, 024720 (2015).
  * Ferrari and Becca (2020) Francesco Ferrari and Federico Becca, “Gapless spin liquid and valence-bond solid in the $J_{1}$–$J_{2}$ Heisenberg model on the square lattice: Insights from singlet and triplet excitations,” Phys. Rev. B 102, 014417 (2020).
  * Liu  _et al._ (2020) Wen-Yuan Liu, Shou-Shu Gong, Yu-Bin Li, Didier Poilblanc, Wei-Qiang Chen, and Zheng-Cheng Gu, “Gapless quantum spin liquid and global phase diagram of the spin-1/2 $J_{1}$–$J_{2}$ square antiferromagnetic Heisenberg model,” (2020), arXiv:2009.01821 [cond-mat.str-el] .
  * Nomura and Imada (2020) Yusuke Nomura and Masatoshi Imada, “Dirac-type nodal spin liquid revealed by machine learning,” (2020), arXiv:2005.14142 [cond-mat.str-el] .
  * Nori and Zimanyi (1991) F Nori and G. T Zimanyi, “Comparison of Frustration and Doping Effects in Generalized Flux Phases,” Europhysics Letters (EPL) 16, 397–403 (1991).
  * Nori  _et al._ (1992) Franco Nori, Eduardo Gagliano, and Silvia Bacci, “Does frustration describe doping in models for high-temperature superconductors?” Phys. Rev. Lett. 68, 240–243 (1992).
  * Nori  _et al._ (1995) Franco Nori, R. Merlin, Stephan Haas, Anders W. Sandvik, and Elbio Dagotto, “Magnetic Raman Scattering in Two-Dimensional Spin-1/2 Heisenberg Antiferromagnets: Spectral Shape Anomaly and Magnetostrictive Effects,” Phys. Rev. Lett. 75, 553–556 (1995).
  * Liang  _et al._ (2018) Xiao Liang, Wen-Yuan Liu, Pei-Ze Lin, Guang-Can Guo, Yong-Sheng Zhang, and Lixin He, “Solving frustrated quantum many-particle models with convolutional neural networks,” Phys. Rev. B 98, 104426 (2018).
  * Choo  _et al._ (2019) Kenny Choo, Titus Neupert, and Giuseppe Carleo, “Two-dimensional frustrated ${J}_{1}\text{$-$}{J}_{2}$ model studied with neural network quantum states,” Phys. Rev. B 100, 125124 (2019).
  * Ferrari  _et al._ (2019) Francesco Ferrari, Federico Becca, and Juan Carrasquilla, “Neural Gutzwiller-projected variational wave functions,” Phys. Rev. B 100, 125131 (2019).
  * Westerhout  _et al._ (2020) Tom Westerhout, Nikita Astrakhantsev, Konstantin S. Tikhonov, Mikhail I. Katsnelson, and Andrey A. Bagrov, “Generalization properties of neural network approximations to frustrated magnet ground states,” Nature Communications 11, 1593 (2020).
  * Szabó and Castelnovo (2020) Attila Szabó and Claudio Castelnovo, “Neural network wave functions and the sign problem,” Phys. Rev. Research 2, 033075 (2020).
  * Nomura (2021) Yusuke Nomura, “Helping restricted boltzmann machines with quantum-state representation by restoring symmetry,” J. Phys.: Condens. Matter 33, 174003 (2021).
  * Carleo  _et al._ (2019b) Giuseppe Carleo, Kenny Choo, Damian Hofmann, James E.T. Smith, Tom Westerhout, Fabien Alet, Emily J. Davis, Stavros Efthymiou, Ivan Glasser, Sheng-Hsuan Lin, Marta Mauri, Guglielmo Mazzola, Christian B. Mendl, Evert van Nieuwenburg, Ossian O’Reilly, Hugo Theveniaut, Giacomo Torlai, Filippo Vicentini, and Alexander Wietek, “Netket: A machine learning toolkit for many-body quantum systems,” SoftwareX 10, 100311 (2019b).

Supplemental Materials for: Purifying Deep Boltzmann Machines for Thermal
Quantum States

##  Appendix A Practical details for Method (I)

###  A.1 Analytical imaginary time evolution for transverse-field Ising model

In Method (I), we construct a DBM such that the imaginary-time evolution is
realized analytically within the Trotter error. As a concrete example, let us
consider the transverse-field Ising (TFI) model on the $N_{\rm site}$ spin
chain with periodic boundary condition,

$\displaystyle\mathcal{H}=\mathcal{H}_{1}+\mathcal{H}_{2},$ (S1)

where ${\mathcal{H}_{1}}=-J\sum_{i}\sigma^{z}_{i}\sigma^{z}_{i+1}$ is the
Ising-type interaction with amplitude $J$ and
${\mathcal{H}_{2}}=-\Gamma\sum_{i}\sigma^{x}_{i}$ is the transverse magnetic
field with amplitude $\Gamma$, where $\sigma_{i}^{a}\ (a=x,y,z)$ denotes the
Pauli matrix operating on the $i$-th site. In the following, we provide the
solutions to realize a small-time propagation $\delta_{\tau}$ by updating the
DBM parameters from ${\mathcal{W}}$ to $\bar{\mathcal{W}}$, i.e.,
$e^{-\delta_{\tau}\mathcal{H}_{\nu}}|\Psi_{\mathcal{W}}\rangle=C|\Psi_{\bar{\mathcal{W}}}\rangle$
($C$ is a constant) Carleo  _et al._ (2018).

Figure S1:  Analytical construction of a DBM which realizes the imaginary-time
evolution of the transverse-field Ising model. (a) Starting from the DBM
representing the infinite-temperature state, we first encode the interaction
propagator ($\nu=1$) by introducing a single hidden spin for every interacting
visible spins. Then, for the magnetic-field propagator ($\nu=2$), we add deep-
and hidden-spin layers in between the visible layer and neighboring hidden
layer. (b) As we repetitively encode the propagators, the network structure
grows vertically, with the depth increasing in proportion to the number of
Trotter steps.

For $\nu=1$ (interaction propagator), a solution is to add one hidden spin per
bond and put the couplings with strength $\frac{1}{2}{\rm
arcosh}(e^{2J\delta_{\tau}})$ to the visible spins on each bond [See Fig
S1(a)].

This is followed by the propagator for $\nu=2$ (transverse-field propagator),
for which we perform the following for each site: (1) cut the existing
couplings between the hidden and visible spins and (2) add a new deep spin
$d_{\rm new}$ that is coupled with the hidden spins.

After applying (1) and (2) for each site, we furthermore introduce a new
hidden spin which is coupled to the visible and the newly added deep spins
$\\{d_{\rm new}\\}$ with the interaction $\frac{1}{2}{\rm
arcosh}\left(\frac{1}{{\rm tanh}(\Gamma\delta_{\tau})}\right)$.

By repeating these operations until the desired inverse temperature is
reached, we obtain a DBM that exactly encodes the Gibbs state via the quantum-
to-classical correspondence [See Fig S1(b)]. Here, the number of the hidden
and deep spins scales as $O(N_{\tau}N_{\rm site})$, where $N_{\tau}$ is the
number of Trotter steps and $N_{\rm site}$ is the number of physical spins.

While we have described the case for the lowest-order Suzuki-Trotter
decomposition, we may alternatively consider higher-order ones so that the
errors are suppressed. In the actual numerical implementation for the 1D TFI
model, we employ the 2nd-order expansion

$\displaystyle{\rm
exp}\left[-\delta_{\tau}(\mathcal{H}_{1}+\mathcal{H}_{2})\right]$
$\displaystyle=$ $\displaystyle{\rm exp}(-\delta_{\tau}\mathcal{H}_{1}/2)\
\\!{\rm exp}(-\delta_{\tau}\mathcal{H}_{2})\ \\!{\rm
exp}(-\delta_{\tau}\mathcal{H}_{1}/2)+\mathcal{O}(\delta_{\tau}^{3}).$ (S2)

###  A.2 Monte Carlo sampling scheme for Method (I)

Here, we describe how we calculate the expectation value of a physical
observable $\langle{\mathcal{O}}\rangle$ from the constructed DBM states.
Considering that the norm of the DBM states is given by

$\displaystyle\langle\Psi|\Psi\rangle=\sum_{\sigma,\sigma^{\prime}}\bigl{|}\Psi(\sigma,\sigma^{\prime})\bigr{|}^{2}=\sum_{\sigma,\sigma^{\prime}}\sum_{h_{1},d_{1},h_{2},d_{2}}\phi^{\ast}(\sigma,\sigma^{\prime};h_{1},d_{1})\phi(\sigma,\sigma^{\prime};h_{2},d_{2})=\sum_{\sigma,\sigma^{\prime}}\sum_{h,d}w(\sigma,\sigma^{\prime};h,d)$
(S3)

with
$w(\sigma,\sigma^{\prime};h,d)\equiv\phi^{*}(\sigma,\sigma^{\prime};h_{1},d_{1})\phi(\sigma,\sigma^{\prime};h_{2},d_{2})$,
the expectation value $\langle{\mathcal{O}}\rangle$ can be evaluated as

$\displaystyle\langle{\mathcal{O}}\rangle=\frac{\langle\Psi(T)|{\mathcal{O}}\otimes{\mathbb{1}}^{\prime}|\Psi(T)\rangle}{\langle\Psi(T)|\Psi(T)\rangle}=\frac{\sum_{\sigma,\sigma^{\prime}}\sum_{h,d}w(\sigma,\sigma^{\prime};h,d)O_{\rm
loc}(\sigma,\sigma^{\prime};h,d)}{\sum_{\sigma,\sigma^{\prime}}\sum_{h,d}w(\sigma,\sigma^{\prime};h,d)},$
(S4)

where the sum over $\sigma$, $\sigma^{\prime}$, $h$, and $d$ is numerically
approximated by MC sampling with weight $w(\sigma,\sigma^{\prime};h,d)$. Here,
the “local” observable $O_{\rm loc}(\sigma,\sigma^{\prime};h,d)$ reads

$\displaystyle O_{\rm
loc}(\sigma,\sigma^{\prime};h,d)=\frac{1}{2}\sum_{\varsigma}\left(\left\langle\varsigma\right|\mathcal{O}\left|\sigma\right\rangle\frac{\phi^{*}(\varsigma,\sigma^{\prime};h_{1},d_{1})}{\phi^{*}(\sigma,\sigma^{\prime};h_{1},d_{1})}+\left\langle\sigma\right|\mathcal{O}\left|\varsigma\right\rangle\frac{\phi(\varsigma,\sigma^{\prime};h_{2},d_{2})}{\phi(\sigma,\sigma^{\prime};h_{2},d_{2})}\right),$
(S5)

where the sum over $\varsigma=(\varsigma_{1},...,\varsigma_{N_{\rm site}})$ is
taken over all configurations. Note that non-diagonal elements
$\langle\varsigma|\mathcal{O}|\sigma\rangle$ ($\varsigma\neq\sigma$) are
mostly zero when the size of the support of $\mathcal{O}$ is finite. For
instance, if $\mathcal{O}$ is a product of Pauli operators with the total
number of bit-flipping operators (i.e., $\sigma_{x}$ or $\sigma_{y}$) given as
$k$, the number of non-zero elements contributing in Eq. (S5) is $2^{k}$.

Alternatively, we can use a marginal probability as the weight for the Monte
Carlo (MC) method by tracing out one of $h$ and $d$ degrees of freedom. Here,
as an example, we show the case where the $h$ spins are traced out (we can
also trace out the $d$ spins and sample the $h$ spins). In this case, the
marginal probability is given by

$\displaystyle\tilde{w}(\sigma,\sigma^{\prime};d)\equiv\sum_{h_{1},h_{2}}\phi^{*}(\sigma,\sigma^{\prime};h_{1},d_{1})\phi(\sigma,\sigma^{\prime};h_{2},d_{2})=\tilde{\phi}^{*}(\sigma,\sigma^{\prime};d_{1})\tilde{\phi}(\sigma,\sigma^{\prime};d_{2}),$
(S6)

with

$\displaystyle\tilde{\phi}(\sigma,\sigma^{\prime};d)=\prod_{j}2\cosh\left[b_{j}+\sum_{i}(W_{ji}\sigma_{i}+W^{\prime}_{ji}\sigma_{i}^{\prime})+\sum_{k}W^{\prime}_{jk}d_{k}\right].$
(S7)

The formula for the expectation value $\langle{\mathcal{O}}\rangle$ is recast
as

$\displaystyle\langle{\mathcal{O}}\rangle=\frac{\sum_{\sigma,\sigma^{\prime}}\sum_{d}\tilde{w}(\sigma,\sigma^{\prime};d)\tilde{O}_{\rm
loc}(\sigma,\sigma^{\prime};d)}{\sum_{\sigma,\sigma^{\prime}}\sum_{d}\tilde{w}(\sigma,\sigma^{\prime};d)},$
(S8)

where the local observable $\tilde{O}_{\rm loc}(\sigma,\sigma^{\prime};d)$ is
given by

$\displaystyle\tilde{O}_{\rm
loc}(\sigma,\sigma^{\prime};d)=\frac{1}{2}\sum_{\varsigma}\left(\left\langle\varsigma\right|\mathcal{O}\left|\sigma\right\rangle\frac{\tilde{\phi}^{*}(\varsigma,\sigma^{\prime};d_{1})}{\tilde{\phi}^{*}(\sigma,\sigma^{\prime};d_{1})}+\left\langle\sigma\right|\mathcal{O}\left|\varsigma\right\rangle\frac{\tilde{\phi}(\varsigma,\sigma^{\prime};d_{2})}{\tilde{\phi}(\sigma,\sigma^{\prime};d_{2})}\right).$
(S9)

In the present study, we evaluated the physical observables based on Eq. (S8),
since the degrees of freedom participating in the sampling are mitigated.

##  Appendix B Practical details of Method (II)

###  B.1 Stochastic Reconfiguration method

The Stochastic Reconfiguration (SR) method Sorella (1998, 2001) employs a
variational principle such that the exact imaginary-time evolution is
approximated within the expressive power of the wave function ansatz. In other
words, the SR method provides an update rule of variational parameters so that
the distance, whose definition relies on the employed variational principle,
between states following exact and approximate imaginary-time evolution is
minimized. In the following, we provide a concise derivation of the parameter
update rule used in the SR method.

Let $\ket{\Psi_{\theta}}$ be a variational wave function with a set of complex
variational parameters $\theta$. In the SR method, we update the variational
parameters as $\theta\leftarrow\theta+\delta\widetilde{\theta}$ according to a
variational principle based on the Fubini-Study metric $\mathcal{F}$:

$\displaystyle\delta\widetilde{\theta}$ $\displaystyle=$
$\displaystyle\mathop{\rm
arg~{}min}\limits_{\delta\theta}\mathcal{F}[e^{-\delta_{\tau}\mathcal{H}}\ket{\Psi_{\theta}},\ket{\Psi_{\theta+\delta\theta}}],$
(S10) $\displaystyle\mathcal{F}[\ket{\psi},\ket{\phi}]:$ $\displaystyle=$
$\displaystyle\arccos\sqrt{\frac{\bra{\psi}\ket{\phi}\bra{\phi}\ket{\psi}}{\bra{\psi}\ket{\psi}\bra{\phi}\ket{\phi}}}.$
(S11)

Here, $\mathcal{H}$ is the Hamiltonian of the system and $\delta_{\tau}$ is a
small step of imaginary-time evolution. Note that the metric $\mathcal{F}$
measures the distance between two quantum states $\ket{\psi}$ and $\ket{\phi}$
after imposing an appropriate normalization, and therefore it is closely
related with the fidelity between pure states $F$ as
$\mathcal{F}=\arccos\sqrt{F^{2}}$, where the fidelity is defined as
$F^{2}[\ket{\psi},\ket{\phi}]=\frac{\bra{\psi}\ket{\phi}\bra{\phi}\ket{\psi}}{\bra{\psi}\ket{\psi}\bra{\phi}\ket{\phi}}.$
Consequently, the parameter update $\delta\widetilde{\theta}$ also satisfies

$\displaystyle\delta\widetilde{\theta}=\mathop{\rm
arg~{}max}\limits_{\delta\theta}F^{2}[e^{-\delta_{\tau}\mathcal{H}}\ket{\Psi_{\theta}},\ket{\Psi_{\theta+\delta\theta}}].$
(S12)

After some algebraic calculations, we obtain the explicit expression, up to
the second order, of the squared fidelity $F^{2}$ as

$\displaystyle F^{2}$ $\displaystyle=$ $\displaystyle
1-\left(\sum_{k,l}\delta\theta_{k}^{*}S_{kl}\delta\theta_{l}+\delta_{\tau}\sum_{k}(f_{k}\delta\theta_{k}^{*}+{\rm
c.c.})+\delta_{\tau}^{2}(\langle\mathcal{H}^{2}\rangle-\langle\mathcal{H}\rangle^{2})\right),$
(S13)

where the variance of the energy, contributing to the intrinsic algorithm
error, is estimated using the MC sampling as
$\langle\cdot\rangle:=\frac{\bra{\Psi_{\theta}}\cdot\ket{\Psi_{\theta}}}{\bra{\Psi_{\theta}}\ket{\Psi_{\theta}}}$.
[See B.3 for further information.] We have introduced an element of the
covariance matrix $S_{kl}$ as

$\displaystyle S_{kl}$ $\displaystyle=$
$\displaystyle\frac{\bra{\partial_{k}\Psi_{\theta}}\ket{\partial_{l}\Psi_{\theta}}}{\bra{\Psi_{\theta}}\ket{\Psi_{\theta}}}-\frac{\bra{\partial_{k}\Psi_{\theta}}\ket{\Psi_{\theta}}}{\bra{\Psi_{\theta}}\ket{\Psi_{\theta}}}\frac{\bra{\Psi_{\theta}}\ket{\partial_{l}\Psi_{\theta}}}{\bra{\Psi_{\theta}}\ket{\Psi_{\theta}}}$
(S14) $\displaystyle=$ $\displaystyle\langle
O^{\dagger}_{k}O_{l}\rangle-\langle O_{k}^{\dagger}\rangle\langle
O_{l}\rangle,$ (S15)

and the generalized force $f_{k}$ as

$\displaystyle f_{k}$ $\displaystyle=$
$\displaystyle\frac{\bra{\partial_{k}\Psi_{\theta}}\mathcal{H}\ket{\Psi_{\theta}}}{\bra{\Psi_{\theta}}\ket{\Psi_{\theta}}}-\frac{\bra{\partial_{k}\Psi_{\theta}}\ket{\Psi_{\theta}}}{\bra{\Psi_{\theta}}\ket{\Psi_{\theta}}}\frac{\bra{\Psi_{\theta}}\mathcal{H}\ket{\Psi_{\theta}}}{\bra{\Psi_{\theta}}\ket{\Psi_{\theta}}}$
(S16) $\displaystyle=$ $\displaystyle\langle
O^{\dagger}_{k}\mathcal{H}\rangle-\langle
O_{k}^{\dagger}\rangle\langle\mathcal{H}\rangle,$ (S17)

where $O_{k}$ denotes a diagonal matrix with each element given as
$O_{k}(\sigma)=\frac{\partial_{k}\Psi_{\theta}(\sigma)}{\Psi_{\theta}(\sigma)}$.
Equations. (S15) and (S17) indicates that all the elements $S_{kl}$ and
$f_{k}$ can be estimated by MC sampling. Finally, by using the stationary
condition for the optimal solution of Eq. (S13), the expression for the
parameter update can be obtained as

$\displaystyle\delta\widetilde{\theta}_{k}=-\delta_{\tau}\sum_{l}S_{kl}^{-1}f_{l}.$
(S18)

In the actual calculation, we add a stabilization factor to the diagonal
elements of $S$ as $S_{kk}+\epsilon_{k}$, where $\epsilon_{k}$ is typically
taken uniformly as $\sim 10^{-4}\bar{S}_{\rm diag}$, where $\bar{S}_{\rm
diag}$ is the mean value of diagonal elements. We observe that the calculation
becomes more reliable when the time step $\delta_{\tau}$ is initially taken to
be small; $O(10^{-4})$. After tens to hundreds of iterations, $\delta_{\tau}$
is increased gradually up to $O(10^{-2})$.

###  B.2 Symmetrization

In Method (II), we numerically optimize a purified DBM wave function with the
form:

$\displaystyle\Psi(\sigma,\sigma^{\prime})=\prod_{j}2\cosh\Bigl{[}b_{j}\\!+\\!\sum_{i}(W_{ji}\sigma_{i}\\!+\\!W^{\prime}_{ji}\sigma_{i}^{\prime})\Bigr{]}.$
(S19)

Here, $b$, $W$, and $W^{\prime}$ are variational parameters. We assign the
site index $i$ for the pair of physical $\sigma$ and ancilla $\sigma^{\prime}$
spins.

To improve the quality of the calculation, we consider the symmetry of the
system. In the case of the 2D $J_{1}$-$J_{2}$ Heisenberg model, we encode the
translational and point-group symmetry in the purified DBM wave function. We
observe that the initial state, i.e., the purified infinite-temperature state,
is invariant under certain symmetry operations: (1) translation $T_{\bf R}$
that shifts all the spins by the amount ${\bf R}$ as
$(\sigma,\sigma^{\prime})\mapsto(T_{\bf R}\sigma,T_{\bf R}\sigma^{\prime})$,
and (2) symmetry operation $R$ of $C_{4v}$ point group that maps a spin
configuration as $(\sigma,\sigma^{\prime})\mapsto(R\sigma,R\sigma^{\prime})$.
Alternatively, we can understand that the purified state is in the zero wave-
number sector and belongs to the $A_{1}$ representation of the $C_{4v}$ point
group (see also Table S1) of a bilayer system. Along the imaginary time
evolution, the purified DBM wave function stays within the identical symmetry
sector at arbitrary temperature. Hence, we impose such a condition by
symmetrizing the purified DBM wave function as

$\displaystyle\Psi_{{\rm sym.}}(\sigma,\sigma^{\prime})=\sum_{R,{\bf
R}}\Psi(T_{\bf R}R\sigma,T_{\bf R}R\sigma^{\prime}).$ (S20)

We employ this formula to impose the symmetry on the wave function on the
left-hand side $\Psi_{\rm sym.}(\sigma,\sigma^{\prime})$. Observe that
$\Psi_{\rm sym.}(\sigma,\sigma^{\prime})$ satisfies the symmetry even when the
bare DBM wave function $\Psi(\sigma,\sigma^{\prime})$ on the right-hand side
does not preserve the symmetry. We emphasize that the corresponding Gibbs
state of the target system consists of contributions from all symmetry sectors
of the original Hamiltonian.

Table S1: Character table of the $C_{4v}$ point group. |  E |  $2C_{4}$ |  $C_{2}$ |  $2\sigma_{v}$ | $2\sigma_{d}$
---|---|---|---|---|---
$A_{1}$ |  1 |  1 |  1 |  1 |  1
$A_{2}$ |  1 |  1 |  1 |  $-1$ | $-1$
$B_{1}$ |  1 |  $-1$ |  1 |  1 | $-1$
$B_{2}$ |  1 |  $-1$ |  1 |  $-1$ |  1
$E$ |  2 |  0 |  $-2$ |  0 |  0

###  B.3 Computing expectation values of physical observables

As in Method (I), expectation values are numerically evaluated by the MC
method. The formula for the expectation value $\langle{\mathcal{O}}\rangle$ in
Method (II) reads

$\displaystyle\langle{\mathcal{O}}\rangle=\frac{\sum_{\sigma,\sigma^{\prime}}p(\sigma,\sigma^{\prime})O_{\rm
loc}(\sigma,\sigma^{\prime})}{\sum_{\sigma,\sigma^{\prime}}p(\sigma,\sigma^{\prime})},$
(S21)

where $p(\sigma,\sigma^{\prime})$ is the weight
$p(\sigma,\sigma^{\prime})\equiv|\Psi_{\rm
sym.}(\sigma,\sigma^{\prime})|^{2}$, and the local observable $O_{\rm
loc}(\sigma,\sigma^{\prime})$ is given by

$\displaystyle O_{\rm
loc}(\sigma,\sigma^{\prime})=\frac{1}{2}\sum_{\varsigma}\left(\left\langle\sigma\right|\mathcal{O}\left|\varsigma\right\rangle\frac{\Psi_{{\rm
sym.}}(\varsigma,\sigma^{\prime})}{\Psi_{{\rm
sym.}}(\sigma,\sigma^{\prime})}+c.c.\right).$ (S22)

We perform the Metropolis sampling over the $\sigma$ and $\sigma^{\prime}$
spins with the weight $p(\sigma,\sigma^{\prime})$ to compute expectation
values.

###  B.4 Calculation conditions

In the present calculations for the 2D $J_{1}$–$J_{2}$ Heisenberg model on the
$6\times 6$ lattice (Figure 3 in the paper), we introduce 8$N_{\rm site}$
(=288) hidden spins. We set the $b$ parameters (magnetic field) to zero, so
that the “up” and “down” spins are equivalent, and optimize only the $W$ and
$W^{\prime}$ parameters. With this setting, the purified DBM wave function
becomes even with respect to the global spin inversion.

The initial $W$ and $W^{\prime}$ parameters are prepared to represent the
infinite-temperature state. The infinite-temperature state can be reproduced
exactly by setting
$W_{ji}=W^{\prime}_{ji}=\mathrm{i}\frac{\pi}{4}\delta_{ji}$, for $1\leq j\leq
N_{\rm site}$, and $W_{ji}=W^{\prime}_{ji}=0$, for $N_{\rm site}+1\leq j\leq
8N_{\rm site}$. In the actual calculations, to make the initial gradient of
the parameter optimization finite, we put small perturbations to the above
setting by adding small random numbers.

As we describe in the main text, starting from the initial $W$ and
$W^{\prime}$ values, we optimize the parameters with the SR method Sorella
(2001), which makes it possible to reproduce the imaginary time evolution as
much as possible within the representability of the DBM. To reduce the number
of variational parameters, we take half of the $W$ and $W^{\prime}$ parameters
to be complex and the rest real as in Ref. Nomura, 2021. We do not impose
symmetry constraints on the $W$ and $W^{\prime}$ parameters, but instead, the
symmetry is restored with the projection [Eq. (S20) in the case of the 2D
$J_{1}$–$J_{2}$ Heisenberg model]. We observe that the most time-consuming
part of the calculation is the MC sampling to estimate the expectation values
of the energy and the gradient of the parameter optimization; its time scales
as ${\mathcal{O}}(N_{h}N_{\rm site}^{2})$, where $N_{h}$ is the number of
hidden spins ($N_{h}=8N_{\rm site}$ in the present case).

###  B.5 Benchmark calculations for 1D and 2D Heisenberg models

Figure S2:  Method (II) results (symbols) for the specific heat $C$ for (a)
the 1D Heisenberg model on the 16-site chain and (b,c) the 2D Heisenberg model
on the 4x4 lattice. In all cases, we take the periodic boundary condition, and
the number of sites is 16 ($N_{\rm site}=16$). In (a) and (b), we impose
translational symmetry on the purified DBM wave function for the extended
system. In (c), we use both translational and $C_{4v}$ point-group symmetry.
The solid black curves are obtained by exact diagonalization.

Here, we show the benchmark results of Method (II) for the 1D and 2D
Heisenberg models. The Hamiltonian reads ${\mathcal{H}}=J_{1}\sum_{\langle
i,j\rangle}{\bf S}_{i}\cdot{\bf S}_{j}$, where $J_{1}=1$ and $\langle
i,j\rangle$ denotes a pair of neighboring sites. In the 2D case, the
Hamiltonian corresponds to the $J_{2}=0$ case for the 2D $J_{1}$–$J_{2}$
Heisenberg model. We take the 16-site chain for the 1D model, and the $4\times
4$ square lattice for the 2D model. In both cases, periodic boundary
conditions are assumed and the number of sites is 16 ($N_{\rm site}=16$). We
consider the zero magnetization sector ($\sum_{i}S_{i}^{z}=0$) as in the main
text.

As we have discussed above and in the main text, the quality of the Method
(II) calculations can be improved by increasing the number of hidden units
$N_{h}$ and/or by imposing symmetry. For a numerical demonstration, we compare
$N_{h}=N_{\rm site}$ and $N_{h}=8N_{\rm site}$ cases. We follow the conditions
described above (Sec. B.4), except that all the parameters are taken to be
complex in the $N_{h}=N_{\rm site}$ case (note that when $N_{h}=N_{\rm site}$,
all the hidden units need to have complex couplings to prepare infinite-
temperature states).

Figure S2 shows the Method (II) results for the specific heat $C$ for (a) the
1D model when imposing translational symmetry, (b) the 2D model when imposing
translational symmetry, and (c) the 2D model when imposing translational and
point-group symmetries. We see that in the case of the 1D model, the Method
(II) result agrees well with the exact result already at $N_{h}=N_{\rm site}$
[Figure S2(a)]. On the other hand, in the case of the 2D model, when only the
translational symmetry is imposed as in the 1D model case, the result with
$N_{h}=N_{\rm site}$ shows a visible (but small) deviation from the exact
result [Figure S2(b)]. The result improves by increasing $N_{h}$ [Figure
S2(b)] or imposing an additional symmetry (point-group symmetry) [Figure
S2(c)]. From this benchmark, we indeed observe that the number of hidden units
$N_{h}$ and the symmetrization are important factors to ensure the quality of
the Method (II) calculations.