Constructing exact representations of quantum many-body systems with deep neural networks
Giuseppe Carleo Yusuke Nomura Masatoshi Imada






Center for Computational Quantum Physics, Flatiron Institute, 162
5th Avenue, New York, NY 10010, USA

Institute for Theoretical Physics, ETH Zurich, Wolfgang-Pauli-Str. 27,
8093 Zurich, Switzerland



Department of Applied Physics, The University of Tokyo, 7-3-1 Hongo,
Bunkyo-ku, Tokyo 113-8656, Japan



Department of Applied Physics, The University of Tokyo, 7-3-1 Hongo,
Bunkyo-ku, Tokyo 113-8656, Japan

We develop a constructive approach to generate artificial neural networks
representing the exact ground states of a large class of many-body
lattice Hamiltonians. It is based on the deep Boltzmann
machine architecture, in which two layers of hidden neurons mediate
quantum correlations among physical degrees of freedom in the visible
layer. The approach reproduces the
exact imaginary-time Hamiltonian evolution, and is completely deterministic. 
In turn, compact and exact network representations for the ground states are obtained without 
stochastic optimization of the network parameters.
The number of neurons grows linearly with the system size and total
imaginary time, respectively. Physical quantities can be measured
by sampling configurations of both physical and neuron degrees of
freedom. 
We provide specific examples for
the transverse-field Ising and Heisenberg models by implementing efficient sampling. 

As a compact, classical representation for many-body
quantum systems, our approach is an alternative to the standard path integral, 
and it is potentially useful also to
systematically improve on numerical approaches based on the restricted Boltzmann machine architecture. 

Constructing exact representations of quantum

 many-body systems with deep neural networks
    Masatoshi Imada
    December 30, 2023
===========================================================================================





§ INTRODUCTION


A tremendous amount of successful developments in quantum physics
builds upon the mapping between many-body quantum systems and effective
classical theories. The probably most well known mapping is due to
Feynman, who introduced an exact representation of many-body quantum
systems in terms of statistical summations over classical particles
trajectories <cit.>. Effective classical representations
of quantum many-body systems are however not unique, and other approaches
rely on different inspiring principles, such as perturbative expansions
<cit.>, or decomposition of interactions with auxiliary
degrees of freedom <cit.>.
The classical representations of quantum states allow both for novel
conceptual developments and efficient numerical simulations. On one
hand, perturbative approaches based on the graphical resummation of
classes of diagrams are at the heart of many-body analytical approaches
in various fields of research, ranging from particle to condensed-matter
physics <cit.>. On the other hand, several non-perturbative
numerical methods for many-body quantum systems are also based on
these mappings. Quantum Monte Carlo (QMC) methods are among the most
successful numerical techniques, relying on continuos-space polymer
representations <cit.>,
world-line lattice path integrals <cit.>, continuous
time algorithm <cit.>, summation of perturbative diagrams <cit.>.
Effective classical representations are also the building block of
variational methods based on correlated many-body wave-functions <cit.>.
Several successful variational techniques make extensive use of parametric
representations of quantum states, where the effective parameters
are determined by means of the variational principle <cit.>.
In matrix-product and tensor-network-states the ground-state is expressed
as a classical network <cit.>. In general, finding alternative,
efficient classical representations of quantum states can help establishing
novel numerical and analytical techniques to study challenging open
issues.

Recently, an efficient variational representation of many-body systems
in terms of artificial neural networks, which consists of classical
degrees of freedom, has been introduced <cit.>.
Numerical results have shown that artificial neural networks can represent
many-body states with high accuracy <cit.>.
The majority of the variational approaches adopted so-far are based
on shallow neural networks, called Restricted Boltzmann Machines (RBM),
in which the physical degrees of freedom interact with an ensemble
of hidden degrees of freedom (neurons). While shallow RBM states have
promising features in terms of entanglement capacity <cit.>,
only deep networks are guaranteed to provide a complete and efficient
description of the most general quantum states <cit.>.

In this Paper we introduce a constructive approach to explicitly generate
deep network structures corresponding to exact quantum many-body ground
states. We demonstrate this construction for interacting lattice spin
models, including the transverse-field Ising and Heisenberg models.
Our constructions are fully deterministic, in stark contrast to the
shallow RBM case, in which the numerical optimization of the network
parameters is inevitable. The number of neurons required in the construction
scales only polynomially with the system size, thus the present approach
constitutes a new family of efficient quantum-to-classical mappings
exhibiting a prominent representational flexibility. Given as a simple
set of iterative rules, these constructions can be used both as a
self-standing tool, or to systematically improve results obtained
with variational shallow networks. The latter improves the efficiency
of the method because the numerically optimized shallow RBM states
are already good approximations for ground states. Finally, we discuss
sampling strategies from the generated deep networks and show numerical
results for one-dimensional spin models.





§ CONSTRUCTION OF DEEP NEURAL STATES


The ground state of a generic Hamiltonian, ℋ, can be
found through imaginary-time evolution, |Ψ(τ)⟩=e^-τℋ|Ψ_0⟩,
for a sufficiently large τ≫Δ E^-1. Here Δ E
is the energy gap between the ground and the first excited state,
and |Ψ_0⟩ is an arbitrary initial state non-orthogonal
to the exact ground state. For a finite system, the energy gap is
typically finite, and the total propagation time needed to reach the
ground state within an arbitrary given accuracy is expected to grow
at most polynomially with the system size (for systems becoming gapless
in the thermodynamic limit). 

Here, we introduce a representation of
the wave-function coefficients in terms of a deep Boltzmann machine
(DBM) <cit.>. For the sake of concreteness,
let us consider the case of N spins, described by the quantum numbers
|σ^z⟩=|σ_1^z…σ_N^z⟩. Then,
we represent generic many-body amplitudes ⟨σ_1^z…σ_N^z|Ψ⟩≡Ψ(σ^z)
in the two-layer DBM form: 

    Ψ_𝒲(σ^z)=∑_{h,d}exp[∑_ia_iσ_i^z+∑_ijσ_i^zW_ijh_j+.
    .+∑_jb_jh_j+∑_jkh_jd_kW_jk^'+∑_kb_k^'d_k]

where we have introduced M hidden units h, M^' deep
units d, and a set of couplings and bias terms 𝒲≡(a,b,b^',W,W^').
A sketch of the DBM architecture is shown in Fig. <ref>.

In the following, we specialize to the case of spin 1/2, thus all
the units are taken to be σ^z,h,d=±1. This representation
is the natural deep-network generalization of the shallow RBM, introduced
as variational ansatz in Ref. <cit.>. As for the
RBM form, also in this case direct connections between variables in
the same layer are not allowed. A crucial difference is however that
the layer of deep variables makes, in general, the evaluation of the
wave-function amplitudes not possible analytically. At variance with
RBM, the DBM form is known to be universal, as proven by Gao and Duan
recently <cit.>. In order to find explicit expressions
for the parameters 𝒲 that represent |Ψ(τ)⟩
for arbitrary imaginary time, we start considering a second-order
Trotter-Suzuki decomposition <cit.>: 

    |Ψ(τ)⟩=𝒢_1(δ_τ/2)𝒢_2(δ_τ)…𝒢_1(δ_τ)𝒢_2(δ_τ)𝒢_1(δ_τ/2)|Ψ_0⟩,

where we have decomposed the Hamiltonian into two non-commuting parts,
ℋ=ℋ_1+ℋ_2, and introduced the
short-time propagators 𝒢_ν(δ_τ)=e^-ℋ_νδ_τ.
The problem of finding an exact representation for |Ψ(τ)⟩
then reduces to finding an exact representation for each of the two
type of propagators. As shown in the following concrete examples for
paradigmatic spin models, thanks to the high representability of DBM,
the imaginary time evolution can be tracked exactly by dynamically
modifying the DBM network structure. In practice, this is achieved
either by changing parameters 𝒲 at each step of the imaginary
time evolution, or by introducing additional parameters in 𝒲,
adding new neurons and creating new connections in the network.





 §.§ Transverse-Field Ising model


We start considering the transverse-field Ising (TFI) model on an
arbitrary interaction graph. In this case, we decompose the Hamiltonian
into two parts: ℋ_1=-∑_lΓ_lσ_l^x,
and ℋ_2=∑_l<mV_lmσ_l^zσ_m^z,
where σ denote Pauli matrices, Γ_l (>0) are site-dependent
transverse fields, and V_lm are arbitrary coupling constants.

In order to implement the mapping to a DBM, we first consider the
action of the diagonal propagator e^-δ_τV_lmσ_l^zσ_m^z,
acting on a bond V_lm. In this case, the goal of finding an exact
DBM representation can be rephrased as finding solutions to 

    ⟨σ^z|e^-δ_τV_lmσ_l^zσ_m^z|Ψ_𝒲⟩    =     CΨ_𝒲̅(σ^z),

i.e. finding a set of new parameters 𝒲̅ that exactly
reproduces the imaginary time evolution on the left hand side. Here
C is an arbitrary finite normalization constant. The diagonal propagator
introduces an interaction between two visible, physical spins, which
is not directly available in the DBM architecture. This interaction
can be mediated by a new hidden unit in the first layer, h_[lm]
which is only connected to the visible spins on that bond, i.e. W̅_l[lm]
and W̅_m[lm] are finite, but W̅_i[lm]=0,∀ i≠ l,m
and W̅_j[lm]^'=0,∀ j [see Fig. <ref>(a)].

More concretely, the new wave function has then the form: 

    Ψ_𝒲̅(σ^z)     =    ∑_h_[lm]e^σ_l^zW_l[lm]h_[lm]+σ_m^zW_m[lm]h_[lm]Ψ_𝒲(σ^z)
        =     2cosh(σ_l^zW_l[lm]+σ_m^zW_m[lm])Ψ_𝒲(σ^z).

Equation (<ref>) is then satisfied if 

    e^-δ_τV_lmσ_l^zσ_m^z=2Ccosh(σ_l^zW_l[lm]+σ_m^zW_m[lm])

for all the possible values of σ_l^z and σ_m^z.
By means of a useful identity [Eq. (<ref>) in Methods],
the new parameters W_l[lm] and W_m[lm] are given by 

    W_l[lm]    =    1/2arcosh(e^2|V_lm|δ_τ)
    
    W_m[lm]    =     -sgn(V_lm)× W_l[lm].

In this way the classical two-body interaction can, in general, be
represented exactly by the shallow RBM.

Next, to exactly represent the off-diagonal propagator e^δ_τΓ_lσ_l^x|Ψ_𝒲⟩,
we must solve: 

    cosh(Γ_lδ_τ)Ψ_𝒲(σ^z)+sinh(Γ_lδ_τ)Ψ_𝒲(σ_l^z→-σ_l^z)=
    
    =CΨ_𝒲̅(σ^z)

for the new weights 𝒲̅, and for an appropriate finite
normalization constant C. In this case, one possible solution is
obtained by adding one deep d_[l] and one hidden h_[l] neurons.
For d_[l], we create new couplings W_j[l]^' to the
existing hidden neurons h_j which are connected to σ_l^z.
We simultaneously allow for changes in the existing parameters. By
the procedure given in Methods, after applying the off-diagonal propagator
for the site l, a solution of Eq.(<ref>) is found
by the matching condition of the hidden unit interactions on the left
and the right hand sides of Eq.(<ref>). Overall, the
solution results in a three-step process [Fig. <ref>(b)]:
First, the hidden units attached to σ_l^z are connected
to the newly introduced deep unit d_[l] as 

    W_j[l]^'    =     -W_lj

(see Eq.(<ref>)). Second, all the hidden units previously
connected to the spin σ_l^z lose their connection, i.e.,
W̅_lj=0,∀ j. Third, the spin σ_l^z and
the deep unit d_[l] are connected to the new hidden unit, h_[l],
through the interactionW_l[l] and W'_[l][l], respectively
as 

    W_l[l]    =    1/2arcosh(1/tanh(Γ_lδ_τ)),
    
    W_[l][l]^'    =     -W_l[l].


Using the given expressions for the parameters 𝒲̅
we can then exactly implement a single step of imaginary-time evolution.
The full imaginary-time evolution is achieved by applying the above
procedure for ℋ_1 and ℋ_2 alternately
and repeatedly. Example applications of these rules, for both the
diagonal and the off-diagonal propagators are shown in Fig. <ref>.






 §.§ Heisenberg model


We now consider the anti-ferromagnetic Heisenberg (AFH) model, on
bipartite lattices. In one dimension, we decompose the Hamiltonian
into odd and even bonds: ℋ_1=∑_⟨ l,m⟩^oddℋ_lm^bond
and ℋ_2=∑_⟨ l,m⟩^evenℋ_lm^bond,
with ℋ_lm^bond=J(σ_l^xσ_m^x+σ_l^yσ_m^y+σ_l^zσ_m^z),
where σ denote Pauli matrices. Because the bond Hamiltonian
ℋ_lm^bond is a building block also in higher
dimensional models, construction of an exact DBM representation of
the ground states can be achieved by finding solutions for the bond-propagator
⟨σ^z|e^-δ_τℋ_lm^bond|Ψ_𝒲⟩=C⟨σ^z|Ψ_𝒲̅⟩,
where the parameters 𝒲̅ are such that the previous
equation is satisfied for all the possible ⟨σ^z|,
and for an arbitrary finite normalization constant C. More explicitly,
we need to satisfy 

    δ_σ_l^z,σ_m^ze^-Jδ_τΨ_𝒲(σ^z)+(1-δ_σ_l^z,σ_m^z)e^Jδ_τcosh(2Jδ_τ)×
    (Ψ_𝒲(σ^z)-tanh(2Jδ_τ)Ψ_𝒲(σ_l^z↔σ_m^z))=CΨ_𝒲̅(σ^z).

The basic strategy of finding a solution for Eq.(<ref>)
is similar to that for Eq.(<ref>) in the transverse
Ising model. Several possibilities arise when looking for solutions
of the bond-propagator equation, Eq. (<ref>). The
existence of non-equivalent solutions prominently shows the non-uniqueness
of DBM structure to represent the very same state and, at the same
time, provides us flexibility in designing DBM architectures. Here,
we show three concrete constructions. See Methods and Supplementary
Information (II-B) for a detailed derivation of the DBM construction
for the Heisenberg model, including anisotropic and bond-disordered
coupling cases.



  §.§.§ 1 deep, 3 hidden


The first construction is dubbed “1 deep, 3 hidden” (1d-3h).
It amounts to adding an extra deep neuron, d_[lm], and three
more hidden neurons to satisfy Eq. (<ref>). A crucial
difference with respect to the TFI model is that the introduced deep
spin d_[lm] has a constraint depending on the state of the spins
on the bond: σ_l^z and σ_m^z. Specifically,
when σ_l^z=σ_m^z the deep spin is constrained
to be d_[lm]=σ_l^z=σ_m^z, whereas when σ_l^z≠σ_m^z,
its value is unconstrained. From a pictorial point of view, the action
of the bond propagator is a four-step process [see Fig. <ref>(a)].
Starting from a given initial network (uppermost structures in Fig.
<ref>), d_[lm] is added and connected,
through W'_j[lm] given in Eq. (<ref>), to the existing
hidden units h_j connected to σ_l^z and σ_m^z.
Second, spin σ_l^z is disconnected to all hidden units
and reconnected to those hidden units the spin σ_m^z is
attached to [see Eq. (<ref>)]. Third, two new hidden
units are introduced. One of the hidden units, h_[lm1], mediates
the interaction between σ_l^z and d_[lm] [Eq. (<ref>)],
and the other hidden unit h_[lm2] mediates a direct spin-spin
interaction between σ_l^z and σ_m^z [Eq.
(<ref>)]. Fourth, a further hidden unit connected
to σ_l^z, σ_m^z and d_[lm] is inserted,
in such a way that the constraint previously described is satisfied.
For all but the last step, the DBM weights are real-valued. In the
last step instead the constraint is enforced by introducing imaginary-valued
interactions (dotted lines in Fig. <ref>),
referred to the “iπ/6` trick, resulting in a sign-problem
free global term cos(π/6(σ_l^z+σ_m^z-d_[lm]))
after the summation over ±1 for the lastly added hidden unit
h_[lm3]: ∑_h_[lm3]=±1exp[iπ/6(σ_l^z+σ_m^z-d_[lm])h_[lm3]].
The constraint mentioned above is assured by this cosine term. 



  §.§.§ 2 deep, 6 hidden


The second construction is dubbed “2 deep, 6 hidden” (2d-6h),
and is more similar to the lattice path-integral formulation. In this
representation, we introduce two auxiliary deep spins per bond, d_[l]
and d_[m] with constraint d_[l]+d_[m]=σ_l^z+σ_m^z,
and six hidden neurons. The action of the bond propagator is schematically
illustrated in Fig. <ref>(b): first, two deep
units d_[l] and d_[m] are introduced, connecting, respectively,
to the hidden units spins σ_l^z and σ_m^z
are attached to [see Eqs. (<ref>), (<ref>)].
Second, all the connections between spins σ_l^z, σ_m^z
and hidden units h_j are cut off [Eqs. (<ref>),
(<ref>)]. Third, four hidden units h_[lm1],…,h_[lm4]
are introduced, to mediate interactions between the two deep units
and the physical spins l,m [Eqs. (<ref>),
(<ref>)]. Finally, two hidden units h_[lm5]
and h_[lm6] are introduced, connecting both to d_[l],d_[m]
and σ_l^z,σ_m^z with imaginary-valued weights.
The last step realizes the constraint d_[l]+d_[m]=σ_l^z+σ_m^z,
through the “iπ/4, iπ/8` trick discussed in Methods
and Supplementary Information (II-B.2). 

In this representation, if
the hidden neurons are traced out, the imaginary-time evolution becomes
equivalent to that of the path-integral Monte Carlo method. More specifically,
the number of deep neurons introduced at each time slice is exactly
the same as the number of visible spins, and the deep neurons at each
time slice can be regarded as additional classical spin degrees of
freedom in the path-integral. Moreover, the constraint d_[l]+d_[m]=σ_l^z+σ_m^z
ensures that the total magnetization is conserved at each time slice.
Finally, the W and W' interactions reproduce the matrix element
of exp(-δ_τℋ_lm^ bond) between
neighboring time slices. See Supplementary Information (II-B.2) for
more detail on this point.



  §.§.§ 2 deep, 4 hidden


A further possible solution to Eq. (<ref>) is dubbed
“2 deep, 4 hidden” (2d-4h) construction. In this case,
we introduce two auxiliary deep variables d_[l] and d_[lm].
We also introduce four hidden units h_[l], h_[m], h_[lm1],
and h_[lm2]. Before the imaginary time evolution, e^-δ_τℋ_lm^ bond,
the physical variables σ_n^z (n=l or m) are already
coupled to each hidden variable h_j with a coupling W_nj.
After the time evolution e^-δ_τℋ_lm^ bond,
as shown schematically in Fig. <ref>(c), the
coupling parameters are updated in the following way based on the
old W_nj: First, the first deep unit d_[l] becomes coupled
to the already existing hidden variables h_j through the coupling
W'_j[l] given in Eq. (<ref>). The second deep unit d_[lm]
becomes similarly coupled to h_j through a term Z_lmj given
in Eq. (<ref>). Second, W_nj is updated to W̅_nj=W_nj+Δ W_nj
[see Eq. (<ref>)]. Third, newly introduced h_[n] (n=l
or m) gets coupled to d_[l] through W'_[n][l], and also
to σ_n^z through W_n[n] [Eqs. (<ref>), (<ref>)].

Within this construction, and as clarified in Methods, we also need to satisfy the constraint
d_[l]d_[lm]=σ_l^zσ_m^z. Such a constraint
is represented in DBM form as 

    ∑_h_[lm1],h_[lm2]exp[iπ/4(h_[lm1]+h_[lm2])(σ_l^z+σ_m^z+d_[l]+d_[lm])],

which ensures d_[l]d_[lm]=σ_l^zσ_m^z after
explicit summation of h_[lm1] and h_[lm2]. 

Finally, we remark that the three constructions presented here have different intrinsic
network topologies. In particular, 2d-6h gives rise to a local topology
(because of the equivalence with the path-integral contruction), 1d-3h
has a local structure in the first layer and non-local in the second
one, and 2d-4h is purely non-local in both layers (see Supplementary Information II.B). 




§ SAMPLING STRATEGIES


With network structures explicitly determined, we now focus on the
problem of extracting meaningful physical quantities from them. To
this end, it is convenient to decompose the DBM weight into two parts,
such that 

    Ψ_𝒲(σ^z)=∑_{h,d}P_1(σ^z,h)P_2(h,d),

where P_1(σ^z,h)=e^σ^z· a+σ^z· W· h+h· b,
and P_2(h,d)=e^h· W^'· d+d· b^'.
The expectation value of an arbitrary (few-body) operator 𝒪
can then be computed through the expression 

    ⟨𝒪⟩    =    ∑_{σ^z,h,h^'d,d^'}Π(σ^z,h,h^',d,d^')O_loc(σ^z,h,h^')/∑_{σ^z,h,h^'d,d^'}Π(σ^z,h,h^',d,d^'),

where we have introduced the pseudo-probability density Π(σ^z,h,h^',d,d^')≡ P_1(σ^z,h)P_2(h,d)P_1^⋆(σ^z,h^')P_2^⋆(h^',d^'),
and the “local” estimator O_loc(σ^z,h,h^')=1/2∑_σ^' z⟨σ^z|𝒪|σ^' z⟩(P_1(σ^' z,h)/P_1(σ^z,h)+P_1(σ^' z,h^')^⋆/P_1(σ^z,h^')^⋆).

For the sampling over the Π distribution, a block Gibbs sampling
analogous to what performed in standard DBM architectures can be performed
<cit.>.
Alternatively, it is possible to devise a set of Metropolis local
updates sampling the exactly known marginals Π̃(σ^z,h,h^')=∑_{d,d^'}Π(σ^z,h,h^',d,d^')
or Π̃'(σ^z,d,d^')=∑_{h,h^'}Π(σ^z,h,h^',d,d^').
We can also employ efficient cluster updates. Sampling is discussed
more in detail in the Supplementary Information (III).



§ NUMERICAL RESULTS


We have implemented numerical algorithms to sample and obtain physical
properties from the DBM previously derived. In Fig. <ref>
(a) we show results for the one-dimensional TFI model. Specifically,
we show the expectation value of the energy following the imaginary-time
evolution starting from a fully polarized (in the x direction)
initial state. The initial state corresponds to an empty network,
where all the DBM parameters are set to zero. The DBM results closely
match the exact imaginary-time evolution, thus verifying the correctness
of our construction. 

Numerical results for the one-dimensional Heisenberg
model are shown in Figure <ref> (b-c). Specifically,
<ref>(b) shows the numerical check for the DBM (construction
2d-6h) time evolution for one-dimensional Heisenberg model for N=16.
As expected, the DBM results also in this case follow the exact time
evolution. Figure <ref>(c) shows the dependence of the
energy from the initial state, for N=80 case. Specifically, by
taking a pre-optimized variational RBM as an initial state as an initial
state, we can significantly decrease the time τ needed to reach
the ground state.

In the case of the TFI model, sampling from the DBM is realized through
the Gibbs scheme previously sketched, in conjunction with a parallel
tempering scheme, to improve ergodicity in the sampling. In the AFH
model with 2d-6h representation, we employ loop update <cit.>
used in the path-integral QMC method, because the imaginary-time evolution
in the 2d-6h representation has a direct correspondence to the path-integral
formulation, allowing for an efficient handling of the constraint
d_[l]+d_[m]=σ_l^z+σ_m^z.



§ DISCUSSION


We have shown how exact ground states of interacting spin Hamiltonians
can be explicitly constructed using artificial neural networks comprising
only two layers of hidden variables. In contrast to approaches based
on one-layer RBMs, the constructions we have derived here do not require
further variational optimization of the network parameters. In the
case of the Heisenberg model, all of the explicit algorithms presented
here give rise to sign-problem-free representations, if the lattice
is bipartite. The DBM representation has an intrinsic conceptual value,
as an alternative to the path-integral representation. Notably, the
additional deep hidden layer in the DBM plays a similar role as an
additional dimension in statistical mechanics. Whereas a single layer
(RBM) is enough to describe exactly the state of a classical system
[see Eq. (<ref>)], a second layer
is necessary to describe exactly quantum mechanical states. 

DBM-based
schemes can be further used to systematically improve upon existing
RBM variational results. More generally, the initial state for the
present DBM scheme can be generic variational states or even combinations
of RBMs and more conventional wave functions <cit.>.
We have shown that, by starting the DBM construction from a pre-optimized
variational state, a fast convergence to the exact ground state is
observed. In conjunction with very accurate initial RBM states, this
kind of scheme opens the possibility of characterizing the ground
state even in the case of non-bipartite lattices with frustration
effects, exploiting the transient regime in which the sign problem
can be still efficiently handled numerically, as for example discussed
in Ref. <cit.>.



§ METHODS




 §.§ Useful identities


It is useful to introduce several identities, which can be
used when more complicated interactions between the visible spins
σ^z, hidden variables h and deep variables d beyond
the standard form Eq. (<ref>) are needed. The first identity
reads 

    e^s_1s_2V=C∑_s_3=±1e^s_1s_3Ṽ_1+s_2s_3Ṽ_2=2Ccosh(s_1Ṽ_1+s_2Ṽ_2).

with 

    C     =    1/2e^-|V|
    Ṽ_1    =    1/2 arcosh(e^2|V|)
    Ṽ_2    =     sgn(V)×Ṽ_1

for Ising variables s_1 and s_2, and a real interaction
V. This is a gadget for decomposing two-body interactions, and
can be proven by examining all the cases of s_1 and s_2.

By taking s_1 and s_2 as visible (physical) variables
σ^z and s_3 as a hidden variable h, the direct classical
two-body interaction between physical variables [the leftmost part
in Eq. (<ref>)] is cut and instead mediated by the hidden
neuron h. Furthermore, a direct interaction between σ^z
and d can also be decomposed: In the following derivations for
the DBM wave constructions, for convenience, we sometimes introduce
the direct interaction between σ^z and d, which is not
allowed in the DBM structure. However, by taking s_1 as a visible
spin σ^z, s_2 as a deep variable d, and s_3
as a hidden variable h in Eq. (<ref>), one can eliminate
the direct interaction between σ^z and d and decompose
it into the interaction mediated only by h with trade-off of the
summation over the hidden variable h. With this trick, one can
recover the standard DBM form in Eq. (<ref>).

Another identity (decomposition of four-body interaction)
is 

    e^s_1s_2s_3s_4V    =    1/4∑_s_5,s_6,s_7exp[iπ/4(s_5+s_6)(s_1+s_2+s_3+s_7)]
       ×   exp(s_4s_7V)
        =    ∑_s_7cos^2[π/4(s_1+s_2+s_3+s_7)]exp(s_4s_7V)

for Ising variables s_i with i=1,⋯,4. Although we have
introduced complex couplings in the first line, each term in the summation
in the second line of Eq. (<ref>) is positive definite if V
is real. The second line remains nonzero only if s_1s_2=s_3s_7,
which proves the identity. This identity with s_1 and s_2
as physical variables, s_4, s_5, and s_6 as hidden
variables, and s_3 and s_7 as deep variables, reads 

    e^σ_1σ_2d_1h_1V    =    1/4∑_h_2,h_3,d_2exp[iπ/4(h_2+h_3)(σ_1+σ_2+d_1+d_2)]
          ×exp(h_1d_2V),

Note that the right hand side fits the DBM structure.

General three-body and two-body interactions can also be
represented by the two-body form just by putting some of s_1,⋯ s_4
as constants in Eq.(<ref>). These could be used instead of Eq.
(<ref>), although we employ Eq. (<ref>) in the formalism
below for the decoupling of the two-body interaction.

Finally, we discuss the gadgets for decomposing general N-body
classical interactions using complex bias term b_j in addition
to the couplings W and W', whereas the gadgets Eqs. (<ref>)
and (<ref>) are represented only by W and W'
interactions. The gadget reads 

    e^σ_1σ_2…σ_NV    =     Ccos^2(b+π/4∑_i=1^Nσ_i)
        =    C/4∑_h_1,h_2e^ib(h_1+h_2)e^iπ/4(h_1+h_2)(σ_1+σ_2+…+σ_N)

with 

    b     =    arctan(e^-V)-π/4mod(N,4),
    
    C     =    1/cos(arctan(e^-V))×sin(arctan(e^-V)).

This fact suggests that any classical partition function defined for
Ising spins can be written exactly in terms of an RBM. Although the
RBM is shown to be powerful in representing also the quantum states,
there is no analytical way to map quantum states to the RBM and one
must rely on numerical optimizations to get the RBM parameters. In
the present study, we show analytical mappings from quantum states
to the DBM, which has additional hidden layer. In the statistical
mechanics, it is known that quantum systems with D dimension can
be mapped on (D+1)-dimensional classical systems. Therefore, having
additional hidden layer in neural network language is equivalent to
acquiring additional dimension in statistical mechanics.



 §.§ Transverse-Field Ising model


The solution of Eq.(<ref>) is found in the
following way. The left hand side of Eq.(<ref>) can
be rewritten by using the notation Eq.(<ref>) as 

    ∑_{h,d}P_1(σ^z,h)P_2(h,d)[1+tanh(Γ_lδ_τ)e^-2σ_l^z∑_jh_jW_lj]
        =     CΨ_𝒲̅(σ^z).

We look for a solution by adding one deep neuron d_[l] and creating
new couplings W_j[l]^' to the existing hidden neurons
h_j which are connected to σ_l^z. We also allow for
changes in the existing interaction parameters. In particular we set
the new couplings to be W̅_lj=W_lj+Δ W_lj, (with
Δ W_lj to be determined). Moreover, we introduce one hidden
neuron h_[l] coupled to σ_l^z and d_[l] through
the interactions W_l[l] and W'_[l][l], respectively. If
we trace out h_[l], the hidden neuron h_[l] mediates the
interaction between σ_l^z and d_[l] (denoted as W”_l[l]).

With this choice, we have (in the representation where h_[l]
is traced out): 

    Ψ_𝒲̅(σ^z)     =    ∑_{h,d}∑_d_[l]P_1(σ^z,h)P_2(h,d)
           e^σ_l^z∑_jΔ W_ljh_j+d_[l]∑_jh_jW_j[l]^'+σ_l^zd_[l]W_l[l]^''.

The equations to be verified are obtained considering the two possible
values of σ_l^z=±1: 

    e^∑_jh_j(Δ W_lj+W_j[l]^')+W_l[l]^''+e^∑_jh_j(Δ W_lj-W_j[l]^')-W_l[l]^''
        =     C×(1+tanh(Γ_lδ_τ)e^-2∑_jh_jW_lj)
           e^∑_jh_j(-Δ W_lj+W_j[l]^')-W_l[l]^''+e^∑_jh_j(-Δ W_lj-W_j[l]^')+W_l[l]^''
        =     C×(1+tanh(Γ_lδ_τ)e^2∑_jh_jW_lj).

This equation has a solution from the requirement that the hidden
unit interactions on the left and right hand sides match, thus we
require 

    Δ W_lj+W_j[l]^'    =     -2W_lj
    Δ W_lj-W_j[l]^'    =     0,

and 

    W_l[l]^''    =    logtanh(Γ_lδ_τ)/2.

Notice that when Γ_l>0, W_l[l]^'' is also
real. By using Eq. (<ref>) with the following replacement s_1→σ_l^z,
s_2→ d_[l], s_3→ h_[l], V→ W_l[l]^”,
Ṽ_1→ W_l[l] and Ṽ_2→ W_[l][l]^',
the last condition determines the real couplings W_l[l] and W_[l][l]^'
as Eqs.(<ref>)
and (<ref>).



 §.§ Heisenberg model


Here, we show the derivation for the general form of bond
Hamiltonian allowing anisotropy and bond-disorder: ℋ_lm^bond=J_lm^xy(σ_l^xσ_m^x+σ_l^yσ_m^y)+J_lm^zσ_l^zσ_m^z.
In the case of the bipartite lattice and the antiferromagnetic exchange
J_lm^z,J_lm^xy>0, we further apply a local gauge transformation
by a π rotation around the z axis in the spin space as σ^x→-σ^x
and σ^y→-σ^y on one of the sublattices,
which gives a - sign for σ_l^xσ_m^x and σ_l^yσ_m^y
interactions. This transformation is equivalent to taking 

    J_lm^xy→-J_lm^xy.

The gauge transformation enables to design a DBM neural network with
real couplings {W,W'} except for those to put “constraint`
on the values of deep neuron spins (see more detail about the constraint
in the following sections). It ensures that the DBM algorithm has
no negative sign problems.

In the case of the antiferromagnetic Heisenberg model after
the gauge transformation on the bipartite lattice, we must solve,
for each bond, 

    δ_σ_l^z,σ_m^ze^-δ_τJ_lm^zΨ_𝒲(σ^z)+(1-δ_σ_l^z,σ_m^z)e^δ_τJ_lm^z
    (Ψ_𝒲(σ^z)cosh(2J_lm^xyδ_τ)+Ψ_𝒲(σ_l^z↔σ_m^z)sinh(2J_lm^xyδ_τ))
    
    =C⟨σ^z|Ψ_𝒲̅⟩.

It is also useful to explicitly write the expression for the exchange
term in the second line above: 

    Ψ_𝒲(σ^z)cosh(2J_lm^xyδ_τ)+Ψ_𝒲(σ_l^z↔σ_m^z)sinh(2J_lm^xyδ_τ)
    
    =∑_{h,d}P_1(σ^z,h)P_2(h,d)[cosh(2J_lm^xyδ_τ)+
    sinh(2J_lm^xyδ_τ)e^(σ_m^z-σ_l^z)∑_jh_j(W_lj-W_mj)].

In the following derivations, for the antiferromagnetic Hamiltonian
(J_lm^z,J_lm^xy>0) after the gauge transformation, we look
for a solution with zero bias terms (a_i, b_j, b'_k=0,
∀ i,j,k). We can also derive a sign-problem free solution
for the imaginary time evolution in the absence of the explicit gauge
transformation by introducing a complex bias term a_i. Indeed,
in the “2 deep, 4 hidden` representation,
we will explicitly show that taking a specific set of complex bias
term a_i on physical spins is equivalent to the gauge transformation,
making a solution free from the sign problem.

In a way similar to the TFI model, solutions of Eq. (<ref>)
can be found by specifying the structure of the deep Boltzmann machine
and the three examples are the following.



  §.§.§ 1d-3h construction


We assume the structure of the updated wave function (corresponding
to Eq. (<ref>) for the TFI model) to be 

    Ψ_𝒲̅(σ^z)=∑_{h,d}∑_d_[lm]=±1
    
    d_[lm]=σ_l^z  if σ_l^z=σ_m^zP_1(σ^z,h)P_2(h,d)
    
    e^σ_l^z∑_jΔ W_ljh_j+d_[lm]∑_jh_jW_j[lm]^'+d_[lm]σ_l^zW_l[lm]^''+V_[lm]σ_l^zσ_m^z.

Similarly to the case of the TFI model, a solution of Eq. (<ref>)
is given by 

    Δ W_lj    =     -W_lj+W_mj
    
    W_j[lm]^'    =     W_lj-W_mj.

and 

    W_l[lm]^''    =     -(logtanh(2J_lm^xyδ_τ))/2
    
    V_[lm]    =     -(logcosh(2J_lm^xyδ_τ))/2-J_lm^zδ_τ

Notice that the first condition is equivalent to cutting all connections
from spin l to the hidden units and attaching the spin l to
all the hidden units connected to spin m, with an interaction W_mj.

Although the terms proportional to W_l[lm]^” and
V_lm do not satisfy the standard DBM form, they can be transformed
to the DBM form by introducing new hidden neurons h_[lm1] and
h_[lm2] [see the gadget Eq. (<ref>)]: 

    e^σ_l^zd_[lm]W_l[lm]^''=C_[lm1]∑_h_[lm1]e^σ_l^zh_[lm1]W_l[lm1]+h_[lm1]d_[lm]W'_[lm1][lm],

with 

    W_l[lm1]=W'_[lm1][lm]=1/2 arcosh(1/tanh(2J_lm^xyδ_τ)).

Similarly, the coupling V_[lm] is decomposed as 

    e^σ_l^zσ_m^zV_[lm]=C_[lm2]∑_h_[lm2]e^σ_l^zh_[lm2]W_l[lm2]+σ_m^zh_[lm2]W_m[lm2],

with 

    W_l[lm2]=-W_m[lm2]=1/2 arcosh(cosh(2J_lm^xyδ_τ)e^2J_lm^zδ_τ).

Finally, as discussed in the main text, the constraint d_[lm]=σ_l^z
when σ_l^z=σ_m^z can be satisfied by adding the
third neuron h_[lm3], introducing pure complex iπ/6 couplings.



  §.§.§ 2d-6h construction


In this case, the form of the new wave function reads

    Ψ_𝒲̅(σ^z)=∑_{h,d}∑_d_[l],d_[m]
    
    d_[l]+d_[m]=σ_l^z+σ_m^zP_1(σ^z,h)P_2(h,d)
    
    e^∑_j∑_n=l,mh_j(Δ W_njσ_n^z+W_j[n]^'d_[n])
    +∑_n=l,mσ_n^z(W”_n[l]d_[l]+W”_n[m]d_[m]).

A solution of Eq. (<ref>) is given by 

    W_j[l]^'    =     W_lj,
    
    W_j[m]^'    =     W_mj,
    Δ W_lj    =     -W_lj,
    Δ W_mj    =     -W_mj,

and 

    W_l[l]^''=W_m[m]^''    =     -J_lm^zδ_τ/2-1/4logsinh(2J_lm^xyδ_τ),
    
    W_l[m]^''=W_m[l]^''    =     -J_lm^zδ_τ/2-1/4logcosh(2J_lm^xyδ_τ).



The direct interactions between (σ_l^z,d_[l]),
(σ_m^z,d_[m]), (σ_l^z,d_[m]), and (σ_m^z,d_[l]),
are mediated by h_[lm1], h_[lm2], h_[lm3], and h_[lm4],
respectively, as follows 

    e^σ_l^zd_[l]W_l[l]^''    =     C_[lm1]∑_h_[lm1]e^σ_l^zh_[lm1]W_l[lm1]+h_[lm1]d_[l]W'_[lm1][l],
    
    e^σ_m^zd_[m]W_m[m]^''    =     C_[lm2]∑_h_[lm2]e^σ_m^zh_[lm2]W_m[lm2]+h_[lm2]d_[m]W'_[lm2][m],
    
    e^σ_l^zd_[m]W_l[m]^''    =     C_[lm3]∑_h_[lm3]e^σ_l^zh_[lm3]W_l[lm3]+h_[lm3]d_[m]W'_[lm3][m],
    
    e^σ_m^zd_[l]W_m[l]^''    =     C_[lm4]∑_h_[lm4]e^σ_m^zh_[lm4]W_m[lm4]+h_[lm4]d_[l]W'_[lm4][l].

By applying the gadget Eq. (<ref>), the new W and W'
interactions are given by, for small δ_τ (such that e^-J_lm^zδ_τ/√(sinh(2J_lm^xyδ_τ))>1),

    W_l[lm1]=W'_[lm1][l]=W_m[lm2]=W'_[lm2][m]
    
    =1/2arcosh(e^-J_lm^zδ_τ/√(sinh(2J_lm^xyδ_τ)))

and 

    W_l[lm3]=-W'_[lm3][m]=W_m[lm4]=-W'_[lm4][l]
    
    =1/2arcosh(√(cosh(2J_lm^xyδ_τ))× e^J_lm^zδ_τ).

Finally, the constraint d_[l]+d_[m]=σ_l^z+σ_m^z
can be put by introducing additionally two hidden neurons h_[lm5]
and h_[lm6], and by introducing complex couplings 

    ∑_h_[lm5],h_[lm6]e^iπ/4((σ_l^z+σ_m^z)h_[lm5]-h_[lm5](d_[l]+d_[m]))
    × e^iπ/8((σ_l^z+σ_m^z)h_[lm6]-h_[lm6](d_[l]+d_[m]))

This term gives interactions among d_[l], d_[m], σ_l^z
and σ_m^z: 4cos(π/4(σ_l^z+σ_m^z-d_[l]-d_[m]))cos(π/8(σ_l^z+σ_m^z-d_[l]-d_[m])),
which realize the constraint.



  §.§.§ 2d-4h construction


For this construction, we assume the following structure
for the wave-function after the propagator: 

    Ψ_𝒲̅(σ^z)=∑_{h,d}∑_d_[l]P_1(σ^z,h)P_2(h,d)e^∑_j,n=l,mσ_n^zh_jΔ W_nj
    × e^∑_jh_jd_[l]W_j[l]^'+∑_n=l,mσ_n^zd_[l]W_n[l]^”+∑_jσ_l^zσ_m^zh_jd_[l]Z_lmj.

In this case, we also look for a solution for the bond operator without
the gauge transformation. This shows that the introduction of a complex
bias term a_i can play the same role as the gauge transformation.
Then, we need to solve: 

    δ_σ_l^z,σ_m^ze^-δ_τJ_lm^zΨ_𝒲(σ^z)+(1-δ_σ_l^z,σ_m^z)e^δ_τJ_lm^z
    (Ψ_𝒲(σ^z)cosh(2J_lm^xyδ_τ)-Ψ_𝒲(σ_l^z↔σ_m^z)sinh(2J_lm^xyδ_τ))
    
    =C⟨σ^z|Ψ_𝒲̅⟩.

Note that the sign for Ψ_𝒲(σ_l^z↔σ_m^z)sinh(2J_lm^xyδ_τ)
term is different from that in Eq. (<ref>).

A solution of Eq. (<ref>) is obtained
as 

    Δ W_lj    =     -Δ W_mj=-1/2(W_lj-W_mj),

where W_nj (n=l,m) is updated to W̅_nj with the
increment Δ W_nj as W̅_nj=W_nj+Δ W_nj.
The new couplings W_j[l]^', Z_lmj and W”_n[l]
are also given by 

    W_j[l]^'=-Z_lmj=-1/2(W_lj-W_mj)

and 

    W”_l[l]    =    1/4[log[-e^-2a_l-mtanh(2J_lm^xyδ_τ)]
        +     2 arcosh[e^-2J_lm^zδ_τ/√(-2e^-2a_l-msinh(4J_lm^xyδ_τ))]]
    
    W”_m[l]    =    1/4[-log[-e^
    -2a_l-mtanh(2J_lm^xyδ_τ)]
        +     2 arcosh[e^-2J_lm^zδ_τ/√(-2e^-2a_l-msinh(4J_lm^xyδ_τ))]]

with a_l-m=a_l-a_m. On a bipartite lattice, to avoid the
negative sign (or complex phase) problem we need to keep W”_l[l]
and W”_m[l] real. This can be achieved by choosing a_l=0
for any l if J_lm<0 (ferromagnetic case). For J_lm>0
(antiferromagnetic case), a_l=nπ i with an arbitrary integer
n if the site l belongs to the sub-lattice A and a_l=(n+1/2)π i
if l belongs to the sub-lattice B. This local gauge for J_lm>0
is equivalent to the transformation J_lm^xy→-J_lm^xy
and a_l=0 for any site l. We further notice that W”_m[l]
can be taken positive if we take a sufficiently small δ_τ
in Eq (<ref>), with the leading order term -log(2J_lm^xyδ_τ)/2.
On the other hand, in Eq. (<ref>), the leading order term is
negative (=-J_lmδ_τ).

To recover the original form of the DBM, we first use Eq.
(<ref>) with the replacement s_1→σ_n^z,
s_2→ d_[l], s_3→ h_[n], C→ D_n,
V→ W”_n[l] Ṽ_1→ W_n[n] and
Ṽ_2→ W'_[n][l] for n=l,m. Then a solution
for D_n, W_n[n], and W'_[n][l] are represented by using
W_n[l]^” as 

    D_n    =    1/2exp[-W_n[l]^”]
    
    W_n[n]    =     W'_[n][l]=1/2 arcosh(exp[2W_n[l]^”]),

for positive W_n[l]^” and 

    D_n    =    1/2exp[W_n[l]^”]
    
    W_n[n]    =     -W'_[n][l]=1/2 arcosh(exp[-2W_n[l]^”]),

for negative W_n[l]^” to give real W_n[n] and W'_[n][l].

To completely recover the original DBM form, we next use
Eq. (<ref>) by replacing σ_1 with σ_l^z,
σ_2 with σ_m^z, d_1 with d_[l], d_2
with d_[lm], h_1 with h_j, h_2 with h_[lm1],
h_3 with h_[lm2], and V with Z_lmj.

With these solutions, by ignoring the trivial constant factors
including D_l and D_m, the evolution is described by introducing
two deep and four hidden additional variables d_[l], d_[lm],
h_[l], h_[m], h_[lm1], and h_[lm2] as 

    Ψ_𝒲̅(σ^z)=∑_{h̅,d̅}P_1(σ^z,h)P_2(h,d)exp[∑_j,n=l,mσ_n^zh_jΔ W_nj
    
    +∑_jh_jd_[l]W_j[l]^'+∑_n=l,mh_[n](σ_n^zW_n[n]+d_[l]W_[n][l]^')
    
    +d_[lm]∑_jh_jZ_lmj+iπ/4(h_[lm1]+h_[lm2])(σ_l^z+σ_m^z+d_[l]+d_[lm])],

where {h̅,d̅} is a set consisting of the existing and
new neurons.




§ ACKNOWLEDGEMENTS


G.C. acknowledges useful discussions with Xun Gao, and Markus
Heyl. Y.N. and M.I. are grateful for the useful discussions with Youhei
Yamaji and Andrew S. Darmawan. Y.N. was financially supported by Grant-in-Aids
for Scientific Research (JSPS KAKENHI) (No. 17K14336). M.I. and Y.N.
were financially supported by a Grant-in- Aid for Scientific Research
(No. 16H06345) from Ministry of Education, Culture, Sports, Science
and Technology, Japan. Part of the calculations were done at Supercomputer
Center, Institute for Solid State Physics, University of Tokyo. This
work was also supported in part by MEXT as a social and scientific
priority issue (Creation of new functional devices and high-performance
materials to support next-generation industries CDMSI) to be tackled
by using post-K computer. We also thank the support provided by the
RIKEN Advanced Institute for Computational Science through the HPCI
System Research project (hp170263) supported by Ministry of Education,
Culture, Sports, Science, and Technology, Japan. 

10
url<#>1urlprefixURL 


feynman1948spacetime
authorFeynman, R. P.
titleSpace-Time Approach to Non-Relativistic
  Quantum Mechanics.
journalReviews of Modern Physics
  volume20, pages367–387
  (year1948).
<https://link.aps.org/doi/10.1103/RevModPhys.20.367>.

dyson1949thes
authorDyson, F. J.
titleThe S Matrix in Quantum Electrodynamics.
journalPhysical Review
  volume75, pages1736–1755
  (year1949).
<https://link.aps.org/doi/10.1103/PhysRev.75.1736>.

hubbard1959calculation
authorHubbard, J.
titleCalculation of Partition Functions.
journalPhysical Review Letters
  volume3, pages77–78 (year1959).
<https://link.aps.org/doi/10.1103/PhysRevLett.3.77>.

stratonovich1957ona
authorStratonovich, R. L.
titleOn a Method of Calculating Quantum
  Distribution Functions.
journalSoviet Physics Doklady
  volume2, pages416 (year1957).
<http://adsabs.harvard.edu/abs/1957SPhD....2..416S>.

abrikosov1975methods
authorAbrikosov, A. A.
titleMethods of Quantum Field Theory in
  Statistical Physics (publisherDover Publications,
  addressNew York, year1975), editionrevised
  edition edn.

binder
authorBinder, K.
titleApplications of the Monte Carlo Method in
  Statistical Physics (publisherSpringer Verlag,
  addressBerlin, year1984).

takahashi1984
authorTakahashi, M. & authorImada, M.
titleMonte carlo calculation of quantum systems.
journalJ. Phys. Soc. Jpn.
  volume53, pages963 (year1984).

takahashi1984-2
authorTakahashi, M. & authorImada, M.
titleMonte carlo calculation of quantum systems. ii.
  higher order correction.
journalJ. Phys. Soc. Jpn.
  volume53, pages3765 (year1984).

ceperley1995pathintegrals
authorCeperley, D.
titlePath-Integrals in the Theory of Condensed
  Helium.
journalReviews of Modern Physics
  volume67, pages279–355
  (year1995).

suzuki
authorSuzuki, M.
titleRelationship 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.
journalProg. Theor. Phys.
  volume56, pages1454 (year1976).

hirsch
authorHirsch, J. E., authorSugar, R.,
  authorScalapino, D. & authorBlankenbecler, R.
titleMonte carlo simulations of one-dimensional fermion
  systems.
journalPhys. Rev. B volume26,
  pages5033 (year1982).

wiese1996
authorBeard, B. & authorWiese, U.-J.
titleSimulations of discrete quantum systems in continuous
  euclidean time.
journalJ.Phys. Rev. Lett.
  volume77, pages5130 (year1996).

sandvik1999stochastic
authorSandvik, A. W.
titleStochastic series expansion method with operator-loop
  update.
journalPhysical Review B
  volume59, pagesR14157–R14160
  (year1999).
<http://arxiv.org/abs/cond-mat/9902226>.
noteArXiv: cond-mat/9902226.

prokofev2007bolddiagrammatic
authorProkof'ev, N. & authorSvistunov, B.
titleBold Diagrammatic Monte Carlo: When Sign
  Problem is Welcome.
journalPhysical Review Letters
  volume99 (year2007).
<http://arxiv.org/abs/cond-mat/0702555>.
noteArXiv: cond-mat/0702555.

feynman
authorFeynman, R. P.
titleAtomic theory of the two-fluid model of liquid
  helium.
journalPhysical Review
  volume94, pages262 (year1954).

gros
authorGros, C.
titlePhysics of projected wavefunctions.
journalAnn. Phys. volume189,
  pages53 (year1989).

kashima2001
authorKashima, T. & authorImada, M.
titlePath-integral renormalization group method for
  numerical study on ground states of strongly correlated electronic systems.
journalJ. Phys. Soc. Jpn.
  volume70, pages2287 (year2001).

tahara
authorTahara, D. & authorImada, M.
titleVariational monte carlo method combined with
  quantum-number projection and multi-variable optimization.
journalJ. Phys. Soc. Jpn.
  volume77, pages114701
  (year2008).

becca2017quantum
authorBecca, F. & authorSorella, S.
titleQuantum Monte Carlo Approaches for
  Correlated Systems (publisherCambridge University Press,
  addressCambridge, United Kingdom ; New York, NY,
  year2017).

white
authorWhite, S. R.
titleDensity-matrix algorithms for quantum renormalization
  groups.
journalPhysical Review B
  volume48, pages10345 (year1993).

orus
authorOrús, R.
titleA practical introduction to tensor networks: Matrix
  product states and projected entangled pair states.
journalAnn. Phys. volume349,
  pages117 (year2014).

carleo2017solving
authorCarleo, G. & authorTroyer, M.
titleSolving the quantum many-body problem with artificial
  neural networks.
journalScience volume355,
  pages602–606 (year2017).
<http://science.sciencemag.org/content/355/6325/602>.

torlai2017manybody
authorTorlai, G. et al.
titleMany-body quantum state tomography with neural
  networks.
journalarXiv:1703.05334  (year2017).
<http://arxiv.org/abs/1703.05334>.
noteArXiv: 1703.05334.

nomura2017restrictedboltzmannmachine
authorNomura, Y., authorDarmawan, A. S.,
  authorYamaji, Y. & authorImada, M.
titleRestricted boltzmann machine learning for solving
  strongly correlated quantum systems.
journalPhys. Rev. B volume96,
  pages205152 (year2017).

deng2017quantum
authorDeng, D.-L., authorLi, X. &
  authorDas Sarma, S.
titleQuantum Entanglement in Neural Network
  States.
journalPhysical Review X
  volume7, pages021021 (year2017).
<https://link.aps.org/doi/10.1103/PhysRevX.7.021021>.

rocchetto2017learning
authorRocchetto, A., authorGrant, E.,
  authorStrelchuk, S., authorCarleo, G. &
  authorSeverini, S.
titleLearning hard quantum distributions with variational
  autoencoders.
journalarXiv:1710.00725 [quant-ph, stat]
  (year2017).
<http://arxiv.org/abs/1710.00725>.
noteArXiv: 1710.00725.

glasser2017neuralnetworks
authorGlasser, I., authorPancotti, N.,
  authorAugust, M., authorRodriguez, I. D. &
  authorCirac, J. I.
titleNeural Networks Quantum States, String-Bond
  States and chiral topological states.
journalarXiv:1710.04045 [cond-mat, physics:quant-ph,
  stat]  (year2017).
<http://arxiv.org/abs/1710.04045>.
noteArXiv: 1710.04045.

kaubruegger2017chiraltopological
authorKaubruegger, R., authorPastori, L. &
  authorBudich, J. C.
titleChiral Topological Phases from Artificial
  Neural Networks.
journalarXiv:1710.04713 [cond-mat,
  physics:quant-ph]  (year2017).
<http://arxiv.org/abs/1710.04713>.
noteArXiv: 1710.04713.

cai2017approximating
authorCai, Z.
titleApproximating quantum many-body wave-functions using
  artificial neural networks.
journalarXiv:1704.05148 [cond-mat]
  (year2017).
<http://arxiv.org/abs/1704.05148>.
noteArXiv: 1704.05148.

saito2017machine
authorSaito, H. & authorKato, M.
titleMachine Learning Technique to Find Quantum
  Many-Body Ground States of Bosons on a Lattice.
journalJournal of the Physical Society of Japan
  volume87, pages014001
  (year2017).
<http://journals.jps.jp/doi/10.7566/JPSJ.87.014001>.

saito2017solving
authorSaito, H.
titleSolving the Bose–Hubbard Model with
  Machine Learning.
journalJournal of the Physical Society of Japan
  volume86, pages093001
  (year2017).
<http://journals.jps.jp/doi/10.7566/JPSJ.86.093001>.

chen2018equivalence
authorChen, J., authorCheng, S., authorXie,
  H., authorWang, L. & authorXiang, T.
titleEquivalence of restricted Boltzmann machines and
  tensor network states.
journalPhysical Review B
  volume97, pages085104
  (year2018).
<https://link.aps.org/doi/10.1103/PhysRevB.97.085104>.

clark2017unifying
authorClark, S. R.
titleUnifying Neural-network Quantum States and
  Correlator Product States via Tensor Networks.
journalarXiv:1710.03545 [cond-mat,
  physics:quant-ph]  (year2017).
<http://arxiv.org/abs/1710.03545>.
noteArXiv: 1710.03545.

deng2017machine
authorDeng, D.-L., authorLi, X. &
  authorDas Sarma, S.
titleMachine learning topological states.
journalPhysical Review B
  volume96, pages195145
  (year2017).
<https://link.aps.org/doi/10.1103/PhysRevB.96.195145>.

gao_efficient_2017
authorGao, X. & authorDuan, L.-M.
titleEfficient representation of quantum many-body states
  with deep neural networks.
journalNature Communications
  volume8, pages662 (year2017).
<https://www.nature.com/articles/s41467-017-00705-2>.

huang2017neuralnetwork
authorHuang, Y. & authorMoore, J. E.
titleNeural network representation of tensor network and
  chiral states.
journalarXiv:1701.06246  (year2017).
<http://arxiv.org/abs/1701.06246>.
noteArXiv: 1701.06246.

salakhutdinov2009deepboltzmann
authorSalakhutdinov, R. & authorHinton, G.
titleDeep Boltzmann Machines.
In booktitlePMLR, pages448–455
  (year2009).
<http://proceedings.mlr.press/v5/salakhutdinov09a.html>.

trotter
authorTrotter, H. F.
titleOn the product of semi-groups of operators.
journalProc. Amer. Math. Soc.
  volume10, pages545 (year1959).

salakhutdinov2012anefficient
authorSalakhutdinov, R. & authorHinton, G.
titleAn efficient learning procedure for deep Boltzmann
  machines.
journalNeural Computation
  volume24, pages1967–2006
  (year2012).

PhysRevLett.70.875
authorEvertz, H. G., authorLana, G. &
  authorMarcu, M.
titleCluster algorithm for vertex models.
journalPhys. Rev. Lett.
  volume70, pages875–879
  (year1993).

ceperley-alder
authorCeperley, D. M. & authorAlder, J.
titleGround state of the electron gas by a stochastic
  method.
journalPhys. Rev. Lett.
  volume45, pages566 (year1980).

doi:10.1143/PTP.58.1377
authorSuzuki, M., authorMiyashita, S. &
  authorKuroda, A.
titleMonte carlo simulation of quantum spin systems. i.
journalProg. Theor. Phys.
  volume58, pages1377–1387
  (year1977).




1






PART: 
*Supplementary Information








§ DEEP BOLTZMANN MACHINES


Deep Boltzmann machine representation of quantum states.
In the main text we have considered a representation of the many-body
wave-function in terms of a two-layers deep Boltzmann Machine (DBM).
In the following we specialize to the case of N spin 1/2 particles,
described by the quantum numbers |σ^z ⟩=|σ_1^z…σ_N^z⟩
with σ_i^z=±1.
Then, we represent the amplitudes ⟨σ_1^z…σ_N^z|Ψ⟩≡Ψ(σ^z)
in the DBM form:

    Ψ_𝒲(σ^z)     =    ∑_{h}e^∑_ia_iσ_i^ze^∑_ijσ_i^zh_jW_ij+∑_jb_jh_j∑_{d}e^∑_jkh_jd_kW_jk^'+∑_kb'_kd_k.

Here, we have introduced M hidden units h_j, M^'
deep units d_k, and a set of couplings and bias terms 𝒲≡(a,b,b',W,W^').
Hereafter, we call the neurons in the 1st hidden layer just hidden neurons
 and distinguish them from the neurons in the 2nd hidden layer, which are called deep neurons.

All those parameters, in general, must be taken complex-valued to
represent a generic many-body state. The hidden and deep units are
taken here to be of spin 1/2, i.e. h_j=±1, d_k=±1,
and the summations are over all the possible values of those variables.
From a pictorial point of view, the DBM architecture features direct
connections (interactions) between nearest-neighboring layers. In
particular, the visible layer of physical degrees of freedom (σ_1^z…σ_N^z)
is connected only to the first layer of hidden variables (h_1… h_M),
whereas the first layer is connected both to the visible spins and
to the deep spins (d_1… d_M^'). 

For the following derivations, it is useful to write the DBM amplitudes
as: 

    Ψ_𝒲(σ^z)     =    ∑_{h,d}P_1(σ^z,h)P_2(h,d),

where we have introduced the two quantities: 

    P_1(σ^z,h)     =     e^∑_ia_iσ_i^ze^∑_ijσ_i^zh_jW_ij+∑_jb_jh_j
    
    P_2(h,d)     =     e^∑_jkh_jd_kW_jk^'+∑_kb'_kd_k.

Notice that, in general, those weights are complex-valued, and cannot
be interpreted as genuine Boltzmann weights. From these expressions,
it is also straightforward to see that the Restricted Boltzmann Machine
(RBM) expression for the wave-function is recovered when M^'=0,
i.e. taking 

    Ψ_𝒲^RBM(σ^z)     =    ∑_{h}P_1(σ^z,h)
        =     e^∑_ia_iσ_i^zΠ_j^M2cosh(∑_i^Nσ_i^zh_jW_ij+b_j),

where we have explicitly performed the summation of the hidden variables.
At variance with the RBM case, in the more general case when M^'>0, it is not possible to analytically obtain the DBM amplitudes.




Useful gadgets in constructing DBM neural network.
In the Methods we have discussed several useful identities to decompose spin interactions. In particular, those identities are very useful if we need more complicated interactions between the visible spins σ^z, hidden variables h and deep variables d beyond the standard form Eq. (<ref>).
For the sake of completeness of this Supplementary Information, we reproduce here the identities for decomposing two-body, three-body, and four-body interactions. 


The first identity reads 

    e^s_1 s_2 V     =  C∑_s_3  = ±1  e^s_1  s_3  Ṽ_1  + s_2 s_3 Ṽ_2 = 2C cosh (s_1 Ṽ_1 + s_2Ṽ_2   ) .

with 

    C    =   1/2 e^ - | V | 
    Ṽ_1    =   1/2 arcosh ( e^2|V|) 
    Ṽ_2    =    sgn (V ) ×Ṽ_1

for Ising variables s_1 and s_2, and a real interaction V.
This is the gadget for decomposing two-body interactions discussed in Methods. 
In the following, we will use this identity to decompose either interactions between visible (physical spins) (in that case s_1 and s_2 are both σ^z variables), or to decompose direct interactions between a σ^z spin and a deep unit d. 


Another identity (decomposition of four-body interaction) is

    e^s_1s_2s_3 s_4 V   =   1/4∑_s_5,s_6,s_7exp[i π/4( s_5+s_6 )(s_1+s_2+s_3 + s_7)  ] exp  (s_4 s_7 V) 
     
       =   ∑_s_7cos^2  [π/4(s_1+s_2+s_3+s_7)  ]  exp  (s_4 s_7 V)

for Ising variables s_i with i=1, ⋯, 4.
Although we have introduced complex couplings in the first line, 
each term in the summation in the second line of Eq. (<ref>) is positive definite if V is real. 
The second line remains nonzero only for 
s_7=1 if s_1s_2s_3=1 and only for s_7=-1 if s_1s_2s_3=-1, 
which proves the identity. 
This identity with s_1 and s_2 as physical variables, s_4, s_5, and s_6 as hidden variables, and s_3 and s_7 as deep variables, which reads

    e^σ_1 σ_2 d_1 h_1 V   =   1/4∑_h_2,h_3, d_2exp[i π/4( h_2+h_3 )(σ_1+σ_2+ d_1 + d_2)  ] exp  (h_1 d_2 V),

will be used in Sec. <ref>. 
Note that the right hand side fits the DBM structure. 


Although identities for decomposing three-body interactions
are not used in the following derivation, it is nonetheless useful to show them: 

    e^s_1s_2s_3 V   =   1/4∑_s_4,s_5,s_6exp[i π/4( s_4+s_5 )(s_1+s_2+s_3 + s_6)  ] exp  ( s_6 V) 
     
       =   ∑_s_6cos^2  [π/4(s_1+s_2+s_3+s_6)  ]  exp  (s_6 V).

This gadget for three-body interactions is obtained by fixing s_4=1 in Eq. (<ref>) (and changing variables).
 Alternative form is obtained by replacing s_3 with 1 in 
Eq. (<ref>), which gives, 


    e ^ s_1 s_2 s_3 V   =   1/4∑_s_4, s_5, s_6exp [ i π/4  (s_4 + s_5)    ( s_1 + s_2 + s_6 +1 )  ]    exp (  s_3 s_6 V )  
       
       =   ∑_s_6cos^2  [ π/4 ( s_1 + s_2 + s_6 +1 )  ]    exp(  s_3 s_6 V ).

As we see, the gadgets for three-body interactions [Eqs. (<ref>) and (<ref>)] have been derived from the gadget for four-body interactions [Eq. (<ref>)] trivially.  

Gadgets for two-body interactions which are different from Eq. (<ref>) can also be obtained from Eq. (<ref>)
by fixing two variables out of s_1, s_2, s_3, s_4 to be 1. 
These could be used instead of (<ref>), 
although we employ (<ref>) in the formalism below for the decoupling of the two-body interaction.




§ REPRESENTING GROUND-STATES


As discussed in the main text, our goal is to construct explicit DBM
representations of ground-states of local Hamiltonians. This goal
is achieved 
by finding a representation of the imaginary-time evolved
state: 

    |Ψ(τ)⟩    =     e^-τℋ|Ψ_0⟩,

where |Ψ_0⟩ is empty RBM (⟨σ^z |Ψ_0⟩=const.) or 
pre-optimized RBM state, converging
to the exact ground-state for large enough τ. To achieve this
goal we first consider a second-order Trotter-Suzuki decomposition:

    |Ψ (τ) ⟩    =    𝒢_1(δ_τ/2)𝒢_2(δ_τ)…𝒢_1(δ_τ)𝒢_2(δ_τ)𝒢_1(δ_τ/2)|Ψ_0⟩,

where δ_τ is a small time step, the Hamiltonian is decomposed
into two non-commuting parts, ℋ=ℋ_1+ℋ_2,
and 𝒢_ν(δ_τ)=e^-ℋ_νδ_τ are
short-time propagators. For given Hamiltonian, we then need to find
specific rules to apply the short-time propagators to a generic DBM,
and obtain a new (time-evolved) DBM, possibly with a larger total
number of hidden and deep neurons. 
In the following, we show concrete examples for the transverse-field Ising and Heisenberg models. 



 §.§ Transverse-Field Ising model


Let us start with the case of the transverse-field Ising model. We
consider a Trotter-Suzuki decomposition of the imaginary-time propagator,
into two parts: ℋ_1=-∑_iΓ_iσ_i^x, and ℋ_2=∑_l<mV_lmσ_l^zσ_m^z.
In the following derivation, we assume that Γ_i is positive (Γ_i >0).
In this case, we look for a solution with zero bias terms: a_i = b_j = b'_k = 0, ∀ i,j,k. 
The case of negative Γ_i can also be treated, and is discussed more in detail at the end of this section. 





Interaction propagator.
The interaction propagator e^-δ_τV_lmσ_l^zσ_m^z

is diagonal in the σ^z basis, and applying it to

 a DBM will lead to a modification in the
DBM parameters. 
In particular, the goal is to satisfy the equation:

    ⟨σ^z  |e^-δ_τV_lmσ_l^zσ_m^z|Ψ_𝒲⟩    =     C⟨σ^z |Ψ_𝒲̅⟩,

i.e. to explicitly find a set of parameters 𝒲̅ that
satisfies the previous equation for all the possible ⟨σ^z |,
and for an arbitrary constant C. 

We can achieve this goal adding a hidden unit in the first layer,
h_[lm] such that it is only connected to the visible spins: W_[lm]k^'=0, ∀ k.
The new wave function has then the form: 

    Ψ_𝒲̅(σ^z)     =    ∑_{h,d}∑_ h_[lm] P_1(σ^z,h)P_2(h,d)e^σ_l^z W_l[lm]h_[lm]+σ_m^z W_m[lm]h_[lm]
        =     2cosh(σ_l^z W_l[lm]+σ_m^z W_m[lm])Ψ_𝒲(σ^z).

Equation (<ref>) is then satisfied if

    e^-δ_τV_lmσ_l^zσ_m^z = 2C
    cosh(σ_l^z W_l[lm]+σ_m^z W_m[lm])

for all the possible values of σ_l^z and σ_m^z.
By using the gadget Eq. (<ref>), the new parameters W_l[lm] and W_m[lm] are given by 

    W_l[lm]    =    1/2arcosh(e^2|V_lm|δ_τ)
    
    W_m[lm]    =     -sgn(V_lm) × W_l[lm].





Transverse-field propagator.
The propagator involving the transverse-field e^δ_τΓ_lσ_l^x
is off-diagonal in σ^z basis.
For this off-diagonal part, we must solve a slightly more involved
equation:

    ⟨σ^z |e^δ_τΓ_lσ_l^x|Ψ_𝒲⟩    =    Ψ_𝒲(σ^z )×cosh(Γ_lδ_τ)+Ψ_𝒲(σ_1^z,…-σ_l^z,…, σ_N^z)×sinh(Γ_lδ_τ)
        =     C⟨σ^z |Ψ_𝒲̅⟩,

for the new parameters 𝒲̅, and for an arbitrary finite
normalization constant C. In turn, this equation is equivalent to: 

    ∑_{h,d}P_1(σ^z,h)P_2(h,d)[1+tanh(Γ_lδ_τ)e^-2σ_l^z∑_jh_jW_lj]     =     CΨ_𝒲̅(σ^z).

We look for a solution
by adding one deep neuron d_[l] 
and creating new couplings W_j[l]^' to the existing hidden neurons h_j 
which are connected to σ_l^z. 
We also allow for 
changes in the existing interaction parameters. 
In particular we set the new couplings to be W̅_lj=W_lj+Δ W_lj,
(with Δ W_lj to be determined).
 

Moreover, we introduce one hidden neuron h_[l] coupled to σ_l^z and d_[l]
through the interactions W_l[l] and W'_[l][l], respectively. 
If we trace out h_[l], the hidden neuron h_[l] mediates the interaction between σ_l^z and d_[l] (denoted as W”_l[l]).   

With this choice, we have (in the representation where h_[l] is traced out): 

    Ψ_𝒲̅(σ^z)     =    ∑_{h,d}∑_ d_[l]P_1(σ^z,h)P_2(h,d)e^σ^z_l∑_jΔ W_ljh_j+d_[l]∑_jh_jW_j[l]^'+σ^z_ld_[l]W_l[l]^''.


The equations to be verified are obtained considering the two
possible values of σ^z_l=±1: 

    e^∑_jh_j(Δ W_lj+W_j[l]^')+W_l[l]^''+e^∑_jh_j(Δ W_lj-W_j[l]^')-W_l[l]^''    =     C×(1+tanh(Γ_lδ_τ)e^-2∑_jh_jW_lj)
    
    e^∑_jh_j(-Δ W_lj+W_j[l]^')-W_l[l]^''+e^∑_jh_j(-Δ W_lj-W_j[l]^')+W_l[l]^''    =     C×(1+tanh(Γ_lδ_τ)e^2∑_jh_jW_lj).

This equation has a solution if the hidden unit interactions on the
l.h.s. and on the r.h.s match, i.e. when: 

    Δ W_lj+W_j[l]^'    =     -2W_lj
    Δ W_lj-W_j[l]^'    =     0,

which in turn are verified when 

    W_j[l]^'    =     -W_lj
    Δ W_lj    =     -W_lj,

and if 

    W_l[l]^''    =    logtanh(Γ_lδ_τ)/2.


When Γ_l > 0, W_l[l]^'' is real. 
By using Eq. (<ref>) with the following replacement s_1 →σ^z_l, s_2  → d_[l], s_3  → h_[l], V  → W_l[l]^”, Ṽ_1  → W_l[l] and Ṽ_2  → W_[l][l]^',

the last condition determines the real couplings W_l[l] and
W_[l][l]^',  
which read 

    W_l[l]    =    1/2arcosh(1/tanh(Γ_lδ_τ))
    
    W_[l][l]^'    =     -W_l[l].

Notice that because of condition (<ref>), after applying
the off-diagonal propagator all the interactions W_lj between
spin l and hidden units h_j are set to zero. However, because
of condition (<ref>), the spin l is reconnected
to the new hidden unit h_[l] with the W_l[l] interaction. 




Negative transverse field.
When Γ_i < 0, it is still possible to recover a DBM representation with purely real interaction weights W and W^'. In order to do so, we apply the gauge transformation σ_i^x → -σ_i^x 
and σ_i^y → -σ_i^y (π spin rotation around the z axis), which maps onto the Hamiltonian with positive Γ_i. 
This gauge transformation can be achieved by taking a finite bias terms a_i in Eq. (<ref>) as a_i = iπ/2
and fix them during the imaginary time evolution. 
With this complex bias term a_i = iπ/2, | ↑⟩ (|↓⟩) state at the ith site 
acquires a phase as follows | ↑⟩→ e^iπ/2| ↑⟩ = i | ↑⟩
(| ↓⟩→ e^-iπ/2| ↓⟩ = -i | ↓⟩), 
which is equivalent to a π spin rotation around the z axis.
In the case when Γ_i is originally positive, we can set all the bias terms { a, b, b' } to be zero. 



 §.§ Heisenberg Model

 

We now consider the case of the Heisenberg model, whose Hamiltonian reads

    ℋ   =   ∑_⟨ lm ⟩ℋ_lm
    ℋ_lm   =   ℋ_lm^z+ℋ_lm^xy
    ℋ_lm^z   =      J_lm^zσ_l^zσ_m^z
    ℋ_lm^xy   =    J_lm^xy (σ_l^xσ_m^x + σ_l^yσ_m^y) =  2J_lm^xy(σ_l^+σ_m^-+σ_l^-σ_m^+)
 
with J_lm^z=J_lm^xy=J. 
We write the Hamiltonian in a general form because the following DBM algorithm can be straightforwardly extended to the more general case of
anisotropic/disordered bonds.
As a starting point for our construction, we decompose the Hamiltonian into pieces by a Trotter-Suzuki decomposition of the imaginary-time
propagator: e^-δ_τℋ∼∏_⟨ lm ⟩e^-δ_τℋ_lm+O(δ_τ^2).
Then in this Section, we represent e^-δ_τℋ_lm by using the DBM in three different forms, which are all exact. By taking δ_τ small enough and operating e^-δ_τℋ_lm many times, those constructions ensure that the ground state is obtained with any controlled accuracy.


For e^-δ_τℋ_lm, 
and the antiferromagnetic exchange J_lm^z,J_lm^xy>0,
if the lattice is bipartite, we further apply a local gauge transformation by π rotation around z axis in the spin space as

    σ^x → -σ^x
         and    σ^y → -σ^y

 on one of the sublattices, which gives a - sign for the σ_l^xσ_m^x and σ_l^yσ_m^y interactions. It is equivalent to the following transformation in the couplings:

    J_lm^xy→ -J_lm^xy.

The gauge transformation enables to design a DBM neural network with real couplings { W, W' } 
except for those necessary to enforce local constraints on the values of deep neuron spins (see more detail about the constraint in the following sections). Overall, we show in the following that the 3 different DBM constructions have no negative sign problem.  

On the bipartite lattice, the Suzuki-Trotter decomposition is frequently expressed by decomposing the Hamiltonian 
ℋ into several groups. For instance, on the one dimensional chain, if it is natural to decompose it into odd and even bonds: 

    ℋ_1=∑_⟨ l,m⟩∈ odd  bondℋ_lm
    ,    ℋ_2=∑_⟨ l,m⟩∈ even  bondℋ_lm,

further decompositions  
e^-δ_τℋ_1=∏_⟨ l,m ⟩∈  odd  bonde^-δ_τℋ_lm
and 
e^-δ_τℋ_2=∏_⟨ l,m ⟩∈  even  bonde^-δ_τℋ_lm
contain commuting elements and are therefore exact. For the square lattice, a similar procedure requires the decomposition of the Hamiltonian into 4 parts, in a checkerboard fashion.
In all cases, the fundamental ingredient to represent the ground-state as a DBM is to find an exact expression for the bond propagator, e^-δ_τℋ_lm, when applied to an existing DBM state. 

In the case of antiferromagnetic Heisenberg model after the gauge transformation on the bipartite lattice, we must solve, 
for each bond, 

    ⟨σ^z | e^δ_τJ_lm^xy(σ_l^xσ_m^x+σ_l^yσ_m^y)
    -δ_τJ_lm^zσ_l^zσ_m^z
    | Ψ_𝒲⟩
    
    = δ_σ_l^z,σ_m^ze^-δ_τJ_lm^zΨ_𝒲(σ^z) 
    +(1-δ_σ_l^z,σ_m^z)e^δ_τJ_lm^z(Ψ_𝒲(σ^z)×cosh(2J_lm^xyδ_τ)+Ψ_𝒲(σ_l^z↔σ_m^z)×sinh(2J_lm^xyδ_τ))
    
    = C⟨σ^z|Ψ_𝒲̅⟩.

It is also useful to explicitly write the expression for the exchange
term in the second line above: 

    Ψ_𝒲(σ^z )×cosh(2J_lm^xyδ_τ)+Ψ_𝒲(σ_l^z↔σ_m^z)×sinh(2J_lm^xyδ_τ)
    
    =∑_{h,d}P_1(σ^z,h)P_2(h,d)[cosh(2J_lm^xyδ_τ)+sinh(2J_lm^xyδ_τ)e^(σ_m^z-σ_l^z)∑_jh_j(W_lj-W_mj)].

In the following derivations, for the antiferromagnetic Hamilonian (J_lm^z, J_lm^xy >0)
after the gauge transformation, we look for a solution with zero bias terms (a_i,  b_j,  b'_k = 0, ∀ i, j, k). 
We can also derive a sign-problem free solution for the imaginary time evolution in the absence of the explicit gauge transformation by introducing complex bias term a_i. 
Indeed, in the “2 deep, 4 hidden" representation in Sec. <ref>, 
we will explicitly show that taking a specific set of complex bias term 
a_i on physical spins is equivalent to 
the gauge transformation, making a solution free from the sign problem. 



  §.§.§ 1 deep, 3 hidden (1d-3h) representation


Strategy.
The first representation we propose is obtained adding one deep neuron
d_[lm], 
which gives new couplings W'_j[lm] to the hidden units h_j connected to σ_l^z and σ_m^z. 
We also allow for changes
in the existing DBM parameters. 
In particular we set the new couplings to be
W̅_lj=W_lj+Δ W_lj, (with Δ W_lj to be
determined). 
We introduce a coupling W_l[lm]^''
between
σ_l^z
and 
d_[lm],  
and a coupling V_[lm]
between σ_l^z and σ_m^z, which are both not allowed in the DBM architecture. 
By using the gadget Eq. (<ref>), 
these interactions can be mediated by hidden neurons 
h_[lm1] and h_[lm2], respectively, and the DBM form is recovered. 
Furthermore, we look for a solution with a constraint: 
d_[lm]=σ^z_l  when σ^z_l=σ^z_m
(when σ^z_l≠σ^z_m, the d_[lm] value is not constrained).
Imposing the constraint on the value of the deep unit 
is a crucial difference from the DBM solution for the TFI model.
We will show that this constraint can be achieved by adding additional hidden neuron h_[lm3]
and introducing complex couplings (“iπ/6" trick).  
We discuss this trick in more detail later. 

In total, we introduce one deep and three hidden neurons. 
After tracing out the three hidden neurons h_[lm1], h_[lm2], and h_[lm3], 
the new wave function reads 

    Ψ_𝒲̅(σ^z)     =    ∑_{h,d}∑_d_[lm]  = ±1 
     d_[lm]= σ_l^z   if  σ_l^z = σ_m^z P_1(σ^z,h)P_2(h,d)e^σ^z_l∑_jΔ W_ljh_j+d_[lm]∑_jh_jW_j[lm]^'+d_[lm]σ^z_lW_l[lm]^''+V_[lm]σ^z_lσ^z_m.




Derivation for the update of parameters.
The equations to be verified are then obtained considering 
all the
possible values of σ^z_l=±1 and σ^z_m=±1, in addition
to the constraints on d_[lm] previously introduced. 
We then have two equations for 

σ^z_l = σ^z_m = ±1: 

    e^∑_j(Δ W_lj+W_j[lm]^')h_j+W_l[lm]^''+V_[lm]    =     C×exp(-J_lm^zδ_τ)
    
    e^∑_j(-Δ W_lj-W_j[lm]^')h_j+W_l[lm]^''+V_[lm]    =     C×exp(-J_lm^zδ_τ),

and the other two equations for 
σ^z_l = -σ^z_m = ±1:

    e^∑_j(Δ W_lj+W_j[lm]^')h_j+W_l[lm]^''-V_[lm]+e^∑_j(Δ W_lj-W_j[lm]^')h_j-W_l[lm]^''-V_[lm]
    
    =C×exp(J_lm^zδ_τ)(cosh(2J_lm^xyδ_τ)+sinh(2J_lm^xyδ_τ)e^-2∑_jh_j(W_lj-W_mj)),


    e^∑_j(-Δ W_lj+W_j[lm]^')h_j-W_l[lm]^''-V_[lm]+e^∑_j(-Δ W_lj-W_j[lm]^')h_j+W_l[lm]^''-V_[lm]
    
    =C×exp(J_lm^zδ_τ)(cosh(2J_lm^xyδ_τ)+sinh(2J_lm^xyδ_τ)e^2∑_jh_j(W_lj-W_mj)).

These equations have a solution if the hidden unit interactions on
the l.h.s. and on the r.h.s match, i.e. when:

    Δ W_lj+W_j[lm]^'    =     0
    Δ W_lj-W_j[lm]^'    =     -2(W_lj-W_mj)

which implies 

    Δ W_lj    =     -W_lj+W_mj
    
    W_j[lm]^'    =     W_lj-W_mj.

Notice that the first condition 
gives W̅_lj = W_lj + Δ W_lj = W_mj, 
which is equivalent to cutting all connections
from spin l to the hidden units and attaching the spin l to
all the hidden units connected to spin m, with an interaction W_mj. 

In order to match the coefficients we must also have: 

    W_l[lm]^''+V_[lm]    =    log C-J_lm^zδ_τ
    
    W_l[lm]^''-V_[lm]    =    log C+logcosh(2J_lm^xyδ_τ)+J_lm^zδ_τ
    
    -W_l[lm]^''-V_[lm]    =    log C+logsinh(2J_lm^xyδ_τ)+J_lm^zδ_τ,

which has the solution:

    W_l[lm]^''    =     -(logtanh(2J_lm^xyδ_τ))/2
    
    V_[lm]    =     -(logcosh(2J_lm^xyδ_τ))/2-J_lm^zδ_τ





Recovery of standard DBM.
The coupling W_l[lm]^'' between the deep unit d_[lm] 
and the visible spin σ_l^z is mediated by the hidden unit h_[lm1] coupled to 
σ_l^z by W_l[lm1] and 
d_[lm] by W'_[lm1][lm]: 


    exp( σ_l^z d_[lm] W_l[lm]^'' ) = 
    C_[lm1]∑_h_[lm1]exp( σ_l^z h_[lm1] W_l[lm1] + h_[lm1]  d_[lm]W'_[lm1][lm] ).

By using Eq. (<ref>) with the following replacement s_1 →σ^z_l, s_2  → d_[lm], s_3  → h_[lm1], V  → W”_l[lm], Ṽ_1  → W_l[lm1] and Ṽ_2  → W'_[lm1][lm],
W_l[lm1] and W'_[lm1][lm] are given by 

    W_l[lm1] = W'_[lm1][lm] = 1/2 arcosh(  1/tanh(2J_lm^xyδ_τ)).


Similarly, the coupling V_[lm] between visible spins σ_l^z and σ_m^z 
is mediated by the hidden unit h_[lm2] coupled to 
σ_l^z by W_l[lm2] and 
σ_m^z by W_m[lm2]: 


    exp( σ_l^z σ_m^z V_[lm] ) = 
     C_[lm2]∑_h_[lm2]exp( σ_l^z h_[lm2] W_l[lm2] + σ_m^z h_[lm2] W_m[lm2] ).

By using Eq. (<ref>) with the following replacement s_1 →σ^z_l, s_2  →σ^z_m, s_3  → h_[lm2], V  → V_lm, Ṽ_1  → W_l[lm2] and Ṽ_2  → W_m[lm2],
W_l[lm2] and W_m[lm2] are given by 

    W_l[lm2] =  - W_m[lm2] = 1/2 arcosh( cosh(2J_lm^xyδ_τ)  e^2J_lm^z δ_τ).





How to enforce the constraint  when 
 (“" trick).
The constraint d_[lm] = σ^z_l when σ^z_l = σ^z_m
can be exactly satisfied by introducing pure complex connections. 
We can replace the sum with the constraint in Eq. (<ref>) as follows
(we ignore trivial constant factor): 

    ∑_d_[lm]  = ±1 
     d_[lm]= σ_l^z   if  σ_l^z = σ_m^z ⟶  ∑_ d_[lm]∑_h_[lm3] 
    e^  i π/6 ( (σ^z_l + σ^z_m )   h_[lm3]   -   h_[lm3]  d_[lm]) 
    = ∑_ d_[lm] 
    2 cos ( π/6 (σ^z_l + σ^z_m -  d_[lm] )    )

One can easily see that the cosine term in the rightmost part gives nonzero value only when d_[lm] = σ^z_l if σ^z_l = σ^z_m.
On the other hand, if σ^z_l ≠σ^z_m, both d_[lm]± 1 contributions survive. 










Proof of no negative sign.
Here, we show that the marginal probability density  
Π̃'(σ^z,d,d^')=∑_h,h^'Π(σ^z,h,h^',d,d^')
obtained by tracing out the hidden unit 
is non-negative definite. 
Therefore, we can perform the Metropolis sampling using Π̃' density without suffering from the negative signs
(see more detail on the sampling scheme in Sec. <ref>).
To prove this, it is sufficient to show 

    ∑_{h }  P_1 (σ^z, h) P_2 (h, d) ≥ 0.
 
for all possible σ^z and d configurations. 

In the 1d-3h representation, i π / 6 complex couplings are originally introduced to put the constraint locally. 
However, as time evolves, these complex couplings become non-local
(see Fig. <ref>).
Because the pure complex couplings give cosine terms after tracing out hidden variables, 
they have a potential to give negative signs. 
Here, we prove that this is not the case. 

We assume that Eq. (<ref>) is satisfied for all possible σ^z and d after several steps of the imaginary time evolution.
Then, we apply the bond propagator e^-ℋ_lmδ_τ to obtain the new wave function. 
In the case when σ^z_l = -σ^z_m = 1, the solution in the 1d-3h representation can be rewritten as  

    ∑_ h_[lm1] , h_[lm2] , h_[lm3]P̅_1 (σ^z,  h̅ ) P̅_2 (  h̅  , d, d_[lm] = 1)     =     P_1 (σ^z, h) P_2 (h, d ) × ( positive  constant) 
    ∑_   h_[lm1] , h_[lm2] , h_[lm3]P̅_1 (σ^z,   h̅) P̅_2 ( h̅ , d, d_[lm] = -1 )     =     P_1 (σ_l^z↔σ_m^z, h) P_2 (h, d ) ×  ( positive  constant)

where P̅_1 ×P̅_2 on the left hand side is the new weight after the imaginary time evolution, 
and {h̅} consists of the existing hidden neurons {h } and 
the newly introduced hidden neurons h_[lm1], h_[lm2], and h_[lm3]. 
By taking the summation on the existing hidden variables on both sides, we get 

    ∑_{h̅}P̅_1 (σ^z,  h̅) P̅_2 ( h̅, d, d_[lm] = 1)     =   ∑_{ h } 
         P_1 (σ^z,  h) P_2 (h, d ) × ( positive  constant) ≥  0  
    ∑_{h̅}P̅_1 (σ^z, h̅ ) P̅_2 ( h̅ , d, d_[lm] = -1)     =   ∑_{ h } P_1 (σ_l^z↔σ_m^z, h) P_2 (h, d ) × ( positive  constant)   ≥ 0

Here, we used Eq. (<ref>) to obtain the rightmost inequality. 





It proves that the new weight with the hidden variables being traced out is also non-negative. 
In the same way, we can show the non-negativeness of the new weight for σ_l = -σ_m = -1.

Next we consider the case σ^z_l = σ^z_m = 1. 
In this case, 

    ∑_   h_[lm1] , h_[lm2] , h_[lm3]P̅_1 (σ^z,  h̅) P̅_2 ( h̅, d, d_[lm] = 1)     =     P_1 (σ^z, h) P_2 (h, d ) × ( positive  constant), 
    ∑_   h_[lm1] , h_[lm2] , h_[lm3]P̅_1 (σ^z,  h̅ ) P̅_2 ( h̅, d, d_[lm] = -1)     =    0.

By taking the summation on the existing hidden variables on both sides, we obtain

    ∑_{h̅}P̅_1 (σ^z,  h̅ ) P̅_2 ( h̅ , d, d_[lm] = 1)     =   ∑_h P_1 (σ^z, h) P_2 (h, d ) × ( positive  constant) ≥  0  
    ∑_{h̅}P̅_1 (σ^z,  h̅ ) P̅_2 ( h̅ , d, d_[lm] = -1)     =     0.

Therefore, the non-negativeness of the weight is ensured. 
The proof for σ^z_l = σ^z_m = -1 case can be done in an analogous way.  

We have proven that the new weight after applying the bond propagator e^-ℋ_lmδ_τ 
is non negative for all the possible σ^z and d̅ configurations:

    ∑_{h̅}P̅_1 (σ^z,  h̅ ) P̅_2 ( h̅, d̅ ) ≥ 0

with {d̅} consisting of { d } and d_[lm]. 
It ensures the non-negativeness of the weight at any time during the imaginary time evolution.   




Summary of 1d-3h representation.
The action the bond propagator is summarized as follows. 
First, the new deep neuron d_[lm] is attached to the existing hidden neurons connected to σ_l^z and σ_m^z. 
Second, σ_l^z is disconnected to all hidden units 
and reconnected to the hidden units having finite couplings to σ_m^z 
(W̅_lj =  W_mj).
Third, four couplings are inserted, involving new hidden neurons h_[lm1] and 
h_[lm2]: σ_l^z ↔ h_[lm1],  
h_[lm1]↔ d_[l], σ_l^z ↔ h_[lm2] and σ_m^z ↔ h_[lm2],
Finally, the new hidden neuron h_[lm3] puts the constraint on the d_[l] sum by the imaginary couplings to 
σ_l^z, σ_m^z, and d_[l].
 

By successively applying the imaginary-time evolutions, the W' couplings become nonlocal or long ranged.  
On the other hand, the W couplings stay local (see Fig. <ref>). 

 



  §.§.§ 2 deep, 6 hidden (2d-6h) representation



Strategy.
We look for a solution where we add two deep neurons d_[l] and d_[m],
giving new couplings W'_j[l], W'_j[m] to the existing hidden spins h_j connected to σ_l^z and σ_m^z. 
We also allow for changes in the existing  
W parameters:  
We set the new couplings to be W̅_lj =W_lj+Δ W_lj and 
W̅_mj =W_mj+Δ W_mj
(with Δ W_lj,   Δ W_mj to be determined).



Furthermore, we add four hidden neurons h_[lm1], h_[lm2], h_[lm3], and h_[lm4] to mediate the interactions between (σ^z_l, d_[l]), (σ^z_m, d_[m]), (σ^z_l, d_[m]), and (σ^z_m, d_[l]), respectively.   
We solve the equation with the constraint σ^z_l + σ^z_m = d_[l] + d_[m]. 
This constraint can be achieved, for example, by adding two further hidden neurons (h_[lm5] and h_[lm6], respectively) and introducing complex connections  (“i π/4, iπ/8" trick).
This trick will be discussed in detail later.  

In total, we add two deep neurons (d_[l] and d_[m]) and six hidden neurons (h_[lm1], …, h_[lm6]).
In the following, to make equations simple, we employ a representation in which the new hidden neurons are analytically traced out. 
The interactions between (σ_l^z, d_[l]), (σ_m^z, d_[m]), (σ_l^z, d_[m]), and (σ_m^z, d_[l]), 
which are mediated by 1st to 4th hidden neurons, will be denoted as W”_l[l], W”_m[m], W”_l[m], and W”_m[l], respectively. 
The 5th and 6th hidden neurons filter out σ^z_l + σ^z_m ≠ d_[l] + d_[m] contributions. 
With this setting, the new wave function is represented as 

    Ψ_𝒲̅(σ^z) =  ∑_{h,d}∑_d_[l], d_[m]
     d_[l] + d_[m] = σ^z_l + σ^z_m        
    P_1(σ^z,h)P_2(h,d) e^∑_j   h_j  (  Δ W_ljσ^z_l   + W_j[l]^' d_[l] ) 
    + ∑_j   h_j  (   Δ W_mjσ^z_m  + W_j[m]^' d_[m] ) 
          × e^σ^z_l (W”_l[l] d_[l] + W”_l[m] d_[m]) 
    +   σ^z_m (W”_m[l] d_[l] +  W”_m[m] d_[m])  .





Derivation for the update of parameters.
When the lth and mth physical spins are anti-parallel (σ^z_l  = - σ^z_m = ±1), 
d_[l]  =  - d_[m]  = ± 1 contributions survive in the sum over d_[l] and d_[m] variables
in Eq. (<ref>), and thus the equations to be satisfied are


    e^∑_jh_j(Δ W_l-m,j+W_j[l]^' - W_j[m]^') +  W_l-m,[l]^'' - W_l-m,[m]^''  
     + e^∑_jh_j(Δ W_l-m,j-W_j[l]^' + W_j[m]^')-W_l-m,[l]^'' + W_l-m,[m]^''
      =   
     C e^J_lm^z δ_τ( cosh( 2 J_lm^xyδ_τ )  + sinh( 2 J_lm^xyδ_τ )  e^-2∑_jh_jW_l-m,j)

for σ^z_l  = - σ^z_m = 1 and 

    e^∑_jh_j( - Δ W_l-m,j+W_j[l]^' - W_j[m]^') -   W_l-m,[l]^'' + W_l-m,[m]^''  
     + e^∑_jh_j(- Δ W_l-m,j-W_j[l]^' + W_j[m]^') + W_l-m,[l]^'' - W_l-m,[m]^''
     
     = 
    
    C e^J_lm^z δ_τ( cosh( 2 J_lm^xyδ_τ )  + sinh( 2 J_lm^xyδ_τ )  e^2∑_jh_jW_l-m,j)

for σ^z_l  = - σ^z_m = -1, respectively. 
Here, W_l-m,[α] = W_l[α] - W_m[α],  Δ W_l-m,[α] = Δ W_l[α] - Δ W_m[α], W_l-m,[α]^'' = W_l[α]^'' - W_m[α]^'' with α = l, m.

When the lth and mth physical spins are parallel (σ^z_l  = σ^z_m = ±1), 
only d_[l] =  d_[m] = σ^z_l = σ^z_m contribution survives in the sum over d_[l] and d_[m] variables in Eq. (<ref>), and thus the equations to be satisfied are

    e^∑_jh_j(Δ W_l+m,j+W_j[l]^' + W_j[m]^') +  W_l+m,[l]^'' + W_l+m,[m]^''  =  C e^- J_lm^z  δ_τ

 for σ^z_l  = σ^z_m =  1
 
    e^∑_jh_j(-Δ W_l+m,j - W_j[l]^'  - W_j[m]^') +  W_l+m,[l]^'' + W_l+m,[m]^''
     =  C e^- J_lm^z  δ_τ

 for σ^z_l  = σ^z_m =  -1, respectively. 
 Here, W_l+m,[α] = W_l[α] + W_m[α ],  Δ W_l+m,[α] = Δ W_l[α] + Δ W_m[α], W_l+m,[α]^'' = W_l[α]^'' + W_m[α]^'' with α =l, m. 

The equations (<ref>),  (<ref>),  (<ref>), and (<ref>) are satisfied if 

    Δ W_l-m,j - W_j[l]^' + W_j[m]^'    =     -2W_l-m,j,
    Δ W_l-m,j + W_j[l]^'  - W_j[m]^'    =     0, 
    Δ W_l+m,j + W_j[l]^' + W_j[m]^'    =     0,

and  

    W_l-m,[l]^''   - W_l-m,[m]^''    =    log C +  J_lm^z δ_τ +  logcosh(2J_lm^xyδ_τ  ), 
     
     - W_l-m,[l]^''  + W_l-m,[m]^''    =    log C +  J_lm^z δ_τ +  logsinh(2J_lm^xyδ_τ) , 
     
    W_l+m,[l]^'' + W_l+m,[m]^''    =    log C -J_lm^z  δ_τ.

These conditions give

    W_j[l]^'    =     W_lj, 
    
    W_j[m]^'    =     W_mj, 
    Δ W_lj    =     - W_lj, 
    Δ W_mj    =      - W_mj,

and 

    W_l [l]^''  =    W_m [m] ^''   =    - J_lm^z δ_τ/2   -  1/4logsinh(2J_lm^xyδ_τ ),
    
     
    W_l [m]^''  =    W_m[l] ^''   =      -J_lm^z  δ_τ/2 - 1/4logcosh(2J_lm^xyδ_τ ).





Recovery of standard DBM. The direct interactions between (σ_l^z, d_[l]), (σ_m^z, d_[m]), (σ_l^z, d_[m]), and (σ_m^z, d_[l]), 
are mediated by h_[lm1],  h_[lm2],  h_[lm3],  and h_[lm4], respectively, as follows 

    exp( σ_l^z d_[l] W_l [l]^''  )    =    
    C_[lm1]∑_h_[lm1]exp( σ_l^z h_[lm1] W_l[lm1] + h_[lm1]  d_[l]W'_[lm1][l] ), 
    exp( σ_m^z d_[m] W_m [m]^''  )    =    
    C_[lm2]∑_h_[lm2]exp( σ_m^z h_[lm2] W_m[lm2] + h_[lm2]  d_[m]W'_[lm2][m] ), 
    exp( σ_l^z d_[m] W_l [m]^''  )    =    
    C_[lm3]∑_h_[lm3]exp( σ_l^z h_[lm3] W_l[lm3] + h_[lm3]  d_[m]W'_[lm3][m] ), 
    exp( σ_m^z d_[l] W_m [l]^''  )    =    
    C_[lm4]∑_h_[lm4]exp( σ_m^z h_[lm4] W_m[lm4] + h_[lm4]  d_[l]W'_[lm4][l] ).

By applying the gadget Eq. (<ref>), the new W, W' interactions are given by, for small δ_τ
(such that e^-J_lm^z δ_τ/√(sinh(2J_lm^xyδ_τ) ) > 1): 

    W_l[lm1] = W'_[lm1][l] = W_m[lm2] = W'_[lm2][m] = 1/2 arcosh( e^-J_lm^z δ_τ/√(sinh(2J_lm^xyδ_τ) )), 
          W_l[lm3] = - W'_[lm3][m] = W_m[lm4] = - W'_[lm4][l] = 1/2 arcosh( √(cosh(2J_lm^xyδ_τ)  )×  e^J_lm^z δ_τ).




How to enforce the constraint   (“" trick). 
Here, we discuss how to design the network to satisfy the constraint σ^z_l + σ^z_m = d_[l] + d_[m]. 

We rewrite the sum with the constraint in Eq. (<ref>) as follows (we ignore trivial constant factor):


    ∑_d_[l], d_[m]
     d_[l] + d_[m] = σ^z_l + σ^z_m  ⟶∑_  d_[l] , d_[m]∑_  h_[lm5] , h_[lm6] 
       e^  i π/4 ( (σ^z_l + σ^z_m )   h_[lm5]   -   h_[lm5] ( d_[l] +  d_[m] )   )× e^  i π/8 ( (σ^z_l + σ^z_m )   h_[lm6]   -   h_[lm6] ( d_[l] +  d_[m] )   )
    
     = ∑_  d_[l] , d_[m] 2 cos (  π/4    (σ^z_l + σ^z_m -  d_[l]  - d_[m] )  ) × 
     2cos (  π/8    (σ^z_l + σ^z_m -  d_[l]  - d_[m] )  )

One can easily see that the second line of the equation 
gives nonzero contribution only when  d_[l] + d_[m] = σ^z_l + σ^z_m. 








Summary of the 2d-6h  representation.
The network changes induced by the bond propagator at each imaginary time step are summarized as follows. 
Eqs. (<ref>) and (<ref>) imply that W̅_lj = W_lj + Δ W_lj = 0
and  W̅_mj  = W_mj + Δ W_mj = 0, i.e., all the existing connections between physical spins and hidden neurons vanish.
Then, the lth and mth physical spins will be connected to the new hidden neurons h_[lm1], …, h_[lm6], 
The new deep neurons d_[l] and d_[m] are also connected to 
h_[lm1], …, h_[lm6].

In total, we have 16 new connections in the deep Boltzmann network. 

By continuing the imaginary time evolution, the neural network grows as in Fig. <ref>. 
The number of neurons increases linearly with the number N_ slice of Suzuki-Trotter time slice. 
For example, in the case of the one-dimensional Heisenberg model, the total number of deep and hidden neurons are N_ site (2  N_ slice  +1) and 3 N_ site (2 N_ slice +1), respectively. 
The number of nonzero connections in the network is 8 N_ site (2  N_ slice  +1). 
The origin of 2 N_ slice +1 is coming from the fact that we apply 𝒢 propagators 2 N_ slice +1 times when we apply the second-order Suzuki-Trotter decomposition.  
The “iπ/4, iπ/8" trick plays a role to preserve the total magnetization for deep spins at each imaginary-time step, i.e.,   ∑_k d_k(t +1) = ∑_k d_k(t),  
where d(t +1) [d(t)] are the deep neurons introduced at (t+1)-th [t-th] step.  






Relationship between the 2d-6h representation and the path-integral quantum Monte Carlo method. 
In the final part of this section, we discuss the similarity between the 2d-6h representation and the imaginary-time path-integral quantum Monte Carlo method <cit.>.  
We will show that, in the 2d-6h representation, the deep neurons can be regard as the additional degrees of freedom along the imaginary time in the path-integral formulation. 

In the quantum Monte Carlo simulations using Suzuki-Trotter decomposition <cit.>,  the partition function Z is evaluated as 

    Z     =   ⟨σ^z (0) |  e^- βℋ  |  σ^z(0) ⟩
       ≃   ∑_σ^z(0),…, σ^z(2N_ slice-1) ⟨σ^z(0) |   e^-ℋ_2 δ_τ   | σ^z(2N_ slice -  1)  ⟩⟨σ^z(2N_ slice- 1) |  e^-ℋ_1 δ_τ   | σ^z (2N_ slice-2) ⟩
    
     
     …
          ∑_σ^z(0),…, σ^z(2N_ slice-1) …⟨σ^z (4)  |  e^-ℋ_2 δ_τ  | σ^z(3)  ⟩⟨σ(3)  |
     e^-ℋ_1 δ_τ   | σ^z(2)  ⟩⟨σ^z (2) | 
     e^-ℋ_2 δ_τ
      | σ^z (1) ⟩⟨σ^z (1) | 
     e^-ℋ_1 δ_τ 
     | σ^z (0) ⟩

In the evaluation of the matrix element of  ⟨σ^z (t + 1) |  e^-ℋ_νδ_τ  | σ^z (t) ⟩ (ν=1 or 2), 
in the case of one-dimensional Heisenberg model, it is sufficient to consider one specific bond, 
⟨σ^z_l (t + 1) σ^z_m (t + 1) |  e^-ℋ_lmδ_τ  | σ^z_l(t) σ^z_m (t) ⟩. The matrix elements are given by

    ⟨σ^z_l (t + 1) σ^z_m (t + 1) |  e^-ℋ_lmδ_τ  | σ^z_l(t) σ^z_m (t) ⟩ = 
    e^ J_lm^z δ_τ[     e^ - 2J_lm^z δ_τ                    0                    0                    0;                    0  cosh(2  J_lm^xyδ_τ) sinh (2  J_lm^xyδ_τ)                    0;                    0 sinh (2  J_lm^xyδ_τ)  cosh(2  J_lm^xyδ_τ)                    0;                    0                    0                    0     e^ - 2J_lm^z δ_τ;                      ]

in the basis {  | ↑↑⟩,    | ↑↓⟩,   |  ↓↑⟩,   | ↓↓⟩}. 

On the other hand,

the imaginary time evolution in Eq. (2) in the main text [or equivalently, Eq. (<ref>)] can be rewritten as 

    ⟨σ^z |Ψ (τ)⟩    =    ∑_σ^z(1),…, σ^z(2N_ slice +1) ⟨σ^z |   e^-  ℋ_1 δ_τ/2   |  σ^z(2N_ slice+1) ⟩⟨σ^z(2N_ slice+ 1) |  e^-ℋ_2 δ_τ  | σ^z (2N_ slice) ⟩
    
     
     …
          ∑_d(1),…, d(2N_ slice+1) …⟨σ^z(4)  | 
     e^-ℋ_1 δ_τ
      | σ^z(3) ⟩⟨σ^z (3)  | 
     e^-ℋ_2 δ_τ
      | σ^z(2 ) ⟩⟨σ^z  (2)  | 
     e^-ℋ_1 δ_τ/2 
     | σ^z (1) ⟩⟨σ^z (1)  | Ψ_0⟩

by inserting complete basis sets at each time slice.
The matrix element used here is exactly the same as that of QMC in Eq. (<ref>).
Here, the D dimensional quantum spin system is mapped on the D+1 dimensional classical system as in the case of the 
path integral quantum Monte Carlo method.  
Because the neuron spins are defined as the classical Ising-type spins, we can represent the summation over σ^z(1),…, σ^z(2N_ slice+1) by the summation over N_ site (2  N_ slice  +1) neuron spins.  
Assuming that these N_ site (2  N_ slice  +1) neuron spins are in the deep layer,

the imaginary time evolution in Eq. (<ref>) reads

    ⟨σ^z |Ψ (τ)⟩    =    ∑_d(1),…, d(2N_ slice + 1) ⟨σ^z |   e^-ℋ_1 δ_τ/2   | d(2N_ slice+1) ⟩⟨  d(2N_ slice+ 1) |  e^-ℋ_2 δ_τ  | d(2N_ slice) ⟩
    
     
     …
          ∑_d(1),…, d(2N_ slice+1) …⟨  d(4)  | 
     e^-ℋ_1 δ_τ
      | d (3) ⟩⟨ d(3)  | 
     e^-ℋ_2 δ_τ
      | d (2 ) ⟩⟨ d (2)  | 
     e^-ℋ_1 δ_τ/2 
     | d(1) ⟩⟨  d(1)  | Ψ_0⟩.

The matrix element ⟨  d_l (t + 1) d_m (t + 1) |  e^-ℋ_lmδ_τ  | d_l(t) d_m(t) ⟩ can be reproduced, for example, by 
the following interaction 


    e^W”_1(  d_l (t + 1)  d_l (t) +  d_m (t + 1) d_m (t) )
    e^W”_2   (  d_l (t + 1)  d_m (t) +  d_m (t + 1) d_l (t) ) 
          ×cos (    π/4 ( d_l (t   +  1)  + d_m (t  + 1) -     d_l (t) - d_m (t)   )  )  
    cos (   π/8 ( d_l (t  + 1)  + d_m (t  + 1) -     d_l (t) - d_m (t)    )  )

with 

    W”_1     =    - J_lm^zδ_τ/2   -  1/4logsinh(2J_lm^xyδ_τ ) , 
    
     
    W”_2     =      -J_lm^z δ_τ/2 - 1/4logcosh(2J_lm^xyδ_τ ),

This interaction can be mediated by adding hidden neurons and mediating the interactions between d (t + 1) and d(t). 
Then, Eq. (<ref>)  can be mapped onto the DBM representation. 

Indeed, the 2d-6h representation presented in this section correspond to this specific DBM construction: 
In the 2d-6h representation, two deep neurons are introduced for each bond at each imaginary time evolution. 
Because each imaginary time evolution acts on either even or odd bonds, the number of deep neurons introduced at one step is exactly same as the number of physical spins. 
In this case, the deep neurons can be considered as the spin degrees of freedom in the imaginary time layers d(1),…, d(2N_ slice+1). 
The interactions in Eqs. (<ref>) and (<ref>) are equivalent to those in Eqs. (<ref>) and (<ref>).
The “iπ/4, iπ/8" trick appears to put constraint to conserve the total magnetization at each layer. 
Therefore, the 2d-6h representation is equivalent to the path-integral formulation.
Indeed, if we rearrange the neurons in this DBM construction (Fig. <ref>), one can see a clear correspondence between the DBM network and the path-integral formulation. 
The extended systems including physical spins and deep neurons can be regard as the D+1 dimensional classical spin systems mapped from D dimensional quantum systems. 




  §.§.§ 2 deep, 4 hidden (2d-4h) representation



Strategy. 
We first extend DBM in the following way:

    Ψ_𝒲̅(σ^z)
        =    ∑_{h ,  d }∑_d_[l] P_1(σ^z,h)P_2(h,d)
    e^∑_j,n=l,mσ^z_nh_jΔ W_nj + ∑_jh_jd_[l]W_j[l]^' 
      + ∑_n=l,mσ^z_nd_[l]W^”_n[l]+∑_jσ^z_l σ^z_mh_jd_[l]Z_lmj.

Here, we have introduced terms which break the standard DBM form, in particular the terms proportional to W^”_n[l] and Z_lmj with n=l,m. Those are essential for this construction, and their reduction to the pure DBM will be shown later.  
Also notice that the sum over j runs through all the hidden neuron sites coupled to σ_l^z and σ_m^z,
 thus it incorporates nonlocal couplings between hidden variables (h), physical (σ^z) and deep (d) variables.
The term proportional to a_i in P_1(σ^z,h) is a local site-dependent magnetic-field term in the DBM acting on the physical variables σ^z, which can also flexibly represent any local gauge transformation, if a_i is taken complex. Here we fix a_i to be site-dependent constants, which stay unchanged through the imaginary time evolution.
We later use the fact that the gauge transformation σ^x → -σ^x and σ^y → -σ^y on one of the sublattices (or J^xy→ -J^xy) on a bipartite lattice as in Eq. (<ref>) is equivalent to a_i= iπ/2 if i is on this sublattice and a_i=0 on the other sublattice as a special choice of a_i.

In the imaginary time evolution of ℋ_lm, we update 
W̅_nj (n=l,m) 
with the increment Δ W_nj, in such a way that W̅_nj=W_nj+Δ W_nj. In addition to the deep variable d_[l], we further introduce one additional deep variable d_[lm] to recover the standard DBM by transforming the term proportional to W^” and Z, with supplementary four hidden variables. 




Derivation for the update of parameters. 
For σ^z_lσ^z_m=-1, the imaginary time evolution of the bond ℋ_lm is given as

    ⟨σ^z | e^-δ_τ(J_lm^zσ_l^zσ_m^z+2J_lm^xy(σ_l^+σ_m^-+σ_l^-σ_m^+))|Ψ_𝒲⟩    =    Ψ_𝒲(σ^z) e^J_lm^zδ_τcosh(2J_lm^xyδ_τ) 
       -   Ψ_𝒲(σ_1^z,…-σ_l^z,⋯ -σ_m^z…) e^J_lm^zδ_τsinh(2J_lm^xyδ_τ)
        =     C'⟨σ^z |Ψ_𝒲̅⟩,

which is equivalent to 

    ∑_{ h,d }Ψ_𝒲[1-tanh(2J_lm^xyδ_τ)e^-2∑_n=l,m (σ_n^z∑_jh_jW_nj+a_nσ_n^z ) ]     =     CΨ_𝒲̅

and C=(e^-J_lm^zδ_τ/cosh( 2J_lm^xyδ_τ) )C'.
Notice that, here, we keep the bias term a_n in Eq. (<ref>) instead of 
applying the gauge transformation in Eq. (<ref>). 

For  σ^z_lσ^z_m=1, we obtain

    ∑_{ h,d }Ψ_𝒲e^-2J_lm^zδ_τ /cosh ( 2J_lm^xyδ_τ )    =     CΨ_𝒲̅.



To make these imaginary time evolutions exact, 
W_nj (n=l,m)
is updated to W̅_nj with the increment Δ W_nj as W̅_nj=W_nj+Δ W_nj 
with

    Δ W_lj   =   -Δ W_mj= -1/2(W_lj-W_mj).

The new couplings W_j[l]^', Z_lmj and W”_n[l] are also
given by

    W_j[l]^'  = -Z_lmj= -1/2(W_lj-W_mj)

and from 

    2 ( W_l[l]^''  -   W_m[l]^''  )
        =    log[-e^
    -2a_l-mtanh (2J_lm^xyδ_τ)]

and

    2cosh(W^”_l[l] + W^”_m[l]) =e^-2J_lm^zδ_τ-W^”_l-m/cosh(2J_lm^xyδ_τ),

we obtain

    W”_l[l]   =   1/4 [log [-e^-2a_l-mtanh (2J_lm^xyδ_τ)  ]
    +2 arcosh [e^-2J_lm^zδ_τ/√(-2e^-2a_l-msinh (4J_lm^xyδ_τ))]  ] 
    
    W”_m[l]   =   1/4 [-log [-e^
    -2a_l-mtanh (2J_lm^xyδ_τ)  ]
    +2 arcosh [e^-2J_lm^zδ_τ/√(-2e^-2a_l-msinh (4J_lm^xyδ_τ))]  ]

with a_l-m = a_l -a_m.
On a bipartite lattice, to avoid the negative sign (or complex phase) problem  we need to keep  
W”_l[l] and W”_m[l] real.

This can be achieved by choosing a_l=0 for any l if J_lm<0 (ferromagnetic case). For J_lm>0 (antiferromagnetic case), a_l =nπ i with an arbitrary integer n if the site l belongs to the sublattice A and a_l =(n+1/2)π i if l belongs to the sublattice B.
This local gauge for J_lm>0 is equivalent to take J_lm^xy→ -J_lm^xy and 
a_l=0 
for any site l as is formulated in Eq.(<ref>).
We further note that W”_m[l] can be taken positive if we take sufficiently small δ_τ in Eq. (<ref>), with the leading order term -log (2J^xy_lmδ_τ )/2. On the other hand, 
in Eq. (<ref>), the leading order term is negative (=-J_lmδ_τ).




Recovery of the standard DBM form. 
To recover the original form of the DBM, we first use Eq. (<ref>)
with the replacement s_1 →σ^z_n, s_2 →  d_[l], 
s_3 →  h_[n], 
C →  D_n,
V →  W”_n[l]
Ṽ_1 →  W_n[n] and Ṽ_2 →  W'_[n][l]
for n=l,m. 
We have added here two hidden variables h_[l] and h_[m].
Then a solution for D_n,
W_n[n], and W'_[n][l] are represented by using W^”_n[l] as 

    D_n    =   1/2exp[-W^”_n[l]] 
    
    W_n[n]    =     W'_[n][l]=1/2 arcosh(exp[2W^”_n[l]]),

if W^”_n[l] is positive (as in the case of W^”_m[l] for small δ_τ), which gives real W_n[n] and W'_[n][l]. On the other hand, if  W^”_n[l] is negative (as in the case of W^”_l[l] for small δ_τ),
we should take 

    D_n    =   1/2exp[W^”_n[l]] 
    
    W_n[n]    =     -W'_[n][l]=1/2 arcosh(exp[-2W^”_n[l]]),

to give real W_n[n] and W_[n][l].

To completely recover the original DBM form, we next use Eq. (<ref>) by replacing
σ_1 with σ^z_l, σ_2 with σ^z_m, d_1 with d_[l], d_2 with d_[lm], h_1 with h_j, h_2 with h_[lm1], h_3 with h_[lm2], and V with Z_lmj.


With these solutions, 
by ignoring the trivial constant factors including D_l and D_m, 
the evolution is described by introducing two deep and four hidden additional variables d_[l], d_[lm], h_[l], h_[m], h_[lm1], and h_[lm2] as

    Ψ_𝒲̅(σ^z)
        =    ∑_{h̅ ,  d̅} P_1(σ^z,h)P_2(h,d) 
    exp [
    ∑_j,n=l,mσ^z_nh_jΔ W_nj
    +∑_jh_jd_[l]W_j[l]^'
       +   ∑_n=l,mh_[n](σ^z_nW_n[n]+d_[l]W^'_[n][l]) +d_[lm]∑_jh_jZ_lmj
       +   iπ/4(h_[lm1]+h_[lm2])(σ^z_l+σ^z_m+d_[l]+d_[lm])   ],

where {h̅ , d̅} is a set consisting of the existing and new neurons.  
Equation (<ref>) recovers the standard form of deep Boltzmann machine, where the physical spins σ^z as well as the deep variables d are not interacting each other and couples only to the hidden variables h. 






Summary.
After summing over {h̅},
we reach

    Ψ_𝒲̅(σ)
        =    ∑_{d̅}exp  [ ∑_n=l,ma_nσ^z_n]∏_j [2cosh[∑_iσ^z_iW_ij+∑_kW^'_jkd_k
    +d_[l]W_j[l]^'+d_[lm]Z_lmj ]  ]
    
       ×   
    ∏_n=l,m ( 2cosh [σ^z_nW_n[n]+d_[l]W^'_[n][l]]   ) 
    
     [2 cos[π/4(σ^z_l+σ^z_m+d_[lm]+d_[lm])]  ]^2

where the parameters W, W' and Z are given 
in Eqs. (<ref>), (<ref>), (<ref>), (<ref>), and (<ref>) (or (<ref>)). 

We have introduced 2 deep and 4 hidden variables. Among them, h_[lm1] and h_[lm2] are simply to relate  d_[l] and d_[lm] to σ^z_l and σ^z_m. 
With this trick, one can constrain d_[l]=d_[lm] for σ^z_l=σ^z_m and d_[l]=-d_[lm] for σ^z_l=-σ^z_m.
After repeatedly operating Eq. (<ref>), for all the combinations of l,m, 
the DBM structure becomes nonlocal as we see in Fig. <ref>.
After the sufficiently long imaginary-time evolution,  with the analytical sum on { h } and the Monte sampling over { d }, 
 
 one can obtain the ground state wave function. 





§ SAMPLING


Once we have determined specific rules to obtain the parameters of the
DBM, the remaining question to be addressed is how to compute expectation
values of physical quantities. Consider a quantum operator 𝒪,
then its expectation value over the DBM is given by the expression

    ⟨𝒪⟩    =    ∑_{σ^z,h,h^'d,d^'}Π(σ^z,h,h^',d,d^')O_loc(σ^z,h,h^')/∑_{σ^z,h,h^'d,d^'}Π(σ^z,h,h^',d,d^'),

where we have introduced the pseudo-probability density Π(σ^z,h,h^',d,d^')≡ P_1(σ^z,h)P_2(h,d)P_1^⋆(σ^z,h^')P_2^⋆(h^',d^'),
and the “local” estimator O_loc(σ^z,h,h^')=1/2∑_{σ^' z}⟨σ^z|𝒪|σ^' z⟩(P_1(σ^' z,h)/P_1(σ^z ,h)+P_1(σ^' z,h^')^⋆/P_1(σ^z ,h^')^⋆).
For a large number of spins and hidden/deep units, it is not possible
to compute those sums numerically, because of the exponential number
of terms involved. However, there are specific cases in which efficient
sampling strategies can be devised, allowing to stochastically compute
the quantum expectation values. In general, when the DBM weights are
all real Π(σ^z,h,h^',d,d^')≥0, and it can
be interpreted as an (unnormalized) probability density. Thus, Markov-chain
sampling techniques can be applied, similarly to the case of applications
in standard machine learning. In the case of complex-valued weights,
the straightforward probabilistic interpretation breaks down, and
a sign (phase) problem arises. However, there are specific cases in
which one can still recover a properly defined probability density,
and efficiently sample from it. In the following we describe two main
sampling methods based on Markov chain techniques. First, Gibbs sampling,
then Metropolis-Hastings sampling. In both cases we discuss when the
sign problem can be circumvented. 



 §.§ Gibbs sampling 


We start discussing a strategy which is the natural generalization
of what traditionally used in most applications of DBM in machine
learning. The approach is based on Gibbs sampling, a strategy which
amounts to generate samples using the exact conditional probabilities
for block of variables. In practice, we introduce three kind of moves,
which allow to generate a Markov chain of visible, hidden, and deep
variables distributed according to Π(σ^z,h,h^',d,d^').



  §.§.§ Sampling visible spins


The first kind of move consists in freezing all the hidden and deep
variables, and sampling the visible spins σ^z. Specifically,
we generate new visible spin configurations according to the conditional
probability: 

    Π(σ^z  |  h,h^',d,d^')     =    P_1(σ^z , h )P_2(h,d)P_1^⋆(σ^z ,h^')P_2^⋆(h^',d^')/∑_{σ̃^z }P_1(σ̃^z,h)P_2(h,d)P_1^⋆(σ̃^z,h^') P_2^⋆(h^',d^') 
        =    P_1(σ^z,h)P_1^⋆(σ^z,h^')/∑_{σ̃^z }P_1(σ̃^z,h)P_1^⋆(σ̃^z,h^')
        =    Π_i^Nexp{σ_i^z [∑_j(h_jW_ij+h_j^'W_ij^⋆)+2a_i^r]}/Π_i^N2cosh(∑_j(h_jW_ij+h_j^'W_ij^⋆)+2a_i^r).

Here, a_i^r is a real part of a_i. 
A particularly appealing aspect of this transition probability is
that each visible spin can be treated independently from the others,
thus we can update in parallel all visible spins at once. The probability
of a given spin to be up for example is: 

    P(σ^z_i=1|h,h^',d,d^')     =    Logistic (2λ_i^[σ^z]),

with λ_i^[σ^z]=∑_j(h_jW_ij+h_j^'W_ij^⋆)+2a_i^r,
and Logistic(x)=1/1+exp(-x). Thus, during this
phase we generate N random numbers η_i uniformly distributed
in [0,1), and set the spin σ_i^z=1 if η_i<Logistic(2λ_i^[σ^z]).
For this approach to be feasible, we must have that the λ_i^[σ^z]
are real. 
Necessary conditions for this condition to be satisfied
are discussed at the end of this section.




  §.§.§ Sampling hidden spins


The second type of move consists in freezing visible and deep spins,
and sampling hidden variables h and h^'. For example,
to sample h the transition probability reads: 

    Π(h| σ^z,h^',d,d^')     =    P_1(σ^z ,h)P_2(h,d)/∑_{h̃}P_1(σ^z,h̃)P_2(h̃,d)
        =    Π_j^Mexp[h_j(∑_iσ^z_iW_ij+b_j+∑_kd_kW_jk^')]/Π_j^M2cosh(∑_iσ^z_iW_ij+b_j+∑_kd_kW_jk^').

The probability of having h_j=1 is then: 

    P(h_j=1|σ^z,h^',d,d^')     =    Logistic(2λ_j^[h]),

with λ_j^[h]=∑_iσ^z_iW_ij+b_j+∑_kd_kW_jk^'.
Again, one can therefore efficiently update all the M hidden spins
at once, without rejection. Analogously, for h^' we have
λ_j^[h^']=∑_iσ^z_iW_ij^⋆+b_j^⋆+∑_kd_k^'W_jk^'^⋆. 



  §.§.§ Sampling deep spins


The final set of moves consists in freezing visible and hidden spins,
and sample from deep variables d and d^'. For example,
to sample d the transition probability is: 

    Π(d|σ^z,h,h^',d^')     =    P_2(h,d)/∑_{d̃}P_2(h,d̃)
        =    Π_k^M^'exp[d_k(∑_jh_jW_jk^'+c_k)]/Π_k^M^'2cosh(∑_jh_jW_jk^'+c_k).

The probability of having d_k=1 is then: 

    P(d_k=1|σ^z,h,h^',d^')     =    Logistic(2λ_k^[d]),

with λ_k^[d]=∑_jh_jW_jk^'+c_k. Analogously,
we have λ_k^[d^']=∑_jh_j^'W_jk^'⋆+c_k^⋆. 



  §.§.§ Overall scheme: alternate block sampling 


The overall sampling scheme is therefore realized putting together
all those individual Gibbs samplings. In particular, we can devise
a two-step block sampling, which takes into account the conditional
dependence of all the probabilities previously derived. 

The overall sampling scheme then works as follow: 


  * Sample h and h^', fixing all the other variables. This
is realized using the probabilities (<ref>) for all
the hidden spins.

  * Sample σ^z,d,d^' fixing the values of h and
h^'.This is realized using the probabilities (<ref>)
and (<ref>) for all the visible and deep spins, respectively.

  * Cycle between 1 and 2.




  §.§.§ Phase problem in the Gibbs scheme


In order to get a consistent sampling scheme, we must have that all
the quantities λ_i^[σ^z],λ_j^[h],λ_j^[h^'],λ_k^[d],λ_k^[d^']
are real valued. In the absence of this condition, we have a phase
problem, and we cannot directly use a stochastic approach to sample
from the DBM. Looking more closely at what conditions are needed,
we start noticing that the visible bias can take arbitrary (complex)
values, since only the real parts, a_i^r, enter λ_i^[σ^z].
In general, there might be specific choices of the DBM parameters which still
guarantee absence of phase problem. One possibility is realized, for
example, when fixing the total magnetizations in the three layers,
i.e. the constraints ∑_iσ^z_i=σ^z_tot,
 ∑_jh_j=h_tot, ∑_kd_k=d_tot.
We further assume that Im(W_ij)=W^I, a constant,
as well as Im(W_jk^')=W^^'I.
Then, it is easy to see that the phase problem is avoided when b_j^I=-σ^z_totW^I-d_totW^^'I
and c_k^I=-h_totW^'I. Notice
that those are just a specific set of conditions, and less stringent
ones can be found using other sampling schemes. 



When each sample has the imaginary part or negative signs, another possibility of avoiding the phase problem is to take the partial trace summation explicitly so that such partial sum gives always a real nonnegative value.  
We will discuss this point in more detail in the next section. 





 §.§ Metropolis sampling





  §.§.§ Marginal probability density

Because there are no intralayer interactions in the DBM architecture, one can analytically trace out either one of h, h' and d,d'. 
Then we get marginal probability density: 
Π̃(σ^z,h,h^')=∑_{ d,d^'}Π(σ^z,h,h^',d,d^')
or Π̃'(σ^z,d,d^')=∑_{ h,h^'}Π(σ^z,h,h^',d,d^'). 
Defining P̃ ( σ^z, h ) and P̃' (σ^z , d) as 

    P̃ ( σ^z, h )    =   ∑_{ d  } P_1(σ^z , h) P_2(h ,d ) 
        = 
         e^∑_i   a_i σ_i^z + ∑_ijσ_i^z  h_j W_ij  +∑_j  b_j h_j ×∏_k   2 cosh (c_k  + ∑_k h_j W'_jk )

and 

    P̃ '( σ^z, d )    =   ∑_{ h  } P_1(σ^z , h) P_2(h ,d ) 
        = 
        ∏_j   2 cosh (b_j  + ∑_i σ^z_i W_ij + ∑_k d_k W'_jk )   
       ×    e^∑_i a_i σ_i^z + ∑_k c_k d_k ,

respectively, the marginal probability densities are given by 

    Π̃(σ^z,h,h^')    =   ∑_{ d,d^'}Π(σ^z,h,h^',d,d^') =  
     P̃ ( σ^z, h ) P̃^⋆  ( σ^z, h' ),   
    Π̃'(σ^z,d,d^')    =   ∑_{ h,h^'}Π(σ^z,h,h^',d,d^') =  
      P̃^' ( σ^z, d ) P̃^'⋆  ( σ^z, d' ).


With these marginal probability densities, we perform the Metropolis sampling to measure physical quantities. 
The expectation value of a quantum operator 𝒪 is given by 

    ⟨𝒪⟩  =  ∑_{σ^z,h,h^'}Π̃ (σ^z,h,h^') Õ_loc(σ^z,h,h^')/∑_{σ^z,h,h^'}Π̃ (σ^z,h,h^') 
    =  ∑_{σ^z,d,d^'}Π̃ '(σ^z,d,d^') Õ'_loc(σ^z,d,d^')/∑_{σ^z,d,d^'}Π̃' (σ^z,d,d^')

with 

    Õ_loc(σ^z,h,h^')   =   1/2∑_{σ^' z}⟨σ^z|𝒪|σ^' z⟩(P̃ (σ^' z,h)/P̃(σ^z ,h)+P̃(σ^' z,h^')^⋆/P̃ (σ^z ,h^')^⋆),  
    Õ'_loc(σ^z,d,d^')    =   1/2∑_{σ^' z}⟨σ^z|𝒪|σ^' z⟩(P̃ '(σ^' z,d)/P̃'(σ^z ,d)+P̃'(σ^' z,d^')^⋆/P̃' (σ^z ,d^')^⋆).





  §.§.§ Phase problem in the Metropolis scheme


An advantage of choosing the marginal probability density is that 
by taking the summation over h and d, the sign problem can sometimes be avoided
even if the DBM has complex parameters. 
An example is to take the summation over the hidden variables h analytically 
in the three DBM constructions for the Heisenberg models presented in Sec.  <ref>. 
In all the three cases, only those W and W' couplings used to enforce the constraints are complex-valued, and the summation over h eliminates the negative weight. 
For example, in the case of the 2d-4h representation in Sec. <ref>,  
though each sample may have a finite imaginary part as in each term of Eq. (<ref>), 
the total weight becomes real and nonnegative, after the explicit summation over the h degrees of freedom is performed as in Eq. (<ref>).

When the lattice is not bipartite, we can still write down the DBM solutions to exactly follow the imaginary time evolutions. 
However, in this case, we will have imaginary W and W' parameters even for the units not involved in enforcing the constraints. 
In this case, the sampling may suffer from sign problem. 
However, as we discuss in the main text, in contrast to the conventional quantum Monte Carlo simulations, we can make the number of imaginary time step to reach the ground state short by starting the analytical DBM time evolution [Eq. (<ref>)]
from a good stating point | Ψ_0 ⟩. 
For example, numerically optimized RBM wave functions can be used for | Ψ_0 ⟩, or more generally, 
| Ψ_0 ⟩ can be wave functions used in the conventional wave function techniques. 
In this case, before we suffer from a severe sign problems, we might be able to reach the ground state with good statistical accuracy. 




  §.§.§ Overall scheme

We sample over σ^z, h, h' [or σ^z, d, d'] with the marginal probability density Π̃(σ^z,h,h^')
[ Π̃'(σ^z,d,d^') ]. 
The physical quantities are measured following Eq. (<ref>).
In the case of Heisenberg model, after tracing out the h spins, 
we have constraints over the values of σ^z, d, d^'.  
In that case, a cluster update rather than a local update will be more efficient. 
In particular, in the 2d-6h representation,
 since the imaginary-time evolution of the DBM is equivalent to the path-integral formalism, 
we can apply an efficient cluster update used in the conventional quantum Monte Carlo method, such as so called loop update <cit.>. 


