# Scaling of neural-network quantum states for time evolution

Sheng-Hsuan Lin  Department of Physics, TFK, Technische Universität München,
James-Franck-Straße 1, 85748 Garching, Germany  Frank Pollmann  Department of
Physics, TFK, Technische Universität München, James-Franck-Straße 1, 85748
Garching, Germany

###### Abstract

Simulating quantum many-body dynamics on classical computers is a challenging
problem due to the exponential growth of the Hilbert space. Artificial neural
networks have recently been introduced as a new tool to approximate quantum-
many body states. We benchmark the variational power of the restricted
Boltzmann machine quantum states and different shallow and deep neural
autoregressive quantum states to simulate global quench dynamics of a non-
integrable quantum Ising chain. We find that the number of parameters required
to represent the quantum state at a given accuracy increases exponentially in
time. The growth rate is only slightly affected by the network architecture
over a wide range of different design choices: shallow and deep networks,
small and large filter sizes, dilated and normal convolutions, with and
without shortcut connections.

††preprint: APS/123-QED

##  I  Introduction

Numerical simulations are pivotal for the understanding of the emergent
physics of quantum-many body systems. Of particular current interest are
quantum many-body systems out of equilibrium, where several exciting phenomena
have recently been discovered such as many-body localization [1, 2, 3, 4, 5]
and quantum many-body scars [6, 7, 8]. The interest in such systems is driven
by rapid experimental progress, e.g., cold atoms and trapped ions allow for an
unprecedented level of control and long coherence times [9, 10, 11, 12].

Over the past decades, a large arsenal of numerical tools has been introduced
to tackle such systems. Widely used are exact diagonalization (ED) and Krylov-
subspace methods, which simulate time evolution exactly. The system size,
however, is strongly limited by the exponential growth of the many-body
Hilbert space. In many cases, the size limitations can be remedied by using
tensor network (TN) methods [13, 14, 15, 16, 17]. These include matrix-product
states (MPS) [18, 19, 20] methods for one-dimensional (1D) systems, which
represent powerful numerical methods that are not limited by the system size
but instead by the entanglement. While ground states of local gapped
Hamiltonians [21] satisfy the entanglement area-law [22, 23] and can thus be
efficiently simulated, the application of MPS to study non-equilibrium
dynamics [24, 25, 26, 27, 28, 29] is strongly limited by the fast growth of
entanglement. In particular, the entanglement generically grows linearly and
hence leads to an exponential growth of the required parameters.

A promising new approach to simulate the dynamics of many-body systems is
based on a representation of quantum states using artificial neural networks
(ANNs) [30]. Recent works have shown not only promising results for simulating
non-equilibrium dynamics [31], but also in describing two-dimensional (2D)
critical ground states [32] and states with chiral topological order [33].
Importantly, states based on ANNs can efficiently encode volume law
entanglement [34, 35] and are thus per se not limited by the entanglement
growth in non-equilibrium systems. It remains an open question what limits the
representation based on ANNs and how the number of parameters generically
scales.

Figure 1:  Illustration of a neural autoregressive quantum state. The neural
network approximates the conditional probability amplitudes
$\Psi_{1},\Psi_{2},\ldots$ in the log space. Each of the conditional
probability amplitudes are normalized to give a normalized wavefunction. The
network connectivity is restricted (masked) for causality, such that
$\Psi_{j}$ only depends on $\sigma_{1}\ldots\sigma_{j-1}$. The network itself
could be either a fully-connected neural network or a convolutional neural
network.

In this paper, we benchmark the variational power of a number of different
ANNs approaches to simulate quench dynamics in 1D and compare to MPS based
approximations. We find numerical evidence that states based on ANNs can
describe quantum many-body states under quenched time evolution of local
Hamiltonian with a number of parameters that grows exponentially in time—which
in fact resembles the scaling of MPS based simulations. The remainder of the
paper is as follows: We first describe the variational wavefunctions based on
ANNs and the variational method to approximate quantum states using supervised
learning in Sec. II. To study the time scales accessible by different
variational wavefunctions, we variationally optimize wavefunctions to
approximate the exact time evolution in Section III. In Section IV, We
conclude by discussing the implication on numerical algorithms and possible
future directions.

##  II  Methods

In this section, we review variational wavefunctions based on general neural
networks and autoregressive neural networks. Moreover, we introduce a novel
cost function to approximate complex-valued quantum states with autoregressive
neural networks in the supervised learning framework [36, 37].

###  II.1  Neural-Network Quantum States

A pure quantum state $\ket{\Psi}$ describing a quantum many-body system in the
Hilbert space $\mathcal{H}=\otimes_{i}^{N}\mathcal{H}_{i}$ can be represented
as

$\displaystyle\ket{\Psi}$
$\displaystyle=\sum_{\sigma_{1},\sigma_{2}\cdots,\sigma_{N}}\Psi_{\sigma_{1}\sigma_{2}\cdots\sigma_{N}}\ket{\sigma_{1}}\otimes\ket{\sigma_{2}}\otimes\cdots\otimes\ket{\sigma_{N}},$

where $\\{\mathcal{H}_{i}\\}$ and $\ket{\sigma_{i}}$ are the local Hilbert
spaces and the corresponding local basis states, respectively. The exact
representation $\Psi_{\sigma_{1}\sigma_{2}\cdots\sigma_{N}}$ requires a number
of parameters that grows exponentially with the system size $N$. Neural-
network quantum states (NQS) [30] are variational wavefunctions that
approximate the amplitudes by neural networks $f$ with a polynomial number of
parameters $\text{poly}(N)$,

$\Psi_{\sigma_{1}\sigma_{2}\cdots\sigma_{N}}\approx
f(\sigma_{1},\sigma_{2},\cdots,\sigma_{N};\bm{w})\coloneqq\Psi^{\text{NN}}_{\sigma_{1}\sigma_{2}\cdots\sigma_{N}},$
(1)

where $\bm{w}$ are the parameters. One of the motivations to consider neural
networks is the universal approximation theorem [38, 39, 40, 41, 42], which
guarantees an arbitrarily close approximation provided a sufficiently large
number of hidden neurons is included (a review of neural networks can be found
in Appendix A). Various different types of neural networks have been
considered for approximating quantum states, including restricted Boltzmann
machines (RBMs) [43, 44] in Ref. 30, fully-connected feed-forward neural
networks (NNs) in Ref. 45, convolutional neural networks (CNNs) in Ref. 46,
and recurrent neural networks (RNNs) [47, 48] in Refs. 49, 50. NQS found
applications in studying ground states and low-energy states [51, 52, 53, 54,
55], quantum state tomography [56], simulating open quantum systems [57, 58,
50, 59, 60, 61], and real-time evolution [30, 31, 62, 63, 64, 65].

In general, NQS have the following properties: The unnormalized probability
amplitude and its gradient $\nabla_{\bm{w}}\Psi^{NN}$ can be efficiently
evaluated at a cost of $\mathcal{O}(N^{2})$. This property allows one to
combine NQS with conventional Monte Carlo methods. For example, the
expectation value of few-body observables $\hat{O}$ can be computed as

$\frac{\braket{\Psi^{NN}}{\hat{O}}{\Psi^{NN}}}{\braket{\Psi^{NN}}{\Psi^{NN}}}=\sum_{\bm{\sigma}}\frac{\lvert\Psi^{NN}(\bm{\sigma})\rvert^{2}}{\braket{\Psi^{NN}}{\Psi^{NN}}}O_{\text{local}}(\bm{\sigma}),$
(2)

where
$O_{\text{local}}(\bm{\sigma})=\sum_{\bm{\sigma}^{\prime}}\braket{\bm{\sigma}}{\hat{O}}{\bm{\sigma}^{\prime}}\Psi^{NN}(\bm{\sigma})/\Psi^{NN}(\bm{\sigma}^{\prime})$.
This expectation value can be estimated using Markov Chain Monte Carlo (MCMC)
sampling from the normalized probability $P_{\Psi^{NN}}(\bm{\sigma})$. One can
then minimize the energy expectation value by stochastic gradient descent
(SGD). This method is also known as variational Monte Carlo (VMC) [66, 67, 68,
69] and is often combined with the second-order optimization method,
stochastic reconfiguration [70, 71, 72], a.k.a. imaginary time-dependent
variational principle (iTDVP) or natural gradient descent [73]. While being a
versatile method, the NQS-based method has potential drawbacks including the
intractable normalization and the long auto-correlation time in MCMC, which
can lead to difficulties in training large deep neural networks.

###  II.2  Neural Autoregressive Quantum States

A proposed solution to resolve the potential drawbacks mentioned above is the
neural autoregressive quantum states (NAQS) [32, 49, 50, 74]. NAQS are NQS
with causal restrictions in the connectivity of the networks. Concretely, NAQS
approximate the conditional probability amplitudes of a given ordering of
sites $\\{1,\ldots,N\\}$ by a neural network,

$\displaystyle\Psi_{\sigma_{1}\sigma_{2}\cdots\sigma_{N}}=\Psi_{\sigma_{1}}\Psi_{\sigma_{2}|\sigma_{1}}\cdots\Psi_{\sigma_{N}|\sigma_{1}\cdots\sigma_{N-1}}$
(3) $\displaystyle\qquad\approx
f_{\text{NN}}^{[1]}(\sigma_{1})f_{\text{NN}}^{[2]}(\sigma_{2}|\sigma_{1})\cdots
f_{\text{NN}}^{[N]}(\sigma_{N}|\sigma_{1}\cdots\sigma_{N-1})$ (4)

where the network has $N$ outputs. Each output represents a conditional
probability amplitude on each site
$\Psi_{j}\coloneqq\Psi_{\sigma_{j}|\sigma_{i<j}}$. The fact that $\Psi_{j}$
can only depend on $\sigma_{i<j}$ leads to a causal structure. Since every
single site conditional probability amplitude $\Psi_{j}$ can be easily
normalized, the multiplication of them gives a normalized variational quantum
state (see Fig. 1).

As a subclass of NQS, NAQS share the same properties of NQS where one can
efficiently evaluate the probability amplitudes and the gradients.
Furthermore, a NAQS represents an autoregressive model. An autoregressive
model can be sampled by the ancestral sampling algorithm which allows sampling
directly from the probability
$P_{\Psi^{NN}}(\bm{\sigma})=\lvert\Psi^{NN}(\bm{\sigma})\rvert^{2}$ (see
Appendix B for further details). As a result, quantum states represented by a
deep neural network with millions of parameters can be trained efficiently
using the SGD algorithm for studying ground states [32].

###  II.3  Supervised learning for Quantum State Approximation

We now describe the variational problem of approximating complex-valued
wavefunctions in the supervised learning setting and propose a new cost
function for normalized variational wavefunctions, e.g. NAQS, to jointly learn
the probabilities and phases.

Given a target wave function $|\Phi\rangle$, we would like to find the optimal
NAQS $|\Psi^{\text{NN}}\rangle$ approximating the target wave function by
minimizing a cost function $\mathcal{R}(\Phi,\Psi^{\text{NN}})$. In the
following, we consider the setup that we are given the samples
$\\{\bm{\sigma},\Phi(\bm{\sigma})\\}$ according to the probability
distribution $P_{\Phi}(\bm{\sigma})=\lvert\Phi(\bm{\sigma})\rvert^{2}$ instead
of the full wave function. It is a practical setup for the exponentially large
Hilbert space. In the supervised learning setup, a loss function $\mathcal{L}$
measures the difference between the data points
$\braket{\bm{\sigma}}{\Phi},\braket{\bm{\sigma}}{\Psi^{\text{NN}}}$. A cost
function $\mathcal{R}$ is the expectation value of the loss function
$\mathcal{L}$ with respect to the probability $P_{\Phi}(\bm{\sigma})$ and is
therefore a measure of the difference between the quantum states
$\ket{\Phi},\ket{\Psi^{\text{NN}}}$. As a result, the optimization problem can
be solved using the SGD method, where the stochastic gradient is computed from
the sampled cost function $\hat{R}$ (see Appendix C for a detailed
discussion).

Here we present a joint learning scheme for complex-valued wave function where
the cost function $\mathcal{R}_{\text{joint}}$ is the sum of the so-called
Kullback–Leibler (KL) divergence of the probabilities
$\mathcal{R}_{\text{KL}}$ and the weighted $\mathcal{L}_{2}$ distance of the
phases $\mathcal{R}_{\theta}$,

$\mathcal{R}_{\text{joint}}=\mathcal{R}_{\text{KL}}+\mathcal{R}_{\theta}.$ (5)

In order to learn the magnitude of the probability amplitude
$\text{Re}\left[\log\Psi^{\text{NN}}(\bm{\sigma})\right]$, we consider the
forward KL-divergence $\mathcal{R}_{\text{KL}}$, which is defined as

$\displaystyle\begin{split}\mathcal{R}_{\text{KL}}&=\sum_{i}^{D_{\mathcal{H}}}P_{\Phi}(\bm{\sigma}_{i})\log{\left(\frac{P_{\Phi}(\bm{\sigma}_{i})}{P_{\Psi^{\text{NN}}}(\bm{\sigma}_{i})}\right)}\\\
&\approx\frac{\sum_{\bm{\sigma}_{i}\sim|\Phi(\bm{\sigma}_{i})|^{2}}2\text{Re}\left[\log\Phi(\bm{\sigma}_{i})-\log\Psi^{\text{NN}}(\bm{\sigma}_{i})\right].}{N_{\text{samples}}}\end{split}$
(6)

It is a common cost function for learning a probability distribution when one
can sample from the exact distribution, which is equivalent to maximal
likelihood learning.

For learning the phase of the probability amplitude
$\text{Im}\left[\log\Psi^{\text{NN}}(\bm{\sigma})\right]$, we consider the
weighted $\mathcal{L}_{2}$ distance of the phases embedded on a unit circle,

$\displaystyle\begin{split}\mathcal{R}_{\theta}&=\sum_{i}^{N}|\Phi(\bm{\sigma}_{i})|^{2}\
\text{dist}(\text{Im}\left[\log\Phi(\bm{\sigma}_{i})\right],\text{Im}\left[\log\Psi^{\text{NN}}(\bm{\sigma}_{i})\right])\\\
&=\sum_{i}^{N}|\Phi(\bm{\sigma}_{i})|^{2}\
\text{dist}(\theta^{\Phi}(\bm{\sigma}_{i}),\
\theta^{\Psi^{\text{NN}}}(\bm{\sigma}_{i}))\\\
&\approx\frac{1}{N_{\text{samples}}}\sum_{\bm{\sigma}_{i}\sim|\Phi(\bm{\sigma}_{i})|^{2}}\text{dist}(\theta^{\Phi}(\bm{\sigma}_{i}),\
\theta^{\Psi^{\text{NN}}}(\bm{\sigma}_{i}))\\\ \end{split},$ (7)

with

$\text{dist}(\theta_{1},\
\theta_{2})=\big{(}\cos\left(\theta_{1}\right)-\cos\left(\theta_{2}\right)\big{)}^{2}+\big{(}\sin\left(\theta_{1}\right)-\sin\left(\theta_{2}\right)\big{)}^{2}.$

The stochastic gradient can be computed from the sampling estimate of the cost
functions in Eq. (6) and Eq. (7). During the optimization procedure, the
magnitudes and phases are learned at the same time. This cost function is only
valid for normalized wavefunction because the KL-divergence is unbounded from
below if $P_{\Psi^{\text{NN}}}(\bm{\sigma}_{i})$ is an unnormalized
probability distribution. As a result, we consider NAQS which ensures the
normalization without additional cost.

The advantage of choosing the cost function $\mathcal{R}_{\text{joint}}$ is
that the resulting stochastic gradient is unbiased and with low variance. The
optimization problem can be solved by SGD with a small mini-batch size or even
a single data input at a time. It is important to have this property for
efficiently solving the variational state approximation problem with a large
neural network and when the Hilbert space is large.

Another cost function, which also can be rewritten as an expectation value and
gives a physical interpretation, is the real part of the negative overlap,

$\displaystyle\begin{split}\mathcal{R}_{\text{neg.
overlap}}&=-\text{Re}\left[\sum_{i}^{D_{\mathcal{H}}}{\Phi}^{*}(\bm{\sigma}_{i})\Psi^{\text{NN}}(\bm{\sigma}_{i})\right]\\\
&=-\sum_{i}^{D_{\mathcal{H}}}P_{\Phi}(\bm{\sigma}_{i})\text{Re}\left[\frac{\Psi^{\text{NN}}(\bm{\sigma}_{i})}{\Phi(\bm{\sigma}_{i})}\right]\end{split}$
(8)

We observed that it has similar performance as the joint cost function. In
Appendix C, we show the two cost functions are indeed equivalent to first
order in error. In principle, one should include a constant term
$c_{\text{ratio}}$ in the equation
$\mathcal{R}_{\text{joint}}=\mathcal{R}_{\text{KL}}+c_{\text{ratio}}\mathcal{R}_{\theta}$,
where the $c_{\text{ratio}}$ is a hyper-parameter to be tuned controlling the
ratio between the two cost functions. We consider $c_{\text{ratio}}=1$ because
the equivalence to $\mathcal{R}_{\text{neg. overlap}}$ in the small error
limit.

In passing we mention, that several different proposals for the cost functions
have been discussed in the literature. This includes the $\mathcal{L}_{2}$
distance of the probability amplitudes [45, 75], the weighted
$\mathcal{L}_{2}$ distance of log of the probability amplitudes [76], and
negative log fidelity [77]. We discuss the detail of the difference between
different cost functions in Appendix C

Note that the variational algorithm to approximate quantum states using neural
networks has found applications in imaginary-time evolution [75], real-time
evolution [62], simulating open quantum system [50], classical simulation of
quantum computation [77] and the study of the expressivity and generalization
properties [45, 78, 79, 76].

##  III Numerical Results

To test the variational power of different networks, we consider far-from-
equilibrium dynamics resulting from global quenches for a paradigmatic 1D
model. The fast (ballistic) growth of entanglement prohibits the classical
simulation based on MPS at long times. Recent work with NQS suggests that NQS
could reach time scales comparable with or exceeding those of the state-of-
the-art tensor network methods [31]. Here, we consider a similar setup with a
$x$-polarized initial product state
$|\psi_{0}\rangle=\prod_{i}|\rightarrow\rangle_{i}$. We then evolve the state
using the Hamiltonian of the quantum Ising model with transverse $(g)$,
longitudinal $(h)$ fields and an interaction term $(k)$:

$\hat{H}=-J\left[\sum_{j=1}^{N-1}\left(\hat{\sigma}^{z}_{j}\hat{\sigma}^{z}_{j+1}+k\hat{\sigma}^{x}_{j}\hat{\sigma}^{x}_{j+1}\right)+\sum_{j=1}^{N}\left(g\hat{\sigma}^{x}_{j}+h\hat{\sigma}^{z}_{j}\right)\right].$
(9)

We consider two cases that exhibit different growth of entanglement: a weak
quench to $g=3$, $h=0.1$, $k=0$ with a slow growth of entanglement and a
strong quench near the critical point $g=1$, $h=0$, $k=0.25$ with fast growth.
Note that the model is non-integrable for both parameter sets. In the
following, we consider open boundary condition (OBC) and provide additional
data for periodic boundary condition (PBC) in Appendix E.

###  III.1 Approximation with MPS

Figure 2: MPS with bond dimension $\chi$ approximating the exact time-evolved
states following a quantum quench from an initial paramagnetic product state.
The MPS are obtained by supervised learning with SGD. The quenched Hamiltonian
is given in Eq. (9) for a chain of length $N=20$ with (a) quantum Ising model
in the paramagnetic phase with weak longitudinal field ($g=3$, $h=0.1$, $k=0$)
and (b) close to critical point with interacting term ($g=1$, $h=0$,
$k=0.25$).

We first variationally approximate the exact time-evolved state using the MPS
ansatz with small bond dimension $\chi$ and reproduce the well known result
that the number of parameters scales exponentially in time. In particular, we
perform the exact simulation of time evolution following the quench to obtain
the target states $|\Phi_{\text{exact}}(t)\rangle$ at certain discrete times
using ED. We then find the MPS with given bond dimension $\chi$ approximating
the target state at each time $t$ by SGD over the cost function
$\mathcal{R}_{\text{joint}}$. Note that the MPS found by supervised learning
with SGD method is close to the optimal result obtained by iteratively
minimizing the $\mathcal{L}_{2}$ distance of the quantum states [14]. We
include this additional result in Appendix E for comparison. The details of
the parameters used in the optimization are shown in Appendix D.

In Fig. 2 we show the error in fidelity after the quantum state approximation
for the two parameter sets of the Ising model. We also show the expectation
value $\langle\hat{\sigma}_{x}\rangle$ in the middle of the chain and the
half-chain von Neumann entanglement entropy $S_{vN}$. For a fixed bond
dimension $\chi$, an MPS can accurately approximate the state under time
evolution to a certain timescale after which the error in fidelity grows
exponentially to a saturation value. This timescale coincides with the
deviation in the local observable $\langle\hat{\sigma}_{x}\rangle$ and the
saturation of entanglement.

To quantitatively access this time scale, we determine the reachable time
$t^{*}$ by the time when the error in fidelity exceeds a threshold value
$1-\mathcal{F}=10^{-3}$. This is indicated by the grey dashed lines in Fig. 2.
We then plot the number of parameters of the MPS as a function of the time
$t^{*}$ that can be reached in Fig. 4. The data shows a clear exponential
growth of the number of parameters as a function of $t^{*}$. This is directly
related to the linear growth of the entanglement $S$ entropy [80], which
yields an exponentially growing bond dimension $\chi\sim e^{S}$ to correctly
capture the underlying states [81, 82].

Note that the exponential scaling of the number of parameters with $t^{*}$ is
not affected by the threshold value chosen as long as the threshold value for
the error in fidelity $1-\mathcal{F}$ is chosen below the value where the
curves flatten due to finite size effects.

###  III.2 Approximation with shallow NAQS

We now investigate the variational power of NAQS in a similar setup as before.
Namely, given the exact time-evolved state $\ket{\Phi(t)}$ at the selected
time $t$, we perform supervised learning to approximate the state
variationally with NAQS. As a first example, we consider a one-hidden layer
convolutional (CNN-1) NAQS. The ratio of the number of hidden neurons to the
system size $\alpha=N_{\text{hidden}}/N$ controls the width of the network.
For CNNs, $\alpha$ is the number of channels. Increasing the width of the
network increase the number of parameters and the variational power of the
network. We optimize the cost function described in Section II.3 using SGD
with the Adam [83] optimizer. We begin with a batch size of $512$ samples and
a learning rate of $10^{-3}$ which is then decreased to $10^{-4}$ once the
sampled cost converges. Technical details are described in Appendix A.

Figure 3: One-hidden layer convolutional NAQS of different width $\alpha$
approximating the exact time-evolved states following a quantum quench from an
initial paramagnetic product state. The quenched Hamiltonian is given in Eq.
(9) for a chain of length $N=20$ with (a) quantum Ising model in the
paramagnetic phase with weak longitudinal field ($g=3$, $h=0.1$, $k=0$) and
(b) close to critical point with interacting term ($g=1$, $h=0$, $k=0.25$).

In Fig. 3, we show the results obtained for the quantum Ising model described
in Eq. (9). We observe for a fixed network that the accuracy of the
approximation decreases with time. The error increases exponentially and
saturates at a size dependent final value as in the case of MPS. Moreover, we
find that the accuracy improves when increasing the width $\alpha$ of the
networks. Identically to the procedure for MPS, we determine the reachable
time $t^{*}$ by the time when the error in fidelity exceeds the threshold
value $1-\mathcal{F}=10^{-3}$ indicated by the grey dashed lines. We notice
that the entanglement entropy also deviates from the exact value but does not
saturate even at long times after the accessible time $t^{*}$. This result is
plausible since NAQS can indeed represent volume law states.

Figure 4:  Effect of increasing network width: The number of parameters
required to reach time $t^{*}$ for following a quantum quench for different
approximation schemes. The quenched Hamiltonian is given in Eq. (9) for a
chain of length $N=20$ with (a) quantum Ising model in paramagnetic phase with
weak longitudinal field ($g=3$, $h=0.1$, $k=0$) and (b) close to critical
point with interacting term ($g=1$, $h=0$, $k=0.25$). The black line with star
(cross) marker shows the results for MPS (RBM). For NAQS, we denote each
combination as (network type)-(number of hidden layers)-(filter size). The
circles, squares, triangles represent 1,2,3-layer networks respectively. The
gradient colors describe same networks of different number of hidden layers.
NNs (Blue); CNNs with global connectivity (Green-Blue); CNNs with local
connectivity (Red).

In Fig. 4, we plot the number of parameters against this reachable time
$t^{*}$ with various widths $\alpha$. The result suggests an exponential
scaling of the number of parameters with the reachable time $t^{*}$. This is
identical to the scaling behavior of MPS but the growth constant is in fact
larger.

A big advantage of NAQS is the flexibility in the design of network
architectures. Essentially, NAQS are broad classes of wavefunctions based on
different network architectures. With a fixed network architecture, two common
factors to control the model complexity are the network width and network
depth [84, 85, 86]. An important question is how these choices affect the
network expressivity, hence the quality of the approximation.

To understand the effect of the network width on the expressivity, we focus on
shallow NNs with at most 3-hidden layers and consider three types of network
architectures: (i) fully-connected neural networks (NNs), (ii) CNNs with a
large filter size equaling to the system size, and (iii) CNNs with small
filter sizes but have receptive fields (causal cone) covering the full system.
For simplicity, we fix the same network width $\alpha$ over the hidden layers.
For a given combination of the network type and the number of hidden layers,
we repeat the procedure of state approximation described above with a
increasing network width $\alpha$. We plot the number of parameters of the
network against the obtained accessible times $t^{*}$ in Fig. 4. We denote
each combination as (network type)-(number of hidden layers) and with (filter
size) if a CNN is considered.

From the data, we observe consistent exponential scaling of the number of
parameters with the reachable time $t^{*}$ by increasing the width $\alpha$
for all network architectures and depths considered. Changing the network
architecture results in a constant shift and a slight change in the growth
rate in the scaling. There is a significant difference between NNs and CNNs,
where we see an order of magnitude improvement likely due to the sharing of
parameters. The difference between CNNs with small and large filter sizes is,
however, not obvious. While here we focus on shallow neural networks, the data
still shows a slight improvement in the growth rate when increasing the
network depth. Notably, the scaling of NNs deviates more from that of MPS for
parameter set (b). It is different from the naive expectation that MPS would
perform less efficiently when the entanglement grows rapidly.

We would like to point out that the scaling does not depend on boundary
conditions—the same exponential scaling is observed for systems with PBC (see
Appendix E). Due to the causality requirement of the NAQS, we have zero-
padding instead of periodic padding for CNNs. This implies that the
wavefunction is not translational invariant; hence, it can be applied to both
OBC and PBC systems. While we found a significant constant shift in scaling
between CNNs and NNs due to parameter sharing, we expect an even more
significant shift for other types of NQS where one enforces exact
translational symmetry.

###  III.3 Approximation with deep NAQS

To study the effect of network depth on expressivity, we consider two types of
deep network architectures: (i) CNNs with gated activation unit and residual
connection (Gated-CNNs) (ii) WaveNets [87]. Both networks consist of modular
building blocks. One can systematically increase the network depth by stacking
more blocks. The building block of Gated-CNNs includes causal convolution
following a gated activation unit and a residual connection [84]. The building
block of WaveNets replaces the standard causal convolution with dilated causal
convolution, which allows the size of the receptive field (causal cone) to
grow exponentially with the network depth, and it includes an additional
parameterized skip connection. We show the building block and the network
architecture in Fig 7 and describe the detail of the network architectures in
Appendix A. For a fixed network architecture and width $\alpha$, we increase
the network depth to increase the number of parameters and repeat the
procedure of state approximation. We vary the number of blocks between $6$ to
$14$ for WaveNets and $10$ to $20$ for Gated-CNNs.

In Fig. 5, we show the results for Gated-CNNs of width $\alpha=12$. We find a
similar result as in Fig. 3. With a larger network, it is now obvious that the
entanglement entropy deviates but does not saturate long after the accessible
time. We even observe slightly larger entanglement entropy at the deviation
(See also data for WaveNets in Appendix E).

In Fig. 6, we plot the number of parameters of networks with different depths
against the accessible times $t^{*}$ and denote each combination as (network
type)-(network width $\alpha$). We observe again consistent exponential
scaling of the number of parameters with the reachable time $t^{*}$ by
increasing the network depth for the two networks and widths considered.
Similarly, we see a small constant shift for different network width favoring
“narrower” networks. More importantly, data for both network architectures
roughly fall on the same line despite the drastically different design choice
of using dilated or normal convolution. Compared with the data from shallow
CNNs where we vary the network width, all data again fall roughly on the same
line. This suggests that increasing the number of parameters by increasing
network width has the same effect on expressivity as increasing the network
depth for the cases considered. Even more surprisingly, such results are
independent of the network architectures of whether one considers global or
local, dilated or normal convolution, with or without skip and residual
connections. Note that, however, the optimization problem may be strongly
affected by the choice of different network architectures. For example, a deep
neural network without shortcut connections may be hard to optimize.

Figure 5: Gated-CNNs NAQS of different depth, i.e. number of blocks $N_{b}$,
approximating the exact time-evolved states following a quantum quench from an
initial paramagnetic product state. The network is of width $\alpha=12$. The
quenched Hamiltonian is given in Eq. (9) for a chain of length $N=20$ with (a)
quantum Ising model in the paramagnetic phase with weak longitudinal field
($g=3$, $h=0.1$, $k=0$) and (b) close to critical point with interacting term
($g=1$, $h=0$, $k=0.25$).  Figure 6:  Effect of increasing network depth: The
number of parameters required to reach time $t^{*}$ for following a quantum
quench for following a quantum quench for different approximation schemes. The
quenched Hamiltonian is given in Eq. (9) for a chain of length $N=20$ with (a)
quantum Ising model in paramagnetic phase with weak longitudinal field ($g=3$,
$h=0.1$, $k=0$) and (b) close to critical point with interacting term ($g=1$,
$h=0$, $k=0.25$). The black line with star (cross) marker shows the results
for MPS (RBM). For NAQS, we denote each combination as (network type)-(network
width $\alpha$). The diamonds and hexagons represent network of width
$\alpha=8,12$ respectively. The gradient colors describe same networks of
different network width. Gated-CNNs (Red); WaveNets (Green); we keep the
shallow NAQS data (light grey) from Fig. 4 for comparison.

We further note that there are two main differences in the network structure
considered here and in previous works [30, 65, 31, 62]. The first difference
is the autoregressive property. It is unknown how the autoregressive
constraint affects the expressivity of the neural networks. The second
difference is the activation function. In Ref. 31, 62, the choice of
polynomial activation with complex-valued weight is crucial to the result.
Notice that such networks [42] do not have the universal approximation
property. In contrast, we consider real-valued weights with rectified linear
unit (ReLU) [88], tanh, and sigmoid activation functions. For shallow NNs, we
observe that changing activation function does not affect the result
significantly. The difference in the setup calls into question whether the
good result before [31] is not a general property of neural networks.

###  III.4  Approximation with RBMQS

As the last example, we study the expressivity of the restricted Boltzmann
machine quantum states (RBMQS) [30] with increasing width. One prominent
difference between RBMQS and NAQS is that RBMQS belongs to general (non-
autoregressive) neural network, hence is not normalized. The normalization
constant is in general intractable. Nevertheless, it is possible to estimate
the fidelity to another wavefunction by stochastic sampling. The gradient of
such stochastic estimate gives an biased stochastic gradient, which works well
in practice for supervised learning with stochastic gradient methods.

We repeat the similar benchmark of state approximation on networks of
increasing width $\alpha=N_{\text{hidden}}/N=\\{1,2,4,8\\}$. We minimize the
negative fidelity between the RBMQS and the target state by SGD. We confirm
our implementation yields similar result to the NetKet [78] implementation
using stochastic reconfiguration (SR) method. For more detail, see the
discussion of the cost function in Appendix C and the additional data for
consistency check in Appendix E.

The result for state approximation using RBMQS is shown in Fig. 12 in Appendix
E, which resembles the results of NAQS shown in Fig. 3 and Fig. 5. We plot the
reachable time $t^{*}$ in Fig. 4 and Fig. 6. We observe again consistent
exponential scaling of the number of parameters with the reachable time
$t^{*}$ by increasing the network width.

We would like to point out our observation in the reachable time is consistent
with previous works [30, 65] which show similar difficulties around $Jt=0.5$
for representing states quenched by Hamiltonian near the critical point using
RBM-based NQS. It is also observed in [62] that a lack of expressibility from
the RBM-based NQS leads to a rapid increase in the error in fidelity in a
short time scale for a quench across critical point.

Lastly, similar scaling plots as in Fig. 4 and Fig. 6 could also be obtained
in a more realistic setup when exact data is not known. For MPS, one could
perform the time-evolving block decimation (TEBD) algorithm [26, 27], and only
increase the bond dimension $\chi$ when the approximation (truncation) error
reaches the threshold value $10^{-3}$. For NQS, one could perform a time
evolution algorithm based on supervised learning [62, 50] and increase the
size of the neural network when the approximation error reaches the threshold
value $10^{-3}$. If we plot the number of parameters of the variational
wavefunctions to the time steps, we could obtain a similar scaling plot.

##  IV  Conclusion

We provide numerical evidence that the required number of parameters to
represent quantum states following a global quench by neural networks grows
exponentially in time. Thus, despite the ability to represent highly-entangled
states, the neural network quantum states (NQS), including neural
autoregressive quantum states (NAQS) and restricted Boltzmann machine quantum
states (RBMQS), considered in this work cannot represent time-evolved states
efficiently. An important remaining question is to understand this scaling.
The scaling is agnostic to the design of network architectures we considered
but is only affected by the difference in architecture related to symmetry,
i.e. convolution. The scaling of NQS resembles that of MPS. However, in the
one-dimensional models considered, NQS shows a larger growth rate of
parameters in time while its entanglement entropy saturates slower. Explaining
this difference and understanding the limitation of NQS may require a new
measure for the complexity of states different from the entanglement entropy.

The proposed cost function for joint learning magnitudes and phases with NAQS
ensures unbiased and low-variance stochastic gradient. It could be applied to
the real-time evolution for pure states and mix states [62, 50] and the
imaginary-time evolution [75]. We expect it to speed up the state
approximation at each time step and to be suitable for learning with deep
neural networks. Monitoring the training and convergence of the two costs
$\mathcal{R}_{\text{KL}}$ and $\mathcal{R}_{\theta}$ gives information about
the hardness of learning probability distributions and the phases. We observe
similar difficulty in learning both the magnitudes and the phases for time-
evolved states. This may provide additional insight for learning frustrated
ground states, which previous works [89, 76] show the main difficulty comes
from learning the phases.

Finally, our findings do not invalidate the simulation for time evolution with
NQS in general but rather provides insight on the choice of networks to reach
the desired time scale. It also suggests the importance to consider
dynamically increasing the network size for the time evolution algorithm based
on supervised learning. It is still important to search for different network
architectures [84, 90, 91, 92, 93, 87], which might be more efficient in
optimization or better in the scaling for different Hamiltonians. Network
architectures incorporating the symmetry of the states have the potential to
give better scaling. Moreover, the approach based on neural networks is less
affected by the spatial dimension of the systems than that based on tensor
networks. This suggests the potential advantage of NQS over tensor network
states in two and higher dimensions.

While understanding the expressivity and scaling of NQS is the essential step
for designing practical algorithms, the recent work [94] studies the numerical
instability in obtaining the NQS using time-dependent variational Monte Carlo
when the time-evolved states can be expressed as NQS accurately. The article
[94] and this work complement each other and provide insights for simulating
time evolution using NQS. The code implementation of this work is based on
TensorFlow [95] and is available at [96].

###### Acknowledgements.

S.L. would like to thank Markus Schmitt, Rohit Dilip, Adam Smith, and Jin-Guo
Liu for helpful discussion. This work is supported by the European Research
Council (ERC) under the European Union’s Horizon 2020 research and innovation
program (grant agreement No. 771537). S.L. and F.P. are supported by the
Deutsche Forschungsgemeinschaft (project number 107745057, TRR 80). F.P.
acknowledges the support of the Deutsche Forschungsgemeinschaft (DFG, German
Research Foundation) under Germany’s Excellence Strategy EXC-2111-390814868.

##  Appendix A Network Architectures

Feed-forward neural networks are functions consisting of alternating affine
maps and elementwise non-linear functions. The input of the network is the
one-hot encoding of the computational basis. That is

$\displaystyle X:\bm{\sigma}\in\\{1,2,\ldots,M\\}^{N}\longrightarrow
X(\bm{\sigma})\in\mathbb{R}^{M\times N}$ (10) $\displaystyle
X(\bm{\sigma})_{i,j}=\begin{cases}1,\qquad\text{if}\ i=\sigma_{j}\\\
0,\qquad\text{else}\end{cases}$ (11)

Fully-connected neural networks (NNs) are defined by a recursion relation at
each layer

$\textbf{z}^{(l)}=g({W}^{(l)}\textbf{z}^{(l-1)}+\textbf{b}^{(l)}),$ (12)

where the weight matrix $W^{(l)}$ and bias vector $\textbf{b}^{(l)}$ at each
layer $(l)$ together are the variational parameters. The $g$ is the
elementwise non-linear function, a.k.a activation function. At the first
layer, the input is first flattened as a vector, i.e.
$z^{(0)}=\text{Vec}(X(\bm{\sigma}))$.

We consider NNs with one- and two-hidden layers and denote them as NN-1 and
NN-2. The number of neurons are $(N,\alpha N,4N)$ for NN-1 and $(N,\alpha
N,\alpha N,4N)$ for NN-2 respectively. The real-valued output is of dimension
$4N$ to encode the $M=2$ different complex-valued conditional probability
amplitude on $N$ sites.

Similarly, convolutional neural networks (CNNs) are defined by a recursion
relation at each layer

${z}^{(l)}=g({W}^{(l)}\ast{z}^{(l-1)}+{b}^{(l)}),$ (13)

where now we replace the matrix multiplication with the convolution operation
$*$. To be explicit, we can write it as

${z}^{(l)}_{c,j}=g(\sum_{c^{\prime}=1}^{\alpha
N}\sum_{k=1}^{N_{F}}{W}^{(l)}_{c,c^{\prime},k}{z}^{(l-1)}_{c^{\prime},(j+k-N_{F})}+{b}^{(l)}_{c,j}),$
(14)

where $N_{F}$ denotes the filter size. We consider two different setups of
global connectivity $N_{F}=N$ and local connectivity $N_{F}<N$. We pad the
input with zeros to keep the width of the network fixed for simplicity and
only change the number of channels at each layer. For shallow networks, we
consider CNNs up to three hidden layers and denote them as CNN-[depth]. The
number of channels are $(2,\alpha,4)$ for CNN-1, $(2,\alpha,\alpha,4)$ for
CNN-2, and $(2,\alpha,\alpha,\alpha,4)$ for CNN-3. The output layer is with 4
channels to represent the $M=2$ different complex-valued conditional
probability amplitudes. All shallow NNs and CNNs have ReLU as activation
function, i.e. $g(z)=\text{max}(0,z)$.

To turn the networks into networks satisfying autoregressive properties, the
connectivity is restricted by masking the weights. For example, the type-A
masking makes sure the $z^{(l)}_{i}$ does not depend on $z^{(l-1)}_{i}$ and
component afterwards by restricting the sum to
${z}^{(l)}_{i}=g(\sum_{j<i}{W}^{(l)}_{i,j}{z}^{(l-1)}_{j}+{b}^{(l)}_{i})$. The
type-B masking makes sure the $z^{(l)}_{i}$ does not depend on
$z^{(l-1)}_{i+1}$ and component afterwards by restricting the sum to
${z}^{(l)}_{i}=g(\sum_{j\leq i}{W}^{(l)}_{i,j}{z}^{(l-1)}_{j}+{b}^{(l)}_{i})$.
Combining the type-A masking at the first layer and type-B masking at all
layers afterward, we ensure the autoregressive properties of the network. For
more detail about autoregressive models see Appendix B and see Fig. 1 for a
graphical illustration.

Figure 7:  The building block and network architecture of WaveNet-like network
adapted from [87]. The building block is indicated by the dotted box. We
repeat the building block to increase the depth. Except for the first and the
last block, the output of the residual connection is the input of the next
following block. The Gated-CNN is a simplification of the WaveNet without
skip-connections and the following layers. Additionally, it uses convolution
instead of dilated convolution.

We consider two different deep neural networks. In Fig. 7, we show the network
adapted from the WaveNet [87]. The WaveNet uses dilated convolution, which is
the convolution applying over a larger area by skipping input values. The
stacked dilated convolutions with exponential growing dilation factors enable
exponential growing receptive field (causal cone). Similar to the WaveNet, we
double the dilation factors at each layer and repeat again once it reach the
system size, e.g. $1,2,4,8,16,1,2,4,\cdots$ for $N=20$.

Following the dilated convolution is the gated activation unit [92],

$z=\text{tanh}(W^{f}*x)\odot\sigma(W^{g}*x)$ (15)

where $\odot$ is the element-wise multiplication. The $\sigma(\cdot)$ is the
sigmoid activation function and should not be confused with the local
computation basis. The $W^{f}$ and $W^{g}$ are the weight for filter and gate
respectively.

After the gated activation unit, there are two special connections. The
residual connections [84] are implemented by the $1\times 1$ convolution and
the summation with the shortcut from the input of the dilated convolution. The
output of the residual connection is the input of the next block. The
parameterized skip connections are implemented by the $1\times 1$ convolution
to form shortcuts from each block directly to the end of repeating blocks,
which are then summed together with the final output. These “shortcut”
connections are the key to training deep neural networks.

In addition to the WaveNet, we consider a simple deep CNN, Gated-CNN, where we
remove the skip connections and use only convolution instead of the dilated
convolution. Without the dilated convolution, the receptive field of a Gated-
CNN grows linear with depth. For the receptive field to cover the system size
$N=20$, we consider Gated-CNNs with $10$ to $20$ blocks. As for WaveNets, we
consider having $6$ to $14$ blocks. While WaveNets include one block of
$1\times 1$ convolution before the final causal convolution, the output of the
Gated-CNN blocks is followed directly by the causal convolution. The design
choice is arbitrary and turns out not to affect the scaling result. Similar to
shallow CNNs, the final output of both networks has 4 channels to represent
the $M=2$ different complex-valued conditional probability amplitudes.

The code implementation based on TensorFlow [95] is available [96] and the
hyperparameters and general setup of the optimization can be found in Appendix
D.

##  Appendix B Review of Neural Autoregressive Quantum States

A joint probability distribution $P$ over $N$ discrete random variables
$x_{i}$ can be factorized as the multiplication of conditional probabilities.

$\displaystyle P(x_{1},x_{2},\ldots,x_{N})$ $\displaystyle=P(x_{1})\times
P(x_{2}|x_{1})\times P(x_{3}|x_{1},x_{2})\ldots$
$\displaystyle=\prod_{i=1}^{N}p(x_{i}|x_{1},x_{2},\ldots,x_{i-1})$
$\displaystyle=\prod_{i=1}^{N}p(x_{i}|{\bf x}_{<i})$

The decomposition imposes an ordering over the random variables and allows for
exact sampling. However, the exact representation of conditional probabilities
also scales exponentially with the number of random variables. Autoregressive
models are parameterized functions approximating the conditional
probabilities.

One important class are neural autoregressive models [97, 98]. By utilizing
neural networks to approximate the conditional probabilities, it allows exact
and efficient sampling and inference. It has been applied to studying
statistical physics in the pioneering work [99] and in [100].

In principle, one can approximate each conditional probability with a
different neural network,

$P_{\text{AR}}=f(x_{1})\times g(x_{2}|x_{1})\times
h(x_{3}|x_{1},x_{2})\times\ldots$

A simpler and common approach is to approximate all $N$ conditional
probability by one neural network with $N$ inputs and $N$ outputs,

$P_{\text{AR}}=f^{[1]}(x_{1})\times f^{[2]}(x_{2}|x_{1})\times
f^{[3]}(x_{3}|x_{1},x_{2})\times\ldots$

The network can no longer be fully-connected, but with connectivity satisfying
the autoregressive property. That is $f^{[2]}$ can only depend on input
$x_{1}$ and $f^{[3]}$ only on inputs $x_{1},x_{2}$ and so on. It is common to
enforce this by zero-masking the original weights $w$. Equivalently, we say we
have masked weight multiplication or masked convolution (see Fig. 1).

Following the same principle, a probability amplitude can be rewritten as a
product of conditional probability amplitudes [101].

$\Psi_{\sigma_{1}\sigma_{2}\cdots\sigma_{N}}=\Psi_{\sigma_{1}}\Psi_{\sigma_{2}|\sigma_{1}}\cdots\Psi_{\sigma_{N}|\sigma_{1}\cdots\sigma_{N-1}}$
(16)

Neural autoregressive quantum states (NAQS) [32] approximate the conditional
probability amplitudes by neural networks.

$\Psi^{NAQS}_{\sigma_{1}\sigma_{2}\cdots\sigma_{N}}=f_{\text{NN}}^{[1]}(\sigma_{1})f_{\text{NN}}^{[2]}(\sigma_{2}|\sigma_{1})\cdots
f_{\text{NN}}^{[N]}(\sigma_{N}|\sigma_{1}\cdots\sigma_{N-1})$ (17)

In the following, we denote $f^{[i]}_{\text{NN}}$ as $\psi_{i}$.

It is shown that if each conditional probability amplitude is normalized
$\sum_{\sigma_{i}}|\Psi_{\sigma_{i}|\sigma_{j<i}}|^{2}=1$, then the full
wavefunction is normalized
$|\Psi_{\sigma_{1}\sigma_{2}\cdots\sigma_{N}}|^{2}=1$. Therefore, we can
normalize the full wavefunction by normalizing each conditional probability
amplitude. Suppose the local basis is of dimension $M$, each of the $N$
outputs is actually a complex vector
$\bm{v_{i}}=(v_{i,1},v_{i,2},\ldots,v_{i,M})$, which represents the logarithm
of the conditional probability amplitude at site $i$.

$\psi_{i}(\sigma_{i}|\sigma_{1},\ldots,\sigma_{i-1})=\frac{e^{v_{i,\sigma_{i}}}}{\sum_{\sigma_{i}^{\prime}}|e^{v_{i,\sigma_{i}^{\prime}}}|^{2}}$
(18)

gives the normalized conditional probability amplitudes. Furthermore, each
$|\psi_{i}|^{2}$ gives a valid conditional probability induced by Born’s rule.
This means a NAQS contains a neural autoregressive model, which permits
efficient and exact sampling.

##  Appendix C Cost Functions and Stochastic Gradient Descent

In this section, we give a brief review of the setup of supervised learning
problems and discuss the difference between it and the general optimization
problem. We further give an overview of different cost functions for
variational states approximation and the resulting properties. The discussion
is general and applies to general variational wavefunctions.

Supervised learning is a minimization problem to learn a function mapping from
$X$ to $Y$. Given the data $\\{x,y\\}^{N}$, where the data pair $\\{x,y\\}\sim
P(x,y)$ are sampled from an unknown probability distribution $P(x,y)$, and the
parameterized function $f:x\in X\rightarrow f(x)\in Y$ with parameters $w$, we
would like to minimize the cost function $\mathcal{R}(f)$ with respect to
parameters $w$. $\mathcal{R}(f)$ is defined as

$\mathcal{R}(f)=\int\mathcal{L}(f(x;w),y)dP(x,y)$ (19)

where $\mathcal{L}(f(x;w),y)$ is the loss function measuring the difference
between $f(x;w)$ and $y$. In general, the expression cannot be evaluated.
Instead, the supervised learning framework, a.k.a. empirical risk
minimization, relies on minimizing the cost from the sampling (empirical risk)
$\hat{\mathcal{R}}(f)$,

$\hat{\mathcal{R}}(f)=\sum_{\\{x,y\\}\sim P(x,y)}\mathcal{L}(f(x;w),y).$ (20)

One can compute the full batch gradient based on the empirical risk and
perform a gradient descent algorithm. In modern applications, often even the
empirical risk itself is too expensive to evaluate at all steps given the
large size of the dataset. One turns to compute the mini-batch stochastic
gradient
$\widehat{\nabla\mathcal{R}}=\sum_{i}\nabla_{w}\mathcal{L}(f(x_{i};w),y_{i})$
where the sum is taken over the mini batch samples. In terms of deep learning,
the supervised learning framework fits nicely with the stochastic gradient
descent algorithm.

Another type of optimization problem is to minimize the cost function
$\mathcal{R}$, which is not in the form as in Eq. (19), i.e. an integral or a
sum of the loss function $\mathcal{L}$. The cost function $\mathcal{R}$ itself
is the thing we want to optimize. In this case, gradient-based algorithms
would require computing the gradient $\nabla_{w}\mathcal{R}$. In certain
cases, the gradient $\nabla_{w}\mathcal{R}$ itself cannot be computed exactly
and approximation is required. Moreover, the gradient $\nabla_{w}\mathcal{R}$
might contain intractable terms as integral or summation, which can be
approximated by sampling. In this case, the sampled “stochastic gradient” has
quite different properties comparing to the previous one in the setup of
supervised learning. More precisely, _the stochastic gradient is unbiased when
the expectation value over the sampling of this stochastic gradient equal to
the true gradient_ , i.e.

$\nabla\mathcal{R}=\mathbb{E}[\widehat{\nabla\mathcal{R}}]$

The stochastic gradient is _biased_ if the expectation value of the sampled
stochastic gradient does not equal to the true gradient. Stochastic gradients
in supervised learning setup are unbiased while in the general optimization
problem they are not. With this in mind, we can distinguish the cost function
for variational state approximation by whether it has the form of Eq. (19).

In terms of the variational state approximation, given the target wavefunction
$|\Phi\rangle$, the variational wavefunction $|\Psi\rangle$ and the dimension
of the many-body Hilbert space $D_{\mathcal{H}}$. here we review different
possible cost functions that work for both unnormalized wavefunctions and
normalized wavefunctions introduced in previous works.

In [45, 75], the $\mathcal{L}_{2}$ distance of the probability amplitudes
$\mathcal{R}_{\text{MSE}}$ are considered,

$\displaystyle\mathcal{R}_{\text{MSE}}$
$\displaystyle=\frac{1}{D_{\mathcal{H}}}\sum_{i}^{D_{\mathcal{H}}}|\Psi(X_{i})-\Phi(X_{i})|^{2}$
$\displaystyle=\frac{1}{D_{\mathcal{H}}}\sum_{i}^{D_{\mathcal{H}}}P_{\Phi}(X_{i})\frac{|\Psi(X_{i})-\Phi(X_{i})|^{2}}{\lvert\Phi(X_{i})\rvert^{2}}.$

Similarly, in Ref. 76, the weighted $\mathcal{L}_{2}$ distance of log of the
probability amplitudes is considered,

$\displaystyle\mathcal{R}_{\text{weighted MSE\ log}}$
$\displaystyle=\frac{1}{D_{\mathcal{H}}}\sum_{i}^{D_{\mathcal{H}}}P_{\Phi}(X_{i})\Big{|}\log\Psi(X_{i})-\log\Phi(X_{i})\Big{|}^{2}$

while in [45] a uniform sampling scheme is considered, and in Ref. 75, 76 the
samples are taken with respect to $P_{\Phi}$. These cost functions are
unbiased but with high variance. Furthermore, the problem is posed as a
regression problem instead of taking into account the probability nature of
the problem.

A different setup is considered in Ref. 77 with the negative log fidelity as
the cost function.

$\displaystyle\mathcal{R}_{\text{neg. log
fidelity}}=-\log\Big{[}\frac{\braket{\Psi}{\Phi}\braket{\Phi}{\Psi}}{\braket{\Psi}{\Psi}\braket{\Phi}{\Phi}}\Big{]}$
$\displaystyle=-\log\Bigg{[}\frac{\Big{(}\sum_{i}^{D_{\mathcal{H}}}\Psi^{*}(X_{i})\Phi(X_{i})\Big{)}\Big{(}\sum_{i}^{D_{\mathcal{H}}}\Phi^{*}(X_{i})\Psi(X_{i})\Big{)}}{\Big{(}\sum_{i}^{D_{\mathcal{H}}}\Psi^{*}(X_{i})\Psi(X_{i})\Big{)}\Big{(}\sum_{i}^{D_{\mathcal{H}}}\Phi^{*}(X_{i})\Phi(X_{i})\Big{)}}\Bigg{]}$

This cost function and the corresponding gradient are in the form that could
be estimated by stochastic sampling. However, the resulting stochastic
gradient is biased. In practice, this cost function works well as long as the
batch size is taken large enough. In this work, we use directly the negative
fidelity as the cost function when optimizing the RBMQS.

$\displaystyle\mathcal{R}_{\text{neg. fidelity}}$
$\displaystyle=-\frac{(\sum_{i}\Psi_{i}^{*}\Phi_{i})(\sum_{i}\Psi_{i}\Phi^{*}_{i})}{(\sum_{i}\Psi_{i}^{*}\Psi_{i})(\sum_{i}\Phi_{i}^{*}\Phi_{i})}$
$\displaystyle=-\frac{\mathbb{E}_{i\sim
P_{\Phi}}[\Psi^{*}_{i}/\Phi_{i}^{*}]\mathbb{E}_{i\sim
P_{\Phi}}[\Psi_{i}/\Phi_{i}]}{\mathbb{E}_{i\sim
P_{\Phi}}[|\Psi_{i}|^{2}/|\Phi_{i}|^{2}]}$

With the stochastic estimate of the fidelity, the biased stochastic gradient
is then obtained from the gradient evaluated by automatic differentiation
framework.

In this work, we consider two cost functions that work only for normalized
wavefunction, which give unbiased stochastic gradient with low variance. In
Section II.3, we mainly focus on the cost function

$\mathcal{R}_{\text{joint}}=\mathcal{R}_{\text{KL}}+\mathcal{R}_{\theta}.$
(21)

An alternative cost function is the real part of the negative overlap,

$\displaystyle\mathcal{R}_{\text{neg. overlap}}$
$\displaystyle=-\sum_{i}^{D_{\mathcal{H}}}\text{Re}\left[\braket{\Phi}{X_{i}}\braket{X_{i}}{\Psi}\right]$
$\displaystyle=-\sum_{i}P_{\Phi}(X_{i})\text{Re}\left[\frac{\braket{X_{i}}{\Psi}}{\braket{X_{i}}{\Phi}}\right]$

Minimizing the real part $\text{Re}\left[\braket{\Phi}{\Psi}\right]$ is
equivalent to minimizing $\lvert\braket{\Phi}{\Psi}\rvert$ because we can
always absorb a phase factor into our complex-valued wavefunction. Minimizing
only the real part gives the form of expectation value and again we can apply
the stochastic gradient descent method. This, however, cannot apply to
unnormalized wavefunction because of the intractable normalization constant.

The fidelity $\mathcal{F}$ is the square of the overlap,
$\mathcal{F}=\lvert\braket{\Phi}{\Psi}\rvert^{2}$. It may be tempting to
maximize the fidelity $\mathcal{F}$ instead by the sampling estimate of the
overlap. However, again we encounter the problem of a biased gradient.
Although, the fidelity estimate cannot serve as a cost function. We can still
evaluate it as a check of the validity of the final result.

Finally, we show that these two cost functions are equivalent to the first
order.

$\displaystyle\mathcal{R}_{\text{neg. overlap}}$
$\displaystyle=-\sum_{i}P_{\Phi}(X_{i})\text{Re}\left[\frac{\braket{X_{i}}{\Psi}}{\braket{X_{i}}{\Phi}}\right]$
$\displaystyle=-\sum_{i}P_{\Phi}(X_{i})\text{Re}\left[\frac{e^{r^{\Psi}_{i}+i\theta^{\Psi}_{i}}}{e^{r^{\Phi}_{i}+i\theta^{\Phi}_{i}}}\right]$
$\displaystyle=-\sum_{i}P_{\Phi}(X_{i})\text{Re}\left[e^{(r^{\Psi}_{i}-r^{\Phi}_{i})+i(\theta^{\Psi}_{i}-\theta^{\Phi}_{i})}\right]$
$\displaystyle\approx-\sum_{i}P_{\Phi}(X_{i})\left[(1+(r^{\Psi}_{i}-r^{\Phi}_{i}))\right(1+\frac{(\Delta\theta_{i})^{2}}{2}\left)\right]$
$\displaystyle\approx\sum_{i}P_{\Phi}(X_{i})\left[-1+(r^{\Phi}_{i}-r^{\Psi}_{i})-\frac{(\Delta\theta_{i})^{2}}{2}\right]$
$\displaystyle\approx-1+\frac{1}{2}(\mathcal{R}_{\text{KL}}+\mathcal{R}_{\theta})$

In practice, we observe slightly different dynamics at the beginning of the
learning process.

##  Appendix D Setup of Optimizations

Below, we first discuss the common setup of optimizations for both NAQS and
MPS.

  * •

_Optimizer and batch size:_ We consider Adam optimizer and a batch size of
$512$.

  * •

_Convergence criterion:_ The convergence criterion is based on monitoring the
error in fidelity. We set the criterion to be either if the error were below
$10^{-4}$ or if average error of the latest $500$ batches were higher than the
average error of the latest $5000$ to $2500$ batches. In addition, there is a
break condition set by the maximal number of steps.

  * •

_Maximal number of steps:_ We set the maximal number of steps to be $3\times
10^{5}$, which is almost never reached. The typical number of steps are below
$10^{5}$.

  * •

_Number of runs_ : For each data point, we run the optimization with random
initialization $4$ times and report the best result.

_Matrix Product States:_ The MPS considered are complex-valued. We start out
with a fix learning rate $10^{-2}$ until convergence. Then, we repeat the
procedure with the learning rate decreased to $10^{-3}$ and $10^{-4}$.

Notice that the canonical (isometric) form or normalization of MPS is not
enforced during the optimization. In fact, we observed that the result is
better than the result enforcing the isometric form using Riemannian
optimization [102, 103].

_Neural Network:_ All the networks considered are parameterized by real-valued
weights and biases. We start out with a fix learning rate $10^{-3}$ until the
convergence criterion is met. Then, we decrease the learning rate to $10^{-4}$
and continue again the optimization until the convergence.

Here, we provide the setup of using NetKet to perform state approximation with
RBMQS.

_NetKet:_

  * •

_Version_ : 2.1

  * •

_Optimizer and batch size_ : stochastic reconfiguration and stochastic
gradient descent. The batch size is $1000$. We take the best result from
different step sizes: $1\times 10^{-1},3\times 10^{-2},1\times 10^{-2}$.

  * •

_Convergence criterion:_ The convergence criterion is based on monitoring the
average error in fidelity. We set the criterion to be either if the average
error were below $10^{-4}$ or if average error of the latest $50$ batches were
close to the average error of the latest $100$ to $50$ batches with a relative
error $3\times 10^{-6}$ and absolute error $3\times 10^{-7}$. In addition,
there is a break condition set by the maximal number of steps.

  * •

_Maximal number of steps:_ We start checking for convergence after 5000
iterations, and set a maximum iteration of $10000$ which is almost never
reached.

  * •

_Number of runs_ : For each data point, we run the optimization with random
initialization $4$ times and report the best result.

##  Appendix E Additional Data

In Fig. 8, we include the result of MPS approximating the target state by the
alternating least square algorithm, which iteratively minimizes the
$\mathcal{L}_{2}$ distance of the quantum states [14] and yield slightly
better result than SGD optimization with supervised learning as shown in Fig.
2.

Figure 8: MPS with bond dimension $\chi$ approximating the exact time-evolved
states following a quantum quench from an initial paramagnetic product state.
The MPS are obtained by iteratively minimizing the $\mathcal{L}_{2}$ distance.
The quenched Hamiltonian is given in Eq. (9) for a chain of length $N=20$ with
(a) quantum Ising model in the paramagnetic phase with weak longitudinal field
($g=3$, $h=0.1$, $k=0$) and (b) close to critical point with interacting term
($g=1$, $h=0$, $k=0.25$).

In Fig. 9, we include the results with shallow NAQS for periodic boundary
condition (PBC) and plot the number of parameters of the network against the
obtained accessible time $t^{*}$.

Figure 9:  The number of parameters required to reach time $t^{*}$ for
following a quantum quench for different approximation schemes. The quenched
Hamiltonian is similar to Eq. (9) but with PBC on a chain of length $N=20$
with (a) quantum Ising model in paramagnetic phase with weak longitudinal
field ($g=3$, $h=0.1$, $k=0$) and (b) close to critical point with interacting
term ($g=1$, $h=0$, $k=0.25$). The black line shows the results for MPS. For
NAQS, we denote each combination as (network type)-(number of hidden
layers)-(filter size). The circles, squares, triangles represents 1,2,3-layer
networks respectively. The gradient colors describe same networks of different
number of layers. NNs (Blue); CNNs with global connectivity (Green-Blue); CNNs
with local connectivity (Red).

In Fig. 10, we show the results for WaveNets of width $\alpha=8$.

In Fig. 11, we plot the errors of the results for Gated-CNNs of width
$\alpha=12$ and different depths. We see the errors in
$\langle\hat{\sigma}^{x}_{L/2}\rangle$ expectation value and half-chain
entanglement entropy $S_{\text{vN}}$ follows similar tendency as of the error
in fidelity. Interestingly, we observe a peak in both errors at a relatively
short time before the network reaches its accessible time $t^{*}$.

Figure 10: WaveNet NAQS of different depths, i.e. number of blocks $N_{b}$,
approximating the exact time-evolved states following a quantum quench from an
initial paramagnetic product state. The network is of width $\alpha=8$. The
quenched Hamiltonian is given in Eq. (9) for a chain of length $N=20$ with (a)
quantum Ising model in the paramagnetic phase with weak longitudinal field
($g=3$, $h=0.1$, $k=0$) and (b) close to critical point with interacting term
($g=1$, $h=0$, $k=0.25$).  Figure 11: The errors of Gated-CNNs NAQS of
different depth, i.e. number of blocks $N_{b}$, approximating the exact time-
evolved states following a quantum quench from an initial paramagnetic product
state. The network is of width $\alpha=12$. The quenched Hamiltonian is given
in Eq. (9) for a chain of length $N=20$ with (a) quantum Ising model in the
paramagnetic phase with weak longitudinal field ($g=3$, $h=0.1$, $k=0$) and
(b) close to critical point with interacting term ($g=1$, $h=0$, $k=0.25$).

In Fig. 12, we show the results for RBMQS.

In Fig, 13, we compare the results for state approximation by RBMQS using SGD
and exact gradient optimization and the SGD with stochastic reconfiguration
from the NetKet [78] implementation.

In Fig, 14, we compare the results for state approximation by one-hidden layer
NAQS using SGD and GD over negative fidelity and using SGD over joint cost
function.

Figure 12:  RBMQS of different widths, i.e. $\alpha=N_{\text{hidden}}/N$,
approximating the exact time-evolved states following a quantum quench from an
initial paramagnetic product state. The quenched Hamiltonian is given in Eq.
(9) for a chain of length $N=20$ with (a) quantum Ising model in the
paramagnetic phase with weak longitudinal field ($g=3$, $h=0.1$, $k=0$) and
(b) close to critical point with interacting term ($g=1$, $h=0$, $k=0.25$).
Figure 13:  Comparison of results obtained from different optimizers and
implementations for RBMQS of two different widths $\alpha=2,4$ approximating
the exact time-evolved states following a quantum quench from an initial
paramagnetic product state. The quenched Hamiltonian is given in Eq. (9) for a
chain of length $N=12,16,20$ with quantum Ising model close to critical point
with interacting term ($g=1$, $h=0$, $k=0.25$).  Figure 14:  Comparison of
results obtained by different optimizers and cost functions for one-hidden
layer NAQS of two different widths $\alpha=2,4$ approximating the exact time-
evolved states following a quantum quench from an initial paramagnetic product
state. The quenched Hamiltonian is given in Eq. (9) for a chain of length
$N=12,16,20$ with quantum Ising model close to critical point with interacting
term ($g=1$, $h=0$, $k=0.25$).

## References

  * Basko  _et al._ [2006] D. M. Basko, I. L. Aleiner, and B. L. Altshuler, Metal–insulator transition in a weakly interacting many-electron system with localized single-ptticle states, Annals of physics 321, 1126 (2006).
  * Pal and Huse [2010] A. Pal and D. A. Huse, Many-body localization phase transition, Physical review b 82, 174411 (2010).
  * Schreiber  _et al._ [2015] M. Schreiber, S. S. Hodgman, P. Bordia, H. P. Lüschen, M. H. Fischer, R. Vosk, E. Altman, U. Schneider, and I. Bloch, Observation of many-body localization of interacting fermions in a quasirandom optical lattice, Science 349, 842 (2015).
  * Choi  _et al._ [2016] J.-y. Choi, S. Hild, J. Zeiher, P. Schauß, A. Rubio-Abadal, T. Yefsah, V. Khemani, D. A. Huse, I. Bloch, and C. Gross, Exploring the many-body localization transition in two dimensions, Science 352, 1547 (2016).
  * Smith  _et al._ [2016] J. Smith, A. Lee, P. Richerme, B. Neyenhuis, P. W. Hess, P. Hauke, M. Heyl, D. A. Huse, and C. Monroe, Many-body localization in a quantum simulator with programmable random disorder, Nature Physics 12, 907 (2016).
  * Bernien  _et al._ [2017] H. Bernien, S. Schwartz, A. Keesling, H. Levine, A. Omran, H. Pichler, S. Choi, A. S. Zibrov, M. Endres, M. Greiner, _et al._ , Probing many-body dynamics on a 51-atom quantum simulator, Nature 551, 579 (2017).
  * Turner  _et al._ [2018] C. J. Turner, A. A. Michailidis, D. A. Abanin, M. Serbyn, and Z. Papić, Weak ergodicity breaking from quantum many-body scars, Nature Physics 14, 745 (2018).
  * Choi  _et al._ [2019] S. Choi, C. J. Turner, H. Pichler, W. W. Ho, A. A. Michailidis, Z. Papić, M. Serbyn, M. D. Lukin, and D. A. Abanin, Emergent su (2) dynamics and perfect quantum many-body scars, Physical review letters 122, 220603 (2019).
  * Greiner  _et al._ [2002] M. Greiner, O. Mandel, T. W. Hänsch, and I. Bloch, Collapse and revival of the matter wave field of a bose–einstein condensate, Nature 419, 51 (2002).
  * Bloch  _et al._ [2008] I. Bloch, J. Dalibard, and W. Zwerger, Many-body physics with ultracold gases, Rev. Mod. Phys. 80, 885 (2008).
  * Polkovnikov  _et al._ [2011] A. Polkovnikov, K. Sengupta, A. Silva, and M. Vengalattore, Colloquium: Nonequilibrium dynamics of closed interacting quantum systems, Rev. Mod. Phys. 83, 863 (2011).
  * Gärttner  _et al._ [2017] M. Gärttner, J. G. Bohnet, A. Safavi-Naini, M. L. Wall, J. J. Bollinger, and A. M. Rey, Measuring out-of-time-order correlations and multiple quantum spectra in a trapped-ion quantum magnet, Nature Physics 13, 781 (2017).
  * Verstraete  _et al._ [2008] F. Verstraete, V. Murg, and J. I. Cirac, Matrix product states, projected entangled pair states, and variational renormalization group methods for quantum spin systems, Advances in Physics 57, 143 (2008).
  * Schollwöck [2011] U. Schollwöck, The density-matrix renormalization group in the age of matrix product states, Annals of Physics 326, 96 (2011).
  * Bridgeman and Chubb [2017] J. C. Bridgeman and C. T. Chubb, Hand-waving and interpretive dance: an introductory course on tensor networks, Journal of Physics A: Mathematical and Theoretical 50, 223001 (2017).
  * Hauschild and Pollmann [2018] J. Hauschild and F. Pollmann, Efficient numerical simulations with tensor networks: Tensor network python (tenpy), SciPost Physics Lecture Notes (2018).
  * Hackbusch [2012] W. Hackbusch,  _Tensor spaces and numerical tensor calculus_ , Vol. 42 (Springer, 2012).
  * Fannes  _et al._ [1992] M. Fannes, B. Nachtergaele, and R. Werner, Abundance of translation invariant pure states on quantum spin chains, letters in mathematical physics 25, 249 (1992).
  * Perez-Garcia  _et al._ [2006] D. Perez-Garcia, F. Verstraete, M. M. Wolf, and J. I. Cirac, Matrix product state representations, arXiv preprint quant-ph/0608197 (2006).
  * Perez-Garcia  _et al._ [2007] D. Perez-Garcia, F. Verstraete, M. Wolf, and J. Cirac, Matrix product state representations, Quantum Information & Computation 7, 401 (2007).
  * White [1992] S. R. White, Density matrix formulation for quantum renormalization groups, Physical review letters 69, 2863 (1992).
  * Eisert  _et al._ [2010] J. Eisert, M. Cramer, and M. B. Plenio, Colloquium: Area laws for the entanglement entropy, Reviews of Modern Physics 82, 277 (2010).
  * Hastings [2007] M. B. Hastings, An area law for one-dimensional quantum systems, Journal of Statistical Mechanics: Theory and Experiment 2007, P08024 (2007).
  * White and Feiguin [2004] S. R. White and A. E. Feiguin, Real-time evolution using the density matrix renormalization group, Physical review letters 93, 076401 (2004).
  * Daley  _et al._ [2004] A. J. Daley, C. Kollath, U. Schollwöck, and G. Vidal, Time-dependent density-matrix renormalization-group using adaptive effective hilbert spaces, Journal of Statistical Mechanics: Theory and Experiment 2004, P04005 (2004).
  * Vidal [2003] G. Vidal, Efficient classical simulation of slightly entangled quantum computations, Physical review letters 91, 147902 (2003).
  * Vidal [2004] G. Vidal, Efficient simulation of one-dimensional quantum many-body systems, Physical review letters 93, 040502 (2004).
  * Hubig  _et al._ [2020] C. Hubig, A. Bohrdt, M. Knap, F. Grusdt, and J. I. Cirac, Evaluation of time-dependent correlators after a local quench in iPEPS: hole motion in the t-J model, SciPost Phys. 8, 21 (2020).
  * Czarnik and Dziarmaga [2018] P. Czarnik and J. Dziarmaga, Time evolution of an infinite projected entangled pair state: An algorithm from first principles, Phys. Rev. B 98, 045110 (2018).
  * Carleo and Troyer [2017] G. Carleo and M. Troyer, Solving the quantum many-body problem with artificial neural networks, Science 355, 602 (2017).
  * Schmitt and Heyl [2020] M. Schmitt and M. Heyl, Quantum many-body dynamics in two dimensions with artificial neural networks, Physical Review Letters 125, 100503 (2020).
  * Sharir  _et al._ [2020] O. Sharir, Y. Levine, N. Wies, G. Carleo, and A. Shashua, Deep autoregressive models for the efficient variational simulation of many-body quantum systems, Physical review letters 124, 020503 (2020).
  * Glasser  _et al._ [2018] I. Glasser, N. Pancotti, M. August, I. D. Rodriguez, and J. I. Cirac, Neural-network quantum states, string-bond states, and chiral topological states, Phys. Rev. X 8, 011006 (2018).
  * Deng  _et al._ [2017] D.-L. Deng, X. Li, and S. D. Sarma, Quantum entanglement in neural network states, Physical Review X 7, 021021 (2017).
  * Levine  _et al._ [2019] Y. Levine, O. Sharir, N. Cohen, and A. Shashua, Quantum entanglement in deep learning architectures, Physical review letters 122, 065301 (2019).
  * LeCun  _et al._ [2015] Y. LeCun, Y. Bengio, and G. Hinton, Deep learning, nature 521, 436 (2015).
  * Goodfellow  _et al._ [2016] I. Goodfellow, Y. Bengio, A. Courville, and Y. Bengio,  _Deep learning_ , Vol. 1 (MIT press Cambridge, 2016).
  * Hornik  _et al._ [1989] K. Hornik, M. Stinchcombe, and H. White, Multilayer feedforward networks are universal approximators, Neural networks 2, 359 (1989).
  * Cybenko [1989] G. Cybenko, Approximation by superpositions of a sigmoidal function, Mathematics of control, signals and systems 2, 303 (1989).
  * Leshno  _et al._ [1993] M. Leshno, V. Y. Lin, A. Pinkus, and S. Schocken, Multilayer feedforward networks with a nonpolynomial activation function can approximate any function, Neural networks 6, 861 (1993).
  * Kim and Adali [2002] T. Kim and T. Adali, Universal approximation of fully complex feed-forward neural networks, in  _2002 IEEE International Conference on Acoustics, Speech, and Signal Processing_ , Vol. 1 (IEEE, 2002) pp. I–973.
  * Voigtlaender [2020] F. Voigtlaender, The universal approximation theorem for complex-valued neural networks, arXiv preprint arXiv:2012.03351 (2020).
  * Hinton and Salakhutdinov [2006] G. E. Hinton and R. R. Salakhutdinov, Reducing the dimensionality of data with neural networks, science 313, 504 (2006).
  * Salakhutdinov  _et al._ [2007] R. Salakhutdinov, A. Mnih, and G. Hinton, Restricted boltzmann machines for collaborative filtering, in  _Proceedings of the 24th international conference on Machine learning_ (2007) pp. 791–798.
  * Cai and Liu [2018] Z. Cai and J. Liu, Approximating quantum many-body wave functions using artificial neural networks, Physical Review B 97, 035116 (2018).
  * Choo  _et al._ [2019] K. Choo, T. Neupert, and G. Carleo, Two-dimensional frustrated ${J}_{1}\text{$-$}{J}_{2}$ model studied with neural network quantum states, Phys. Rev. B 100, 125124 (2019).
  * Rumelhart  _et al._ [1986] D. E. Rumelhart, G. E. Hinton, and R. J. Williams, Learning representations by back-propagating errors, nature 323, 533 (1986).
  * Hochreiter and Schmidhuber [1997] S. Hochreiter and J. Schmidhuber, Long short-term memory, Neural computation 9, 1735 (1997).
  * Hibat-Allah  _et al._ [2020] M. Hibat-Allah, M. Ganahl, L. E. Hayward, R. G. Melko, and J. Carrasquilla, Recurrent neural network wave functions, Physical Review Research 2, 023358 (2020).
  * Luo  _et al._ [2020] D. Luo, Z. Chen, J. Carrasquilla, and B. K. Clark, Probabilistic formulation of open quantum many-body systems dynamics with autoregressive models, arXiv preprint arXiv:2009.05580 (2020).
  * Choo  _et al._ [2018] K. Choo, G. Carleo, N. Regnault, and T. Neupert, Symmetries and many-body excitations with neural-network quantum states, Physical review letters 121, 167204 (2018).
  * Choo  _et al._ [2020] K. Choo, A. Mezzacapo, and G. Carleo, Fermionic neural-network states for ab-initio electronic structure, Nature communications 11, 1 (2020).
  * Nomura [2021] Y. Nomura, Helping restricted boltzmann machine with quantum-state representation by restoring symmetry, Journal of Physics: Condensed Matter (2021).
  * Vieijra  _et al._ [2020] T. Vieijra, C. Casert, J. Nys, W. De Neve, J. Haegeman, J. Ryckebusch, and F. Verstraete, Restricted boltzmann machines for quantum states with non-abelian or anyonic symmetries, Phys. Rev. Lett. 124, 097201 (2020).
  * Vieijra and Nys [2021] T. Vieijra and J. Nys, Many-body quantum states with exact conservation of non-abelian and lattice symmetries through variational monte carlo, arXiv e-prints , arXiv (2021).
  * Carrasquilla  _et al._ [2019] J. Carrasquilla, G. Torlai, R. G. Melko, and L. Aolita, Reconstructing quantum states with generative models, Nature Machine Intelligence 1, 155 (2019).
  * Hartmann and Carleo [2019] M. J. Hartmann and G. Carleo, Neural-network approach to dissipative quantum many-body dynamics, Phys. Rev. Lett. 122, 250502 (2019).
  * Vicentini  _et al._ [2019] F. Vicentini, A. Biella, N. Regnault, and C. Ciuti, Variational neural-network ansatz for steady states in open quantum systems, Phys. Rev. Lett. 122, 250503 (2019).
  * Reh  _et al._ [2021] M. Reh, M. Schmitt, and M. Gärttner, Time-dependent variational principle for open quantum systems with artificial neural networks, arXiv preprint arXiv:2104.00013 (2021).
  * Nagy and Savona [2019] A. Nagy and V. Savona, Variational quantum monte carlo method with a neural-network ansatz for open quantum systems, Phys. Rev. Lett. 122, 250501 (2019).
  * Yoshioka and Hamazaki [2019] N. Yoshioka and R. Hamazaki, Constructing neural stationary states for open quantum many-body systems, Phys. Rev. B 99, 214306 (2019).
  * Gutiérrez and Mendl [2019] I. L. Gutiérrez and C. B. Mendl, Real time evolution with neural-network quantum states, arXiv preprint arXiv:1912.08831 (2019).
  * Fabiani and Mentink [2019] G. Fabiani and J. H. Mentink, Investigating ultrafast quantum magnetism with machine learning, SciPost Phys. 7, 4 (2019).
  * Fabiani  _et al._ [2021] G. Fabiani, M. Bouman, and J. Mentink, Supermagnonic propagation in two-dimensional antiferromagnets, arXiv preprint arXiv:2101.10945 (2021).
  * Czischek  _et al._ [2018] S. Czischek, M. Gärttner, and T. Gasenzer, Quenches near ising quantum criticality as a challenge for artificial neural networks, Phys. Rev. B 98, 024311 (2018).
  * McMillan [1965] W. L. McMillan, Ground state of liquid he 4, Physical Review 138, A442 (1965).
  * Foulkes  _et al._ [2001] W. Foulkes, L. Mitas, R. Needs, and G. Rajagopal, Quantum monte carlo simulations of solids, Reviews of Modern Physics 73, 33 (2001).
  * Yokoyama and Shiba [1987] H. Yokoyama and H. Shiba, Variational monte-carlo studies of hubbard model. i, Journal of the Physical Society of Japan 56, 1490 (1987).
  * Sorella [2005a] S. Sorella, Wave function optimization in the variational monte carlo method, Physical Review B 71, 241103 (2005a).
  * Sorella [2001] S. Sorella, Generalized lanczos algorithm for variational quantum monte carlo, Physical Review B 64, 024512 (2001).
  * Sorella [2005b] S. Sorella, Wave function optimization in the variational monte carlo method, Phys. Rev. B 71, 241103 (2005b).
  * Neuscamman  _et al._ [2012] E. Neuscamman, C. Umrigar, and G. K.-L. Chan, Optimizing large parameter sets in variational quantum monte carlo, Physical Review B 85, 045103 (2012).
  * Amari [1998] S.-I. Amari, Natural gradient works efficiently in learning, Neural computation 10, 251 (1998).
  * Luo  _et al._ [2021] D. Luo, Z. Chen, K. Hu, Z. Zhao, V. M. Hur, and B. K. Clark, Gauge invariant autoregressive neural networks for quantum lattice models, arXiv preprint arXiv:2101.07243 (2021).
  * Kochkov and Clark [2018] D. Kochkov and B. K. Clark, Variational optimization in the ai era: Computational graph states and supervised wave-function optimization, arXiv preprint arXiv:1811.12423 (2018).
  * Westerhout  _et al._ [2020] T. Westerhout, N. Astrakhantsev, K. S. Tikhonov, M. I. Katsnelson, and A. A. Bagrov, Generalization properties of neural network approximations to frustrated magnet ground states, Nature communications 11, 1 (2020).
  * Jónsson  _et al._ [2018] B. Jónsson, B. Bauer, and G. Carleo, Neural-network states for the classical simulation of quantum computing, arXiv preprint arXiv:1808.05232 (2018).
  * Carleo  _et al._ [2019] G. Carleo, K. Choo, D. Hofmann, J. E. Smith, T. Westerhout, F. Alet, E. J. Davis, S. Efthymiou, I. Glasser, S.-H. Lin, _et al._ , Netket: A machine learning toolkit for many-body quantum systems, SoftwareX 10, 100311 (2019).
  * Borin and Abanin [2020] A. Borin and D. A. Abanin, Approximating power of machine-learning ansatz for quantum many-body states, Physical Review B 101, 195141 (2020).
  * Kim and Huse [2013] H. Kim and D. A. Huse, Ballistic spreading of entanglement in a diffusive nonintegrable system, Phys. Rev. Lett. 111, 127205 (2013).
  * Verstraete and Cirac [2006] F. Verstraete and J. I. Cirac, Matrix product states represent ground states faithfully, Phys. Rev. B 73, 094423 (2006).
  * Eisert [2013] J. Eisert, Entanglement and tensor network states, arXiv preprint arXiv:1308.3318 (2013).
  * Kingma and Ba [2014] D. P. Kingma and J. Ba, Adam: A method for stochastic optimization, arXiv preprint arXiv:1412.6980 (2014).
  * He  _et al._ [2016] K. He, X. Zhang, S. Ren, and J. Sun, Deep residual learning for image recognition, in  _Proceedings of the IEEE conference on computer vision and pattern recognition_ (2016) pp. 770–778.
  * Zagoruyko and Komodakis [2016] S. Zagoruyko and N. Komodakis, Wide residual networks, arXiv preprint arXiv:1605.07146 (2016).
  * Xie  _et al._ [2017] S. Xie, R. Girshick, P. Dollár, Z. Tu, and K. He, Aggregated residual transformations for deep neural networks, in  _Proceedings of the IEEE conference on computer vision and pattern recognition_ (2017) pp. 1492–1500.
  * Oord  _et al._ [2016a] A. v. d. Oord, S. Dieleman, H. Zen, K. Simonyan, O. Vinyals, A. Graves, N. Kalchbrenner, A. Senior, and K. Kavukcuoglu, Wavenet: A generative model for raw audio, arXiv preprint arXiv:1609.03499 (2016a).
  * Nair and Hinton [2010] V. Nair and G. E. Hinton, Rectified linear units improve restricted boltzmann machines, in  _Icml_ (2010).
  * Szabó and Castelnovo [2020] A. Szabó and C. Castelnovo, Neural network wave functions and the sign problem, Physical Review Research 2, 033075 (2020).
  * Huang  _et al._ [2017] G. Huang, Z. Liu, L. Van Der Maaten, and K. Q. Weinberger, Densely connected convolutional networks, in  _Proceedings of the IEEE conference on computer vision and pattern recognition_ (2017) pp. 4700–4708.
  * Vaswani  _et al._ [2017] A. Vaswani, N. Shazeer, N. Parmar, J. Uszkoreit, L. Jones, A. N. Gomez, L. Kaiser, and I. Polosukhin, Attention is all you need, arXiv preprint arXiv:1706.03762 (2017).
  * Oord  _et al._ [2016b] A. v. d. Oord, N. Kalchbrenner, O. Vinyals, L. Espeholt, A. Graves, and K. Kavukcuoglu, Conditional image generation with pixelcnn decoders, arXiv preprint arXiv:1606.05328 (2016b).
  * Salimans  _et al._ [2017] T. Salimans, A. Karpathy, X. Chen, and D. P. Kingma, Pixelcnn++: Improving the pixelcnn with discretized logistic mixture likelihood and other modifications, arXiv preprint arXiv:1701.05517 (2017).
  * Hofmann  _et al._ [2021] D. Hofmann, G. Fabiani, J. H. Mentink, G. Carleo, and M. A. Sentef, Role of stochastic noise and generalization error in the time propagation of neural-network quantum states, arXiv preprint arXiv:2105.01054 (2021).
  * Abadi  _et al._ [2015] M. Abadi, A. Agarwal, P. Barham, E. Brevdo, Z. Chen, C. Citro, G. S. Corrado, A. Davis, J. Dean, M. Devin, S. Ghemawat, I. Goodfellow, A. Harp, G. Irving, M. Isard, Y. Jia, R. Jozefowicz, L. Kaiser, M. Kudlur, J. Levenberg, D. Mané, R. Monga, S. Moore, D. Murray, C. Olah, M. Schuster, J. Shlens, B. Steiner, I. Sutskever, K. Talwar, P. Tucker, V. Vanhoucke, V. Vasudevan, F. Viégas, O. Vinyals, P. Warden, M. Wattenberg, M. Wicke, Y. Yu, and X. Zheng, TensorFlow: Large-scale machine learning on heterogeneous systems (2015), software available from tensorflow.org.
  * Lin [2021] S.-H. Lin, Neural-network-quantum-state, https://github.com/ShHsLin/Neural-Network-Quantum-State (2021).
  * Bengio and Bengio [2000] Y. Bengio and S. Bengio, Modeling high-dimensional discrete data with multi-layer neural networks, Advances in Neural Information Processing Systems 12, 400 (2000).
  * Uria  _et al._ [2016] B. Uria, M.-A. Côté, K. Gregor, I. Murray, and H. Larochelle, Neural autoregressive distribution estimation, The Journal of Machine Learning Research 17, 7184 (2016).
  * Wu  _et al._ [2019] D. Wu, L. Wang, and P. Zhang, Solving statistical mechanics using variational autoregressive networks, Physical review letters 122, 080602 (2019).
  * Liu  _et al._ [2021] J.-G. Liu, L. Mao, P. Zhang, and L. Wang, Solving quantum statistical mechanics with variational autoregressive networks and quantum circuits, Machine Learning: Science and Technology 2, 025011 (2021).
  * Albareda  _et al._ [2014] G. Albareda, H. Appel, I. Franco, A. Abedi, and A. Rubio, Correlated electron-nuclear dynamics with conditional wave functions, Physical review letters 113, 083003 (2014).
  * Hauru  _et al._ [2021] M. Hauru, M. Van Damme, and J. Haegeman, Riemannian optimization of isometric tensor networks, SciPost Phys 10, 040 (2021).
  * Luchnikov  _et al._ [2020] I. A. Luchnikov, A. Ryzhov, S. N. Filippov, and H. Ouerdane, Qgopt: Riemannian optimization for quantum technologies, arXiv preprint arXiv:2011.01894 (2020).