# Predicting quantum dynamical cost landscapes with deep learning

Mogens Dalgaard  Department of Physics and Astronomy, Aarhus University, Ny
Munkegade 120, 8000 Arhus C, Denmark  Felix Motzoi<EMAIL_ADDRESS>Forschungszentrum Jülich, Institute of Quantum Control (PGI-8), D-52425
Jülich, Germany  Jacob Sherson  Department of Physics and Astronomy, Aarhus
University, Ny Munkegade 120, 8000 Arhus C, Denmark

(June 28, 2021)

###### Abstract

State-of-the-art quantum algorithms routinely tune dynamically parametrized
cost functionals for combinatorics, machine learning, equation-solving, or
energy minimization. However, large search complexity often demands many
(noisy) quantum measurements, leading to the increasing use of classical
probability models to estimate which areas in the cost functional landscape
are of highest interest. Introducing deep learning based modelling of the
landscape, we demonstrate an order of magnitude increases in accuracy and
speed over state-of-the-art Bayesian methods. Moreover, once trained the deep
neural network enables the extraction of information at a much faster rate
than conventional numerical simulation. This allows for on-the-fly
experimental optimizations and detailed classification of complexity and
navigability throughout the phase diagram of the landscape.

##  I Introduction

Figure 1: Machine learning the cost functional landscape structure of a spin
chain system. The diagrams represent two-dimensional cuts through the control
landscape, chosen to include two or three local optima. (a) An illustration of
how distinct optima may be connected. They can either be bridged together or
constitute separate peaks. In the latter case, they may also be of different
heights, with the lower one termed a local trap. (b) The actual fidelity
landscape evaluated at all points between three distinct optimized solutions
(taken randomly from the validation set) compared to the deep neural network
predictions with $g=0$ and $TJ=4.0$ in Eq. (2). The parentheses denote the
computational method and wall time consumption (i.e. the computational time).
We select three distinct solutions $v_{1}$, $v_{2}$, and $v_{3}$, which
depicts (top) a bridge and (bottom) a trap in the lower left corner.

Many quantum algorithms, notably noisy intermediate scale quantum (NISQ)
methods Arute  _et al._ (2019); Zhong  _et al._ (2020); Preskill (2018);
Peruzzo  _et al._ (2014); Kandala  _et al._ (2017); Kokail  _et al._ (2019);
Farhi  _et al._ (2000); Omran  _et al._ (2019); Farhi  _et al._ (2014);
Harrigan  _et al._ (2021), encode problem instances onto the parameter space
of a quantum device, improving the classical sampling cost. The device acts as
a parameter-tunable black box whose output can be optimized via adaptive
measurements. Determining the parameter regime to probe typically relies on
various heuristics for qualified initial ansätze and subsequent optimal
control Glaser  _et al._ (2015); Theis  _et al._ (2018); Müller  _et al._
(2021) or discrete optimization Maslov  _et al._ (2008); Motzoi  _et al._
(2017); Yao  _et al._ (2020); Zhu  _et al._ (2020). Popular heuristic-based
algorithms include variational Peruzzo  _et al._ (2014); Kandala  _et al._
(2017); Kokail  _et al._ (2019), adiabatic Farhi  _et al._ (2000); Omran  _et
al._ (2019), and parametric Farhi  _et al._ (2014); Harrigan  _et al._ (2021)
quantum circuits.

In this context, parameter initialization and optimization may benefit from
classical modelling. For example, machine learning can provide a powerful
ansatz for the many-body Schrödinger equation given by neural-network encoded
quantum states Carleo and Troyer (2017); Schmitt and Heyl (2020), be trained
to correct for measured quantum error syndromes Liu and Poulin (2019), design
improved experiments Goerz  _et al._ (2017); Menke  _et al._ (2021), and
obtain optimal values for quantum dynamics and circuit parametrizations using
reinforcement learning Bukov  _et al._ (2018); Niu  _et al._ (2019); Dalgaard
_et al._ (2020a). These modelling methods focus on parametrizing the behaviour
of one or several high-precision candidate solutions.

In contrast, the black box behaviour of the quantum device (encoding the
problem) may itself be learned, or metamodeled, to some helpful degree Wigley
_et al._ (2016); Bentley and Eisfeld (2018); Kokail  _et al._ (2019); Sauvage
and Mintert (2020); Paulson  _et al._ (2020); Koczor and Benjamin (2019).
Here, the modeling task is more challenging since it requires learning the
complete (many-to-many) dynamics mapping of the cost functional. Metamodeling
the quantum physical process can reduce ad hoc assumptions about good starting
guesses, problem difficulties, and algorithmic hyperparameters. It may also
avoid undersampling in areas of interest while oversampling or getting stuck
in ‘barren plateaus’ McClean  _et al._ (2018). Learning a metamodel may also
be helpful to transfer knowledge between instances of problems or devices.
Perhaps most importantly, offloading the bulk of the modeling cost to the
classical co-processor via a metamodel could greatly speed up NISQ devices.

Discussion of metamodeling is especially relevant to research in the last two
decades about cost functional landscapes, particularly for quantum optimal
control. Here the most salient open question has been the existence and
preponderance of suboptimal solutions (traps) Rabitz  _et al._ (2004); Pechen
and Tannor (2011); Werschnik and Gross (2007); Brif  _et al._ (2010); Caneva
_et al._ (2011); Larocca  _et al._ (2018); Bukov  _et al._ (2018), and other
measures of problem difficulties Shen  _et al._ (2006); Hsieh  _et al._
(2009); Nanduri  _et al._ (2013); Larocca  _et al._ (2018), important when
selecting suitable optimization algorithms and adjusting their parametrical
settings. For instance, discrete optimization methods are well known to
exhibit phase transitions between regimes of different difficulty Smith and
Dyer (1996); Xu and Li (2000); Gent and Walsh (1996), and recent work Bukov
_et al._ (2018); Day  _et al._ (2019) has shown that the quantum cost
functional landscape of binary-valued (i.e. bang-bang) control can be mapped
onto spin glass physics. This has lead to the proposal of universal behavior
of the control landscape under strongly constrained conditions. However, these
ideas have yet to be considered in the far more common setting of continuous-
value controlled Schrödinger evolution.

Recent results in metamodeling have focused on Bayesian estimation, using for
example Gaussian Process (GP) regression Wigley  _et al._ (2016); Bentley and
Eisfeld (2018); Kokail  _et al._ (2019); Sauvage and Mintert (2020); Paulson
_et al._ (2020), and trigonometric expansion of products of Pauli strings
Koczor and Benjamin (2019). The main bottleneck of these approaches is the low
precision of Bayesian estimation due to a large computational overhead, in
particular when large experimental data sets are needed Paulson  _et al._
(2020). In addition, the user-selected covariance function used for fitting
may also by restrictive, limiting generalizability and transferability across
problem choices.

In this work, we show that very high precision metamodeling of the complete
cost functional landscape can be attained using a deep neural network model.
The parametric tunability of the model allows the handling of a much larger
data throughput than earlier probabilistic methods, which, as we demonstrate
is generally needed to reach error rates at least as low as $10^{-3}$ for the
Ising spin chains considered. Even higher throughputs are expected when
sampling from quantum devices such as superconducting processors Walter  _et
al._ (2017), allowing prediction of the complete landscape to an even higher
precision. We also analyze the underlying structure of the cost functional
landscape for the Ising chains considered, and extract key measures that can
be used to select suitable optimization algorithms and adjust their parametric
setting. In the course of our study, we also identify continuous control
phases in the landscapes, where transitions between the phases can be
identified as the most challenging regimes, in contrast to earlier studies
Tibbetts  _et al._ (2012); Bukov  _et al._ (2018), and generalizing from the
binary control case Bukov  _et al._ (2018); Day  _et al._ (2019).

The paper is organized as follows: in Section II we demonstrate high precision
learning of a dynamical quantum cost functional landscape related to the
control of a many-body Ising spin chain. In Section III we consider the
scaling properties of various learning methods and compare to solving the
associated equations with matrix exponentiation. In Section IV we analyze the
underlying functional landscape structure and in Section V we conclude the
paper.

Figure 2: Illustration of a spin chain which consists of a series of coupled
two-level systems (spins). Here we model the spin-spin interactions as $ZZ$
terms, where adjacent spins are coupled with strength $J$ and next-nearest
spins are coupled with strength $g$. In our simulations, we further assume
open boundary conditions.

##  II Learning spin-chain dynamics

We study the quantum control cost landscape for a state-to-state transfer of a
spin-chain. Here, the objective is to manipulate a quantum system from an
initial $\ket{\psi_{0}}$ to a target state $\ket{\psi_{t}}$, up to an
inconsequential global phase, which is achieved by maximizing the fidelity

$\displaystyle F\big{[}\ket{\psi(T)}\big{]}=|\braket{\psi_{t}}{\psi(T)}|^{2},$
(1)

where $\ket{\psi(T)}$ denotes the solution to the Schrödinger equation at
final time $T$ starting from the initial state. Here $F=1$ implies a perfect
transfer; and the smallest amount of time where this is possible (at least to
a satisfactory degree) is called the quantum speed limit (QSL) Caneva  _et
al._ (2009).

We may induce the state-to-state transfer by controlling the quantum system
through externally applied pulses, whose shape can be optimized through
quantum control methods in order to reach local optima of the fidelity
function Eq. (1). When considering distinct optima, different situation can
occur where they may be bridged together as illustrated in Fig. 1(a, top) or
separated peaks as illustrated in Fig. 1(a, middle). The literature has
previously identified the existence of both Larocca  _et al._ (2018); Heck
_et al._ (2018); Larocca  _et al._ (2020). Note that a bridge could either be
a plateau of optima or consist of a set of near-optimal solutions connecting
two or more optima. Fig. 1(a, top) illustrates the latter. Moreover, separated
peaks may be of the same height as in Fig. 1(a, middle) or different heights
as in Fig. 1(a, bottom). In this context, solutions that are not global optima
are referred to as traps Rabitz  _et al._ (2004); Pechen and Tannor (2011);
Caneva  _et al._ (2011); De Fouquieres and Schirmer (2013); Riviello  _et al._
(2015); Rach  _et al._ (2015).

In this work we study a one-dimensional spin-chain as illustrated in Fig. 2.
The system is governed by the Hamiltonian

$\displaystyle
H(t)=-J\sum_{j}\sigma_{j}^{z}\sigma_{j+1}^{z}-g\sum_{j}\sigma_{j}^{z}\sigma_{j+2}^{z}+u(t)\sum_{j}\sigma_{j}^{x},$
(2)

where $J$ denotes the nearest spin-spin interaction, $g$ the next-nearest
interaction, and $u(t)$ is proportional to the amplitude of a global
transversal magnetic field in the $x$-direction, which we may control. The
time evolution we consider occurs via piecewise constant control fields for
$N$ equidistant time steps ($\Delta t=T/N$), such that $\ket{\psi(t+\Delta
t)}=\exp{(-iH(t)\Delta t)}\ket{\psi(t)}$ where $\hbar\equiv 1$. With this
choice, the fidelity functional, Eq. (1), becomes a function of the control
amplitudes $F(u_{1},u_{2},\ldots,u_{N})$ and thereby defines the multivariate
quantum dynamics landscape Rabitz  _et al._ (2004). We model five spins, with
$N=20$ time steps, periodic boundary conditions, and limit the time-dependent
control $u(t)/J\in[-1,+1]$. The initial and target states are respectively
chosen as the two degenerate ground states $\ket{00000}$ and $\ket{11111}$ of
Eq. (2) in the absence of control. In this case, the state transfer can only
be completed by populating a series of intermediate excited states. In the
following, we will investigate the cost functional landscape for two different
regimes: without ($g=0$) and with ($g=J/10$) next-nearest neighbor spin-spin
interactions.

Figure 3: Landscape prediction error as a function of control duration with
$N=20$ time steps and $L=5$ spins for (a) nearest neighbor interaction only
($g=0$). Here the parentheses denote the data size, and the shaded region the
quantum speed limit (QSL), which we estimate from the collected data. (b) We
also plot the predicted versus the actual infidelity at selected durations for
the deep neural network for 1000 points drawn at random from the validation
set, where the specified numbers denote the total control duration $TJ$. We
repeat the analysis for next-nearest neighbor interaction ($g=J/10$) in (c)
and (d).

The cost functional landscape is not simple enough to fully explore with
random sampling alone. Therefore, we utilize the second-order GRAPE
optimization algorithm Khaneja  _et al._ (2005); Motzoi  _et al._ (2011); De
Fouquieres  _et al._ (2011); Dalgaard  _et al._ (2020b), which seeks to
minimize the infidelity $\mathcal{C}=1-F$ by performing gradual updates
$\mathcal{C}(\mathbf{u}_{n+1})\leq\mathcal{C}(\mathbf{u}_{n})$ via gradient or
Hessian based optimization. We save the infidelity and pulse of each gradual
update $(\mathbf{u_{n}},\mathcal{C}(\mathbf{u_{n}}))$ made by GRAPE. This is
done for 10,000 randomly drawn pulses at each control duration, which
typically results in about one million collected data points. We split the
collected data into $80\%$ for training and $20\%$ for validation. The actual
structure of the cost functional landscape $\mathcal{C}(\mathbf{u})$ will be
discussed further in section IV.

We seek to learn the fidelity landscape given by Eq. (1) via a deep neural
network for which we use a feed forward network with several hidden layers
(for implementation details see Appendix A). The landscape is a very high
dimensional (20-dimensional) manifold that cannot easily be visualized.
However, we may still depict low-dimensional slices through this landscape. In
Fig. 1(b, top) we compare the predictions of a deep neural network trained on
the training set to the true landscape obtained by solving the Schrödinger
equation via matrix exponentiation, which is a standard integration technique
Blanes  _et al._ (2009). For comparison, we select three distinct,
representative optima $v_{1},v_{2},$ and $v_{3}$ from the validation set which
defines a 2D cut through the landscape, parametrized by

$\displaystyle v(\alpha,\beta)=v_{1}+\alpha(v_{2}-v_{1})+\beta(v_{3}-v_{1}).$
(3)

For both methods we use a $100\times 100$ equidistant grid with
$\alpha,\beta\in[-0.2,1.2]$, where we plot everything within the physical
boundaries ($u_{n}/J\in[-1,+1]$).

Fig. 1(b, top) depicts a comparison between the actual solution space and the
one predicted by the deep neural network for the same optima $v_{1},v_{2}$,
and $v_{3}$ chosen from the validation set. From the figure we see that the
deep neural network can very accurately recreate the fidelity landscape,
except for small and subtle differences. From the figure we also see the
existence of a bridge between two of the solutions. The analysis is repeated
for three different solutions in Fig.1(b, bottom), where we see the existence
of a trap, i.e., a solution that is not a global optimum.

Besides deep neural networks, we also assess the predictive power of a
different interpolation method, Gaussian processes regression Rasmussen
(2003). Gaussian processes are fast to train and work very well with smaller
data sets. In contrast, deep neural networks are known to take longer to
train, but can handle much larger amounts of data. We refer the reader to
Appendix A for a technical explanation on how these were implemented. To
evaluate the performances, we define the prediction error as the mean over
absolute differences $|\mathcal{C}^{\text{pred}}-\mathcal{C}^{\text{actual}}|$
between the predicted infidelity $\mathcal{C}^{\text{pred}}$ and the actual
infidelity $\mathcal{C}^{\text{actual}}$ in the validation set. As a
comparison baseline, we also use a linear model
$\mathcal{C}^{\text{pred}}(\mathbf{u})=\mathbf{w}^{T}\mathbf{u}+b$, where the
model parameters $(\mathbf{w},b)$ are found by linear least-square regression.

The collected data constitutes roughly one million pulses and infidelities per
control duration. Since Gaussian processes cannot handle such large quantities
of data, it is only tested on a subset of 1000 pulses and infidelities. The
results with nearest neighbor interactions only ($g=0$) is depicted in Fig.
3(a). Gaussian processes perform significantly better that the linear model
with around an order of magnitude improvement across different control
durations. The deep neural network performs even better, with approximately
one and two order of magnitude improvements over the Gaussian processes and
the linear model respectively. In Fig. 3(a), we show the quantum speed limit
(QSL), which we estimate from the collected data to be $T_{\text{QSL}}J\approx
3.33$. To verify that the deep neural network has learned to predict
infidelities, we plot in Fig. 3(b) the predicted versus the actual infidelity
for 1000 randomly selected pulses in the validation set. For all control
durations, the infidelities generally lie on the diagonal line, but with a few
points off at $TJ=4.0$.

We repeat the analysis with next-nearest neighboring interactions ($g/J=0.1$).
This significantly increases the complexity of the learning task and prolongs
the quantum speed limit, which is now around $T_{\text{QSL}}J\approx 5.0$. We
plot the results in Fig. 3(c), where we see some of the same tendencies as
before: the deep neural network performs significantly better than both the
linear model and the Gaussian process. However, now the prediction errors for
the machine learning algorithms are higher and scale worse with time, although
they still achieve overall very accurate predictions. Again, we plot the
predicted versus the actual infidelity of the deep neural network for 1000
randomly selected pulses in the validation set at selected control times in
Fig. 3(d). Here we observe a region of high (low) predictability at small
(large) control durations, which is directly related to the underlying
structure of the cost functional landscape. This we investigate further in
Section IV, where we relate the landscape structure to Fig. 3.

Figure 4: The scaling of the prediction errors when modeling Eq. 2 with
nearest neighbor interaction only (g=0). The prediction error calculated as
the average over the difference in predicted versus actual infidelities
$|\mathcal{C}^{\text{pred}}-\mathcal{C}^{\text{actual}}|$ from the validation
set. ($a$) Scaling with the number of spins $L$ using $N=20$ time steps at
$TJ=3.33$. Note, the prediction problem for the neural network
counterintuitively becomes simpler when more spins are modelled, unlike the
Gaussian process. $(b)$ Scaling with the number of time steps $N$ using $L=5$
spins also at $TJ=3.33$.

##  III Scaling and speed-up

Two questions naturally arise with the models presented so far: how does each
method scale in performance when the control problem is changed, and are there
potential advantages in their use?

To answer these questions, we focus on the example of modelling nearest spin-
spin interactions only ($g=0$) at the quantum speed limit ($TJ=3.33$) for five
spins. We start by looking at the performance of the various methods by
considering how they scale with the number of spins in the spin chain. This is
depicted in Fig. 4(a) where we have otherwise repeated the procedure from the
previous sections. Again, we see the same tendencies as before with the deep
neural network significantly outperforming the other methods. However, Fig.
4(a) also contains a rather counterintuitive result in that the deep neural
network, unlike the Gaussian process, performs better with increasing size of
the spin-chain, i.e., the control problem seems to become somewhat simpler
with increasing chain-size. We will return to this point later in this
section. In Fig. 4(b) we show the scaling with the number of time steps $N$
using $L=5$ spins. Again, we see the same tendencies as earlier when comparing
the different methods, but now with a slight increase in the deep neural
network prediction error with increasing $N$. In both cases, we see that
increasing the complexity of the dynamics does not significantly affect the
ability of the neural network to accurately predict the cost functional
landscape.

We now compare the wall time, i.e., the computational time for numerically
solving the Schrödinger equation via matrix exponentiation and for evaluating
the neural network. Matrix exponentiation scales with the size of the Hilbert
space ($d=2^{L}$), whereas the neural network size necessary for encapsulating
the dynamics scales with the complexity of the control task. For instance,
adding another qubit to the spin chain doubles the size of the Hilbert space,
but does not make the control task twice as difficult. As a matter of fact,
Fig. 4(a) points towards the complexity scaling favorably with the number of
spins, which we attribute to mean field effects causing the landscape to
become simpler with more spins. Note however that the Gaussian process does
not see an improvement with increasing $L$, indicating that the model
parametrization is key to this effect.

In Fig. 5(a) we compare the wall time consumption for numerically solving the
Schrödinger equation via matrix exponentiation with evaluating the neural
network. As expected, the neural network scales independently of the size of
the Hilbert space, whereas solving the Schrödinger equation via matrix
exponentiation scales very unfavorably with increasing dimensionality, leading
to several order of magnitude improvement with the neural network. There exist
of course other techniques which allow for handling larger Hilbert spaces
Perez-Garcia  _et al._ (2007); De Chiara  _et al._ (2008); Orús (2019), but
these still typically scale with the Hilbert space dimension, due to
calculation of the full state dynamics. Hence, using a deep neural network to
predict cost functional outcomes could lead to larger computational
improvements in applications where many repeated evaluations of the
Schrödinger equation are necessary. This must, of course, be compared with the
additional time of collecting data and training, for example through
experimental cost functional sampling. In Fig. 5(b), we also compare the wall
time with increasing number of time steps $N$. The width of the neural network
scales linearly with $N$ (see Appendix A) leading to a slight increase in wall
time with $N$. However, the neural network still performs much faster than
solving the Schrödinger equation with matrix exponentiation. From the figure
we also see a fluctuation around $N=40$, which we believe is due to some
numerical instabilities of measuring the wall time.

Figure 5: The wall time consumption for numerically solving the Schrödinger
equation via matrix exponentiation and evaluating the deep neural network.
Here each point corresponds to the evaluation of $1,000$ pulses. ($a$) The
scaling with an increasing number of spins $L$. Here we obtain a very
favorable scaling with the deep neural network, since its evaluation does not
scale with the size of the Hilbert space ($d=2^{L}$). ($b$) The scaling with
an increasing number of time steps $N$, which is also favorable to the deep
neural network.

At the quantum speed limit for $L=5$ spins and $N=20$ time steps, collecting
the data took around $22$ hours, while the subsequent training of the neural
network took around $37$ hours. From this data along with the speed-up
depicted in Fig. 5(b), we estimate that the two methods become comparable in
total computational time when the number of evaluations of the Schrödinger
equation exceeds a couple of million, which is the presently considered
regime. Moreover, data collection could be accelerated by performing parallel
simulations on different processors and both training and prediction of the
deep neural network could be accelerated by using modern GPUs rather than CPUs
as used in this work.

In an experimental setting, one can sample the data directly from the quantum
device, which for large Hilbert spaces will be much faster than classical
sampling from simulations. This can even be done on-the-fly for global
optimization, as e.g. in Ref. Kokail  _et al._ (2019), by combining with an
appropriate optimization algorithm. Once an approximate map of the device is
learned, one can also quickly extract general features which may transfer from
one device or problem to another. Such general landscape features are
discussed in the next section. In addition to considering the scaling with
number of time steps and spins, we also investigate the scaling of the deep
neural network with the size of the training set. In Fig. 6 we plot the
prediction error of the network at different sizes $M$ of the training set,
with $L=5$ spins and $N=20$ time steps, and at the quantum speed limit
$TJ=3.33$. We make the empirical observation that all of the data points,
except the first, lie on a straight line in the log-scale plot, which is
consistent with a scaling proportional to $M^{-3/4}$. We depict scaling in the
figure with a black dashed line. Hence, increasing the size of the training
set by a factor of 10 would reduce the prediction error with a factor of
$10^{-3/4}=0.1778\sim 1/5$ for this particular problem. It is especially
noteworthy that the error curve plotted in the figure does not saturate within
the considered sizes of data, indicating that significantly lower prediction
errors would be achievable with access to more data. Although this is outside
the computational scope of our computing resources, in an experimental
setting, it is possible to obtain larger data sets directly from the quantum
device. For example superconducting qubits can have very high data acquisition
rates Walter  _et al._ (2017), even scaling favourably with increasing Hilbert
space sizes, thus potentially going well beyond the data sizes considered
here.

Figure 6: The prediction errors of the deep neural network at different size
$M$ of the training data set. The results are consistent with an overall
scaling proportional to $M^{-3/4}$, implying that much lower prediction errors
are achievable with more training data.

##  IV The landscape structure

Having obtained a highly accurate mapping of the dynamical quantum fidelity
landscape, we turn our attention to extracting some of its most important
characteristics.

###  IV.1 Landscape properties

The structure of the cost landscape not only contributes directly to our
understanding of the underlying dynamics but is also important from an
optimization point of view, since given optimization techniques can have both
advantages and disadvantages for different problems. For instance, in Ref.
Dalgaard  _et al._ (2020a) machine learning-enhanced exploration outperformed
optimization through random seeding on a problem with combinatorially few
global optima, while hill-climbing approaches can suffer barren plateaus and
vanishing gradients when random sampling is used McClean  _et al._ (2018).
Furthermore, in the discrete optimization setting it is well known that
different problem classes can map onto classical phase transitions exhibiting
universal behaviour Smith and Dyer (1996); Xu and Li (2000); Gent and Walsh
(1996). This has recently been extended to the discrete control (binary
constrained) quantum optimization setting where a connection is made to
classical spin glass transitions Bukov  _et al._ (2018); Day  _et al._ (2019).
However, relatively little is known about the continuous control optimization
landscape, which is far more common given current experimental capabilities,
while also benefiting from greater controllability. Thus, it is broad
significance to extend these notions of problem difficulty to a continuous
setting.

Earlier, in Fig. 1(a) we depicted examples of how the fidelity landscape may
look on a local scale. In the following we now seek to investigate the global
structure of the landscape. For this purpose, we define a few different
measures, which are intended to capture the difficulty of the quantum
optimization task as well as the global landscape structure. A schematic of
three important characteristics is given in Fig. 7 and this provides the basis
for our measures, with opposite regimes for each feature being shown.

The first is a characterization of the distance a local optimizer must travel
in the quantum landscape in order to reach an optimimum, which is determined
by the _the density of attractors_. An attractor is a local optimum (or
connected set of optima), towards which local gradient- or Hessian-based
optimization will gradually converge. The landscape will have zero derivatives
at the point of an attractor unless if it lies on the boundary of the
landscape.

Secondly, while the density of attractors determines the average distance a
local optimizer must traverse before convergence, the speed of convergence to
the local attractor may also vary greatly between different landscapes.
Generally, local hill-climbing (such as gradient- or Hessian-based)
algorithms, may perform larger gradual updates on slowly changing landscapes
and therefore converge at a faster rate. Hence, we aim to quantify the
_ruggedness_ of the terrain, i.e., the rate at which the landscape changes
with respect to changes in the control parameters. A landscape that changes
slowly we describe as smooth, while we describe a landscape that changes
rapidly as rugged. A smooth landscape will also have the additional benefit
that functional values near an optimum will not vary greatly, thereby also
providing _robustness_ to experimental imperfections.

Thirdly, the optimization effort is greatly influenced by the existence of
suboptimal local attractors, hence the last characteristic is the _the density
of traps_. A trap is an attractor that is not a global optimum, where local
gradient- or Hessian-based optimization will still converge towards, i.e.
become trapped. Therefore, finding the global optima would be difficult for a
quantum landscape with many traps, often requires a global landscape search
strategy. In the opposite regime where there are no traps, local optimization
would converge to an optimum with unit probability. For this reason, the
existence of traps has attracted much attention in the quantum control
literature Rabitz  _et al._ (2004); Pechen and Tannor (2011); Caneva  _et al._
(2011); De Fouquieres and Schirmer (2013); Riviello  _et al._ (2015); Rach
_et al._ (2015); Jensen  _et al._ (2021).

Figure 7: Illustration of opposite regimes of three principle characteristics
of the landscape structure. These are: (top) the density of attractors,
(middle) the ruggedness, and (bottom) the density of traps. Here the landscape
is depicted as a maximization task.

###  IV.2 Landscape measures

We now relate the landscape characteristics described above and illustrated in
Fig. 7 to different quantitative measures of control optimization difficulty.
The first measure is simply the largest (or best) optimized fidelity
$F^{*}=\max_{j}F_{j}$, where $\\{F_{j}\\}$ denotes the set of optimized
fidelities Eq. (1). This corresponds to the tallest of all the peaks in Fig.
7.

In order to describe the computational cost of local optimization we propose
two different measures. The first is the optimization distance, i.e., the
distance from the initial randomly drawn pulse to the final optimized pulse
using BFGS-based GRAPE

$\displaystyle
D_{\text{attractor}}=\frac{1}{Nu_{\max}|\mathcal{D}|}\sum_{\mathbf{u}\in\mathcal{D}}||\mathbf{u}_{\text{optim}}-\mathbf{u}_{\text{initial}}||.$
(4)

Here $N$ is the number of time steps, $u_{\max}=J$ the amplitude bounds
imposed on the pulses, $\mathcal{D}$ denotes the set of saved optimization
trajectories with size $|\mathcal{D}|$ (i.e. the number of random seeds
optimized), and $||\bullet||$ denotes the Euclidean distance. The right hand
side of the above equation is simply the normalized Euclidean distance from
the initial to the optimized pulse averaged over all iterations. This will be
small for a landscape with many attractors and large for a landscape with few
attractors. This measure has previously been studied in the literature to
assess the straightness of local optimization paths Nanduri  _et al._ (2013);
Larocca  _et al._ (2018).

In addition to this measure, we also consider the fraction of repeated
optimized pulses, $1-|\mathcal{D^{*}}|/|\mathcal{D}|$ as in Ref. Day  _et al._
(2019) with $\mathcal{D}^{*}$ denoting the set of unique optimal pulses. Here
we consider two optimized pulses to be different if their normalized Euclidean
distance is larger than $10^{-9}$.

To describe the ruggedness of the control landscape we propose the mean
diagonal Hessian element at the optima

$\displaystyle\rho_{\text{rugged}}=\frac{1}{N|\mathcal{D}|}\sum_{\mathbf{u}\in\mathcal{D}}\sum_{j=1}^{N}\frac{\partial^{2}\mathcal{C}}{\partial
u(t_{j})^{2}}\bigg{|}_{\text{optima}},$ (5)

where $\mathcal{C}=1-F$ denotes the infidelity, with $F$ given by Eq. (1). The
Hessian diagonal describes the rate of change of the gradient, which is
relatively small for a smooth landscape that changes slowly and in contrast
relatively large for a rugged landscape. The eigenvalue decomposition of the
Hessian has previously been utilized to analyze quantum control landscapes
Shen  _et al._ (2006); Hsieh  _et al._ (2009). For our calculations we use the
analytical form of the Hessian found in Ref. Dalgaard  _et al._ (2020b).

In order to describe the density of traps in the control landscape we use the
variance over optimized fidelities, which is zero for a trap-free landscape.

Figure 8: Phases in the control landscape. The best fidelity together with
different measures that characterizes the trap-density, attractor-density, and
the smoothness for ($a$) nearest spin-spin interaction ($g=0$) and ($b$) with
next nearest spin-spin interaction as well $g=0.1J$. Here the two vertical
dashed lines marks phase transitions i.e., points in the total control
duration where the control landscape smoothly changes character, marked by the
peak of Attractor-distance II and Trap-density measure (see text). The shaded
region marks the quantum speed limit estimated from the data.

We will now use the measures given above to determine which characteristics
(c.f. Fig. 7) describe the landscape and its respective phases and relate this
to the potential difficulties of numerical or experimental pulse optimization.

###  IV.3 Continuous control phases

We first consider the different phases of the continuous control landscape
that emerge as the total control duration is changed. This generalizes on the
well-studied cases of discrete optimization Parker and Rardin (2014), which
can be characterized through several measures of optimization-hardness. First,
there is computational complexity Bernstein and Vazirani (1997), which is a
measure of how the worst-case computational-time effort scales with problem
size when finding a global optimum. In contrast, there is also the phase
transition language from condensed matter physics, which discusses average
properties of (energy) landscapes, and draws conclusions about critical
behaviour between different phases or regimes of the landscape Smith and Dyer
(1996); Xu and Li (2000); Gent and Walsh (1996). In particular, the latter is
more suited for discussion of optimization difficulty since measures of
average effort are of the greatest practical relevance (for example the
computational complexity does not capture the difficulty in finding a local
optimum).

The landscape measures that we use generalize the notion of global vs. local
behaviour, which is a central consideration both for combinatorial
optimization and for minimization in energy landscapes. In particular, Ref.
Day  _et al._ (2019) identifies two measures related to the Edwards-Anderson
order parameter: the first is a local measure that calculates the correlation
between optima where improvements are only possible by flipping a single
binary control at a time; the second is a more global measure that does the
same with two controls at a time. The authors find that these correspond
respectively to phase transitions, from overconstained to correlated, and then
to underconstrained optimization. The latter transition also lines up with the
quantum speed limit in their system.

Because the more general continuous-valued pulses used in our work do not
exhibit this discrete bang-bang nature, the measures outlined in the previous
section are instead used, indicating both local difficulty (distance and
ruggedness) and global difficulty (trapping frequency). These are calculated
for different control durations and shown in Fig. 8 for both (a) nearest-
neighbour interactions ($g=0$) and (b) next-nearest-neighbour interactions
($g=J/10$).

Although our measures are markedly different than the discrete case, we find
similar conclusions in terms of over$-,$ critically$-$, and under-constrained
dynamics. In both panels, we make the empirical observation that there exist
two instances of the total control duration, where the landscapes change
behavior, which are the turning points of the density of optima (fraction of
repeated solutions) and the density of traps (infidelity variance). This gives
us three distinct control regimes, marked with dashed lines, which we can
relate to the notion of phases. Note, that the marked phases occur at
different times on the two figures.

In essence, we see in Fig. 8 that these transitions correspond to two
tradeoffs, one corresponding to local effects and the other to global ones.
Locally, as the evolution time increases, the number of different dynamical
trajectories increases, which can be seen in the increased distance between
optima (green squares). Increasing at the same time is the fraction of
repeated solutions, as the increased controllability translates into less
spatial constraining (green diamonds). These are at odds with each other
eventually, and as the distance continues to grow, we see the fraction of
repeated solutions abruptly start to drop, indicating the attraction to
distinct minima in the landscape. However, these continue to have roughly the
same fidelity, as the controls have limited effect at such short durations.
Note also the difference from the discrete case in Ref. Bukov  _et al._
(2018), where overconstraining leads to a convex landscape around a single
optimum at short times. In the present case, the control space (space of
possible controls) is significantly larger, leading to a significantly larger
set of optimal pulses, even at relatively short control durations.

The second region sees the appearance of a second tradeoff, where, as the
fraction of distinct optima increases, the globally optimal fidelity also
continues to improve (as a result of increased controllability) and so
localization of solutions occurs. Therefore, we see the emergence of traps
(orange triangles). Abstractly, this can be seen as critical constraining,
that is where the search space becomes commensurate with the number of
constraints and so satisfying the maximum number of criteria (for maximal
fidelity) becomes increasingly difficult. Thus, we see the landscape changing
from the situation depicted in Fig. 1(b) in the top panel increasingly to that
of the bottom panel.

In the third region, the trapping begins to subside as we also see the average
distance between optima decreasing from its maximum in Fig. 8. This is the so-
called underconstrained phase, where the number of repeated solutions is very
small but nonetheless the constraints for high fidelity are still easier to
achieve. It is also here that we achieve full controllability, with the
quantum speed limit existing in this phase. In the continuum case, this is
also where we expect to see bridge solutions Shen  _et al._ (2006); Hsieh  _et
al._ (2009), that is to say, solution spaces described by a continuous
symmetry.

Finally, it is worth noting differences with conclusions otherwise drawn in
the literature. One key finding is that the quantum speed limit and the phase
transition need not to occur at the same point. Rather, we see that the
maximum trapping frequency region can actually be quite separate from the
quantum speed limit, in contrast to the common conjecture Bukov  _et al._
(2018); Day  _et al._ (2019); Tibbetts  _et al._ (2012). In particular, this
calls into question the generality of the finding in Tibbetts  _et al._ (2012)
that a superexponential effort is needed near the quantum speed limit, as the
trapping frequency may actually be quite small at this control duration. In
the next subsection we examine more closely the question of problem difficulty
from an optimization and learning perspective.

###  IV.4 Optimization difficulty

The landscape of a difficult optimization task is quantified by having large
attractor distance, high trap density, and/or being rugged. In fact, from
these considerations, it is plain to see in Fig. 8 that the first and second
critical durations correspond to where the computational effort of
optimization is the largest. Moreover, these are the largest local and global
optimization efforts, respectively.

At the first transition, characterized where the relative optimization
distance (II) is at a maximum, the other measures of difficulty remain fairly
low. In particular, trapping probability is quite low and the landscape is
very smooth. We can draw a few conclusions about the computational effort
required for different types of algorithms and problems based on the different
measures. Because the trapping is quite low, the complexity in terms of the
number of controls can be relatively benign. That is, a local search may
suffice in the sense that any random control seeding will converge to a near
optimal solution with high probability. This discounts the need for global
optimizers in this regime. On the other hand, a low number of repeated
solutions can be potentially problematic in terms of the Hilbert space
complexity. The large distance between peaks can lead to the emergence of long
flat portions of the landscape known as barren plateaus, which may lead to
exponential complexity as a function of the number of spins McClean  _et al._
(2018). In particular, this demands the use of good initial guesses or global
machine learning based seeding Dalgaard  _et al._ (2020a). Nonetheless, the
smoothness of the landscape indicates that one can improve on basic hill
climbers such as simplex search and gradient descent by incorporating
information about the curvature of the landscape Dalgaard  _et al._ (2020b),
with potentially near quadratic improvement in the convergence. Finally, we
note that the distance measure can be used to tune common hyperparameters,
such as the size of the initial simplex in Nelder-Mead simplex search Nelder
and Mead (1965).

The second transition point occurs where the trapping frequency is maximal.
This corresponds in the literature respectively to critical constraining in
terms of, e.g., constraint satisfaction problems, and frustration in the
minimization of energy for many body physics. In contrast to the first
transition, we see that it is not certain to find a global optimum with local
optimization, and therefore random sampling may be insufficient to achieve
extremely low infidelity, as needed for certain quantum information tasks. In
addition, the Hilbert space complexity may also play a significant role again,
with optimization distance (I) near maximal in both panels, and hybrid global-
local algorithms Dalgaard  _et al._ (2020a) should be implemented with
carefully designed stopping criteria for the local optimizer to avoid
vanishing gradients. The optimization difficulty may also be linked in this
regime to the landscape learning difficulty itself. In Fig. 3 we see that the
critical points ($TJ\approx 3$) are exactly where local maxima occur in the
landscape learning error (within a larger trend of increasing difficulty with
evolution duration). In this sense, both training of controls and of fidelity
functionals suffer from variance in the peaks and critical constraining. Thus
we may argue that the empirically observed difficulty in learning the
landscape and minimizing infidelity can act as proxies for each other as
difficulty measures. Interestingly, this does not appear to be the case for
the quantum speed limits, where the shaded regions in Fig. 8 correspond to
minima rather than maxima in Fig. 3.

Finally, we comment on the longstanding debate about the trap frequency in
dynamics optimization for typical problems Rabitz  _et al._ (2004); Pechen and
Tannor (2011); Werschnik and Gross (2007); Brif  _et al._ (2010); Caneva  _et
al._ (2011); Larocca  _et al._ (2018); Bukov  _et al._ (2018). It is now
generally accepted that constraints are indeed commonly leading to trapping in
the cost landscape. However, it is remarkable that the ‘easy’ phase of the
landscape attributed to Rabitz  _et al._ (2004) is not only occurring far to
the right of the quantum speed limit but also to its left. Indeed, at shorter
times, decoherence is expected to play a much smaller role, and so it may be
of more practical relevance. Thus it seems that for this fairly standard
problem considered here, the most important region, namely the quantum speed
limit, is actually also in the ‘easy’ regime, and therefore not prone to
trapping while remaining of highest practical interest. Of course, this does
not exclude the evidence that for other problems the quantum speed limit and
critical duration may be much closer to each other. Moreover, although strict
trapping may not be a problem, the increased ruggedness at these times
indicates limited utility for second order methods and the possibility that
gradient descent could still suffer from vanishing gradients.

##  V Conclusion

In this work, we have shown that a highly complicated quantum dynamics cost
landscape, the fidelity landscape, can be learned with very high precision
using deep learning, improving by an order of magnitude over Bayesian
estimation based methods, such as Gaussian process regression. Despite being
able to correctly predict tens of millions of points to high accuracy,
evaluating the neural network remains orders of magnitude faster than
simulating the Schrödinger evolution directly. Thus, when combined with direct
sampling of experiments, a complete mapping of many typical experimental
dynamic protocols should be possible with high throughput.

This mapping may not only enable on-the-fly global optimization towards true
optima in the system, but allows us to faithfully characterize given
experimental systems in terms of their phase diagram and difficulty measures.
These measures generalize notions of difficulty from the discrete optimization
case, and provide key insights into algorithmic choices, e.g. global vs local,
gradient vs. Hessian, and their hyperparameters, e.g. distances between
optima, smoothness, etc. Combining the developed learning and characterization
methods thus provides valuable tools and insights for both specific instances
and wider problem classes of parametric quantum optimization.

##  VI Acknowledgements

The authors would like to thank Carrie Ann Weidner, Marin Bukov and Martino
Calzavara for their help in completing this paper. This work was funded by the
Carlsberg Foundation and by the Deutsche Forschungsgemeinschaft (DFG, German
Research Foundation) under Germany’s Excellence Strategy – Cluster of
Excellence Matter and Light for Quantum Computing (ML4Q) EXC 2004/1 –
390534769. The numerical results presented in this work were obtained at the
Centre for Scientific Computing, Aarhus phys.au.dk/forskning/cscaa.

##  Appendix A Implementation details

We used our own implementation of GRAPE, with exact gradients using a L-BFGS-B
optimizer from Ref. Virtanen  _et al._ (2020). In the following we briefly
comment on Gaussian process regression and deep neural networks.

A Gaussian process regression setup Rasmussen (2003) consists of a previous
set of observed data $(X_{1},\mathbf{y_{1}})$, for which we desire to build a
model $f$ that allows us to make new predictions $\mathbf{y_{2}}=f(X_{2})$,
while still producing the old results $f(X_{1})=\mathbf{y_{1}}$. We treat this
via a multivariate normal distribution

$\displaystyle\begin{bmatrix}\mathbf{y}_{1}\\\
\mathbf{y}_{2}\end{bmatrix}\sim\mathcal{N}\Bigg{(}\begin{bmatrix}\boldsymbol{\mu}_{1}\\\
\boldsymbol{\mu}_{2}\end{bmatrix},\begin{bmatrix}\boldsymbol{\Sigma}_{11}&\boldsymbol{\Sigma}_{12}\\\
\boldsymbol{\Sigma}_{21}&\boldsymbol{\Sigma}_{22}\end{bmatrix}\Bigg{)},$ (6)

where $\boldsymbol{\mu}$ and $\boldsymbol{\Sigma}$ denotes the mean and
covariances respectively. Where we model the covariances via a kernel function
such as a radial basis function $i\neq j$

$\displaystyle\boldsymbol{\Sigma}_{i,j}=k(\mathbf{x_{i}},\mathbf{x_{j}})=\exp\bigg{(}-\frac{||\mathbf{x_{i}}-\mathbf{x_{j}}||^{2}}{2l^{2}}\bigg{)}.$
(7)

Here $l$ denotes a length scale parameter that we along with the mean
$\boldsymbol{\mu}=[\boldsymbol{\mu}_{1},\boldsymbol{\mu}_{2}]$ seek to fit in
order to maximize the marginal likelihood over the data

$\displaystyle
p(\mathbf{y}|\boldsymbol{\mu},\boldsymbol{\Sigma})=\frac{1}{(2\pi^{d}|\boldsymbol{\Sigma}|)}\exp\bigg{(}-\frac{1}{2}(\mathbf{y}-\boldsymbol{\mu})^{T}\boldsymbol{\Sigma}^{-1}(\mathbf{y}-\boldsymbol{\mu})\bigg{)}.$
(8)

Here $\mathbf{y}=[\mathbf{y}_{1},\mathbf{y}_{2}]$ and $\boldsymbol{\Sigma}$
the matrix containing all covariances. For Gaussian processes we use the
implementation from Ref. Pedregosa  _et al._ (2011). We used radial basis
functions as kernels and the L-BFGS-B for optimization, with five restarts per
fit.

A neural network consists of layers, where each layer consists of a set of
artificial neurons. Each neuron in a given layer receive an input signal
$x_{j}$ for each neuron in the previous layer and based on this calculates an
output signal $y$ that is transmitted to each neuron in the next layer and so
on. The output signal is calculated via an activation function $y=a(z)$, where
$z=\sum_{j}w_{j}x_{j}+b$ where $w$ and $b$ denotes the weights and biases of a
given neuron. These constitute the model parameters we seek to fit in order to
minimize a given cost function. A deep neural network is a network that has
many layers (and perhaps several million parameters) and hence is expensive to
evaluate, but with the possible advantage that is may fit highly complicated
models. For the deep neural network, we use the implementation from Ref. Abadi
_et al._ (2016). The input layer was linear $a(z)=z$, followed by five hidden
layers using the ReLu activation function $a(z)=\max(0,z)$, and the width
(neurons per layer) of the neural network scaling with the number of time
steps as $20N$. The output layer was also linear, and we used biases for all
layers. For optimization, we used the Adam optimizer, the loss function was
mean squared error, and we used a batch size of 16 with early stopping and a
learning rate of $10^{-5}$.

## References

  * Arute  _et al._ (2019) F. Arute, K. Arya, R. Babbush, D. Bacon, J. C. Bardin, R. Barends, R. Biswas, S. Boixo, F. G. Brandao, D. A. Buell, _et al._ , Nature 574, 505 (2019).
  * Zhong  _et al._ (2020) H.-S. Zhong, H. Wang, Y.-H. Deng, M.-C. Chen, L.-C. Peng, Y.-H. Luo, J. Qin, D. Wu, X. Ding, Y. Hu, _et al._ , Science 370, 1460 (2020).
  * Preskill (2018) J. Preskill, Quantum 2, 79 (2018).
  * Peruzzo  _et al._ (2014) A. Peruzzo, J. McClean, P. Shadbolt, M.-H. Yung, X.-Q. Zhou, P. J. Love, A. Aspuru-Guzik, and J. L. O’Brien, Nat. Commun. 5, 1 (2014).
  * Kandala  _et al._ (2017) A. Kandala, A. Mezzacapo, K. Temme, M. Takita, M. Brink, J. M. Chow, and J. M. Gambetta, Nature 549, 242 (2017).
  * Kokail  _et al._ (2019) C. Kokail, C. Maier, R. van Bijnen, T. Brydges, M. K. Joshi, P. Jurcevic, C. A. Muschik, P. Silvi, R. Blatt, C. F. Roos, _et al._ , Nature 569, 355 (2019).
  * Farhi  _et al._ (2000) E. Farhi, J. Goldstone, S. Gutmann, and M. Sipser, quant-ph/0001106 (2000).
  * Omran  _et al._ (2019) A. Omran, H. Levine, A. Keesling, G. Semeghini, T. T. Wang, S. Ebadi, H. Bernien, A. S. Zibrov, H. Pichler, S. Choi, _et al._ , Science 365, 570 (2019).
  * Farhi  _et al._ (2014) E. Farhi, J. Goldstone, and S. Gutmann, arXiv:1411.4028 (2014).
  * Harrigan  _et al._ (2021) M. P. Harrigan, K. J. Sung, M. Neeley, K. J. Satzinger, F. Arute, K. Arya, J. Atalaya, J. C. Bardin, R. Barends, S. Boixo, _et al._ , Nat. Phys. 17, 332 (2021).
  * Glaser  _et al._ (2015) S. J. Glaser, U. Boscain, T. Calarco, C. P. Koch, W. Köckenberger, R. Kosloff, I. Kuprov, B. Luy, S. Schirmer, T. Schulte-Herbrüggen, _et al._ , Eur. Phys. J. D 69, 279 (2015).
  * Theis  _et al._ (2018) L. Theis, F. Motzoi, S. Machnes, and F. Wilhelm, EPL 123, 60001 (2018).
  * Müller  _et al._ (2021) M. M. Müller, R. S. Said, F. Jelezko, T. Calarco, and S. Montangero, arXiv:2104.07687 (2021).
  * Maslov  _et al._ (2008) D. Maslov, G. W. Dueck, D. M. Miller, and C. Negrevergne, IEEE Trans. Comput.-Aided Des. Integr. Circuits Syst. 27, 436 (2008).
  * Motzoi  _et al._ (2017) F. Motzoi, M. P. Kaicher, and F. K. Wilhelm, Phys. Rev. Lett. 119, 160503 (2017).
  * Yao  _et al._ (2020) J. Yao, M. Bukov, and L. Lin, “Policy gradient based quantum approximate optimization algorithm,” (2020), arXiv:2002.01068 [quant-ph] .
  * Zhu  _et al._ (2020) L. Zhu, H. L. Tang, G. S. Barron, F. Calderon-Vargas, N. J. Mayhall, E. Barnes, and S. E. Economou, arXiv:2005.10258 (2020).
  * Carleo and Troyer (2017) G. Carleo and M. Troyer, Science 355, 602 (2017).
  * Schmitt and Heyl (2020) M. Schmitt and M. Heyl, Phys. Rev. Lett. 125, 100503 (2020).
  * Liu and Poulin (2019) Y.-H. Liu and D. Poulin, Phys. Rev. Lett. 122, 200501 (2019).
  * Goerz  _et al._ (2017) M. H. Goerz, F. Motzoi, K. B. Whaley, and C. P. Koch, npj Quantum Inf. 3, 1 (2017).
  * Menke  _et al._ (2021) T. Menke, F. Häse, S. Gustavsson, A. J. Kerman, W. D. Oliver, and A. Aspuru-Guzik, npj Quantum Inf. 7, 1 (2021).
  * Bukov  _et al._ (2018) M. Bukov, A. G. Day, D. Sels, P. Weinberg, A. Polkovnikov, and P. Mehta, Phys. Rev. X 8, 031086 (2018).
  * Niu  _et al._ (2019) M. Y. Niu, S. Boixo, V. N. Smelyanskiy, and H. Neven, npj Quantum Inf. 5, 1 (2019).
  * Dalgaard  _et al._ (2020a) M. Dalgaard, F. Motzoi, J. J. Sorensen, and J. Sherson, npj Quantum Inf. 6 (2020a).
  * Wigley  _et al._ (2016) P. B. Wigley, P. J. Everitt, A. van den Hengel, J. W. Bastian, M. A. Sooriyabandara, G. D. McDonald, K. S. Hardman, C. D. Quinlivan, P. Manju, C. C. Kuhn, _et al._ , Sci. Rep. 6, 1 (2016).
  * Bentley and Eisfeld (2018) C. Bentley and A. Eisfeld, Journal of Physics B: Atomic, Molecular and Optical Physics 51, 205003 (2018).
  * Sauvage and Mintert (2020) F. Sauvage and F. Mintert, Phys. Rev. X 1, 020322 (2020).
  * Paulson  _et al._ (2020) D. Paulson, L. Dellantonio, J. F. Haase, A. Celi, A. Kan, A. Jena, C. Kokail, R. van Bijnen, K. Jansen, P. Zoller, _et al._ , arXiv:2008.09252 (2020).
  * Koczor and Benjamin (2019) B. Koczor and S. C. Benjamin, arXiv:1912.08660 (2019).
  * McClean  _et al._ (2018) J. R. McClean, S. Boixo, V. N. Smelyanskiy, R. Babbush, and H. Neven, Nat. Commun. 9, 1 (2018).
  * Rabitz  _et al._ (2004) H. A. Rabitz, M. M. Hsieh, and C. M. Rosenthal, Science 303, 1998 (2004).
  * Pechen and Tannor (2011) A. N. Pechen and D. J. Tannor, Phys. Rev. Lett. 106, 120402 (2011).
  * Werschnik and Gross (2007) J. Werschnik and E. Gross, J. Phys. B 40, R175 (2007).
  * Brif  _et al._ (2010) C. Brif, R. Chakrabarti, and H. Rabitz, New Journal of Physics 12, 075008 (2010).
  * Caneva  _et al._ (2011) T. Caneva, T. Calarco, and S. Montangero, Phys. Rev. A 84, 022326 (2011).
  * Larocca  _et al._ (2018) M. Larocca, P. M. Poggi, and D. A. Wisniacki, J. Phys. A 51, 385305 (2018).
  * Shen  _et al._ (2006) Z. Shen, M. Hsieh, and H. Rabitz, J. Chem. Phys. 124, 204106 (2006).
  * Hsieh  _et al._ (2009) M. Hsieh, R. Wu, and H. Rabitz, J. Chem. Phys. 130, 104109 (2009).
  * Nanduri  _et al._ (2013) A. Nanduri, A. Donovan, T.-S. Ho, and H. Rabitz, Phys. Rev. A 88, 033425 (2013).
  * Smith and Dyer (1996) B. M. Smith and M. E. Dyer, Artif. Intell. 81, 155 (1996).
  * Xu and Li (2000) K. Xu and W. Li, J. Artif. Intell. Res. 12, 93 (2000).
  * Gent and Walsh (1996) I. P. Gent and T. Walsh, Artif. Intell. 88, 349 (1996).
  * Day  _et al._ (2019) A. G. Day, M. Bukov, P. Weinberg, P. Mehta, and D. Sels, Phys. Rev. Lett. 122, 020601 (2019).
  * Walter  _et al._ (2017) T. Walter, P. Kurpiers, S. Gasparinetti, P. Magnard, A. Potočnik, Y. Salathé, M. Pechal, M. Mondal, M. Oppliger, C. Eichler, _et al._ , Phys. Rev. Appl. 7, 054020 (2017).
  * Tibbetts  _et al._ (2012) K. W. M. Tibbetts, C. Brif, M. D. Grace, A. Donovan, D. L. Hocker, T.-S. Ho, R.-B. Wu, and H. Rabitz, Phys. Rev. A 86, 062309 (2012).
  * Caneva  _et al._ (2009) T. Caneva, M. Murphy, T. Calarco, R. Fazio, S. Montangero, V. Giovannetti, and G. E. Santoro, Phys. Rev. Lett. 103, 240501 (2009).
  * Heck  _et al._ (2018) R. Heck, O. Vuculescu, J. J. Sørensen, J. Zoller, M. G. Andreasen, M. G. Bason, P. Ejlertsen, O. Elíasson, P. Haikka, J. S. Laustsen, _et al._ , Proceedings of the National Academy of Sciences 115, E11231 (2018).
  * Larocca  _et al._ (2020) M. Larocca, E. Calzetta, and D. A. Wisniacki, Phys. Rev. A 101, 023410 (2020).
  * De Fouquieres and Schirmer (2013) P. De Fouquieres and S. G. Schirmer, Infin. Dimens. Anal. Quantum Probab. Relat. Top. 16, 1350021 (2013).
  * Riviello  _et al._ (2015) G. Riviello, K. M. Tibbetts, C. Brif, R. Long, R.-B. Wu, T.-S. Ho, and H. Rabitz, Phys. Rev. A 91, 043401 (2015).
  * Rach  _et al._ (2015) N. Rach, M. M. Müller, T. Calarco, and S. Montangero, Phys. Rev. A 92, 062343 (2015).
  * Khaneja  _et al._ (2005) N. Khaneja, T. Reiss, C. Kehlet, T. Schulte-Herbrüggen, and S. J. Glaser, J. Magn. Reson. 172, 296 (2005).
  * Motzoi  _et al._ (2011) F. Motzoi, J. M. Gambetta, S. Merkel, and F. Wilhelm, Phys. Rev. A 84, 022307 (2011).
  * De Fouquieres  _et al._ (2011) P. De Fouquieres, S. Schirmer, S. Glaser, and I. Kuprov, J. Magn. Reson. 212, 412 (2011).
  * Dalgaard  _et al._ (2020b) M. Dalgaard, F. Motzoi, J. H. M. Jensen, and J. Sherson, Phys. Rev. A 102, 042612 (2020b).
  * Blanes  _et al._ (2009) S. Blanes, F. Casas, J.-A. Oteo, and J. Ros, Phys. Rep. 470, 151 (2009).
  * Rasmussen (2003) C. E. Rasmussen, in  _Summer School on Machine Learning_ (Springer, 2003) pp. 63–71.
  * Perez-Garcia  _et al._ (2007) D. Perez-Garcia, F. Verstraete, M. M. Wolf, and J. I. Cirac, Quantum Info. Comput. 7, 401–430 (2007).
  * De Chiara  _et al._ (2008) G. De Chiara, M. Rizzi, D. Rossini, and S. Montangero, J. Comput. Theor. Nanosci. 5, 1277 (2008).
  * Orús (2019) R. Orús, Nat Rev Phys 1, 538 (2019).
  * Jensen  _et al._ (2021) J. H. M. Jensen, M. Gajdacz, S. Z. Ahmed, J. H. Czarkowski, C. Weidner, J. Rafner, J. J. Sørensen, K. Mølmer, and J. F. Sherson, Phys. Rev. Research 3, 013057 (2021).
  * Parker and Rardin (2014) R. G. Parker and R. L. Rardin,  _Discrete optimization_ (Elsevier, 2014).
  * Bernstein and Vazirani (1997) E. Bernstein and U. Vazirani, SIAM J. Sci. Comput. 26, 1411 (1997).
  * Nelder and Mead (1965) J. A. Nelder and R. Mead, Comput. J. 7, 308 (1965).
  * Virtanen  _et al._ (2020) P. Virtanen, R. Gommers, T. E. Oliphant, _et al._ , Nat. Methods 17, 261 (2020).
  * Pedregosa  _et al._ (2011) F. Pedregosa, G. Varoquaux, A. Gramfort, V. Michel, B. Thirion, O. Grisel, M. Blondel, P. Prettenhofer, R. Weiss, V. Dubourg, J. Vanderplas, A. Passos, D. Cournapeau, M. Brucher, M. Perrot, and E. Duchesnay, J. Mach. Learn. Res. 12, 2825 (2011).
  * Abadi  _et al._ (2016) M. Abadi, A. Agarwal, P. Barham, E. Brevdo, Z. Chen, C. Citro, G. S. Corrado, A. Davis, J. Dean, M. Devin, _et al._ , arXiv:1603.04467 (2016).


# Parallelized JUNO simulation software based on SNiPER

Tao Lin1 Jiaheng Zou1 Weidong Li1 Ziyan Deng1 Guofu Cao1 Xingtao Huang2 and
Zhengyun You3 (On Behalf of the JUNO Collaboration)  1Institute of High Energy
Physics, Chinese Academy of Sciences, Beijing, China  2Shandong University,
Jinan, China  3Sun Yat-sen University, Guangzhou, China<EMAIL_ADDRESS><EMAIL_ADDRESS>
###### Abstract

The Jiangmen Underground Neutrino Observatory (JUNO) is a neutrino experiment
to determine neutrino mass hierarchy. It has a central detector used for
neutrino detection, which consists of a spherical acrylic vessel containing 20
kt liquid scintillator (LS) and about 18,000 20-inch photomultiplier tubes
(PMTs) to collect light from LS. Around the central detector, there is a water
pool to shield radioactivities. The outer water pool is also equipped with
about 2000 PMTs to measure cosmic ray muons by detecting Cherenkov light.

As one of the important parts in JUNO offline software, the serial simulation
framework is developed based on SNiPER. It is in charge of physics generator,
detector simulation, event mixing and digitization. However Geant4 based
detector simulation of such a large detector is time-consuming and
challenging. It is necessary to take full advantages of parallel computing to
speedup simulation. Starting from version 10.0, Geant4 supports event-level
parallelism. Even though based on pthread, it could be extended with other
libraries such as Intel TBB and MPI. Therefore it is possible to parallelize
JUNO simulation framework via integrating Geant4 and SNiPER.

In this paper, our progress in developing parallelized simulation software are
presented. The SNiPER framework can run in sequential mode, Intel TBB mode or
other modes. The SNiPER task component is in charge of event loop, which is
like a simplified application manager. Two types of tasks are introduced in
the simulation framework, one is global task and another is worker task. The
global task will run only once to initialize detector geometry and physics
processes before any other tasks spawned. Later it is accessed by other tasks
passively. The worker tasks will be spawned after global task is done. In each
worker task, a Geant4 run manager is invoked to do the real simulation.
Therefore the simulation framework and the underlying TBB have been decoupled.
Finally, the software performance of parallelized JUNO simulation software is
also presented.

##  1 Introduction

JUNO [1, 2] is an underground neutrino experiment to measure neutrino mass
hierarchy and oscillation parameters. It will be located in southern China,
about 53 km away from two nuclear power plants. Figure 1 shows the schematic
view of the JUNO detector. To detect neutrinos, the innermost part, which is
called central detector, is filled with 20 kt liquid scintillator (LS). When
neutrinos react with LS, lights are yielded and then collected by the
surrounding photomultiplier tubes (PMTs). To suppress the radioactivities from
PMTs and rocks, water is filled around the spherical acrylic vessel. The outer
water pool is also equipped with PMTs to veto cosmic ray muon events by
detecting Cherenkov light. On the top of the water pool, a top tracker is used
to measure muons.

Figure 1: Schematic view of JUNO detector

In the JUNO experiment, offline software is an important part. Simulation
software is used for detector design and optimization, algorithms tuning and
physics studies. PMT waveform reconstruction is used to calibrate charge and
time for each PMT. Reconstruction algorithms are used to reconstruct events
and get physics objects for physics analysis. These algorithms are built upon
a software framework called SNiPER [3], and compose a full chain of data
processing [4], as shown in Figure 2.

Figure 2: Offline data processing chain

##  2 Towards parallelized detector simulation framework

A serial detector simulation framework [5] is implemented to integrate SNiPER
and Geant4 [6, 7]. Figure 3 shows the core class diagram. It is lightweight
with only a few classes in the design. Simulation software can be easily
migrated from standalone application into simulation framework. However,
simulation of a large detector is not only time consuming but also memory and
I/Os. After construction of the full geometry, simulation takes about 700 MB
in physical memory. More memory can be occupied during event loop. During
simulating of cosmic ray muons, millions of optical photons are produced and
propagated. It’s a challenge to simulate so many optical photons. A
parallelized simulation software is necessary for Monte-Carlo production in
the future. Both underlying framework and simulation libraries are required to
be thread safe and able to support parallel computing.

Figure 3: Design of detector simulation framework

In the SNiPER framework, a parallelized solution based on Intel Threading
Building Blocks (TBB) is proposed. Intel TBB supports task-based programming,
which simplifies the thinking of parallel computing. In SNiPER, Task is an
important component, which is like a lightweight application manager. In
sequential mode, a Task is configured by user and then invokes registered
algorithms. To support parallel computing, SNiPER Muster (Multiple SNiPER Task
Scheduler) is designed and implemented to integrate SNiPER Task and Intel TBB.
Instead of creating a single instance of Task, multiple instances of Task are
created. However, SNiPER Muster does not execute these Tasks directly. It will
create corresponding TBB-based workers to execute the Task. For better
performance, each worker is created only once and then binds to one thread.
Each worker pulls a Task and execute one event. If there are several workers,
Tasks are dispatched to different workers dynamically. As shown in Figure 4,
an I/O Task and three regular Tasks are running while a regular Task is
inactive. If the I/O Task is done, the inactive Task will continue to run.

Figure 4: SNiPER Muster

Starting from version 10.0, Geant4 supports multi-threaded simulation [8].
Events are simulated in different threads. With the evolution of Geant4, it is
possible to run simulation with Intel TBB, MPI and so on [9]. A parallelized
simulation framework is developed to integrate the latest Geant4 and SNiPER
Muster, which takes full advantages of parallel computing. Based on SNiPER
Muster, the simulation framework is decoupled with underlying Intel TBB, and
becomes more flexible and simple.

##  3 SNiPER Muster based simulation framework

The parallelized simulation framework is developed based on SNiPER Muster.
Several new classes are introduced, compared with the serial simulation
framework. The event loop is still controlled by SNiPER instead of Geant4,
however, the control flow is a bit different from serial simulation framework.
As shown in Figure 5, the simulation framework consists of a global task and
several worker tasks. The global task is in charge of initialization of
detector geometry and physics list. A customized master run manager is derived
from Geant4’s G4MTRunManager to take control of simulation. A corresponding
service is created, so that SNiPER can initialize run manager automatically.

Figure 5: Design of parallelized simulation framework

After the initialization of global task, SNiPER Muster starts the worker
tasks. Each worker task manages a simulation algorithm and a slave run manager
service. The simulation algorithm can invoke the slave run manager via the
corresponding service. The slave run manager is derived from Geant4’s
G4WorkerRunManager, which breaks the original event loop and simulates one
event each call. Instead of executing fixed number of events in each worker
task, the execution of worker task is controlled by SNiPER Muster. If there
are events to be simulated, the worker task will invoke the simulation
algorithm. Otherwise, the worker task will be stopped. Hence, this event
dispatching procedure can make full usage of CPUs.

##  4 Performance measurements

In the performance measurements, 1000 events of single $\gamma$s with energy
of 2.2 MeV are generated at detector center. All are done at a blade server
with Intel Xeon CPU E5-2680 v3 @ 2.5 GHz and 64 GB of memory. The operating
system is Scientific Linux 6.5 with GCC 4.9.4. The Geant4 is 10.03.p01, built
with global-dynamic TLS (Thread Local Storage) model. The physics list
G4EmStandardPhysics and G4OpticalPhysics are used in simulation. The
simulation time are measured three times for each case. In order to eliminate
I/O interference, the events are not saved into ROOT files.

Figure 6 shows speedup ratio versus number of threads. The dotted line
represents the ideal result. The gray line with full square represents the
result without any optimization. Up to four threads, the speedup is good.
However, when increasing the number of threads to 24, the speedup ratio is
stuck around 8.

Figure 6: Speedup versus number of threads

Intel VTune Amplifier profiler is used to analyze this performance issue. We
find that the hotspot comes from a mutex used in Geant4’s
G4MaterialPropertiesTable. This mutex is used when each optical photon
calculates its group velocity. The problem is that there are a lot of optical
photons in JUNO simulation, so this mutex is accessed frequently in order to
get the velocity. When there are multiple threads accessing the mutex, the
performance issue happens. The optimization is moving this part into
initialization stage. When refractive index is set, the group velocity is
precalculated during the initialization. As shown in Figure 6, the blue line
with full circle represents the result with optimization. After this
optimization, a better performance is achieved.

##  5 Conclusions

In this proceeding, we show how JUNO simulation software is parallelized based
on SNiPER Muster. By introducing several new classes, Geant4 10 is integrated
into SNiPER Muster. The global task initializes detector geometry and physics
processes while the worker task is invoked by SNiPER Muster to simulate one
event. Events are dispatched dynamically to different workers, which improves
CPU utilization. The software performance studies show that there are
performance issues due to a mutex in material properties table. After
optimization, simulation software achieves a linear speedup, which fulfills
requirements.

To speedup events such as cosmic ray muons, we are investigating how to use
track-level parallelism in simulation framework. The idea is using MPI to
dispatch tracks from master node to different worker nodes. In the future, MPI
and TBB can be applied simultaneously in different levels of parallelized
simulation software in order to speedup simulation.

## Acknowledgments

This work is supported by Joint Large-Scale Scientific Facility Funds of the
NSFC and CAS (U1532258), the Strategic Priority Research Program of the
Chinese Academy of Sciences, Grant No. XDA10010900, National Natural Science
Foundation of China (11575224, 11405279, 11675275).

## References

## References

  * [1] An F et al. (JUNO) 2016 J. Phys. G43 030401 (Preprint 1507.05613)
  * [2] Djurcic Z et al. (JUNO) 2015 (Preprint 1508.07166)
  * [3] Zou J H, Huang X T, Li W D, Lin T, Li T, Zhang K, Deng Z Y and Cao G F 2015 J. Phys. Conf. Ser. 664 072053
  * [4] Huang X T, Li T, Zou J H, Lin T, Li W D, Deng Z Y and Cao G F 2016 PoS ICHEP2016 1051 URL https://pos.sissa.it/282/1051
  * [5] Lin T, Zou J, Li W, Deng Z, Fang X, Cao G, Huang X and You Z (JUNO) 2017 22nd International Conference on Computing in High Energy and Nuclear Physics (CHEP 2016) San Francisco, CA, October 14-16, 2016 (Preprint 1702.05275) URL http://inspirehep.net/record/1514060/files/arXiv:1702.05275.pdf
  * [6] Agostinelli S et al. (GEANT4) 2003 Nucl. Instrum. Meth. A506 250–303
  * [7] Allison J et al. 2006 IEEE Trans. Nucl. Sci. 53 270
  * [8] Allison J et al. 2016 Nucl. Instrum. Meth. A835 186–225
  * [9] Dotti A, Asai M, Barrand G, Hrivnacova I and Murakami K 2016 Proceedings, 2015 IEEE Nuclear Science Symposium and Medical Imaging Conference (NSS/MIC 2015): San Diego, California, United States p 7581867 (Preprint 1605.01792) URL http://inspirehep.net/record/1456144/files/arXiv:1605.01792.pdf