* (90) S. Kim, S. Kundu, E. Lee, J. Lee, S. Minwalla and C. Patel, _Grey Galaxies’ as an endpoint of the Kerr-AdS superradiant instability_ , _JHEP_ 11 (2023) 024, [2305.08922].
  * (91) N. Engelhardt and A. C. Wall, _Decoding the Apparent Horizon: Coarse-Grained Holographic Entropy_ , _Phys. Rev. Lett._ 121 (2018) 211301, [1706.02038].
  * (92) N. Engelhardt and A. C. Wall, _Coarse Graining Holographic Black Holes_ , _JHEP_ 05 (2019) 160, [1806.01281].
  * (93) S. W. Hawking and G. F. R. Ellis, _The Large Scale Structure of Space-Time_.  Cambridge Monographs on Mathematical Physics. Cambridge University Press, 2, 2023, 10.1017/9781009253161.
  * (94) S. Ryu and T. Takayanagi, _Holographic derivation of entanglement entropy from AdS/CFT_ , _Phys. Rev. Lett._ 96 (2006) 181602, [hep-th/0603001].
  * (95) S. Ryu and T. Takayanagi, _Aspects of Holographic Entanglement Entropy_ , _JHEP_ 08 (2006) 045, [hep-th/0605073].
  * (96) V. E. Hubeny, M. Rangamani and T. Takayanagi, _A Covariant holographic entanglement entropy proposal_ , _JHEP_ 07 (2007) 062, [0705.0016].
  * (97) M. Rangamani and T. Takayanagi, _Holographic Entanglement Entropy_ , vol. 931.  Springer, 2017, 10.1007/978-3-319-52573-0.
  * (98) R. Emparan and M. Martinez, _Exact Event Horizon of a Black Hole Merger_ , _Class. Quant. Grav._ 33 (2016) 155003, [1603.00712].
  * (99) P. Calabrese and A. Lefevre, _Entanglement spectrum in one-dimensional systems_ , _Phys. Rev. A_ 78 (2008) 032329, [0806.3059].
  * (100) A. Bernamonti, F. Galli, R. C. Myers and I. A. Reyes, _in progress_ , .
  * (101) M. Lostaglio and K. Korzekwa, _Continuous thermomajorization and a complete set of laws for Markovian thermal processes_ , _Phys. Rev. A_ 106 (2022) 012426, [2111.12130].
  * (102) K. Korzekwa and M. Lostaglio, _Optimizing Thermalization_ , _Phys. Rev. Lett._ 129 (2022) 040602, [2202.12616].
  * (103) D. Marolf, _Microcanonical Path Integrals and the Holography of small Black Hole Interiors_ , _JHEP_ 09 (2018) 114, [1808.00394].
  * (104) D. Marolf and J. E. Santos, _Stability of the microcanonical ensemble in Euclidean Quantum Gravity_ , _JHEP_ 11 (2022) 046, [2202.12360].


# Numerical differentiation of noisy data:
A unifying multi-objective optimization framework

Floris van Breugel1, J. Nathan Kutz2, Bingni W. Brunton3

###### Abstract

Computing derivatives of noisy measurement data is ubiquitous in the physical,
engineering, and biological sciences, and it is often a critical step in
developing dynamic models or designing control. Unfortunately, the
mathematical formulation of numerical differentiation is typically ill-posed,
and researchers often resort to an ad hoc process for choosing one of many
computational methods and its parameters. In this work, we take a principled
approach and propose a multi-objective optimization framework for choosing
parameters that minimize a loss function to balance the faithfulness and
smoothness of the derivative estimate. Our framework has three significant
advantages. First, the task of selecting multiple parameters is reduced to
choosing a single hyper-parameter. Second, where ground-truth data is unknown,
we provide a heuristic for automatically selecting this hyper-parameter based
on the power spectrum and temporal resolution of the data. Third, the optimal
value of the hyper-parameter is consistent across different differentiation
methods, thus our approach unifies vastly different numerical differentiation
methods and facilitates unbiased comparison of their results. Finally, we
provide an extensive open-source Python library pynumdiff to facilitate easy
application to diverse datasets (https://github.com/florisvb/PyNumDiff).

1Department of Mechanical Engineering, University of Nevada, Reno, NV 89557
(e-mail<EMAIL_ADDRESS>
2Department of Applied Math, University of Washington, Seattle, WA, 98195, USA
(e-mail<EMAIL_ADDRESS>
3Department of Biology, University of Washington, Seattle, WA, 98195 (e-mail:
<EMAIL_ADDRESS>
_Keywords:_ Numerical differentiation, derivatives, optimization, data-driven
modeling.

##  1 Introduction

Derivatives describe many meaningful characteristics of physical and
biological systems, including spatial gradients and time rates-of change.
However, these critical quantities are often not directly measurable by
sensors. Although computing derivatives of analytic equations is
straightforward, estimating derivatives from real sensor data remains a
significant challenge because sensor data is invariably corrupted by noise
[1]. More accurate estimation of derivatives would improve our ability to
produce robust diagnostics, formulate accurate forecasts, build dynamic or
statistical models, implement control protocols, and inform policy making.
There exists a large and diverse set of mathematical tools for estimating
derivatives of noisy data, most of which are formulated as an ill-posed
problem regularized by some appropriate smoothing constraints. However, the
level and type of regularization are typically imposed in an ad hoc fashion,
so that there is currently no consensus “best-method” for producing “best-fit”
derivatives.

One particularly impactful application of estimating derivatives is the use of
time-series data in modeling complex dynamical systems. These models are of
the form $d{\mathbf{x}}/dt=\dot{{\mathbf{x}}}=f({\mathbf{x}})$, where
${\mathbf{x}}$ is the state of the system. Models of this kind have been
integral to much of our understanding across science and engineering [2],
including in classical mechanics [3], electromagnetism [4], quantum mechanics
[5], chemical kinetics [6], ecology [7] epidemiology [8], and neuroscience [9,
10, 11]. In some cases, even higher order time derivatives are also crucial
for understanding the dynamics [12]. A recent innovation in understanding
complex dynamical systems uses data-driven modeling, where the underlying
dynamics are learned directly from sensor data using a variety of modern
methods [13, 14, 15]. For this application in particular, a derivative with
both small and unbiased errors is crucial for learning interpretable dynamics.

In principle, the discrete derivative of position can be estimated as the
finite difference between adjacent measurements. If we write the vector of all
noiseless positions in time measured with timestep $\Delta t$ as $\mathbf{x}$,
then

$\mathbf{\dot{x}}_{k}=\frac{\mathbf{x}_{k+1}-\mathbf{x}_{k}}{\Delta t},$ (1)

where $k$ indexes snapshots in time. In reality, however, only noisy
measurements $\mathbf{y}$ are available,

$\mathbf{y}=\mathbf{x}+\boldsymbol{\eta},$

where $\boldsymbol{\eta}$ represents measurement noise. Here we will assume
$\boldsymbol{\eta}$ is zero-mean Gaussian noise with unknown variance. Even
with noise of moderate amplitude, a naïve application of Eq. (1) produces
derivative estimates that are far too noisy to be useful (Fig. 1A). Thus, more
sophisticated methods for data smoothing and/or differentiation of noisy time
series measurements of position $\mathbf{y}$ are required.

Although smoothing mitigates the errors, it can also introduce biases. Our
goal in this paper is to develop a general approach for methodically choosing
parameters that balance the need to minimize both error and bias. We use
$\mathbf{\hat{x}}$ and $\mathbf{\hat{\dot{x}}}$ to denote the smoothed
estimates of the position and its derivative computed from $\mathbf{y}$,
respectively. To evaluate the quality of these estimates, we compare these
estimates to the true discrete time position and its derivative, $\mathbf{x}$
and $\mathbf{\dot{x}}$. Developing approaches for estimating
$\mathbf{\hat{\dot{x}}}$ from noisy measurements $\mathbf{y}$ has been the
focus of intense research for many decades. Despite the diversity of methods
that have been developed, only a few studies have performed a comprehensive
comparison of their performance on different types of problems [16, 17, 1].

In this paper, we tackle the challenge of parameter selection by developing a
novel, multi-objective optimization framework for choosing parameters to
estimate the derivative of noisy data. Our approach minimizes a loss function
consisting of a weighted sum of two metrics computed from the derivative
estimate: the faithfulness of the integral of the derivative and its
smoothness. We suggest these metrics as proxies for minimizing the error and
bias of the estimated derivative, and we show that sweeping through values of
a single hyper-parameter $\gamma$ produces derivative estimates that generally
trace the Pareto front of solutions that minimize error and bias. Importantly,
this optimization framework assumes no knowledge of the underlying true
derivative and reduces the task of selecting many parameters of any
differentiation algorithm to solving a loss function with a single hyper-
parameter. Furthermore, we show that the value of the hyper-parameter is
nearly universal across four different differentiation methods, making it
possible to compare the results in a fair and unbiased way. For real-world
applications, we provide a simple heuristic for automatically determining a
value of $\gamma$ that is derived from the power spectrum and temporal
resolution of the data. All of the functionality described in this paper is
implemented in an open-source Python toolkit pynumdiff, which is found here:
https://github.com/florisvb/PyNumDiff.

Figure 1:  Choice of parameters leads to a diversity of derivative estimates.
A. Noisy time series data, from a Lorenz system, and the corresponding finite
difference derivative. B. To evaluate the quality of a derivative estimate
relative to the ground truth, we consider two metrics: Root Mean Square Error
(RMSE), and the Pearson’s correlation coefficient ($R^{2}$) between the error
and the true value of the derivative. Gray dots show the values of these
metrics for $5,481$ different sets of parameter choices for a smoothed
Savitzky-Golay filter. The violet line is the result of our multi-objective
optimization framework and nearly traces the Pareto front of the metrics. The
derivative estimates and metrics for the five colored points along the Pareto
front are shown in C and D, respectively.

##  2 Motivation for error metrics

What is a “good” estimate of a derivative? Let us start by considering a toy
system with synthetic measurement noise, where we are able to evaluate the
quality of an estimated derivative by comparing to the true, known derivative.
We consider two metrics for evaluating the quality of a derivative Fig. 1B–D;
later, we use these same metrics to evaluate the performance of our
optimization framework, which does not have access to the ground truth.

First, the most intuitive metric is how faithfully the estimated derivative
$\mathbf{\hat{\dot{x}}}$ approximates the actual derivative
$\mathbf{\dot{x}}$. We can measure this using the root-mean-squared error,

$\text{RMSE}(\mathbf{\hat{\dot{x}}},\mathbf{\dot{x}})=\lVert(\mathbf{\hat{\dot{x}}}-\mathbf{\dot{x}})\rVert_{2},$
(2)

where $\lVert\cdot\rVert_{2}$ is the vector 2-norm. If the data are very
noisy, a small RMSE can only be achieved by applying significant smoothing.
However, smoothing the data often attenuates sharp peaks in the data and
results in underestimating the magnitude of the derivative.

To measure the degree to which the derivative estimate is biased due to
underestimates of the actual derivative, we calculate the square of the
Pearson’s correlation coefficient, $R^{2}$, between the errors
$(\mathbf{\hat{\dot{x}}}-\mathbf{\dot{x}})$ and the actual derivative
$\mathbf{\dot{x}}$. We refer to this metric as the error correlation, which is
bounded between 0 and 1. Small error correlations imply that the imposed
dynamics of the differentiation method (e.g. filtering) minimally influenced
the derivative estimate; therefore, the method of estimating derivatives would
have minimal impact on any models that are constructed using these estimates.
Conversely, large error correlations imply that the estimate is significantly
influenced by the dynamics of the differentiation method and typically
correspond to very smooth estimates. In the limit where the derivative
estimate is a horizontal line, the error correlation takes on a value of
unity. Other metrics that measure the smoothness, for example the total
variation or tortuosity, may be substituted for error correlation [1];
however, these metrics are harder to interpret. For instance, if the true
derivative is very smooth, a low total variation is desired, whereas if the
true derivative is quite variable, a high total variation would correspond to
an accurate derivative. In contrast, a low error correlation is desirable for
any true derivative.

For many datasets, the RMSE and error correlation metrics define a Pareto
front, where no single parameter choice minimizes both values (Fig. 1B).
Furthermore, the minimal RMSE can be achieved with a variety of different
error correlations. The most suitable parameter set depends on the application
of the estimated derivative: is a non-smooth derivative with minimal bias
preferred (Fig. 1B-D: teal), or one that is smooth, but biased (Fig. 1B-D:
brown). We suggest that, for most purposes, the estimated derivative that
balances these metrics (Fig. 1B-D: blue and red) serves as a reasonable
starting point.

##  3 Methods for Numerical Differentiation

A large variety of methods for numerical differentiation exist, and a complete
review of them all is beyond the scope of this paper. Instead, we have
selected four differentiation methods (Table 1), which make different
assumptions and represent different approaches to computing the derivative
including both global and local methods [1], to showcase the universal
application of our optimization framework.

One common approach to manage noisy data is to apply a smoothing filter to the
data itself, followed by a finite difference calculation. In this family of
differentiation methods, we chose to highlight the Butterworth filter [18],
which is a global spectral method with two parameters: filter order and
frequency cutoff.

The second family of methods relies instead on building a local model of the
data through linear regression. A common and effective approach involves
making a sliding polynomial fit of the data [19], often referred to as locally
estimated scatterplot smoothing (LOESS) [20]. An efficient approach for
accomplishing the same calculations is the Savitzky-Golay filter, which builds
the polynomial model in the frequency domain [21, 22]. The Savitzky-Golay
filter has two parameters: window size and polynomial order. By default, a
Savitzky-Golay filter provides a jagged derivative because the polynomial
models can change from one window to the next, so here we also apply some
smoothing by convolving the result with a Gaussian kernel. This smoothing adds
a third parameter: a smoothing window size.

The third family we consider is the Kalman filter [23, 24, 25]. The Kalman
filter is most effective when models of the system and of the noise
characteristics are known. Our focus here is the case where neither is known,
so we chose to highlight a constant acceleration forward-backward Kalman
smoother [26] with two parameters: the model and noise covariances.

Finally, we consider an optimization approach to computing derivatives with
the total variation regularization (TVR) method [27, 28]. One advantage of the
TVR methods is that there is only a single parameter, which corresponds to the
smoothness of the derivative estimate. TVR derivatives are not as widely used
as the other three methods we highlight, so we provide a brief overview here.
Solving for the TVR derivative involves first finding $\mathbf{\hat{x}}$ and
its corresponding finite-difference derivative $\mathbf{\hat{\dot{x}}}$
(calculated according to Eq. 1) that minimize the following loss function,

$L=\lVert\mathbf{y}-\mathbf{\hat{x}}\rVert_{2}+\gamma*TV(\mathbf{\hat{\dot{x}}}).$
(3)

Here $TV$ is the total variation,

$TV(\mathbf{\hat{\dot{x}}})=\frac{1}{m}\left\lVert\mathbf{\hat{\dot{x}}}_{0:m-1}-\mathbf{\hat{\dot{x}}}_{1:m}\right\rVert_{1},$
(4)

where $\lVert\cdot\rVert_{1}$ denotes the $\ell_{1}$ norm and $m$ is the
number of time snapshots in the data. The single parameter for this method is
$\gamma$, and larger values result in smoother derivatives. If $\gamma$ is
zero, this formulation reduces to a finite difference derivative.

Solutions for TVR $\mathbf{\hat{\dot{x}}}$ can be found with an iterative
solver [28]. Because both components of the loss function Eq. (3) are convex,
we can also solve for $\mathbf{\hat{\dot{x}}}$ using convex optimization
tools, such as cvxpy [29], and with a convex solver, such as MOSEK [30]. The
two methods are equivalent, if the iterative solver is repeated sufficiently
many times.

The convex solution to penalizing the first order difference in time, as in
Eq. (4), results in a piece-wise constant derivative estimate. By offloading
the calculations to a convex optimization solver, however, we can easily
penalize higher order derivatives by replacing the 1st order finite difference
derivative $\mathbf{\hat{\dot{x}}}$ in Eq. (3) with a 2nd order
($\mathbf{\hat{\ddot{x}}}$) or 3rd order ($\mathbf{\hat{\dddot{x}}}$) finite
difference derivative. Penalizing higher-order time derivatives results in
smoother derivative estimates. For example, penalizing the 2nd order
derivative results in a piece-wise linear derivative estimate, whereas
penalizing the 3rd order derivative, also known as the jerk, results in a
smooth estimate. In this paper, we will use the total variation regularized on
the jerk (TVRJ). For large datasets, solving for the TVRJ derivative is both
computationally expensive and can accumulate small errors. To manage the size
of the optimization problem, we solve for the TVRJ derivative in sliding
windows.

Full name |  | Abbreviated
---
name
# Parameters |  | Computational
---
cost
References
Butterworth filter followed by finite difference | Butterworth | 2 | low | [18]
Smooth Savitzky-Golay filter | Savitzky-Golay | 3 | low | [21, 22]
Constant acceleration forward-backward Kalman smoother | Kalman smooth | 2 | high | [26]
Total Variation Regularized Jerk | TVRJ | 1 | high | [28]
Table 1: Summary of the four differentiation methods highlighted in this
paper.

##  4 Computing derivatives of noisy data with no ground truth

With noisy data collected in the real world, no ground truth is accessible.
The RMSE and error correlation metrics described in the previous section
cannot be calculated and used to optimize parameter choices, so the parameter
selection is an ill-posed problem. Even so—somehow—parameters must be chosen.
In this section, we propose a general approach for choosing parameters and
show that for a wide range of problems, noise levels, time resolutions, and
methods, our approach yields reasonable derivative estimates without the need
for hyper-parameter turning.

###  4.1 Optimization framework without ground truth derivatives

Given noisy position measurements $\mathbf{y}$, we seek to estimate the
derivative in time of the dynamical system that underlies the measurements
$\mathbf{\hat{\dot{x}}}$. When the ground truth $\mathbf{\dot{x}}$ is unknown,
we propose choosing the set of parameters $\Phi$ (for any given numerical
algorithm, including those enumerated in Table 1) that minimize the following
loss function, which is inspired by Eq. (3),

$L=\mbox{RMSE}\bigg{(}\mbox{trapz}(\mathbf{\hat{\dot{x}}}(\Phi))+\mu,\mathbf{y}\bigg{)}+\gamma\bigg{(}{TV}\big{(}\mathbf{\hat{\dot{x}}}(\Phi)\big{)}\bigg{)},$
(5)

where $\mbox{trapz}(\cdot)$ is the discrete-time trapezoidal numerical
integral, $\mu$ resolves the unknown integration constant,

$\mu=\frac{1}{m}\sum_{k=0}^{m}\bigg{(}\mbox{trapz}(\mathbf{\hat{\dot{x}}}(\Phi))-\mathbf{y}\bigg{)},$
(6)

and $\gamma$ is a hyper-parameter. Note that this formulation has a single
hyper-parameter $\gamma$, and a heuristic for choosing $\gamma$ is introduced
in the following section.

The first term of the loss function in Eq. (5) promotes faithfulness of the
derivative estimate by ensuring that the integral of the derivative estimate
remains similar to the data, whereas the second term encourages smoothness of
the derivative estimate. If $\gamma$ is zero, the loss function simply returns
the finite difference derivative. Larger values of $\gamma$ will result in a
smoother derivative estimate.

This loss function effectively reduces the set of parameters $\Phi$ (which
ranges between 1 and 3 or more, depending on the method) to a single hyper-
parameter $\gamma$. Unfortunately, $L$ is not convex, but tractable
optimization routines can be used to solve for the set of $\Phi$ that minimize
$L$. Here we use the Nelder-Mead method [31], a gradient descent algorithm, as
implemented in SciPy [32], with multiple initial conditions.

###  4.2 Heuristics for automated hyper-parameter tuning of $\gamma$

The advantages of our loss function in Eq. (5) are that it does not require
any ground truth data, and it simplifies the process of choosing parameters by
reducing all the parameters associated with any given method for
differentiation to a single hyper-parameter $\gamma$ corresponding to the how
smooth the resulting derivative should be. To understand the qualities of the
derivative estimates resulting from parameters selected by our loss function,
we begin by analyzing the derivative estimates of noisy sinusoidal curves
using the Savitzky-Golay filter and return to our original metrics, RMSE and
error correlation to evaluate the results.

Interestingly, sweeping through values of $\gamma$ results in derivative
estimates with RMSE and error correlation values that generally follow the
Pareto front defined by all possible derivative estimates for that given
method (Fig. 2A). Which of these derivative estimates is best depends on the
intended use of the derivative; nevertheless, we suggest that a good general
purpose derivative is one that corresponds with the elbow in the lower left
corner of the stereotypical curve traced by a sweep of $\gamma$ in the RMSE
vs. error correlation space (the star-shaped markers in Fig. 2A). This point
often, but not always, corresponds to the lowest RMSE (for example, see Fig.
1). Although in many cases a quantitatively better derivative estimate than
the one found by our loss function does exist (the gray dots in Fig. 2A that
lie left of the star), the qualitative differences between these two
derivative estimates are generally small (Fig. 2A middle row).

Variable | Coeff | P-value
---|---|---
intercept | -5.26 | 0
$\log(freq)$ | -1.55 | 0
$\log(dt)$ | -0.74 | 0
$\log(noise)$ | 0.11 | 0.32
$\log(length)$ | 0.10 | 0.32

Table 2: Optimal $\log(\gamma)$ is correlated with frequency and temporal
resolution, but not the noise or length of the dataset. The table provides the
coefficients and associated p-values for a ordinary least squares model, with
an adjusted $R^{2}=0.78$.

Variable | Coeff | P-value
---|---|---
intercept | -5.1 | 0
$\log(freq)$ | -1.6 | 0
$\log(dt)$ | -0.71 | 0

Table 3: Optimal $\log(\gamma)$ can be determined based on the frequency and
temporal resolution of the data. The table provides the coefficients and
associated p-values for a ordinary least squares model, with an adjusted
$R^{2}=0.78$.

In practice, the need to choose even a single parameter can be time consuming
and arbitrary. To alleviate these issues, we derive an empirical heuristic to
guide the choice of $\gamma$ that corresponds with the elbow of the Pareto
front. We found that the best choice of $\gamma$ is dependent on the frequency
content of the data. To characterize this relationship, we evaluated the
performance of derivative estimates achieved by a Savitzky-Golay filter by
sweeping through different values of $\gamma$ for a suite of sinusoidal data
with various frequencies ($f$), noise levels (additive white (zero-mean)
Gaussian noise with variance $\sigma^{2}$), temporal resolutions ($\Delta t$),
and dataset lengths (in time steps, $L$) (Fig. 2A-B).

To describe this empirical relationship between the optimal choice of $\gamma$
and quantitative features of the data, we first considered an all-inclusive
multivariate log-linear model,

$\log(\gamma)=\alpha_{1}\log(f)+\alpha_{2}\log(\Delta
t)+\alpha_{3}\log(\sigma)+\alpha_{4}\log(L)+\alpha_{5}.$ (7)

Fitting the data (Fig. 2B triangles) to this model with ordinary least squares
resulted in an $R^{2}=0.76$, suggesting that, in many cases, it is feasible to
automatically determine a reasonable guess for $\gamma$. Table 3 provides the
coefficients ($\alpha_{k}$) and associated p-values for each of the four terms
and intercept. From this analysis we can conclude that the magnitude of
measurement noise in the data is not an important predictor of $\gamma$. We
note, however, that here we have assumed that the magnitude of noise does not
change within a time-series dataset.

Eliminating the unnecessary terms from our model results in slightly adjusted
coefficients, provided in Table 3. In short, the optimal choice of $\gamma$,
assuming that both low RMSE and low error correlation are valued, can be found
according to the following relationship:

$\log(\gamma)=-1.6\log(f)-0.71\log(dt)-5.1.$ (8)

Figure 2: Optimal choice of $\gamma$ is a function of frequency and temporal
resolution of the data. A. (Top) Four example sine waves of different
frequencies (note the time scales), temporal resolutions, and noise levels.
(Middle) Comparison of the actual derivative (black dashed) with Savitzky-
Golay estimates: lowest achievable RMSE (gray), and the result from our loss
function with the optimal choice of $\gamma$ defined in the bottom panel.
(Bottom) Trade-off between error correlation and RMSE for 5,481 potential
parameter choices (gray) and the options provided by our loss function for a
sweep through $\gamma$ (colored line). The star indicates the optimal choice
of gamma, corresponding to the shoulder of the colored curve. (B) The optimal
choice of $\gamma$ (defined in A) as a function of frequency (Hz), for
different temporal resolutions of data (0.001, 0.01, 0.1 sec). Also included
in the plot, but not indicated, are different noise levels (0-mean normally
distributed with standard deviations of $0.05\%,0.5\%,5\%,$ and $25\%$ of the
amplitude) and length of the dataset ($1,4,5,25,100,500,1000$ sec). The ”+”
markers indicate results from datasets for which the period was greater than
the length of the time series, which were omitted from the fit. The diagonal
lines indicate the empirical heuristics for choosing $\gamma$ based on a
multivariate ordinary least squares model, provided in Eqn. 8 and Table 3.

We analyze the performance of our loss function and heuristic with respect to
a broad suite of representative synthetic problems. Real world data takes on a
much greater diversity of shapes than the sinusoidal timeseries we used to
derive the heuristic for choosing $\gamma$ given in Eq. (8). Because it is
difficult to define a clear quantitative description of the range of shapes
that real data might take on (such as frequency for a sinusoidal function), we
first examine differentiating one component of a Lorenz system (Fig. 3). From
the power spectra, we select a frequency corresponding to the frequency where
the power begins to decrease and the noise of the spectra increases. Although
somewhat arbitrary, this approach (in conjunction with Eq. (8)) allows us to
use a standard signal processing tool to quickly determine a choice of
$\gamma$. Our method produces reliable derivatives without further tuning in
each case except high noise and low temporal resolution (Fig. 3, fourth row),
which is not surprising considering the low quality of the data.

Next we consider four other synthetic problems, all with similarly effective
results (Fig. 3). For the logistic growth problem, the curve traced by our
loss function takes on a more complicated shape, perhaps because the
characteristics of data vary substantially across time. Still, our heuristic
results in a good choice of parameters that correspond to an accurate
derivative. For the triangle wave, the loss function does a good job of
tracing the Pareto front, and the heuristic selects an appropriate value of
$\gamma$, yet the resulting derivative does show significant errors. This is
likely due to two reasons. First, the Savitzky-Golay filter is designed to
produce a smooth derivative, rather than a piece-wise constant one. Second,
the frequency content of the data varies between two extremes, near-zero, and
near-infinity. For the sum of sines problem, selecting the appropriate
frequency cutoff is more straightforward than the previous problems, as we can
simply choose a frequency shortly after the high frequency spike in the
spectra. The final problem is a time-series resulting from a simulated
dynamical system controlled by a proportional-integral controller subject to
periodic disturbances. This data is a challenging problem for numerical
differentiation, as the position data almost appears to be a straight line but
does contain small variations. Our loss function does an excellent job of
tracing the Pareto front in this case, and our heuristic results in an
appropriate choice of $\gamma$.

Figure 3: Heuristic for choosing $\gamma$ is effective across a broad range of
toy problems, using a Savitzky-Golay filter. The first column shows raw
(synthetic) position data, indicating the shape of the data, degree of noise,
and temporal resolution. Next we evaluate the performance of derivative
estimate using the metrics described in the Fig. 1. Gray dots indicate the
range of outcomes for 5,481 parameter choices, the violet line indicates the
options provided by our loss function, and the red star indicates the
performance using the automatically selected value of $\gamma$ according to
Eqn. 8. Frequency of the data is evaluated by inspecting the power spectra;
the red line indicates the frequency used to determine $\gamma$. The final two
columns compare the ground truth and estimates for position and velocity.
Figure 4: Loss function and heuristic for choosing $\gamma$ is equally
effective for different differentiation methods. A. Synthetic noisy data from
the same Lorenz system as shown in Fig. 3. B. Comparison of metrics, position,
and velocity estimates using four differentiation methods, with the same value
of $\gamma$, as determined through the spectral analysis in Fig. 3. C. Overlay
of the Pareto fronts and velocities for all four methods.  Figure 5: Loss
function and heuristic for choosing $\gamma$ is equally effective for
different differentiation methods across a range of toy problems. Data plotted
as in Fig. 4C, for each of the scenarios presented in Fig. 3.

###  4.3 Direct comparison of differentiation methods

We examine how our loss function and heuristic for choosing $\gamma$ might
perform on other differentiation methods beyond the Savitzky-Golay filter.
Figure 4 shows that for a noisy Lorenz system, the possible solution space is
similar for all four methods we highlighted earlier, and our loss function
achieves a similar Pareto front in each case. Note that although the Savitzky-
Golay and Butterworth filters both operate in the frequency domain, the Kalman
smoother and TVRJ methods do not.

Interestingly, for all four differentiation methods, the possible solutions
(the gray dots), and in particular their Pareto front, are quite similar, with
the exception of the TVRJ method. This deviation may be because the TVRJ
method only contains a single parameter. Our loss function, which defines the
colored curves in the RMSE vs error correlation space, results in similar
curves for each method, each of which follows the Pareto front quite closely.
Although there are some differences in the location along the Pareto front
that our heuristic selects as the optimal choice for each method, the
resulting derivative estimates are qualitatively quite similar. A close
comparison of the curves defined by the loss function, and the points selected
by the heuristic, suggest that the Kalman and TVRJ methods produce slightly
more accurate derivative estimates with a lower error correlation. However,
looking at the resulting derivatives we see that the regions where the
derivative estimates have high errors, all four estimates exhibit similar
errors, suggesting that these errors may be a result of the data, not the
method.

These results suggest that our optimization framework is universal across
different methods, a claim further supported by its performance across a range
of synthetic problems (Fig. 5. The most significant result of this analysis is
that all four methods, despite being very different in their underlying
mathematics, behave similarly under both our loss function and heuristic for
choosing $\gamma$ across a wide range of data. Even in the case where they
disagree on a quantitative level (second row, low temporal resolution Lorenz
data), and the Savitzky-Golay filter appears to provide the estimate with the
lowest error correlation, the resulting derivative estimates are in fact
qualitatively quite similar.

Taking a closer look at the errors in the derivative estimates across the
range of toy problems shown in Fig. 5 reveals a subtle point about the
limitations of the differentiation methods we highlight here. For all four
methods, the errors in the derivative estimates are largest for the triangle
problem, and to a lesser extent the proportional-integral control problem.
These errors likely stem from two particular challenges. First, the frequency
content of the data is very heterogeneous: it is near zero between the peaks
and valleys, and near infinite at the peaks and valleys. Furthermore, the
frequency of the oscillations for the triangle increase with time. Second, all
four of the methods we highlighted here are designed to provide smooth
derivatives, whereas the true derivative for the triangle problem is piece-
wise constant. If this were known from the outset, it might be more effective
to choose a method that is designed to return piece-wise constant derivatives,
such as the total variation regularized on the $1^{st}$ derivative.

##  5 Demonstrations on real-world data

The real value of our multi-objective optimization framework is its
straightforward application to real, noisy data where no ground truth data is
available. Here we provide two such examples: differentiation of the new
confirmed daily cases in the United States of COVID-19, the disease caused by
SARS-CoV-2 (Fig. 6), and differentiation of gyroscope data from a downhill ski
(Fig. 7). In both examples, we examine the power spectra of the data to choose
a cutoff frequency that corresponds to the start of the dropoff in power. This
cutoff frequency, in conjunction with the time resolution of the data, are
then used as inputs to our heuristic described by Eq. (8) to determine an
optimal value of $\gamma$. With $\gamma$ chosen, we minimize our loss function
from Eq. (5) to find the optimal parameters for numerical differentiation.

Figure 6: Numerical differentiation of new confirmed daily cases in the United
States of COVID-19 [33] with no parameter tuning. A. Raw new daily cases. B.
Power spectra of the data, indicating the cutoff frequency (red) used for
selecting $\gamma=4.1$. C. Smoothed new daily cases, and their derivative,
using a Savitzky-Golay filter (violet), a Butterworth filter (green), a
constant acceleration Kalman forward-backward smoother (light blue), and total
variation regularized jerk (orange). Note the similarity between all four
methods except in the very last week, despite the significant differences in
how each method works and the automated parameter selection.

The year 2020 has seen a dramatic growth of the prevalence of a novel
coronavirus, SARS-CoV-2, which causes the disease known as COVID-19.
Estimating and understanding the rate of increase of disease incidence is
important for guiding appropriate epidemiological, health, and economic
policies. In the raw data ( [33], https://github.com/CSSEGISandData/COVID-19)
for the raw new confirmed daily cases of COVID-19, Fig. 6A) there is a clear
oscillation with a period of one week, most likely due to interruptions in
testing and reporting during weekends. As such, we selected a lower cutoff
frequency of 2 months, corresponding to the beginning of the steep drop off in
the power spectra (Fig. 6B). If the weekly oscillations were important, one
could just as easily select a cutoff frequency of 1/week. Our heuristic for
choosing $\gamma$ was based on sinusoidal data with a limited domain of time
resolutions ranging from 0.001 to 0.1 seconds, so we scaled the time step
units of the COVID-19 data to be close to this range, using $dt=$ 1 day,
rather than 86,400 seconds. Our chosen cutoff frequency yielded a value of
$\gamma=4.1$.

Using this same value of $\gamma$ for each of the four differentiation methods
under consideration resulted in very similar smoothed daily case estimates and
derivatives, except during the final 2 weeks (Fig. 6C). This highlights an
important application of our method, which facilitates easy and fair
comparison between different smoothing methods. Where these methods disagree,
it is clear that none of the estimates can be trusted. A more subtle
difference between the methods is that the Butterworth filter appears to
preserve a larger remnant of the weekly oscillations seen in the raw data.

Figure 7: Numerical differentiation of noisy gyroscope data from a downhill
ski during one ski run, with no parameter tuning. A. Data from one axis of a
gyroscope attached to the center of a downhill ski. B. Power spectra of the
data, indicating the cutoff frequency (red) used for selecting $\gamma=11.5$.
C. Zoomed in section of the data from A, which was used to optimize parameter
selection. D. Smoothed angular velocities and angular accelerations,
calculated using a Savitzky-Golay filter and the optimal parameters determined
using our heuristic and loss function. E-F. Zoomed in sections from D.

Finally, we consider angular velocity data collected from a gyroscope attached
to a downhill ski over one minute of descent (Fig. 7A) (ICM-20948, SparkFun;
Wildcat Ski, Moment Skis). This type of data is representative of kinematic
data that might be collected during experiments with robots or animals, which
might be used to construct data-driven models of their dynamics [34]. From the
power spectrum, we chose a cutoff frequency of 0.2 Hz (Fig. 7B). This
selection together with the time resolution of 0.0009 seconds yielded an
optimal value of $\gamma=11.5$ using our heuristic. We calculated the smoothed
angular velocity and acceleration estimates using a Savitzky-Golay filter
(Fig. 7D–F). The other methods showed similar results (not shown for visual
clarity), though the total variation method is not recommended for large
datasets like this one due to the compounding computational costs.

##  6 Discussion

In summary, this paper develops a principled multi-objective optimization
framework to provide clear guidance for solving the ill-posed problem of
numerical differentiation of noisy data, with a particular focus on parameter
selection. We define two independent metrics for quantifying the quality of a
numerical derivative estimate of noisy data: the RMSE and error correlation.
Unfortunately, neither metric can be evaluated without access to ground truth
data. Instead, we show that the total variation of the derivative estimate,
and the RMSE of its integral, serve as effective proxies. We then introduced a
novel loss function that balances these two proxies, reducing the number of
parameters that must be chosen for any given numerical differentiation method
to a single universal hyperparameter, which we call $\gamma$. Importantly, the
derivative estimates resulting from a sweep of $\gamma$ lie close to the
Pareto front of all possible solutions with respect to the true metrics of
interest. Although different applications may require different values of
$\gamma$ to produce more smooth or less biased derivative estimates, we derive
an empirical heuristic for determining a general purpose starting point for
$\gamma$ given two features that can easily be determined from timeseries
data: the cutoff frequency and time step. Our method also makes it possible to
objectively compare the outputs for different methods. We found that for each
problem that we tried, the four differentiation methods we explored in depth,
including both local and global methods, all produce qualitatively similar
results.

In our loss function we chose to use the RMSE of the integral of the
derivative estimate and the total variation of the derivative estimate as our
metrics. However, our loss function can be extended to a more general form,

$L=M_{1}(\mathbf{\hat{\dot{x}}},\mathbf{\dot{x}})+\gamma_{2}M_{2}(\mathbf{\hat{\dot{x}}},\mathbf{\dot{x}})+\cdots+\gamma_{p}M_{p}(\mathbf{\hat{\dot{x}}},\mathbf{\dot{x}}),$
(9)

where $M_{1},M_{2},\cdots,M_{p}$ represent $p$ different metrics that could be
used, balanced by $p-1$ hyper-parameters. Alternative metrics include, for
example, the tortuosity of the derivative estimate, the error correlation
between the data and the integral of the derivative estimate, a metric
describing the distribution of the error between the data and the integral of
the derivative estimate. Depending on the qualities of the data and the
specific application, different sets of metrics may be suitable as terms in
the loss function.

Our loss function makes three important assumptions that future work may aim
to relax. The first is that we assume the data has consistent zero-mean
Gaussian measurement noise. How sensitive the loss function and heuristic are
to outliers and other noise distributions remains an open question. It is
possible that once we include other noise models, we will find differences in
the behavior of differentiation methods. The second major limitation is that
our loss function finds a single set of parameters for a given time series.
For data where the frequency content dramatically shifts over time, it may be
better to use time-varying parameters. Presently, this is limited by our
current implementation, which relies on a computationally expensive
optimization step. Future efforts may focus on ways to improve the efficiency
of these calculations. Finally, we have focused on single dimensional time-
series data. In principle, our proposed loss function can be used with multi-
dimensional data, such as 2- and 3-dimensional spatial data, with only minor
modifications.

By simplifying the process of parameter selection for numerical
differentiation to the selection of a single hyper-parameter, our approach
makes it feasible to directly compare the performance of different methods
within a given application. One particular application of interest is that of
data-driven model discovery. Methods such as sparse identification of
nonlinear dynamics (SINDy) [14], for example, rely directly on numerical
derivative estimates, and the characteristics of these estimates can have an
important impact on the resulting models. Using our method, it is now
tractable to systematically investigate the collection of data-driven models
learned from estimated derivatives of different smoothness and explore their
impact on the models.

## Acknowledgements

We are grateful for helpful discussions with Steve Brunton and Pierre
Karashchuk. FvB acknowledges support from a Moore/Sloan Data Science and
Washington Research Foundation Innovation in Data Science Postdoctoral
Fellowship, a Sackler Scholarship in Biophysics, and the National Institute of
General Medical Sciences of the National Institutes of Health (P20GM103650).
JNK acknowledges support from the Air Force Office of Scientific Research
(FA9550-19-1-0011). BWB acknowledges support from the Washington Research
Foundation, the Air Force Office of Scientific Research (FA9550-19-1-0386),
and the National Institute of Health (1R01MH117777).

## References

  * [1] K. Ahnert and M. Abel, “Numerical differentiation of experimental data: local versus global methods,” Computer Physics Communications, vol. 177, pp. 764–774, Nov. 2007.
  * [2] C.-C. Lin and L. A. Segel, Mathematics applied to deterministic problems in the natural sciences.  SIAM, 1988.
  * [3] H. Goldstein, C. Poole, and J. Safko, Classical mechanics.  Addison-Wesley, 2002.
  * [4] J. D. Jackson, Classical electrodynamics.  John Wiley & Sons, 2007.
  * [5] D. J. Griffiths and D. F. Schroeter, Introduction to quantum mechanics.  Cambridge University Press, 2018.
  * [6] R. I. Masel, Chemical Kinetics and Catalysis.  John Wiley & Sons, 2001.
  * [7] M. Kot, Elements of mathematical ecology.  Cambridge University Press, 2001.
  * [8] K. J. Rothman, S. Greenland, and T. L. Lash, Modern epidemiology.  Lippincott Williams & Wilkins, 2008.
  * [9] E. Roth, S. Sponberg, and N. Cowan, “A comparative approach to closed-loop computation,” Current Opinion in Neurobiology, vol. 25, pp. 54–62, Apr. 2014.
  * [10] M. S. Madhav and N. J. Cowan, “The synergy between neuroscience and control theory: The nervous system as inspiration for hard control challenges,” Annual Review of Control, Robotics, and Autonomous Systems, vol. 3, pp. 243–267, May 2020.
  * [11] P. K. Keel, “Generalized local linear approximation of derivatives from time series,” 2010.
  * [12] D. C. Lin, C. P. McGowan, K. P. Blum, and L. H. Ting, “Yank: the time derivative of force is an important biomechanical variable in sensorimotor systems,” The Journal of Experimental Biology, vol. 222, p. jeb180414, Sept. 2019.
  * [13] M. Schmidt and H. Lipson, “Distilling free-form natural laws from experimental data,” Science, vol. 324, pp. 81–85, Apr. 2009.
  * [14] S. L. Brunton, J. L. Proctor, and J. N. Kutz, “Discovering governing equations from data by sparse identification of nonlinear dynamical systems,” Proceedings of the National Academy of Sciences, vol. 113, no. 15, pp. 3932–3937, 2016.
  * [15] B. C. Daniels and I. Nemenman, “Automated adaptive inference of phenomenological dynamical models,” Nature Communications, vol. 6, Aug. 2015.
  * [16] J. A. Walker, “Estimating velocities and accelerations of animal locomotion: a simulation experiment comparing numerical differentiation algorithms,” The Journal of Experimental Biology, vol. 201, pp. 981–995, 1998.
  * [17] M. M. H.C. Crenshaw, C.N. Ciampaglio, “Analysis of the three-dimensional trajectories of organisms: estimates of velocity, curvature and torsion from positional information,” The Journal of Experimental Biology, vol. 203, pp. 961–982, 1998.
  * [18] S. Butterworth, “On the theory of filter amplifiers,” Experimental Wireless and the Wireless Engineer, vol. 7, pp. 536–541, 1930.
  * [19] T. Belytschko, Y. Krongauz, D. Organ, M. Fleming, and P. Krysl, “Meshless methods: An overview and recent developments,” Computer Methods in Applied Mechanics and Engineering, vol. 139, pp. 3–47, Dec. 1996.
  * [20] F. Harrell, Regression modeling strategies : with applications to linear models, logistic and ordinal regression, and survival analysis.  Cham: Springer, 2015.
  * [21] R. Schafer, “What is a savitzky-golay filter? [lecture notes],” IEEE Signal Processing Magazine, vol. 28, pp. 111–117, July 2011.
  * [22] A. Savitzky and M. J. E. Golay, “Soothing and differentiation of data by simplified least squares procedures,” Anal. Chem., vol. 36, pp. 1627–1639, 1964.
  * [23] R. E. Kalman, “A new approach to linear filtering and prediction problems,” Journal of Basic Engineering, vol. 82, pp. 35–45, Mar. 1960.
  * [24] P. Zarchan, Fundamentals of Kalman filtering : a practical approach.  Reston, VA: American Institute of Aeronautics and Astronautics, Inc, 2015\.
  * [25] A. Aravkin, J. V. Burke, L. Ljung, A. Lozano, and G. Pillonetto, “Generalized kalman smoothing: Modeling and algorithms,” Automatica, vol. 86, pp. 63–86, Dec. 2017.
  * [26] J. L. Crassidis and J. L. Junkins, Optimal Estimation of Dynamic Systems, Second Edition (Chapman & Hall/CRC Applied Mathematics & Nonlinear Science).  Chapman & Hall/CRC, 2nd ed., 2011.
  * [27] E. F. Leonid I. Rudin, Stanley Osher, “Nonlinear total variation based noise removal algorithms,” Physica D: Nonlinear Phenomena, vol. 60, pp. 259–268, 1992.
  * [28] R. Chartrand, “Numerical differentiation of noisy, nonsmooth data,” ISRN Applied Mathematics, p. 164564, 2011.
  * [29] S. Diamond and S. Boyd, “CVXPY: A Python-embedded modeling language for convex optimization,” Journal of Machine Learning Research, vol. 17, no. 83, pp. 1–5, 2016.
  * [30] M. ApS, The MOSEK optimization API for Python. Version 8.1., 2018.
  * [31] J. A. Nelder and R. Mead, “A simplex method for function minimization,” The Computer Journal, vol. 7, pp. 308–313, Jan. 1965.
  * [32] P. Virtanen, R. Gommers, T. E. Oliphant, M. Haberland, T. Reddy, D. Cournapeau, E. Burovski, P. Peterson, W. Weckesser, J. Bright, S. J. van der Walt, M. Brett, J. Wilson, K. Jarrod Millman, N. Mayorov, A. R. J. Nelson, E. Jones, R. Kern, E. Larson, C. Carey, İ. Polat, Y. Feng, E. W. Moore, J. Vand erPlas, D. Laxalde, J. Perktold, R. Cimrman, I. Henriksen, E. A. Quintero, C. R. Harris, A. M. Archibald, A. H. Ribeiro, F. Pedregosa, P. van Mulbregt, and S. . . Contributors, “SciPy 1.0: Fundamental Algorithms for Scientific Computing in Python,” Nature Methods, 2020.
  * [33] E. Dong, H. Du, and L. Gardner, “An interactive web-based dashboard to track COVID-19 in real time,” The Lancet Infectious Diseases, vol. 20, pp. 533–534, May 2020.
  * [34] P. Karashchuk, K. L. Rupp, E. S. Dickinson, E. Sanders, E. Azim, B. W. Brunton, and J. C. Tuthill, “Anipose: a toolkit for robust markerless 3d pose estimation,” bioRxiv, 2020.


# Entanglement transition in deep neural quantum states

Giacomo Passetti  Institut für Theorie der Statistischen Physik, RWTH Aachen
University and JARA-Fundamentals of Future Information Technology, 52056
Aachen, Germany  Dante M. Kennes<EMAIL_ADDRESS>Institut für
Theorie der Statistischen Physik, RWTH Aachen University and JARA-Fundamentals
of Future Information Technology, 52056 Aachen, Germany  Max Planck Institute
for the Structure and Dynamics of Matter, Center for Free-Electron Laser
Science, Luruper Chaussee 149, 22761 Hamburg, Germany

###### Abstract

Despite the huge theoretical potential of neural quantum states, their use in
describing generic, highly-correlated quantum many-body systems still often
poses practical difficulties. Customized network architectures are under
active investigation to address these issues. For a guided search of suited
network architectures a deepened understanding of the link between neural
network properties and attributes of the physical system one is trying to
describe, is imperative. Drawing inspiration from the field of machine
learning, in this work we show how information propagation in deep neural
networks impacts the physical entanglement properties of deep neural quantum
states. In fact, we link a previously identified information propagation phase
transition of a neural network to a similar transition of entanglement in
neural quantum states. With this bridge we can identify optimal neural quantum
state hyperparameter regimes for representing area as well as volume law
entangled states. The former are easily accessed by alternative methods, such
as tensor network representations, at least in low physical dimensions, while
the latter are challenging to describe generally due to their extensive
quantum entanglement. This advance of our understanding of network
configurations for accurate quantum state representation helps to develop
effective representations to deal with volume-law quantum states, and we apply
these findings to describe the ground state (area law state) vs. the excited
state (volume law state) properties of the prototypical next-nearest neighbor
spin-1/2 Heisenberg model.

## Introduction

In recent years neural quantum states (NQS) have gained tremendous attention
as a new promising technique for the exact representation of quantum many-body
wave functions 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 18,
19. Initial studies focused on shallow architectures, such as Restricted
Boltzmann Machines (RBM)1, 2. Remarkably, these studies demonstrated how RBMs
can efficiently represent quantum states with volume law entanglement in
certain special cases2, 20, 3, 21, 13, 22. This was a significant advancement
in the research of numerical ansätze for representing complex quantum many-
body states and their unique properties elevated NQSs as a potential
complementary numerical technique to tensor network-related methods or
others13, which are often bound by low-entanglement constraints.

However, there exists also known limitations of shallow network architectures
such as RBMs, which have been rigorously characterized 20, 21. In this context
is was shown that RBMs cannot efficiently represent quantum states generated
by generic polynomial-size quantum circuits. Therefore, it is nowadays well-
established that shallow networks are insufficient to address generic problems
of quantum many-body physics 20, 21. To remedy these limitations of shallow
architectures using deep network architectures was established. Here, it has
been demonstrated that increasing a network’s depth, introducing more hidden
layers, enables a higher efficiency of representation, for generic states in
many cases20, 23, 21. In fact, deep neural networks are known to exhibit an
exponential increase in expressivity as the depth of the network increases24;
a property that is crucial for their ability to model complex, high-
dimensional functions25, 26, 27. To train these deep networks efficiently a
sophisticated network architecture, such as Convolutional Neural Networks
(CNN) or Recurrent Neural Networks (RNN) can be beneficial21.

However, there is a tradeoff. (i) Deep neural networks are computationally
expensive, require large amounts of data for training and are delicate to
handle due to the large number of layers and (hyper-)parameters27, 28. (ii)
Despite the high potential of the method, it has been recently pointed out
that some complex states are hard to learn for NQS 29, 30. (iii) Unlike
traditional algorithms with rather transparent decision-making processes, the
inner workings of deep neural networks are often considered ’black boxes’, and
the characterization of how deep neural networks make decisions remains an
open problem31.

To remedy (iii) a deepened theoretical understanding of neural networks and
their workings needs to be established and has gained intense attention under
the general umbrella of ”understandable AI”. To address this problem a
physics-inspired statistical field theory approach has been successfully
applied32 within a mean field approximation. Within such a mean field analysis
it was understood that tuning hyperparameters can trigger an information phase
transition. This transitions separates an ordered from a disordered phase, in
the sense that two inputs are either correlated or not by the network, in the
ordered and disordered phase, respectively. In a physics-inspired language
this corresponds to the emergence of a mean field order parameter. Deep in
either the ordered or disordered phase the number of layers required to
represent this correlation is exponentially small, while the number of
involved layers diverges at the transition. This is akin to the behavior of
the fluctuations at second order phase transitions in physical systems. This
motivates to use networks tuned to this transition point in order to harness
the full power of deep neural networks 32, 24, 33. Therefore, understanding
these phase transitions and other properties of deep neural networks is key to
comprehending their behavior, to demystify their inner processes and last but
not least to improve their performance32, 33, 34.

As we aim to design state-of-the-art deep learning architectures for the
description of quantum many-body physics, it becomes essential to adapt the
existing knowledge from the broader machine learning community. More
specifically, there is a pressing need to fine-tune neural networks to achieve
efficient and advantageous representations of quantum many-body wave
functions. In computational physics, representability of wavefunctions is
often linked to their entanglement properties. For local, gapped Hamiltonians
it has been shown that ground states exhibit low-entanglement restricted by
the so called area law 35; limiting the entanglement between two subparts of a
system to the surface they share. Excited states generically feature a much
more unfavorable entanglement scaling with system size described a volume
law36.

In this paper, we show that ordered to chaotic deep network phase transitions
in the sense of Refs. 32, 33, 34 transduces a phase transition in the physical
properties of corresponding deep NQSs. In the chaotic vs. ordered phase of the
neural network, the behavior of the entanglement entropy scaling of the
corresponding NQS is markedly different in the sense that the ordered and the
chaotic phase of the network can best represent area-law and volume-law
states, respectively. Our work, thus, provides a reliable prescription on how
to initialize a deep NQS with the desired scaling properties. We put our
deepened understanding of the network’s properties to use by demonstrating
that the tuning of the network properties to either the ordered or chaotic
phase allows us to better approximate the energy of physical ground states or
excited states, respectively. We choose the physically relevant $J_{1}-J_{2}$
Heisenberg model and its ground state (area law) and mid-spectrum excited
state (volume law) as our work horse to illustrate this point.

## Main

### Ordered to Chaotic Phase Transition in Deep Neural Networks

We briefly review some aspects of the mean field formalism developed to study
phase transitions in infinitely deep Feed Forward Neural Network (FFNN) as far
as they are relevant to understanding our new findings. For further details,
we refer to Refs. 32, 33. We consider a deep FFNN defined by a composition of
$\mu\in\mathbb{N}$ layers where the $l$th network layer, $1\leq l\leq\mu$,
applies the transformation

$\displaystyle\begin{split}z_{i}^{(l)}=\sum_{j}W_{ij}^{(l)}y_{j}^{(l)}+b_{i}^{(l)},\qquad
y_{i}^{(l+1)}=\phi(z_{i}^{(l)}),\qquad l=1,\dots,\mu\end{split}$ (1)

from input $y^{(l)}_{i}$ to output signal $y^{(l+1)}_{i}.$ This is an affine
transformation with weight matrix (or kernel) $W^{(l)}$, bias $b^{(l)},$ and a
nonlinear activation function $\phi.$ The networks weights and biases are
drawn from zero mean Gaussian distribution such that $W^{l}_{ij}\sim
N(0,\sigma_{w}^{2}/N_{l-1})$ and $b^{(l)}_{i}\sim N(0,\sigma_{b}^{2})$.

Considering an initial input vector $z^{(0)}_{i;\alpha}$ it has been analyzed
how its magnitude is propagated along the network 32. This, in fact, is
mathematically equivalent to monitoring the second moment of the distribution,
which at each layer $l$ reads
$\mathbb{E}\left[z^{(l)}_{i;\alpha}z^{(l)}_{j;\alpha}\right]=q^{l}_{\alpha}\delta_{ij},$where

$q^{l}_{\alpha}=\sigma_{w}^{2}\int\mathcal{D}z\phi^{2}\left(\sqrt{q_{\alpha}^{l-1}}z\right)+\sigma_{b}^{2}$
(2)

with $\int\mathcal{D}z=\frac{1}{\sqrt{2\pi}}\int dze^{-\frac{1}{2}z^{2}}$. For
any arbitrary choice of $\sigma_{w}^{2}$ and $\sigma_{b}^{2}$ and bounded
activation functions $\phi$ Eq. (2) has a fixed point at
$q^{*}=\lim_{l\rightarrow\infty}q_{\alpha}^{l}$ 32. To quantify how two
different input vectors get correlated across each layer of the network, one
can consider two different initial input vectors $z^{(0)}_{i;\alpha}$ and
$z^{(0)}_{i;\beta}$, and look at the covariance of the layer defined pre-
activation vectors
$\mathbb{E}\left[z^{(l)}_{i;\alpha}z^{(l)}_{j;\beta}\right]=q_{\alpha\beta}^{l}\delta_{ij}$.
The covariance function $q_{\alpha\beta}^{l}$ is then related to the
correlation function
$c^{l}_{\alpha\beta}=q^{l}_{\alpha\beta}/\sqrt{q^{l}_{\alpha}q^{l}_{\beta}}$
and can be evaluated through the recurrence relation

$q^{l}_{\alpha\beta}=\sigma_{w}^{2}\int\mathcal{D}z_{1}\mathcal{D}z_{2}\phi\left(u_{1}\right)\phi\left(u_{2}\right)+\sigma_{b}^{2}$
(3)

with $u_{1}=\sqrt{q_{\alpha}^{l-1}}z_{\alpha}$ and
$u_{2}=\sqrt{q_{\beta}^{l-1}}\left(c_{\alpha\beta}^{l-1}z_{1}+\sqrt{1-(c_{\alpha\beta}^{l-1})^{2}}z_{2}\right)$.
It has been shown how the fix point
$c^{*}=\lim_{l\rightarrow\infty}c^{l}_{\alpha\beta}$ allows to distinguish
between an ordered phase characterized by $c^{*}=1$ and a chaotic phase where
32 $c^{*}=0$. As pointed out by Schoenholz et al.33 correlations decay
exponentially, see Fig. 1 a), through a network to their asymptotic value
according to

$c=|c^{l}_{\alpha\beta}-c^{*}|=e^{-l/\xi_{c}}.$ (4)

The behaviour of $\xi_{c}$ is reported in figure 1 b). We find, in
correspondence to the information phase transition between ordered and chaotic
phase, a divergence of the correlation decay length33.

(a)

(b)

Figure 1: (a) Qualitative behavior for the propagation of the correlation
function $c$ (defined by Eq. (4)) through a deep FFNN, the pink area
illustrates the typical exponential decay characterized by the decay length
$\xi_{c}$. (b) Order parameter $c^{*}$ (red line) and Dependence of the decay
length $\xi_{c}$ (black line) on $\sigma_{w}$, in the asymptotic limit of
infinitely many layers. At the transition point from the ordered phase to the
chaotic phase tracked by $c^{*}$, the decay length $\xi_{c}$ diverges. The
figure shown is obtained by setting $\sigma_{b}=0.01$.

### Deep feed forward Neural Quantum States

We apply the general FFNN architecture presented in the previous paragraph to
define a deep NQS describing spin-1/2 quantum states. The weight matrices,
introduced in Eq.(1), are now complex-valued, with elements still drawn from a
Gaussian distribution. For simplicity, we restrict to the case of
$\sigma_{b}=0$. The input vectors represent configurations of a spin chain,
with a total number of spins $L$. The last layer of the network is transformed
into probability amplitudes using an exponential sum. Details of this
implementation are presented in the Methods section.

As depicted in Fig. 2(a) we then consider a bipartition of the spin chain,
where the subpartition $A$ of the chain contains the first $\frac{L}{2}$ spins
of the chain and the subpartition $B$ is simply the complementary
subpartition. Following the standard prescription, from any given quantum
state represented by $\ket{\psi_{\theta}}$ the associated density matrix
$\rho_{\theta}=\ket{\psi_{\theta}}\bra{\psi_{\theta}}$ allows to write the
bipartite entanglement entropy as

(a)

(b)

Figure 2: (a) Schematics of the setup considered, a deep FFNN represents the
wave function for a spin chain. We consider the bipartition of the A and B
halves of the chain, as sketched in the figure. (b) Bipartite entanglement
entropy $S_{\mathrm{AB}}$ and correlation decay length $\xi_{c}$ varying
$\sigma_{w}$, for different values of network layers $\mu$ (coloured lines,
scale shown by left vertical axis). We see that upon increasing the network
depth, going toward the limit of an infinitely deep network, the profile of
the entanglement entropy converges toward that of $\xi_{c}$ (black line, scale
shown by the right vertical axis). Each data point has been averaged over
$10^{3}$ independent realizations; the error bar corresponds to the standard
deviation.

$S_{\rm AB}=-\text{Tr}\left[\rho_{\rm A}\log\rho_{A}\right],\quad\rho_{\rm
A}=\text{Tr}_{\rm B}\left[\rho\right]$ where $\text{Tr}_{\rm B}$ is the trace
operation with respect to the degrees of freedom of the subpartition B. We
begin our analysis by showing that for a given fixed physical system size $L$
we observe that the profile for the average bipartite entanglement entropy
$S_{\rm AB}$ for the NQS associated to random neural networks, has a
dependence on $\sigma_{W}$ that exhibit a peak at the information phase
transition point of the network, Fig. 2(b). The height of the entanglement
entropy peak increases together with the network depth $\mu$ considered as
expected for a physical phase transition in the entanglement properties. The
resulting profile obtained simulating finite network realizations, approaches
the theoretical prediction of $\xi_{c}$ for increasing numbers of layers.

### Connection between entanglement and information propagation transitions

(a)

(b)

Figure 3: (a) Scaling of the bipartite entanglement entropy $S_{\mathrm{AB}}$
for deep FFNN of $\mu=20$ layers. Labelled different colors correspond to
different value for $\sigma_{w}$, two behaviors are observed: For values of
$\sigma_{W}$ in the ordered phase the bipartite entanglement entropy converges
to a finite value, as for a 1D area law state. At $\sigma_{w}\approx 1.5$ the
bipartite entanglement entropy almost saturates the Page value, and in the
chaotic phase the entanglement displays still a volume law scaling. (b) Top
panel: Bipartite entanglement entropy $S_{\mathrm{AB}}$ for random networks
with $\mu=20$ layers, as function of $\sigma_{w}$. Bottom panel: Considering
the same random network realizations we evaluate the energy expectation with
respect to the hamiltonians $H_{J_{1}-J_{2}}$ (defined by Eq. 6) plotted with
respect to the left axis, and $\left(H_{J_{1}-J_{2}}\right)^{2}$ (right axis).
We set $J_{1}=1$ and $J_{2}=0.2$. Each data point has been averaged over
$10^{3}$ independent realizations, the error bar corresponds to the standard
deviation.

Next, we study how the bipartite entanglement entropy depends on the number of
spins represented by $L$. This scaling analysis allows to determine whether
the quantum states that are represented by random network states satisfy a
volume or an area law. While in volume law states the $S_{\rm AB}$ grows
linearly with the number of spins $L$, for an area law state $S_{\rm AB}$
should in general scale much weaker with systems size (and approach a constant
for example in the one-dimensional case considered later). The result of this
analysis is presented in Fig.3. In the same figure we compare the results with
the Page entropy that corresponds to the entanglement entropy for a completely
random quantum state and in the special case of bipartite entanglement is
given by

$S_{\rm Page}=L\log(2)-\frac{1}{2}.$ (5)

As we can see, for values of the variance smaller than $\sigma_{w}\approx
1.4$, thus in the FFNN ordered phase, $S_{\rm AB}$ converges to a constant
value, corresponding to an area law state. Across the phase transition we see
that in the chaotic case $S_{\rm AB}$ follows the same scaling behaviour as
that of the Page entropy, thus following a volume law scaling. In the chaotic
phase, as the system approaches the phase transition point, $S_{\rm AB}$
increases and almost saturates to the Page value for $\sigma_{w}\approx 1.5$.
Deep FFNNs that are near the critical point, and thus can build up
correlations between two distinct inputs, are then naturally mapped to deep
random NQS which display volume law entanglement entropy. To further emphasize
how the deep network properties are directly translated into physical
properties once the network is cast into the language of a NQS we next
consider the physical Hamiltonian operator

$H_{J_{1}-J_{2}}=J_{1}\sum_{\braket{ij}}\mathbf{S}_{i}\cdot\mathbf{S}_{j}+J_{2}\sum_{\braket{\braket{ij}}}\cdot\mathbf{S}_{i}\cdot\mathbf{S}_{j}$
(6)

with $\mathbf{S}_{i}$ the vector of spin-1/2 Pauli matrices. This corresponds
to the $J_{1}-J_{2}$ model, which is equivalent to the Heisenberg model with
the addition of spin-spin interactions between next-to-nearest neighbour
spins. The ground state of this local Hamiltonian is known to be described by
an area law for the entanglement entropy, up to logarithmic corrections when
the system is tuned to be critical. Highly excited states lying in the middle
of the energy spectrum are instead typically described by a volume law for the
entanglement entropy36. Because of this, the squared Hamiltonian operator
$\left(H_{J_{1}-J_{2}}\right)^{2}$ is characterized by a volume law ground
state. Now we compare the $\sigma_{w}$ dependence of the energy expectation
value for random deep NQS. This is shown in Fig. 3(b), demonstrating a direct
correspondence between the entanglement transition of the random NQS and the
energy associated to the Hamiltonian as well as its square. Importantly, the
energy profile has an opposite behaviour in the two cases: (i) In the ordered
phase, the network is closer in energy to the area law ground state of
$H_{J_{1}-J_{2}}$. (ii) In the chaotic phase, the network is closer in energy
to the volume law ground state of $\left(H_{J_{1}-J_{2}}\right)^{2}$. We
demonstrated that the chaotic to order phase transition of the deep FFNN is
mapped into a transition of the entanglement scaling properties for the
associated NQS architecture, and that this allows to tailor the ansatz wave
functions to be closer in energy to the physical quantum states that share the
same entanglement properties.

## Discussion

We have investigated the entanglement characteristics of NQS constructed as
deep FFNN with weight parameters drawn from a normal distribution. Our studies
have revealed that the information propagation phase transition, originally
explored in the context of random deep networks32, 33, is directly related to
the bipartite entanglement entropy as a function of the standard deviation of
the network weights. Our results provide a direct demonstration that the
correlations established between independent input vectors and propagated
through networks in proximity to the critical point are effectively translated
into the generation of highly entangled quantum states.

This discovery leads to a prescription for the preparation of initial wave
function ansätze that exhibit closer energy proximity to specific target
ground state wave functions. In fact, we demonstrate that transitioning the
network from an ordered phase to a chaotic phase facilitates the description
of states with energy characteristics resembling area law ground states (in
the ordered phase) or volume law states (in the chaotic phase).

Ground state search training protocols of challenging models are sensitive to
the fine-tuning of many hyper-parameters 30. An understandable connection,
such as the one we describe here, is a central step forward in the NQS
architecture development, where currently limits of learnability and
representability are under heavy investigation 23, 30. Our results pave the
way for future investigations on potential boosts in the performance of
optimization protocols that can efficiently target ground and excited states.

##  I Methods

We summarize here the prescriptions defining the NQS architecture studied in
the main part. We start from a conventional FFNN defined by

$\displaystyle\begin{split}z_{i}^{(l)}=\sum_{j}W_{ij}^{(l)}y_{j}^{(l)}+b_{i}^{(l)},\qquad
y_{i}^{(l+1)}=\phi(z_{i}^{(l)}),\qquad l=1,\dots,\mu.\end{split}$ (7)

For this, we consider the weight matrices appearing in Eq. (1) to be complex
valued, such that
$W_{ij}^{(l)}=w_{\mathrm{R};ij}^{(l)}+iw_{\mathrm{I};ij}^{(l)}$ is defined
with $w_{\mathrm{R};ij}^{(l)}$ and $w_{\mathrm{I};ij}^{(l)}$ independently
drawn from Gaussian distribution, setting $\sigma_{\rm R}^{2}=\sigma_{\rm
I}^{2}=\sigma_{w}^{2}/N^{l-1}$. We consider configurations with no biases such
that $b^{(l)}_{i}=0$ and we choose the scaled exponential linear unit (SELU)37
as the activation function $\phi$. The contraction of a network with a
specific vector of the Hilbert space basis is then converted to the associated
probability amplitude by applying an exponential sum, that, together with the
definition in Eq. (7), converts the last layer to the amplitude and is defined
by $\log\braket{x}{\psi_{\theta}}=\log\sum_{i=1}^{\alpha
L}\exp\left[y^{(\mu)}_{i}\right],\quad
z^{(1)}=\sum_{j}W^{(1)}_{ij}x_{i}+b^{(1)}.$ The input vectors represent
configurations of a spin chain, with a total number of spins $L$. Therefore,
each possible spin configuration consists of a vector $x\in\\{0,1\\}^{L}$ ,
where $0$ and $1$ have been arbitrarily associated with spin down and up.
Since we restrict to system sizes $L$ that allow to numerically represent the
base of the full Hilbert space, both the evaluations of the bipartite
entanglement entropy reported in Figs.2(b)-3(a) and the energy expectation
value reported in Fig.3(b) are obtained exactly.

## Data availability

All data supporting the findings of this study are available from the
corresponding authors upon reasonable request.

## Code availability

The code supporting the findings of this study are available from the
corresponding authors upon reasonable request.

## ACKNOWLEDGEMENT

We acknowledge valuable discussions with Javed Lindner, Damian Hofmann,
Claudia Merger and Lukas Grunwald. We further acknowledge support from the Max
Planck-New York City Center for Non-Equilibrium Quantum Phenomena. We
acknowledge funding by the DFG under Germany’s Excellence Strategy - Cluster
of Excellence Matter and Light for Quantum Computing (ML4Q) EXC 2004/1 -
390534769. NQS calculations have been performed using NetKet 3 38, 39 with jax
40. Computations were performed on the HPC system Ada at the Max Planck
Computing and Data Facility (MPCDF), and with computing resources granted by
RWTH Aachen University under project rwth0926.

## References

  * [1] Carleo, G. & Troyer, M.  Solving the quantum many-body problem with artificial neural networks.  _Science_ 355, 602–606 (2017).  URL https://doi.org/10.1126/science.aag2302.
  * [2] Deng, D.-L., Li, X. & Das Sarma, S.  Quantum entanglement in neural network states.  _Phys. Rev. X_ 7, 021021 (2017).  URL https://link.aps.org/doi/10.1103/PhysRevX.7.021021.
  * [3] Glasser, I., Pancotti, N., August, M., Rodriguez, I. D. & Cirac, J. I.  Neural-network quantum states, string-bond states, and chiral topological states.  _Physical Review X_ 8 (2018).  URL https://doi.org/10.1103/physrevx.8.011006.
  * [4] Kaubruegger, R., Pastori, L. & Budich, J. C.  Chiral topological phases from artificial neural networks.  _Physical Review B_ 97 (2018).  URL https://doi.org/10.1103/physrevb.97.195136.
  * [5] Choo, K., Neupert, T. & Carleo, G.  Two-dimensional frustrated j1-j2 model studied with neural network quantum states.  _Physical Review B_ 100 (2019).  URL https://doi.org/10.1103/physrevb.100.125124.
  * [6] Fabiani, G. & Mentink, J. H.  Investigating ultrafast quantum magnetism with machine learning.  _SciPost Phys._ 7, 004 (2019).  URL https://scipost.org/10.21468/SciPostPhys.7.1.004.
  * [7] Schmitt, M. & Heyl, M.  Quantum many-body dynamics in two dimensions with artificial neural networks.  _Physical Review Letters_ 125 (2020).  URL https://doi.org/10.1103/physrevlett.125.100503.
  * [8] Fabiani, G., Bouman, M. & Mentink, J.  Supermagnonic propagation in two-dimensional antiferromagnets.  _Physical Review Letters_ 127 (2021).  URL https://doi.org/10.1103/physrevlett.127.097202.
  * [9] Astrakhantsev, N. _et al._ Broken-symmetry ground states of the heisenberg model on the pyrochlore lattice.  _Physical Review X_ 11 (2021).  URL https://doi.org/10.1103/physrevx.11.041021.
  * [10] Bukov, M., Schmitt, M. & Dupont, M.  Learning the ground state of a non-stoquastic quantum Hamiltonian in a rugged neural network landscape.  _SciPost Phys._ 10, 147 (2021).  URL https://scipost.org/10.21468/SciPostPhys.10.6.147.
  * [11] Roth, C., Szabó, A. & MacDonald, A.  High-accuracy variational monte carlo for frustrated magnets with deep neural networks (2022).  URL https://arxiv.org/abs/2211.07749.
  * [12] Reh, M., Schmitt, M. & Gärttner, M.  Optimizing design choices for neural quantum states.  _Phys. Rev. B_ 107, 195115 (2023).  URL https://link.aps.org/doi/10.1103/PhysRevB.107.195115.
  * [13] Sharir, O., Shashua, A. & Carleo, G.  Neural tensor contractions and the expressive power of deep neural quantum states.  _Phys. Rev. B_ 106, 205136 (2022).  URL https://link.aps.org/doi/10.1103/PhysRevB.106.205136.
  * [14] Çeven, K., Oktel, M. O. & Keleş, A.  Neural-network quantum states for a two-leg bose-hubbard ladder under magnetic flux.  _Phys. Rev. A_ 106, 063320 (2022).  URL https://link.aps.org/doi/10.1103/PhysRevA.106.063320.
  * [15] Zou, E., Long, E. & Zhao, E.  Learning a compass spin model with neural network quantum states.  _Journal of Physics: Condensed Matter_ 34, 125802 (2022).  URL https://dx.doi.org/10.1088/1361-648X/ac43ff.
  * [16] Szoldra, T.  Fidelity and overlap of neural quantum states: Error bounds on the monte carlo estimator (2023).  eprint arXiv:2311.14820.
  * [17] Liu, Z. & Clark, B. K.  A unifying view of fermionic neural network quantum states: From neural network backflow to hidden fermion determinant states (2023).  eprint arXiv:2311.09450.
  * [18] Humeniuk, S., Wan, Y. & Wang, L.  Autoregressive neural Slater-Jastrow ansatz for variational Monte Carlo simulation.  _SciPost Phys._ 14, 171 (2023).  URL https://scipost.org/10.21468/SciPostPhys.14.6.171.
  * [19] Nomura, Y.  Boltzmann machines and quantum many-body problems.  _Journal of Physics: Condensed Matter_ 36, 073001 (2023).  URL https://dx.doi.org/10.1088/1361-648X/ad0916.
  * [20] Gao, X. & Duan, L.-M.  Efficient representation of quantum many-body states with deep neural networks.  _Nature Communications_ 8 (2017).  URL https://doi.org/10.1038/s41467-017-00705-2.
  * [21] Levine, Y., Sharir, O., Cohen, N. & Shashua, A.  Quantum entanglement in deep learning architectures.  _Phys. Rev. Lett._ 122, 065301 (2019).  URL https://link.aps.org/doi/10.1103/PhysRevLett.122.065301.
  * [22] Sun, X.-Q., Nebabu, T., Han, X., Flynn, M. O. & Qi, X.-L.  Entanglement features of random neural network quantum states.  _Phys. Rev. B_ 106, 115138 (2022).  URL https://link.aps.org/doi/10.1103/PhysRevB.106.115138.
  * [23] Carleo, G., Nomura, Y. & Imada, M.  Constructing exact representations of quantum many-body systems with deep neural networks.  _Nature Communications_ 9 (2018).  URL https://doi.org/10.1038/s41467-018-07520-3.
  * [24] Raghu, M., Poole, B., Kleinberg, J., Ganguli, S. & Sohl-Dickstein, J.  On the expressive power of deep neural networks (2016).  eprint arXiv:1606.05336.
  * [25] Georgevici, A. I. & Terblanche, M.  Neural networks and deep learning: a brief introduction.  _Intensive Care Medicine_ 45, 712–714 (2019).  URL https://doi.org/10.1007/s00134-019-05537-w.
  * [26] Sarker, I. H.  Deep learning: A comprehensive overview on techniques, taxonomy, applications and research directions.  _SN Computer Science_ 2 (2021).  URL https://doi.org/10.1007/s42979-021-00815-1.
  * [27] Chun, C. & Lee, D. D.  Sparsity-depth tradeoff in infinitely wide deep neural networks (2023).  eprint arXiv:2305.10550.
  * [28] Raza, S. & Ding, C.  Deep dynamic neural network to trade-off between accuracy and diversity in a news recommender system (2021).  eprint arXiv:2103.08458.
  * [29] Lin, S.-H. & Pollmann, F.  Scaling of neural-network quantum states for time evolution.  _physica status solidi (b)_ 259, 2100172 (2022).  URL https://doi.org/10.1002/pssb.202100172.
  * [30] Passetti, G. _et al._ Can neural quantum states learn volume-law ground states?  _Phys. Rev. Lett._ 131, 036502 (2023).  URL https://link.aps.org/doi/10.1103/PhysRevLett.131.036502.
  * [31] Brocki, L. & Chung, N. C.  Fidelity of interpretability methods and perturbation artifacts in neural networks (2022).  eprint arXiv:2203.02928.
  * [32] Poole, B., Lahiri, S., Raghu, M., Sohl-Dickstein, J. & Ganguli, S.  Exponential expressivity in deep neural networks through transient chaos.  In Lee, D., Sugiyama, M., Luxburg, U., Guyon, I. & Garnett, R. (eds.) _Advances in Neural Information Processing Systems_ , vol. 29 (Curran Associates, Inc., 2016).  URL https://proceedings.neurips.cc/paper/2016/file/148510031349642de5ca0c544f31b2ef-Paper.pdf.
  * [33] Schoenholz, S. S., Gilmer, J., Ganguli, S. & Sohl-Dickstein, J.  Deep information propagation.  _ArXiv_ abs/1611.01232 (2017).
  * [34] Tamai, K., Okubo, T., Duy, T. V. T., Natori, N. & Todo, S.  Absorbing phase transitions in artificial deep neural networks (2023).  eprint arXiv:2307.02284.
  * [35] Eisert, J., Cramer, M. & Plenio, M. B.  Colloquium: Area laws for the entanglement entropy.  _Reviews of Modern Physics_ 82, 277–306 (2010).  URL https://doi.org/10.1103/revmodphys.82.277.
  * [36] Bianchi, E., Hackl, L., Kieburg, M., Rigol, M. & Vidmar, L.  Volume-law entanglement entropy of typical pure quantum states.  _PRX Quantum_ 3, 030201 (2022).  URL https://link.aps.org/doi/10.1103/PRXQuantum.3.030201.
  * [37] Klambauer, G., Unterthiner, T., Mayr, A. & Hochreiter, S.  Self-normalizing neural networks (2017).  URL https://arxiv.org/abs/1706.02515.  eprint arXiv:1706.02515.
  * [38] Vicentini, F. _et al._ NetKet 3: Machine Learning Toolbox for Many-Body Quantum Systems.  _SciPost Phys. Codebases_ 7 (2022).  URL https://scipost.org/10.21468/SciPostPhysCodeb.7.
  * [39] Carleo, G. _et al._ NetKet: A machine learning toolkit for many-body quantum systems.  _SoftwareX_ 10, 100311 (2019).  URL https://doi.org/10.1016/j.softx.2019.100311.
  * [40] Bradbury, J. _et al._ JAX: composable transformations of Python+NumPy programs (2018).  URL http://github.com/google/jax.

## Author contributions

GP performed the calculations and DMK supervised the work. Both authors
collaborated in writing the manuscript

## Competing Interests

The authors declare no competing financial or nonfinancial interests.