11institutetext: Department of Computer Engineering, Boğaziçi University,
34342, Istanbul
11email: {okan.asik, binnur.gorer<EMAIL_ADDRESS>
# End-to-End Deep Imitation Learning: Robot Soccer Case Study

Okan Aşık  Binnur Görer  H. Levent Akın

###### Abstract

In imitation learning, behavior learning is generally done using the features
extracted from the demonstration data. Recent deep learning algorithms enable
the development of machine learning methods that can get high dimensional data
as an input. In this work, we use imitation learning to teach the robot to
dribble the ball to the goal. We use B-Human robot software to collect
demonstration data and a deep convolutional network to represent the policies.
We use top and bottom camera images of the robot as input and speed commands
as outputs. The CNN policy learns the mapping between the series of images and
speed commands. In 3D realistic robotics simulator experiments, we show that
the robot is able to learn to search the ball and dribble the ball, but it
struggles to align to the goal. The best-proposed policy model learns to score
4 goals out of 20 test episodes.

##  1 Introduction

In robot learning, it is a challenging problem to collect training data set in
real-world environment. One possible solution is to isolate the problem from
the entire robotic system. For example, object detection, which is a highly
important problem in robotics vision, can be addressed separately in
development of a home service robot. Image data can be captured from the robot
camera without requiring fully functional robotic system and can be used with
supervised or unsupervised learning methods. In contrast to supervised and
unsupervised learning, reinforcement learning is inherently hard for real-
world robotics as it requires direct interaction with the environment. Let’s
consider the problem of high-level behavior (task) learning such as learning
to play soccer. Reinforcement learning is based on the exploration of the
state-action space and the iterative improvement of the values of experienced
state-action pairs. That exploration takes many time steps and may harm the
robot during the interaction with the environment. Therefore, it is not
feasible to use direct reinforcement learning for behavior learning in
robotics.

One of the most common approaches to solve behavior learning is the imitation
learning. The imitation learning can be classified as one of the supervised
learning approaches since the learning model imitates data provided by
demonstration. A demonstrator controls the robot to complete a specific
behavior(task) using the robot. During the demonstration, data is collected
from the sensors and actuators of the robot. By using many demonstration
sessions, a data set is created that can be used to learn sensor-actuator
tuples. In this study, we use imitation learning to dribble the ball to score
a goal using a humanoid robot.

These behavior learning approaches require a symbolic representation of the
behavior. This representation determines the abstraction of the problem. For
example, in the case of the ball dribbling behavior, we can represent the
world of the robot as the pose of the robot and the position of the ball on
the field. To be able to construct a learning problem and apply it in real-
world, we need to have robot perception algorithm that is able to calculate
the pose of the robot and the position of the ball. Most of the behavior
learning approaches in the literature uses this level of abstraction in
learning [13, 10, 3, 11]. However, this level of abstraction requires solving
complex perception problems such as determining its own pose using the
features extracted from the camera image. In this study, we learn the behavior
using the camera image as an input to our machine learning model. Therefore,
the basic robotics problems such as perception, localization, and planning
would be inherently learned by our model.

We create a training data set that is collected while a teacher realizes the
behavior on the robot. The dataset consists of images taken from the robot
camera and movement commands. Instead of recruiting a human teacher to make
the task to the robot a number of times, we use a robotics software developed
for the robot soccer as a teacher. The simulation environment makes it easy to
code the behavior since we have the pose of the robot and the ball without
running any perception or localization algorithm. We train a deep learning
model using the training data set. The model that can replicate the same speed
commands given the images achieves the imitation learning without requiring
any complex solutions to subproblems in the task.

The robot soccer problem is a partially observable problem because the robot
cannot sense the current state of the environment using its own onboard
sensors. Therefore, a policy that maps the current perception to an action
command may have suboptimal behavior. Therefore, we propose to use recurrent
convolutional neural network since it has the capability of estimating the
current state from a series of perceptions. In addition to the recurrent
convolutional neural network, we also propose to use a hierarchical model to
divide the ball dribbling behavior into sub-skills such as searching for the
ball and aligning to the goal and a more basic convolutional neural network
approach to compare the performances.

We measure the learning performance of our models on the test demonstration
data set. On the test data set, we achieved an average error of $0.16$,
$0.17$, and $0.14$ with the convolutional neural network, hierarchical
convolutional neural network, and recurrent convolutional neural network
policies, respectively. The error is calculated as the average Euclidean
distance between two action vectors. Also, we count the number of goals scored
in 20 test episodes; CNN scores 3 goals, H-CNN scores 3 goals, and R-CNN
scores 4 goals. In this study, we show that we can teach a humanoid robot to
carry out a complex task that requires perception and planning using deep
imitation learning. We test our method on a realistic robotics simulator, but
our future aim is to evaluate the approach using real robots.

##  2 Related Work

Imitation learning is a promising research area since it has the potential of
enabling non-technical domain experts to teach their expertise to a
computational machine. There are many ways of demonstrating a behavior to a
robot such as teleoperation and exoskeleton suit. In this study, a high-level
behavior control software controls the robot in simulation, and raw sensor
readings and mid-level actions are collected as demonstration data set. Our
mid-level actions are 2D speed commands. We use teleoperation based data
demonstration and direct state-action policy representation according to the
categorization of Husseion et al., and Argall et al. since we learn state-
action mapping [7, 2].

This study is based on end to end deep learning. The first study shows that an
agent can learn to play Atari using the screen inputs as the state of the game
[14]. Although end to end learning is a reinforcement learning method, its key
feature is that the policy representation is powerful enough to learn without
abstracting the state or action of the problem. For example, in Atari
learning, state is the screen of the game and action is joystick command so
that no feature engineering is involved. Our study also uses sensor inputs
(camera images) without any feature engineering, but the actions have an
abstraction since humanoid walking problem is still quite complex to be used
in such end to end learning. Also, most of the end to end learning approaches
use deep convolutional neural network to represent the expected reward of
state-action pairs. However, we represent our policy as a simple mapping
function that maps states to actions. We could also use reinforcement
learning, but it would take quite long time to learn to score a goal. We know
that if the reward is deep in reinforcement learning problems, it gets harder
to learn (consider the Atari game Montezuma’s Revenge [14])

We structure our related work in two subsections; end to end deep learning and
ball dribbling behavior. In the first section, we overview the methods using
end to end deep learning. In the second section, we overview the methods that
directly issue the ball dribbling behavior.

###  2.1 End to End Learning

The first study that combines deep learning and reinforcement learning is done
by Mnih et al. [14]. They achieved human-level atari game playing by using a
deep convolutional neural network to represent the policy that evaluates the
expected reward of state-action pairs. Their most important contribution is
using the state representation as the actual screen of the game as if a human
player perceives the game. This approach is also the most important part of
our work where relevant features are automatically discovered by the
convolutional neural network instead of human-engineered features.

Guo et al. use Monte Carlo Planning algorithm to collect data that will be
used to train the convolutional neural network [4]. Although Monte Carlo
Planning has better performance than reinforcement learning algorithms, it is
quite slow compared to neural networks. They show that reinforcement learning
that uses the data generated by Monte Carlo Planning improves the performance.
This is also one of the inspiration for this study where an expert software is
used to teach an end-to-end neural network policy. However, we do not
represent policy as a Q-function [18], but as a state-action mapping function.
Another work that combines the Monte Carlo planning and deep learning is the
amazing Alpha Go that beats the Go champion [16].

One of the most important robotics application of the end-to-end deep learning
method is the study that learns to manipulate objects using the raw camera
images [12]. The robot learns to send joint commands to motors using a
convolutional neural network. In contrast to them, we generate high-level
motion commands to control the 2D motion of the robot using 2D speed commands.
They improve the performance of their method by using many different
preprocessing and pre-learning methods. In contrast to them, we train the
whole network at once using the data generated by the demonstrator.

Husseion et al. combine the deep imitation learning with active learning to
learn to navigate in the 3D maze [6]. The use of imitation learning in 3D
simulation environment to solve such complex tasks are also our aim. In
contrast to their work, we use realistic 3D robotics simulator.

###  2.2 Ball Dribbling Task

One of the early studies that aim to learn to dribble the ball in robot soccer
is done by Latzke et al. [10]. They use imitation learning to improve the
performance of reinforcement learning algorithm. The data generated by the
teacher is used to initialize the reinforcement learning policy. They train a
humanoid robot to dribble the ball to the empty goal. They report a reduction
in training time and increase in the learning performance when imitation data
and function approximation is used. Leottau et al. propose two layer approach
for humanoid robot ball dribbling [11]. They use a fuzzy logic controller for
alignment to the ball and reinforcement learning to push the ball towards the
goal. Although pushing the ball towards the goal after the alignment seems a
simple task, in addition to the dribbling the ball as fast as possible, they
try to keep the ball possession, i.e. keep the ball as close as possible.
Therefore, they use reinforcement learning to optimize for two conflicting
goals; being as fast as possible, but also keeping the ball as close as
possible. Mericli et al. propose to use a corrective human feedback system to
teach the robot to dribble the ball through stationary defender robots [13].
The main contribution of the work is the combination of hand-coded behavior
with the active demonstration of the human. The performance of the ball
dribbling, that is the time to scoring a goal, is improved by the integration
of the demonstrations. All these ball dribbling tasks have lower ambitions
compared to our work. They use the model of the problem, robot perception, and
localization and solve one of the subtasks of the behavior, that is ball
dribbling.

##  3 Methods

This study consists of two parts; dataset creation and the training of deep
neural networks. In imitation learning, the dataset is created by an expert
who knows how to control the robot to carry out a task. In this study, we use
a hand-coded behavior on B-Human simulation environment instead of a human
expert. We define a behavior that searches the ball, goes towards the ball,
aligns to the goal, and dribbles the ball to the goal. The collected data set,
a set of image-speed tuples, is used to train the convolutional neural
networks with three different architectures. These neural networks learn to
predict robot speed commands based on the images taken from the cameras.

###  3.1 Dataset Creation

RoboCup Standard Platform League (SPL) is a competition where robot soccer
teams compete using 5 Nao robots. B-Human is one of the successful teams in
SPL and they make their robot software publicly
available111https://github.com/bhuman/BHumanCodeRelease. They provide all
their software with a 3D robotics simulation software called SimRobot with the
modules that can both run on the robot and in the simulator. A screenshot of
the simulator can be seen in Figure 1(a). In this study, we use the SimRobot
simulation environment and robot software. [15].

(a) A screenshot of SimRobot simulator.

(b) The position of robots (red points) and balls (blue points) of the dataset

Figure 1: (left) SimRobot environment and (right) robot and ball positions

The Nao robot that is used in SPL has two cameras, one views the forward (top
camera) and another one (bottom camera) is aligned to view the feet of the
robot. These cameras do not (almost) overlap. Although we can get the images
from those cameras at different resolutions, most of the teams use 640x480 for
the top camera and 320x240 for the bottom camera. Those cameras provide images
at 33 Hz.

We use the images of these two cameras as input our learning approach. We pre-
process images to lower the dimensionality of images to able to use as input.
We scale and convert images to gray-scale. Scaling is done by averaging the
pixel values corresponding to a particular region by Thumbnail provider of
B-Human software module [15]. We scale top camera images to 160x120 resolution
and bottom camera images to 80x60. A series of sample images can be seen in
Figure 2.

Figure 2: A series of sample pre-processed images. The bigger images on the
left are top camera images, the smaller images on the right are bottom camera
images.

B-Human software components have two main categories; cognition and motion.
The motion commands calculated by the cognition modules are carried out by the
motion modules. The distinction between these two processes is due to the
different operating frequencies of motors and sensors. In this study, due to
the practicality, we developed a simple behavior in B-Human software and
created our data set from this behavior. However, we can safely assume that we
can replace the behavior software with a human operator and collect similar
data. Our behavior has the following components; searching for the ball, going
towards the ball, turning around the ball to align to the closest goal and
dribbling the ball towards the closest goal. By default, the software
calculates a speed command at every time step when a new camera image is
taken. Since we aim to learn a policy that maps image to speed commands, we
create a motion command when images from both top and bottom cameras are
taken. Therefore, at the current time step, if we have not received images
from both cameras, we repeat the last command. Otherwise, if the images
captured from two cameras processed separately, the predicted actions for each
of them may conflict.

We create our dataset by selecting random robot position in 4 meters by 4
meters area on the center of the field, and random ball position in 4 meters
by 2 meters area on either side of the field. We choose the position of the
ball such that the ball is at least one meters away from the center line of
the field. In this way, we aim to choose a side for the ball where the robot
will score a goal. The random positions of the robots and balls can be seen in
Figure 1(b). Data collection starts when the robot starts to move and ends
when the robot scores a goal. We created our dataset using 100 episodes with
the random ball and robot positions. Our dataset has 160x120 and 80x60
grayscale images with speed commands (speed commands consist of forward, left
and turn speeds).

###  3.2 Deep Imitation Learning

There are two important research questions of the imitation learning; handling
the states that are not part of the learning and the generalizing capability
of the model and representation. A teacher demonstrates the expected behavior
by controlling the robot at every time step. That way, we construct state-
action tuples. By using this demonstration data, we learn the policy of the
task. If the teacher’s demonstrations cover very limited part of the state
space or the learning model is not able to generalize well, the performance of
the imitation learning decreases.

In this study, we use convolutional neural networks as our learning models. In
this way, we are able to directly use camera images as input to our learning
model. This is called end-to-end learning since we do not use intermediate
feature extraction approach. In other words, the model itself learns to
extract useful features based on the task. This approach is closer to the
human’s cognition where real neural network processes the image and produces
an action.

When the state is fully observable by the robot, the direct state-action
mapping has the capability of representing the optimal policy taught by the
teacher assuming that the dataset encompasses every state. However, when the
problem is partially observable, where a single observation is not enough to
infer the world state, our learning model needs to learn to infer action from
a series observations. For example, in the robot soccer problem, the robot
perceives the world from its cameras. Therefore, the robot can perceive only a
small part of the field. If there is a ball in the camera, the robot knows the
position of the ball, but if the ball is not in the field of view of the
camera, the robot does not know the current state of the world about the ball.
A convolutional neural network (CNN) can only map the current image to an
action so that previous actions or images does not affect the choice of the
action. However, the robot soccer behavior needs the information of previous
camera images to act properly. Therefore, we extend basic CNN with a
hierarhical neural network (H-CNN) and recurrent convolutional neural network
(Recurrent-CNN).

####  3.2.1 Convolutional Neural Network (CNN)

This is the straightforward neural network to learn state-action mapping. We
assume that the images taken from two cameras of the robot are powerful enough
to learn the behavior operated by the teacher. To be able to process the two
images having different resolutions, we have two different channels of input
layers. The first hidden layer consists of 16 units for top camera image and 8
for bottom camera image 5 by 5 convolutional filters that convolve the images
with stride 2. After the convolution, a 2 by 2 maximum pooling filter is used.
The second hidden layer consists of 32 for top camera image and 16 for bottom
camera image 5 by 5 convolutional filters that convolves the images with
stride 2. After the convolution, a 2 by 2 maximum pooling filter is used. We
also used ReLU [9] activation function after every convolutional layer. After
this four layers, we have a fully connected layer with 128 units and ReLU
activation. The last layer is a fully connected linear output layer with three
units. The layers of the neural network model can be seen in Figure 3.

Figure 3: The Architecture of Deep Convolutional Neural Network

####  3.2.2 Hierarchical Convolutional Neural Network (H-CNN)

We propose a hierarchical approach to learn a behavior that consists of sub-
behaviors or skills. For example, our robot soccer behavior consists of four
sub-behaviors such as searching the ball, going to the ball, aligning to the
goal and dribbling the ball. Our assumption is that our demonstration data
including the actions can be clustered into four classes where each may
correspond to a sub-behavior. A set of example images from different clusters
can be seen in Figure 4. Therefore, we first cluster the whole data into a
predetermined number of classes. Then, we learn a separate convolutional
neural network with the same architecture as Figure 3 for each cluster. During
the test phase, we find the closest cluster to the given images and use the
learned CNN of that cluster.

Figure 4: Sample representative images from different clusters.

For the clustering, we first create histograms of images with 5 bins. We
concatenate the top and bottom histogram features. Then, we use the K-Means
clustering algorithm [1] to cluster the whole dataset. After the clustering,
we obtain a different demonstration dataset for each cluster. Finally, we
train a CNN for each of these datasets.

####  3.2.3 Recurrent Convolutional Neural Network (Recurrent-CNN)

We also propose a recurrent convolutional neural network model (Recurrent-CNN)
to learn the mapping between a series of states and a series of actions. We
augment the CNN model proposed in Section 3.2.1 by adding a Long-Short Term
Memory (LSTM) layer after the convolutional layers as seen in Figure 5. LSTM
cells have the capability to learn to keep which part of the data in its
memory and which part to forget [5]. It is shown that LSTM is the state of the
art method on sequence learning [17].

Figure 5: The Recurrent-CNN with LSTM Layer

The LSTM layer of Recurrent-CNN has 64 units. It is trained using a 100 time
step window of data that corresponds to approximately three seconds of the
data sequence. We empirically determine the window size in order to balance
the training time and the performance. We test and evaluate the model by
sliding a 100 time step window over the data sequence.

##  4 Experiments and Results

Our training data set consists of 100 episodes where every episode includes a
complete robot soccer behavior as explained in Section 3.1. We train all of
our three models using ADAM neural network optimization algorithm [8]. We used
different maximum iterations for different models since we used fixed
processing time. Finally, we measure the overall training error as the root
mean square of the Euclidean distance between the true action vector and
predicted vector. We also measure the performance of the models on 20 test
episodes that is not used in training.

###  4.1 The Learning Performance

One of the most important factors that determine the success of imitation
learning methods is the learning performance. This is based on how well the
model can learn the demonstration data. Also, the model should not memorize
all the demonstrations in order to avoid over-fitting. A robust imitation
learning method is supposed to learn which action to take for the states which
are not in the demonstration. The model should be able to generalize for new
states as well. The learning performance of an imitation learning model can be
measured by the average error on the demonstration data. The generalization
performance can be measured by the average error on the test demonstration
data.

Table 1: The Performance of CNN Models Models | # of Training Iterations | Training Error | Test Error | Goals Scored
---|---|---|---|---
CNN | $10^{6}$ | 0.12 | 0.16 | 3
H-CNN | $15\times 10^{4}$ | 0.15 | 0.17 | 3
Recurrent-CNN | $3\times 10^{3}$ | 0.036 | 0.14 | 4

The overall performance of different models can be seen in Table 1. The R-CNN
has the best performance by having the least training error and least test
error. Test and training errors are reported as the average distance of the
calculated speed vector and the demonstration speed vector. The training
results are provided over all training data and the test results are provided
on the data from 20 new episodes. When we observe the overall learned
behavior, we see that robot struggles to find the direction of the goal and
turns around while dribbling the ball. However, the robot robustly performs
searching for the ball, moving towards the ball and dribbling the ball
actions. When we measure the number of goals scored in test scenarios, the
robot scores 4 goals using R-CNN, 3 goals using CNN and 3 goals using H-CNN
models over 20 test episodes. The trajectories of the robot with R-CNN model
and demonstration trajectories are seen in Figure 6. The left trajectories
present the goal scored ones, the right trajectories present some of the
failed ones.

Figure 6: The comparison of robot trajectories while running test episodes.
The left figures show goal scored episodes. The right figures show some of the
failed test episodes.

###  4.2 Discussion

There are two important research questions about imitation learning; how to
represent and collect demonstration data, and how to represent and optimize
the policy. Our approach generates demonstrations from an expert software such
that the behavior of the robot is consistent. By using a software for the
collection of the demonstration data, we minimize the possibility of
demonstrator errors and instability.

Since demonstrator errors are minimized, the success of our approach basically
depends on the data representation and model selection for the policy. We use
deep convolutional neural networks (CNN) to represent our policies because
they enable us to represent the demonstration data as camera images. It is
known that CNNs perform well on visual classification and detection tasks [9].
The dataset we generated can be viewed as a machine learning dataset for
visual regression task. However, basic use of CNN comes with the assumption
that we can infer the expected action using camera images of the current time
step. Although the problem is not fully observable, we see that the
performance of different policy representations (CNN, H-CNN, Recurrent-CNN)
are close. This might be due to the possible overlap between training data set
and test data set. We argue that learned models are able to memorize some of
the training cases to reduce the error, but not able to generalize or memorize
the whole task.

When we investigate the behaviors generated with different learning models, we
see that all of the models are able to carry out three basic subtasks;
searching the ball, going towards the ball, and dribbling the ball. However,
they generally fail to align with the goal. This is due to the demonstration
behavior. The demonstration behavior decides to align to the goal when the
ball is close to dribbling and orientation of the robot is not towards the
goal. However, the behavior does not use the robot images to carry out
aligning to the goal behavior. It uses the true position of the robot provided
by the simulator. Hence, how to align to the goal using only the camera images
is not explicitly exhibited in the demonstration data.

Although the real world application of this method is limited, the policy
learned in the simulation can be further improved with less amount of real-
world data. Also, we keep the CNN architectures as small as possible to be
able to have real-time performance. On a laptop having 2.8 GHz cpu222Intel
7700HQ 4 cores 2.8 GHz, 32GB 2400MHz RAM, the forward pass of CNN architecture
takes 0.0018 seconds, and Recurrent-CNN architecture takes 0.0021 seconds on
average. Based on these statistics, we expect near-real-time performance on
Nao robot.

##  5 Conclusion

Using the deep imitation learning method, we learned a basic robot soccer
behavior of searching the ball, moving towards the ball, and dribbling the
ball to the goal. Our proposed convolutional neural network uses two images
from two different cameras of the robot as input. We created our dataset using
RoboCup SPL team B-Human software modules and carried out our experiments
using SimRobot 3D realistic robot simulator. We show that using quite a few
samples, we can learn simple robot soccer behaviors using end-to-end training.
In the future, we aim to scale our method to more complex behaviors on real
robots in real environments.

Acknowledgments

This project is supported by Turkey Technology Team Foundation (T3).

## References

  * [1] Alpaydin, E.: Introduction to machine learning. MIT press (2014)
  * [2] Argall, B.D., Chernova, S., Veloso, M., Browning, B.: A survey of robot learning from demonstration. Robotics and autonomous systems 57(5), 469–483 (2009)
  * [3] Aşık, O., Akın, H.L.: Solving multi-agent decision problems modeled as dec-pomdp: A robot soccer case study. In: Chen, X., Stone, P., Sucar, L.E., van der Zant, T. (eds.) RoboCup 2012: Robot Soccer World Cup XVI. pp. 130–140. Springer Berlin Heidelberg, Berlin, Heidelberg (2013)
  * [4] Guo, X., Singh, S., Lee, H., Lewis, R.L., Wang, X.: Deep learning for real-time atari game play using offline monte-carlo tree search planning. In: Advances in neural information processing systems. pp. 3338–3346 (2014)
  * [5] Hochreiter, S., Schmidhuber, J.: Long short-term memory. Neural computation 9(8), 1735–1780 (1997)
  * [6] Hussein, A., Elyan, E., Gaber, M.M., Jayne, C.: Deep imitation learning for 3d navigation tasks. Neural computing and applications pp. 1–16 (2017)
  * [7] Hussein, A., Gaber, M.M., Elyan, E., Jayne, C.: Imitation learning: A survey of learning methods. ACM Computing Surveys (CSUR) 50(2), 21 (2017)
  * [8] Kingma, D.P., Ba, J.L.: Adam: A Method for Stochastic Optimization pp. 132–142 (2014), https://arxiv.org/abs/1412.6980
  * [9] Krizhevsky, A., Sutskever, I., Hinton, G.E.: Imagenet classification with deep convolutional neural networks. In: Advances in neural information processing systems. pp. 1097–1105 (2012)
  * [10] Latzke, T., Behnke, S., Bennewitz, M.: Imitative reinforcement learning for soccer playing robots. In: Robot Soccer World Cup. pp. 47–58. Springer (2006)
  * [11] Leottau, L., Celemin, C., Ruiz-del Solar, J.: Ball dribbling for humanoid biped robots: a reinforcement learning and fuzzy control approach. In: Robot Soccer World Cup. pp. 549–561. Springer (2014)
  * [12] Levine, S., Finn, C., Darrell, T., Abbeel, P.: End-to-end training of deep visuomotor policies. The Journal of Machine Learning Research 17(1), 1334–1373 (2016)
  * [13] Meriçli, Ç., Veloso, M., Akın, H.L.: Task refinement for autonomous robots using complementary corrective human feedback. International Journal of Advanced Robotic Systems 8(2), 16 (2011)
  * [14] Mnih, V., Kavukcuoglu, K., Silver, D., Rusu, A.A., Veness, J., Bellemare, M.G., Graves, A., Riedmiller, M., Fidjeland, A.K., Ostrovski, G., et al.: Human-level control through deep reinforcement learning. Nature 518(7540), 529 (2015)
  * [15] Röfer, T., Laue, T., Bülter, Y., Krause, D., Kuball, J., Mühlenbrock, A., Poppinga, B., Prinzler, M., Post, L., Roehrig, E., Schröder, R., Thielke, F.: B-Human team report and code release 2017 (2017), only available online: http://www.b-human.de/downloads/publications/2017/coderelease2017.pdf
  * [16] Silver, D., Huang, A., Maddison, C.J., Guez, A., Sifre, L., Van Den Driessche, G., Schrittwieser, J., Antonoglou, I., Panneershelvam, V., Lanctot, M., et al.: Mastering the game of go with deep neural networks and tree search. nature 529(7587), 484–489 (2016)
  * [17] Sutskever, I., Vinyals, O., Le, Q.V.: Sequence to sequence learning with neural networks. In: Ghahramani, Z., Welling, M., Cortes, C., Lawrence, N.D., Weinberger, K.Q. (eds.) Advances in Neural Information Processing Systems 27, pp. 3104–3112. Curran Associates, Inc. (2014), http://papers.nips.cc/paper/5346-sequence-to-sequence-learning-with-neural-networks.pdf
  * [18] Sutton, R.S., Barto, A.G.: Introduction to reinforcement learning, vol. 135. MIT press Cambridge (1998)


# CP2K on the road to exascale

Thomas D. Kühne†,§,∗, Christian Plessl‡,§, Robert Schade‡,§, and Ole Schütt¶ †
Dynamics of Condensed Matter and Center for Sustainable Systems Design, Chair
of Theoretical Chemistry, Paderborn University, D-33098 Germany
‡Department of Computer Science, Paderborn University, D-33098 Germany
§ Paderborn Center for Parallel Computing, Paderborn University, D-33098
Germany ¶CP2K Foundation, CH-8006 Zurich Switzerland<EMAIL_ADDRESS>
###### Abstract

The CP2K program package, which can be considered as the swiss army knife of
atomistic simulations, is presented with a special emphasis on ab-initio
molecular dynamics using the second-generation Car-Parrinello method. After
outlining current and near-term development efforts with regards to massively
parallel low-scaling post-Hartree-Fock and eigenvalue solvers, novel
approaches on how we plan to take full advantage of future low-precision
hardware architectures are introduced. Our focus here is on combining our
submatrix method with the approximate computing paradigm to address the
immanent exascale era.

##  1 Background and Current Status

The open-source simulation package CP2K is an extensive quantum chemistry and
condensed matter physics program that comprises a large variety of different
theoretical methods and computational approaches to conduct most diverse
atomistic simulations for large-scale condensed-phase systems, such as
liquids, solids, nanomaterials and molecular structures, to name just a few.
All of this is made possible from the outset by the design of highly efficient
algorithms with a focus on excellent parallel scalability and suitability for
novel high-performance computing architectures, as demonstrated by its
Quickstep electronic structure module. The latter is based on the Gaussian and
plane wave (GPW) approach and its all-electron variant, the Gaussian-augmented
plane wave (GAPW) method, which allows for a particularly efficient treatment
of large-scale orbital-free and Kohn–Sham density functional theory (DFT)
calculations, as well as wavefunction-based correlation methods, such as
Hartree-Fock (HF), second-order Møller–Plesset perturbation theory (MP2),
random phase approximation (RPA), hybrid and double-hybrid DFT and the GW
approximation, all with arbitrary boundary conditions [1].

However, the great appeal of CP2K lies in the possibility to combine all
available total energy and force methods with a wide selection of sampling
techniques, such as Monte Carlo, Ehrenfest/real-time dynamics and most
importantly molecular dynamics, as well as advanced free-energy and rare-event
sampling schemes, to enable realistic simulations at finite-temperature beside
conventional static calculations. On the one hand, this necessitates the
general availability of analytic gradients in particular for periodic boundary
conditions, in order to permit the efficient calculation of nuclear forces. On
the other hand, minimum time to solution is essential to allow for an
extensive sampling via the techniques listed above. In that respect, a unique
selling point of CP2K is the implementation of the second-generation Car-
Parrinello method, which allows to routinely conduct nanosecond long DFT-based
ab-initio molecular dynamics (AIMD) simulations with thousands of atoms. The
superior efficiency of this approach originates from the design of an improved
coupled electron-ion dynamics that keeps the electrons very close to their
instantaneous ground-state using just one preconditioned gradient calculation
per AIMD step, which can thus be thought of as an electronic force to
propagate the electronic degrees of freedom in dimensionless time [2].

##  2 Development Priorities

Beside the implementation of sophisticated spectroscopic properties [3, 4],
which are either based on density functional perturbation theory or time-
dependent DFT within the Tamm-Dancoff approximation [1, 5], the current
development priorities are focused mainly on devising novel low-scaling post-
HF methods including the implementation of analytic nuclear gradients for
arbitrary boundary conditions. As already indicated above, particular
emphasize is on HF and MP2 methods, which are a prerequisite for simulations
based on RPA, hybrid and double-hybrid DFT schemes, as well as GW. The four-
center two-electron repulsion integrals (ERI), which in Mulliken notation
reads as

$(\mu\nu|\lambda\sigma)=\int d\mathbf{r}_{1}\int
d\mathbf{r}_{2}\,\phi_{\mu}^{*}(\mathbf{r}_{1})\phi_{\nu}(\mathbf{r}_{1})\,\frac{1}{|\mathbf{r}_{1}-\mathbf{r}_{2}|}\,\phi_{\lambda}^{*}(\mathbf{r}_{2})\phi_{\sigma}(\mathbf{r}_{2})$
(1)

are of central importance for all wavefunction-based post-HF methods. In
addition to well established integral screening techniques based on the
Schwarz inequality

$|(\mu\nu|\lambda\sigma)|\leq|(\mu\nu|\mu\nu)|^{1/2}\,|(\lambda\sigma|\lambda\sigma)|^{1/2}$
(2)

to reduce the scaling from $\mathcal{O}(N^{4})$ to $\mathcal{O}(N^{2})$, a
similar density matrix screening can also be employed to eventually sustain
linear scaling with respect to the system size $N$. In the latter, the largest
density matrix element
$P_{\max}=\max\\{|P_{\mu\lambda}|,|P_{\mu\sigma}|,|P_{\nu\lambda}|,|P_{\nu\sigma}|\\}$
is used to screen elements smaller than $\epsilon_{\text{Schwarz}}$ via

$P_{\max}\times|(\mu\nu|\mu\nu)|^{1/2}\,|(\lambda\sigma|\lambda\sigma)|^{1/2}\leq\epsilon_{\text{Schwarz}},$
(3)

where $P_{\max}$ is either the density matrix from the previous self-
consistent field iteration, or from a converged semi-local DFT calculation,
but ideally the propagated density matrix via second-generation Car-Parrinello
AIMD of the previous timestep [6]. Yet, at the core of all implemented post-HF
approaches are either the auxiliary density matrix method (ADMM) [6], or the
resolution of identity (RI) approach [7]. The key ingredient of the former is
the use of an auxiliary density matrix $\mathbf{\hat{P}}$, which approximates
the original density matrix $\mathbf{P}$, but is substantially easier to
compute due to being smaller in size, or more rapidly decaying than
$\mathbf{P}$. For the case of computing the Hartree-Fock exchange (HFX), the
exact energy $E_{\text{X}}^{\text{HFX}}[\mathbf{P}]$ is replaced by the
computationally superior expression
$E_{\text{X}}^{\text{HFX}}[\mathbf{\hat{P}}]$, whereas the difference between
the two terms is corrected approximately at the semi-local DFT level. Hence,

$\displaystyle E_{\text{X}}^{\text{HFX}}[\mathbf{P}]$ $\displaystyle=$
$\displaystyle
E_{\text{X}}^{\text{HFX}}[\mathbf{\hat{P}}]+\left(E_{\text{X}}^{\text{HFX}}[\mathbf{P}]-E_{\text{X}}^{\text{HFX}}[\mathbf{\hat{P}}]\right)$
(4) $\displaystyle\approx$ $\displaystyle
E_{\text{X}}^{\text{HFX}}[\mathbf{\hat{P}}]+\left(E_{\text{X}}^{\text{DFT}}[\mathbf{P}]-E_{\text{X}}^{\text{DFT}}[\mathbf{\hat{P}}]\right),$

where $E_{\text{X}}^{\text{DFT}}[\mathbf{P}]$ and
$E_{\text{X}}^{\text{DFT}}[\mathbf{\hat{P}}]$ are the DFT exchange energies
for the primary and auxiliary density matrices, respectively. The RI
approximation, however, allows to substitute the computationally demanding
four-center ERIs by just 2- and 3-center integrals by factorizing them via

$\displaystyle(ia|jb)_{\text{RI}}=\sum_{PQ}\,(ia|P)(P|Q)^{-1}(Q|jb),$ (5)

where $(P|Q)^{-1}$ is the inverse of the Coulomb metric over auxiliary
Gaussian basis functions, i.e.

$\displaystyle(P|Q)=\int d\mathbf{r}_{1}\int
d\mathbf{r}_{2}\,\phi_{P}(\mathbf{r}_{1})\,\frac{1}{|\mathbf{r}_{1}-\mathbf{r}_{2}|}\,\phi_{Q}(\mathbf{r}_{2}).$
(6)

Since the latter is a positive definite matrix, its inverse can be efficiently
obtained by means of the Cholesky decomposition

$\displaystyle(P|Q)=\sum_{R}\,L_{PR}L_{RQ}^{T}$ (7)

followed by an inversion of the triangular matrix $\mathbf{L}$, i.e.

$\displaystyle(P|Q)^{-1}=\sum_{R}\,L_{PR}^{-T}L_{RQ}^{-1}.$ (8)

In this way, the factorization of the integrals $(ia|jb)$ can be written in
compact form as

$\displaystyle(ia|jb)_{\text{RI}}=\sum_{P}\,B_{P}^{ia}B_{P}^{jb},$ (9)

where

$\displaystyle B_{P}^{ia}=\sum_{R}\,(ia|R)L_{PR}^{-1}.$ (10)

Therein, the three-center integrals $(ia|R)$ are computed starting from
integrals over atomic orbitals $(\mu\nu|R)$, so that the final expression for
the elements $B_{P}^{ia}$ reads as

$\displaystyle(ia|P)=\sum_{\nu}\,C_{\nu a}\sum_{\mu}\,C_{\mu
i}\sum_{R}\,(\mu\nu|R)L_{PR}^{-1},$ (11)

where $\mathbf{C}$ is the molecular orbital coefficient matrix, i.e.
$\mathbf{CC}^{T}=\mathbf{P}$.

Another major development direction is the design of novel linear-scaling
algorithms [8]. Within CP2K, the relation

$\displaystyle\text{sign}\begin{pmatrix}\mathbf{0}&\mathbf{A}\\\
\mathbf{I}&\mathbf{0}\end{pmatrix}=\text{sign}\begin{pmatrix}\mathbf{0}&\mathbf{A}^{1/2}\\\
\mathbf{A}^{-1/2}&\mathbf{0}\end{pmatrix}$ (12)

is employed together with various iterative methods to compute the matrix sign
function

$\displaystyle\text{sign}\left(\mathbf{A}\right)=\mathbf{A}\left(\mathbf{A}^{2}\right)^{-1/2}$
(13)

to yield the inverses and (inverse) square roots of large sparse matrices
$\mathbf{A}$ with a computational effort that scales just linearly with system
size. Most importantly, however, the sign function can also be used for the
purification of the Kohn-Sham matrix $\mathbf{H}_{\text{KS}}$ into
$\mathbf{P}$ by the use of

$\displaystyle\mathbf{P}=\frac{1}{2}\left(\mathbf{I}-\text{sign}\left(\mathbf{S}^{-1}\mathbf{H}_{\text{KS}}-\mu\mathbf{I}\right)\right)\mathbf{S}^{-1}.$
(14)

Figure 1: Total walltime for a single AIMD step for bulk water at ambient
conditions as a function of system size on a Cray XT5.

The superior efficiency and linear-scaling potential of this approach has been
demonstrated at the DFT level using a double-$\zeta$ valence polarization
basis set, for up to a million of atoms in a massively parallel
implementation, as shown in Fig. 1.

##  3 Meeting the Exascale Challenges

Standing at the brink of the exascale era, we can see that fundamental
disruption in programming models and processor technology that were predicted
a decade ago did not happen. From the developer’s point of view, three trends
originating in the petascale area are continued and emphasized: 1) the need of
massive parallelism in terms of node and thread counts, 2) an even more
widespread use of GPUs, and 3) the availability of mixed-precision matrix or
tensor operation hardware accelerators.

CP2K is fundamentally well-suited to scale to immense levels of parallelism
because it was designed as a massively parallel MPI application right from the
outset. With the advent of multi-core processors, OpenMP directives were added
to important loops, while leaving the underlying data layout unchanged to
support multi-threading and hybrid parallelism. With that many concurrent
threads the data has to be partitioned to prevent bottlenecks from reductions
or atomic access collisions. The new library for Distributed Block-sparse
Matrices (DBM) is a first step in this direction. It uses a fixed assignment
of matrix block rows to threads, which eliminates the need for synchronization
from most operations. We plan to refactor other primitives in the same way, in
particular the grid data structures that power methods like GPW have a lot to
gain from per-thread partitioning.

The modular structure and use of modern Fortran 2008 allowed the addition of
GPU support at the level of GPU-accelerated libraries (e.g. DBM, COSMA, SpFFT,
SPLA, grid, pw and Sirius), some of which were also spun off from CP2K as
stand-alone libraries such as e.g. DBCSR. The challenge for the exascale era
is the evolution and increasing diversity of GPU architectures. In the early
days we mostly struggled with finding the sweet spot within the tight
constrains set by small register files and scarce shared memory. In present
systems, the bottleneck has now shifted to the PCI bus, where we are often
limited by host-to-device communication. The way forward is to use GPU-aware
MPI, which we are currently adding to DBM and will later roll our to other
parts of the code. The status and results of these efforts can be tracked on
the CP2K GPU dashboard at https://www.cp2k.org/gpu.

Lastly, we expect that the matrix and tensor processing units will have a
profound and long-lasting effect on method and code development. These
computing elements can achieve one to two order of magnitude more FLOPs for
dense linear algebra operations in reduced or mixed precision, e.g. FP16
operations with FP32 accumulation. To exploit this potential, it is key to
develop methods that heavily rely on dense linear algebra on medium sized
local matrices, whose numerical inaccuracies due to mixed- and low-precision
arithmetic can be rigorously compensated by the design of a modified Langevin-
type equation [2].

In recent work, we have developed the submatrix method that is specifically
designed with these design principles in mind. The core idea, as illustrated
in Fig. 2, is to convert the evaluation of a matrix function on a large
distributed sparse matrix into a large-scale parallel evaluation of the matrix
function on many dense, but much smaller matrices.

Figure 2: Schematic representation of the steps of the submatrix method for
the approximate calculation of a matrix function $f(\mathbf{A})$ of a large
sparse matrix $\mathbf{A}$. The first step is the construction of a submatrix
$\mathcal{T}_{i}(\mathbf{A})$ for every column of the matrix $\mathbf{A}$.
Then the matrix function is applied to the dense submatrices, i.e.
$f(\mathcal{T}_{i}(\mathbf{A}))$, and finally the relevant result columns are
inserted into the sparse result matrix. Figure reproduced from [8] under CC-
BY.

In combination with the purification scheme of Eq. 14, we were able to perform
record sized linear-scaling electronic structure computations on systems with
more than 100 million atoms, thereby achieving a sustained performance of 324
PFLOPs with an efficiency of more than $67\%$ [8].

##  4 Concluding Remarks

We conclude by noting that beside exploiting upcoming exascale supercomputers
to facilitate ever more complex and accurate simulations, in the future
alternative approaches such as deep neural networks (DNN), as well as quantum
computing algorithms will find their application within CP2K. On the one hand
this can be so-called learning on-the-fly schemes, in which the decision if
the calculation of an observable can be outsourced to an already sufficiently
accurate DNN or needs to be conducted explicitly, which can then be used as
training data to directly relearn the DNN, is done during the simulation
itself. On the other hand, all sort of hybrid approaches can be imagined,
where the majority of a particular quantity is explicitly calculated using a
computationally simple, but approximate electronic structure method, which is
augmented by a correction term that is represented by a DNN for instance [9].
De facto exact configuration interaction, or reduced density-matrix functional
(RDMF) theory simulations will become feasible by the usage of hybrid quantum-
classical algorithms in which the quantum mechanical expectation values of the
RDMF are evaluated on a quantum computer, whereas the parameters of the trial
states are optimized by a Car-Parrinello-like constrained minimization scheme
within CP2K on a classical computer [10].

##  5 Acknowledgements

The authors would like to thank the whole CP2K development team, who had been
contributing to the code over the last two decades.

## References

## References

  * [1] T. D. Kühne et al., CP2K: An electronic structure and molecular dynamics software package - Quickstep: Efficient and accurate electronic structure calculations, J. Chem. Phys. 152, 194103 (2020).
  * [2] T. D. Kühne, M. Krack, F. R. Mohamed, and M. Parrinello, Efficient and accurate Car-Parrinello-like approach to Born-Oppenheimer molecular dynamics, Phys. Rev. Lett. 98, 066401 (2007).
  * [3] M. Iannuzzi and J. Hutter, Inner-shell spectroscopy by the Gaussian and augmented plane wave method, Phys. Chem. Chem. Phys. 9, 1599 (2007).
  * [4] A. Bussy and J. Hutter, Efficient and low-scaling linear-response time-dependent density functional theory implementation for core-level spectroscopy of large and periodic systems, Phys. Chem. Chem. Phys. 23, 4736 (2021).
  * [5] A.-S. Hehn, B. Sertcan, F. Belleflamme, S. K. Chulkov, M. B. Watkins, and J. Hutter, Excited-state properties for extended systems: Efficient hybrid density functional methods, doi:10.26434/chemrxiv-2022-swkw9.
  * [6] M. Guidon, J. Hutter, and J. VandeVondele, Auxiliary Density Matrix Methods for Hartree-Fock Exchange Calculations, J. Chem. Theory Comput. 6, 2348 (2010).
  * [7] J. Hutter, and J. Wilhelm, V. V. Rybkin, M. Del Ben, and J. VandeVondele, MP2-and RPA-Based Ab Initio Molecular Dynamics and Monte Carlo Sampling, In Handbook of Materials Modeling: Methods: Theory and Modeling (Cham: Springer) p 523.
  * [8] R. Schade et al., Towards electronic structure-based ab-initio molecular dynamics simulations with hundreds of millions of atoms, Parallel Computing 111, 102920 (2022).
  * [9] M. Pauletti, V. V. Rybkin, and M. Iannuzzi, Subsystem Density Functional Theory Augmented by a Delta Learning Approach to Achieve Kohn–Sham Accuracy, J. Chem. Theory Comput. 17, 6423 (2021).
  * [10] R. Schade, C. Bauer, K. Tamoev, L. Mazur, C. Plessl, and T. D. Kühne, Parallel Quantum Chemistry on Noisy Intermediate-Scale Quantum Computers, arXiv:2202.02417.


# Machine learning time-local generators of open quantum dynamics

Paolo P. Mazza  Institut für Theoretische Physik and Center for Quantum
Science, Universität Tübingen, Auf der Morgenstelle 14, 72076 Tübingen,
Germany  Dominik Zietlow  Max Planck Institute for Intelligent Systems, Max-
Planck-Ring 4, 72076 Tübingen, Germany  Federico Carollo  Institut für
Theoretische Physik and Center for Quantum Science, Universität Tübingen, Auf
der Morgenstelle 14, 72076 Tübingen, Germany  Sabine Andergassen  Institut für
Theoretische Physik and Center for Quantum Science, Universität Tübingen, Auf
der Morgenstelle 14, 72076 Tübingen, Germany  Georg Martius  Max Planck
Institute for Intelligent Systems, Max-Planck-Ring 4, 72076 Tübingen, Germany
Igor Lesanovsky  Institut für Theoretische Physik and Center for Quantum
Science, Universität Tübingen, Auf der Morgenstelle 14, 72076 Tübingen,
Germany  School of Physics and Astronomy, University of Nottingham,
Nottingham, NG7 2RD, UK  Centre for the Mathematics and Theoretical Physics of
Quantum Non-Equilibrium Systems, University of Nottingham, Nottingham, NG7
2RD, UK

###### Abstract

In the study of closed many-body quantum systems one is often interested in
the evolution of a subset of degrees of freedom. On many occasions it is
possible to approach the problem by performing an appropriate decomposition
into a bath and a system. In the simplest case the evolution of the reduced
state of the system is governed by a quantum master equation with a time-
independent, i.e. Markovian, generator. Such evolution is typically emerging
under the assumption of a weak coupling between the system and an infinitely
large bath. Here, we are interested in understanding to which extent a neural
network function approximator can predict open quantum dynamics — described by
time-local generators — from an underlying unitary dynamics. We investigate
this question using a class of spin models, which is inspired by recent
experimental setups. We find that indeed time-local generators can be learned.
In certain situations they are even time-independent and allow to extrapolate
the dynamics to unseen times. This might be useful for situations in which
experiments or numerical simulations do not allow to capture long-time
dynamics and for exploring thermalization occurring in closed quantum systems.

_Introduction_ —The investigation of isolated quantum systems out of
equilibrium is a central problem in physics D’Alessio  _et al._ (2016);
Polkovnikov  _et al._ (2011). Often, partial information that concerns a
spatially localized subsystem ${S}$ is sufficient to study dynamical effects
related to relaxation or thermalization. This information is encoded in the
so-called reduced quantum state $\rho_{S}$, [c.f. Fig. 1(a,b)], whose time
dependence is obtained by “integrating out” the remainder of the evolving
many-body system, which is often referred to as bath ${B}$ [see Fig. 1(c)]. In
certain settings, e.g. for systems which are weakly interacting with a large
environment Breuer  _et al._ (2002); Gardiner and Zoller (2004), the dynamics
of the reduced state is effectively described by a quantum master equation
Lindblad (1976); Gorini  _et al._ (1976); Alicki and Lendi (2007). In the
simplest manifestation, the dynamics of the reduced quantum state $\rho_{S}$
is then governed by a time-independent generator $\mathcal{L}$ acting only on
$S$:

$\displaystyle\rho_{S}(t):=\textnormal{Tr}_{B}\left(U_{t}\left(\rho_{S}\otimes\rho_{B}\right)U_{t}^{\dagger}\right)\approx
e^{t\,\mathcal{L}}\left[\rho_{S}\right].$ (1)

Here, $U_{t}$ is the unitary evolution operator of the full many-body system,
the state $\rho_{S}\otimes\rho_{B}$ is the initial state of the many-body
system in product form, and $\textnormal{Tr}_{B}$ indicates the average over
the bath ${B}$. However, the generator does not need to be time-independent
and, in most generic instances Chruściński and Kossakowski (2010), it may
depend locally on time, $\mathcal{L}_{t}$, as sketched in Fig. 1(c).

Figure 1: Reduced local dynamics of a many-body quantum system. The reduced
state $\rho_{S}$ describes local degrees of freedom of a quantum spin system.
The remainder of the system, $B$, takes the role of a “bath”, with initial
state $\rho_{B}$. We consider two scenarios: (a) Model I: spin chain with open
boundary conditions, where the interaction strength $V$ and range $\alpha$ are
varied [see Eq. (2)]; (b) Model II: spin chain with closed boundaries and
nearest-neighbor interaction, where the system-bath coupling $V^{\prime}$ and
the initial temperature $\beta$ of the bath are modified. (c) The dynamics of
the reduced state, $\rho_{S}(t)$, is obtained by tracing out the degrees of
freedom of $B$ and to be learned by a neural network. The learned generator
may be either Markovian (time-independent) or time-dependent.

In this paper, we show how simple neural network architectures can learn time-
local generators that govern the local dynamics of a closed quantum spin
system. The input of the networks consists of the time-dependent average
values of the reduced system observables from which the dynamical generator is
estimated. Firstly, we consider an architecture which provides a time-averaged
generator $\overline{\mathcal{L}}$ corresponding to an effectively Markovian
description of the system dynamics. Our findings indicate that this can yield
indeed a valid approximation for the generator even beyond typical scenarios
justifying a Markovian weak-coupling assumption Breuer  _et al._ (2002);
Gardiner and Zoller (2004). Once the generator is known, the network can
further be exploited to make predictions for times which go beyond the
previously analyzed time frame, as is sketched in Fig. 2. To study settings
with a time-dependent generator, $\mathcal{L}_{t}$, we use a different neural
network architecture which is based on so-called hypermodels and allows to
assess the time-dependence of the generator of the reduced dynamics. Our work
links to recent efforts aiming at understanding quantum dynamics with neural
networks Yoshioka and Hamazaki (2019); Hartmann and Carleo (2019a); Nagy and
Savona (2019); Krastanov  _et al._ (2020); Carleo  _et al._ (2019); Krastanov
_et al._ (2019); Carleo and Troyer (2017); Hartmann and Carleo (2019b); Liu
_et al._ (2020); Luchnikov  _et al._ (2020); Martina  _et al._ (2021); Miles
_et al._ (2020). Our approach, which uses machine learning tools and provides
a directly interpretable object such as the physical generator of the
dynamics, highlights a possible route for the application of neural networks
in the study of the long-time dynamics of local observables in closed quantum
systems. Moreover, it may find applications in the context of quantifying the
degree of non-Markovianity in reduced quantum dynamics Chruściński and
Kossakowski (2010).

Figure 2: Effective time-independent generator: Training and extrapolation. To
obtain a time-independent generator we use a multi-layer perceptron (MLP). The
inputs are the time-dependent expectation values of the reduced system
observables — for a single spin contained in the Bloch vector $\mathbf{v}_{\rm
b}$– in a given training time window (shaded blue region) and for several
initial Bloch vectors (top to bottom). The MLP learns the “propagator”
$M[\theta^{*}]$, which depends on the optimized network parameters
$\theta^{*}$ but not on the actual time. Considering an Euler integration step
of length $dt$, this provides a matrix $\bar{{L}}$ representing the action of
an “averaged” generator on the Bloch vector. When the reduced dynamics is
Markovian, this learned generator allows for the network to make predictions
for unseen future times (red shaded region).

_Many-body spin models_ —For illustrating our ideas, we consider two one-
dimensional systems (in the following referred to as Model I and II)
consisting of $N$ interacting spins, which are inspired by recent state-of-
the-art quantum simulator experiments with Rydberg atoms Bloch  _et al._
(2012); Kim  _et al._ (2018); Ebadi  _et al._ (2020); Browaeys and Lahaye
(2020); Bloch  _et al._ (2008) or trapped ions Islam  _et al._ (2011); Bohnet
_et al._ (2016); Joshi  _et al._ (2020). Model I has open boundary conditions
[see Fig. 1(a)]. It features power-law interactions and is described by the
Hamiltonian

$H_{\mathrm{I}}=\Omega\sum_{i=1}^{N}\sigma_{i}^{x}+V\sum^{N}_{i<j}\frac{n_{i}n_{j}}{|i-j|^{\alpha}}.$
(2)

Here, $n_{i}=(1+\sigma^{z}_{i})/2$ is a projector on the “up”-state of the
$i$-th spin and $\sigma^{x,y,z}_{i}$ are Pauli matrices. The parameter
$\alpha$ accounts for the interaction range, i.e. the power-law decay of the
interaction potential, and $V$ controls the overall interaction strength with
respect to the strength $\Omega$ of a transverse field term. The reduced
system $S$, for this model, is the middle spin of the chain, as shown in Fig.
1(a). Model II is a closed spin chain [sketched in Fig. 1(b)], with nearest-
neighbor interactions and Hamiltonian

$H_{\mathrm{II}}=\Omega\sum_{i=2}^{N}\sigma_{i}^{x}+V\sum_{j=2}^{N-1}n_{j}n_{j+1}+\Omega^{\prime}\sigma_{1}^{x}+V^{{}^{\prime}}n_{1}(n_{2}+n_{N}).$
(3)

Here, the spin with label $1$ is singled out as reference spin (forming the
system $S$) and the interaction strength with its neighbors is given by
$V^{\prime}$, while the interaction strength among all other spins is $V$.
Moreover, we assume that the strength $\Omega^{\prime}$ of the transverse
field acting on the reference spin can be controlled. Both models capture a
whole variety of different scenarios, including short- and long-range
interactions, absence and presence of translational invariance as well as weak
and strong coupling between local and bath degrees of freedom. They moreover
encompass a number of standard scenarios often explored in quantum many-body
physics, such as the so-called PXP-model Sun and Robicheaux (2008); Ates  _et
al._ (2012); Turner  _et al._ (2018), the Ising model in the presence of
longitudinal and transverse fields, and all-to-all connected spin systems with
resemblance to the Dicke model Emary and Brandes (2003) or the Lipkin-Meshkov-
Glick model Dusuel and Vidal (2004).

Our aim is to investigate the reduced dynamics of local degrees of freedom, as
for instance a single spin $S$, as depicted in Fig. 1(a,b). To this end, we
assume that the system is initialized at time $t=0$ in the state
$\rho=\rho_{S}\otimes\rho_{B}$. In our studies concerning Model I, we assume
the bath to be in the infinite temperature state,
$\rho_{B}\propto\mathds{1}_{B}$. When considering Model II, instead, we assume
a finite temperature situation with $\rho_{B}\propto
e^{-\beta\tilde{H}_{\mathrm{II}}}$, where the Hamiltonian
$\tilde{H}_{\mathrm{II}}$ is the one of Eq. (3) with
$\Omega^{\prime}=V^{\prime}=0$, and $\beta$ is the inverse temperature. We let
the state $\rho$ evolve according to the unitary dynamics $\rho(t)=U_{t}\rho
U^{\dagger}_{t}$, with $U_{t}=\textnormal{e}^{-iHt}$ from which we obtain the
full dynamical information about the subsystem $S$. Throughout we will make
use of the Bloch vector representation of the subsystem’s density matrix:
$\rho_{S}(t)=(\mathbf{v}_{\rm b}(t)\cdot\mathbf{\sigma})/2$ with
$\mathbf{\sigma}=(\mathds{1},\sigma^{x},\sigma^{y},\sigma^{z})$. Here the
Bloch vector $\mathbf{v}_{\rm
b}(t)=\left(1,\braket{\sigma_{k}^{x}}_{t},\braket{\sigma_{k}^{y}}_{t},\braket{\sigma_{k}^{z}}_{t}\right)$
is given in terms of the expectation $\braket{\cdot}_{t}$ with respect to the
full quantum state at time $t$, and $k$ labels the system spin [see Fig.
1(a,b)].

In the following, we will characterize the generator of this reduced dynamics;
that is, we want to find a generator acting solely on $S$ that reproduces (or
approximates) the dynamics of $\mathbf{v}_{\rm b}(t)$ for the two models under
consideration [Fig. 1(a,b)]. For Model I we explore for which parameter
combinations $(V,\alpha)$, the local generator, which is in principle time-
dependent $\mathcal{L}_{t}$, can be approximated by an effective “averaged”
time-independent one, $\mathcal{L}_{t}\approx\bar{\mathcal{L}}$. We remark
here that, while $B$ acts as a fictitious bath for subsystem $S$, the
situation described here is far from that of typical weak-coupling limits.
Indeed, $B$ is a finite quantum system and its Hamiltonian has a discrete
spectrum. As such, $B$ can hardly be thought of as an infinite Markovian bath
of bosonic oscillators, whose state is unaffected by the interaction with $S$.
The only aspect in common with a weak-coupling limit is that the interaction
strength can be made small. However, this is not sufficient to argue that the
local dynamics may be described by time-independent generators.

For Model II we consider the case in which an “isolated” spin (on site $1$)
interacts with $B$ with an independently tunable strength $V^{\prime}$.
Experimentally such a situation can be realized, e.g., by exploiting the
distance dependence of dipolar interactions of atoms in Rydberg states
Browaeys and Lahaye (2020). Also for Model II we investigate how well the
generator of the reduced dynamics is approximated by a time-independent model
when the initial state of $B$ (parametrized by the inverse temperature
$\beta$) and the interaction strength $V^{{}^{\prime}}$ is varied. In
addition, for this scenario, we use a hypermodel in order to analyze regimes
that necessitate a time-dependent generator $\mathcal{L}_{t}$.

_Network architectures and training_ —In order to learn dynamical generators
for the reduced dynamics of subsystem $S$, we employ machine learning
algorithms. Our approach is completely data-driven, i.e. the network has no
prior information about the physical system. For learning time-independent
generators, we use a linear multi-layer perceptron (MLP) architecture (see
Fig. 2) which turns out to be the simplest possible one. Every data-point of
our data-set $\mathcal{D}$ contains the triple $(\mathbf{v}_{\rm
b},\mathbf{v}_{\rm b}^{\prime},t)$ where $\mathbf{v}_{\rm b}:=\mathbf{v}_{\rm
b}(t)$ is the value of the Bloch vector at a given time-point $t$, and
$\mathbf{v}_{\rm b}^{\prime}$ is the value of the Bloch vector after a further
discrete time step $dt$, i.e. $\mathbf{v}_{\rm b}^{\prime}:=\mathbf{v}_{\rm
b}(t+dt)$.

The input of our model is $\mathbf{v}_{\rm b}$, and the output is the vector
$\mathbf{o}(t)=M[{\theta}]\mathbf{v}_{\rm b}$, where $M[{\theta}]$ is a
$4\times 4$-matrix (matching the dimension of the Bloch vector), that depends
on the parameters ${\theta}$ of the network. To optimize the network, we
introduce the following loss function, which is given by the norm of the
difference between the next time step $\mathbf{v}_{\rm b}^{\prime}$ and the
output of the network:

$C(\theta)=\mathbb{E}_{(\mathbf{v}_{\rm b},\mathbf{v}_{\rm
b}^{\prime})\sim\mathcal{D}}\left[\left\lVert M[{\theta}]\mathbf{v}_{\rm
b}-\mathbf{v}_{\rm b}^{\prime}\right\rVert\right],$ (4)

where the expectation is taken over the training data $\mathcal{D}$.
Minimizing the above function thus provides a model $M[{\theta}^{\ast}]$ for
the propagation of the Bloch vector over an infinitesimal time step $dt$. This
model is related to the generator of the local dynamics: indeed, we have
$M[{\theta}^{\ast}]=\mathds{1}+\overline{{L}}dt$, where $\bar{L}$ is a
representation of the time-averaged generator $\overline{\mathcal{L}}$ acting
on Bloch vectors. The training data-set is constructed by evolving the Bloch
vector $\mathbf{v}_{\rm b}$ over a time interval $[0,T_{\mathrm{train}}]$. The
trajectories are generated evolving $100$ randomly chosen initial states of
the form

$\rho_{0}=\frac{1}{2}\left(\begin{array}[]{@{}cc@{}}1+z&x-iy\\\
x+iy&1-z\end{array}\right)\otimes\rho_{B},$ (5)

with $x^{2}+y^{2}+z^{2}<1$, i.e. we consider non-pure initial states. We split
our data-set such that we consider the $80\%$ of it for training, i.e. finding
the best parameters of our model. The remaining $20\%$ form the validation set
to assess the performance of the model built in the training phase.

To encode and quantify the time-dependence of dynamical generators,
$\mathcal{L}_{t}$, we use a different architecture and consider a so-called
hypermodel, which computes network weights based on a context input. In our
case, the context input is the time $t$ and a multi-layer perceptron (MLP) $F$
with non-linear activation functions transforms this input to the $4\times 4$
matrix of the generator $M[F(t;\theta)]$.

The data-set is organized in the same way as for the time-independent model.
The only difference with respect to the previous case is that we also pass the
information about the actual time step $t$. The network parameters $\theta$
are optimized by minimizing the loss function

$C_{\mathrm{H}}(\theta)=\mathbb{E}_{(\mathbf{v}_{\rm b},\mathbf{v}_{\rm
b}^{\prime},t)\sim\mathcal{D}}\left[\left\lVert M[F(t;\theta)]\mathbf{v}_{\rm
b}-\mathbf{v}_{\rm b}^{\prime}\right\rVert\right].$ (6)

By training the model in this way we obtain a time-dependent representation of
the propagator of the local dynamics, i.e.
$M[F(t;\theta^{*})]=\mathds{1}+L_{t}dt$. For the training set, we choose
quantum states as in the previous case. For the evaluation data-set, instead,
we consider

$\rho_{0}=\left(\begin{array}[]{@{}cc@{}}1-c&0\\\
0&c\end{array}\right)\otimes\rho_{B},$ (7)

with $c$ being a random number, $c\in(0.01,0.7)$111See supplemental material
for details.

Figure 3: Time-independent generator. (a) Model I: Comparison between the
exact dynamics of the Bloch vector and the one obtained by learning an
effective time-independent generator with neural networks. The density plot
shows the time-averaged norm difference, $\overline{\epsilon}$ [Eq. (8)
averaged over $5$ randomly selected values $c$ of the initial conditions, see
Eq. (7)], for various parameter choices $(V,\alpha)$. The data is calculated
for $N=9$ spins with a training time of $T_{\mathrm{train}}=10/\Omega$,
$dt=0.01/\Omega$ and a total time $T_{\mathrm{tot}}=20/\Omega$ (see main text
for details). We also show two examples for the Bloch vector evolution, where
we choose $c=1$ and $T_{\mathrm{tot}}=20/\Omega$ for illustration. Here solid
curves correspond to the components of the Bloch vector propagated with the
learned time-averaged generator. The dashed lines correspond to the
numerically exact solution. (b) Same as panel (a), but for Model II and the
parameter set $(V^{\prime},\beta)$.

_Time-independent generators_ —With the above numerical approach, we can now
investigate under which circumstances the generator of the dynamics of the
subsystem $S$ is time-independent. We consider first the open spin chain [see
Fig. 1(a) and Eq. (2)] and explore different values of $V$ and $\alpha$. To
test the performance of the model after the training, we study the _time-
averaged_ norm difference between the exact Bloch vector $\mathbf{v}^{\rm
ex}_{\rm b}$ — obtained by simulating the full many-body spin chain dynamics —
and the time-evolution of the same quantity as predicted by our model
$\mathbf{v}^{\rm mod}_{\rm b}$:

$\epsilon=\frac{1}{T_{\rm tot}}\int_{0}^{T_{\rm
tot}}\\!\\!dt\left\lVert\mathbf{v}^{\rm mod}_{\rm b}(t)-\mathbf{v}_{\rm
b}^{\rm ex}(t)\right\rVert\,,$ (8)

computed by numerical integration.

As shown in Fig. 3(a), the error $\epsilon$ is small for short-range (large
$\alpha$) and weakly interacting (with small ratio $V/\Omega$) spin chains.
This means that our time-independent model captures the relevant features of
the reduced dynamics, and suggests that for finite systems with sufficiently
weak interactions, the dynamics of local degrees of freedom can indeed be
effectively described by a _time-independent_ generator $\bar{\mathcal{L}}$.
Furthermore, for this Markovian regime of the dynamics, our model allows to
make predictions for times exceeding the training time.

We now turn to the discussion of Model II [see Fig. 1(b) and Eq. (3)]. Using
the same procedure explained above, we can learn the time-averaged dynamical
generator for subsystem $S$. In Fig. 3(b) we summarize the results of this
analysis. We display the time-averaged error, as defined in Eq. (8), for
different parameter choices, which in this case are the (inverse) temperature
$\beta$ of the initial state $\rho_{B}$ and the interaction strength
$V^{{}^{\prime}}$. We furthermore take $\Omega^{\prime}=\Omega$. As shown in
the figure the dynamics is well approximated by a time-independent generator,
$\bar{\mathcal{L}}$, for weak interactions $V^{\prime}$ and for small $\beta$.
Here also the extrapolation to times that exceed $T_{\mathrm{train}}$ is
possible. The data in Fig. 3 shows that in certain parameter regimes the
approximation of the dynamics through a time-independent generator does not
work well. Here significant deviations are even observable during the training
period.

Figure 4: Hypermodel and time-dependent generator. Hypermodels allow to encode
time-dependent generators $\mathcal{L}_{t}$. The time-dependence is quantified
through the derivative of $L_{t}$ with respect to time, expressed by the
positive quantity $\xi(t)$ (see text for definition). Panel (a) shows three
example curves corresponding to Model II ($N=9$). The function $\xi(t)$
depends only on the parameters of the model and is independent from the
initial condition chosen. Generally, the stronger the interaction
$V^{\prime}$, the more pronounced the time dependence. This trend is clearly
visible in panel (b), where we show the averaged value of $\xi(t)$ in the
interval $0$ to $T=10/\Omega$, denoted as $\Xi$. In order to obtain the
density plot we averaged, also in this case, over $5$ different values of $c$
in Eq. (7), also in this case $dt=0.01\Omega$.

_Time-dependent time-local generators_ —We are now interested in analyzing how
strongly the dynamical generator depends on time, in those instances where the
model introduced in the previous section fails. To this end, we adopt a
hypermodel which can take as input the information about the running time. In
this way, the neural network is able to learn an optimal parametrization of
the “propagator” $M[F(t;\theta^{*})]=\mathds{1}+{{L}_{t}}dt$ for the Bloch
vector dynamics which explicitly depends on time. Here, the matrix $L_{t}$
encodes the action of the generator on the Bloch vector and, essentially,
implements the differential equations obeyed by the entries of the Bloch
vector $\mathbf{v}_{\rm b}$. This architecture allows to learn and accurately
reproduce the dynamics of Model II for all studied parameter regimes (see
Supplemental Material for examples). This is not surprising, but gives us a
handle for analyzing the time-dependence of the dynamical generator $L_{t}$.
To this end we consider the positive quantity $\xi(t)=\sqrt{{\rm
Tr}(\dot{L}_{t}^{\dagger}\dot{L}_{t})}$, where the dot denotes the derivative
with respect to time. This can only be zero when $\dot{L}_{t}=0$, i.e. when
the action of the generator is not depending on the running time. To quantify
an overall time-dependence within a time window $[0,T]$, we define the time-
averaged value $\Xi=\frac{1}{T}\int_{0}^{T}dt\,\xi(t)$. In Fig. 4 we show the
corresponding data. For strong interactions $V^{\prime}$, $\xi(t)$ is non-zero
throughout which indicates a strong explicit time-dependence. This is also
reflected in a large $\Xi$. For weak interactions, on the other hand, $\xi(t)$
remains small for all times considered. The oscillation we attribute to the
finite size of the system. This confirms that, in such parameter regime, it is
indeed possible to approximate the dynamics of the Bloch vector by means of a
time-independent matrix, as considered in the previous section.

_Conclusions_ —We have presented two simple examples of neural network
architectures that can learn the dynamical features of reduced quantum states.
When such time evolution is effectively Markovian, the network can find a
suitable approximation for the generator of the local dynamics. Here one can
extrapolate the dynamics of reduced degrees of freedom to times that were
unexplored during the training procedure. This possibility is particularly
promising for applications in combination with tensor networks, which perform
extremely well for short times. A neural network could learn the time-
independent generator during this time interval and then possibly extrapolate
to long times. In principle, this might enable the exploration of the onset of
stationary or thermalization regimes. When the dynamics is not of Markovian
type, we have shown that hypermodels can recover the generator of the reduced
quantum dynamics. This allows to quantify non-Markovian effects Benatti  _et
al._ (2016), which manifest in an explicit time-dependence of the dynamical
generator Ref. Chruściński and Kossakowski (2010). Furthermore, being able to
reconstruct the generator of the reduced dynamics, as we have done here with
hypermodels, makes it possible to explore different measures of non-
Markovianity based on non-divisibility criteria for quantum dynamical time-
evolutions Breuer  _et al._ (2009); Laine  _et al._ (2010); Breuer (2012);
Rivas  _et al._ (2014).

A possible future development in this regard could be the application of more
advanced machine learning methods for learning time correlations in time
series. This can be achieved with models borrowed from language studies (e.g.
long-short term memory recurrent neural networks (LSTM) Hochreiter and
Schmidhuber (1997) and transformers Vaswani  _et al._ (2017)) or with
algorithms geared towards more interpretable models by either learning
analytical expressions of the differential equation Sahoo  _et al._ (2018) or
by modelling the time-evolution dynamics directly Chen  _et al._ (2018).

###### Acknowledgements.

_Acknowledgments_ —We acknowledge financial support from the Deutsche
Forschungsgemeinschaft (DFG, German Research Foundation) under Germany’s
Excellence Strategy – EXC-Number 2064/1 – Project number 390727645. IL
acknowledges funding from the “Wissenschaftler Rückkehrprogramm GSO/CZS” of
the Carl-Zeiss-Stiftung and the German Scholars Organization e.V., and through
the DFG projects number 428276754 (SPP GiRyd) and 435696605. The authors thank
the International Max Planck Research School for Intelligent Systems (IMPRS-
IS) for supporting DZ.

## References

  * D’Alessio  _et al._ (2016) L. D’Alessio, Y. Kafri, A. Polkovnikov, and M. Rigol, Adv. Phys. 65, 239 (2016).
  * Polkovnikov  _et al._ (2011) A. Polkovnikov, K. Sengupta, A. Silva, and M. Vengalattore, Rev. Mod. Phys. 83, 863 (2011).
  * Breuer  _et al._ (2002) H.-P. Breuer, F. Petruccione, _et al._ ,  _The theory of open quantum systems_ (Oxford University Press on Demand, 2002).
  * Gardiner and Zoller (2004) C. Gardiner and P. Zoller,  _Quantum noise: a handbook of Markovian and non-Markovian quantum stochastic methods with applications to quantum optics_ (Springer Science & Business Media, 2004).
  * Lindblad (1976) G. Lindblad, Commun. Math. Phys. 48, 119 (1976).
  * Gorini  _et al._ (1976) V. Gorini, A. Kossakowski, and E. C. G. Sudarshan, Journal of Mathematical Physics 17, 821 (1976).
  * Alicki and Lendi (2007) R. Alicki and K. Lendi,  _Quantum dynamical semigroups and applications_ , Vol. 717 (Springer, 2007).
  * Chruściński and Kossakowski (2010) D. Chruściński and A. Kossakowski, Phys. Rev. Lett. 104, 070406 (2010).
  * Yoshioka and Hamazaki (2019) N. Yoshioka and R. Hamazaki, Phys. Rev. B 99, 214306 (2019).
  * Hartmann and Carleo (2019a) M. J. Hartmann and G. Carleo, Phys. Rev. Lett. 122, 250502 (2019a).
  * Nagy and Savona (2019) A. Nagy and V. Savona, Phys. Rev. Lett. 122, 250501 (2019).
  * Krastanov  _et al._ (2020) S. Krastanov, K. Head-Marsden, S. Zhou, S. T. Flammia, L. Jiang, and P. Narang, arXiv:2009.03902 (2020).
  * Carleo  _et al._ (2019) G. Carleo, I. Cirac, K. Cranmer, L. Daudet, M. Schuld, N. Tishby, L. Vogt-Maranto, and L. Zdeborová, Rev. Mod. Phys. 91, 045002 (2019).
  * Krastanov  _et al._ (2019) S. Krastanov, S. Zhou, S. T. Flammia, and L. Jiang, Quantum Sci. Technol. 4, 035003 (2019).
  * Carleo and Troyer (2017) G. Carleo and M. Troyer, Science 355, 602 (2017).
  * Hartmann and Carleo (2019b) M. J. Hartmann and G. Carleo, Phys. Rev. Lett. 122, 250502 (2019b).
  * Liu  _et al._ (2020) Z. Liu, L.-M. Duan, and D.-L. Deng, arXiv:2008.05488 (2020).
  * Luchnikov  _et al._ (2020) I. A. Luchnikov, S. V. Vintskevich, D. A. Grigoriev, and S. N. Filippov, Phys. Rev. Lett. 124, 140502 (2020).
  * Martina  _et al._ (2021) S. Martina, S. Gherardini, and F. Caruso, arXiv:2101.03221 (2021).
  * Miles  _et al._ (2020) C. Miles, A. Bohrdt, R. Wu, C. Chiu, M. Xu, G. Ji, M. Greiner, K. Q. Weinberger, E. Demler, and E.-A. Kim, arXiv:2011.03474 (2020).
  * Bloch  _et al._ (2012) I. Bloch, J. Dalibard, and S. Nascimbène, Nat. Phys. 8, 267 (2012).
  * Kim  _et al._ (2018) H. Kim, Y. Park, K. Kim, H.-S. Sim, and J. Ahn, Phys. Rev. Lett. 120, 180502 (2018).
  * Ebadi  _et al._ (2020) S. Ebadi, T. T. Wang, H. Levine, A. Keesling, G. Semeghini, A. Omran, D. Bluvstein, R. Samajdar, H. Pichler, W. W. Ho, _et al._ , arXiv:2012.12281 (2020).
  * Browaeys and Lahaye (2020) A. Browaeys and T. Lahaye, Nat. Phys. 16, 132–142 (2020).
  * Bloch  _et al._ (2008) I. Bloch, J. Dalibard, and W. Zwerger, Rev. Mod. Phys. 80, 885 (2008).
  * Islam  _et al._ (2011) R. Islam, E. Edwards, K. Kim, S. Korenblit, C. Noh, H. Carmichael, G.-D. Lin, L.-M. Duan, C.-C. J. Wang, J. Freericks, _et al._ , Nat. Comm. 2, 1 (2011).
  * Bohnet  _et al._ (2016) J. G. Bohnet, B. C. Sawyer, J. W. Britton, M. L. Wall, A. M. Rey, M. Foss-Feig, and J. J. Bollinger, Science 352, 1297 (2016).
  * Joshi  _et al._ (2020) M. K. Joshi, A. Elben, B. Vermersch, T. Brydges, C. Maier, P. Zoller, R. Blatt, and C. F. Roos, Phys. Rev. Lett. 124, 240505 (2020).
  * Sun and Robicheaux (2008) B. Sun and F. Robicheaux, New J. Phys. 10, 045032 (2008).
  * Ates  _et al._ (2012) C. Ates, J. P. Garrahan, and I. Lesanovsky, Phys. Rev. Lett. 108, 110603 (2012).
  * Turner  _et al._ (2018) C. Turner, A. Michailidis, D. Abanin, M. Serbyn, and Z. Papić, Nat. Phys. 14, 745 (2018).
  * Emary and Brandes (2003) C. Emary and T. Brandes, Phys. Rev. E 67, 066203 (2003).
  * Dusuel and Vidal (2004) S. Dusuel and J. Vidal, Phys. Rev. Lett. 93, 237204 (2004).
  * Note (1) See supplemental material for details.
  * Benatti  _et al._ (2016) F. Benatti, F. Carollo, R. Floreanini, and H. Narnhofer, Phys. Lett. A 380, 381 (2016).
  * Breuer  _et al._ (2009) H.-P. Breuer, E.-M. Laine, and J. Piilo, Phys. Rev. Lett. 103, 210401 (2009).
  * Laine  _et al._ (2010) E.-M. Laine, J. Piilo, and H.-P. Breuer, Phys. Rev. A 81, 062115 (2010).
  * Breuer (2012) H.-P. Breuer, J. Phys. B 45, 154001 (2012).
  * Rivas  _et al._ (2014) Á. Rivas, S. F. Huelga, and M. B. Plenio, Rep. Prog. Phys. 77, 094001 (2014).
  * Hochreiter and Schmidhuber (1997) S. Hochreiter and J. Schmidhuber, Neural Comput. 9, 1735 (1997).
  * Vaswani  _et al._ (2017) A. Vaswani, N. Shazeer, N. Parmar, J. Uszkoreit, L. Jones, A. N. Gomez, L. Kaiser, and I. Polosukhin, arXiv:1706.03762 (2017).
  * Sahoo  _et al._ (2018) S. S. Sahoo, C. H. Lampert, and G. Martius, in  _Proc. 35th International Conference on Machine Learning, ICML 2018, Stockholm, Sweden, 2018_, Vol. 80 (PMLR, 2018) pp. 4442–4450.
  * Chen  _et al._ (2018) R. T. Q. Chen, Y. Rubanova, J. Bettencourt, and D. K. Duvenaud, in  _Advances in Neural Information Processing Systems_, Vol. 31, edited by S. Bengio, H. Wallach, H. Larochelle, K. Grauman, N. Cesa-Bianchi, and R. Garnett (Curran Associates, Inc., 2018) pp. 6571–6583.

Supplemental material: Machine learning time-local generators of open quantum
dynamics

In this supplemental material we report a more detailed analysis concerning
the effectiveness of the hypermodel in predicting the dynamics of the Bloch
vector within the training time-window. In addition, we provide the hyper-
parameters of the used models.

##  I Hypermodel vs. exact dynamics

In the main text we analyze the time-dependence of time-local generators
parametrizing the “propagator” $M[F(t;\theta^{\ast})]=\mathds{1}+dtL_{t}$ for
the Bloch vector, using a hypermodel architecture. The matrix $L_{t}$ encodes
all the information concerning the dynamics of local (one-site) observables.
As stated in the main text, a suitable measure to quantify the dependence on
time of this matrix is given by $\Xi=\frac{1}{T}\int_{0}^{T}dt\,\xi(t)$, where
$\xi(t)=\sqrt{{\rm Tr}(\dot{L}_{t}^{\dagger}\dot{L}_{t})}$. However, this can
be a faithful measure of the time-dependence of the physical generator, only
if $L_{t}$ correctly implements the dynamics of the Bloch vector. It is thus
important to check that the prediction of the hypermodel matches the exact
dynamics, in the time window $[0,T]$ in which we want to analyze the generator
$L_{t}$. To this end, we report in Fig. 1 a detailed comparison between the
predicted time evolution of the Bloch vector components and the exact ones for
six different parameter regimes.

Figure 1: Hypermodel accuracy: Comparison between the exact dynamics and the
dynamics predicted by the hypermodel for initial state $\rho_{0}$ reported in
Eq. (1) and for $N=7$. As it can be observed the agreement between the model
prediction (dotted line) and the exact dynamics (full line) is excellent. This
fact enables us to use the local operator learned using the hypermodel for
studying the time-dependence of $L_{t}$.

We start from a generic initial state of the form

$\rho_{0}=\frac{1}{2}\left(\begin{array}[]{@{}cc@{}}1+z&x-iy\\\
x+iy&1-z\end{array}\right)\otimes\rho_{B},$ (1)

with $x^{2}+y^{2}+z^{2}<1$, in particular we choose $x=0.6$, $y=0.3$ and
$z=0.4$, and we consider different interactions and initial bath temperatures.
In the training time-window $[0,T]$, we observe an excellent agreement between
the two curves for $T=10/\Omega$. In Fig. 2 we report the time evolution of
the quantity $\xi(t)$ for the same parameter regimes.

Figure 2: Time-dependence: Behavior of the quantity $\xi(t)$ for the same
parameter regimes of Fig. 1. The quantity $\Xi$ grows with the interaction,
meaning that the generator of the local dynamics $L_{t}$ becomes more time-
dependent.

As it can be observed for weakly interacting systems the generator $L_{t}$ is
almost time-independent. In this regime, as discussed in the main text, it is
possible to predict the local dynamics using simple time-independent
architectures. Increasing the interaction, the parameter $\Xi$ grows
accordingly signaling that the dependence on time of $L_{t}$ becomes stronger.
The training procedure for the hyper-model is slightly different with respect
what we have done for the time-independent case, in particular it is different
the differentiation between training and evaluation set. In this case, in
fact, we generate two different types of data-sets, the first data-set for
training consists in $100$ trajectories with randomly chosen initial state
$\rho_{0}$ of the the same form reported in Eq. (1). For evaluation we select,
instead, randomly data-points from trajectories generated from initial states
$\rho_{0}^{\mathrm{ev}}$ of type

$\rho_{0}^{\mathrm{ev}}=\left(\begin{array}[]{@{}cc@{}}1-c&0\\\
0&c\end{array}\right)\otimes\rho_{B},$ (2)

with $0.01<c<0.7$. In this way we are performing a separation of the data-sets
in the Hilbert space. This is contrary to what we had for the time-independent
case where the data-set was given by separated time-intervals of the same
trajectory.

##  II Network architectures

Setting |  Architecture |  Optimizer |  Dataset
---|---|---|---
Time-independent |  Linear perceptron (no hidden layers) |  Adam Learning rate: $10^{-3}$ Betas: $0.9$, $0.999$ Epsilon: $10^{-8}$ Batch size: $256$ Batches per epoch: $512$ Epochs: $5$ |  Trajectories (train / val): $100$ / $20$ Total time ($\Omega\,T$): $10$ Time step ($\Omega\,dt$): $0.01$
Time-dependent (hyper-model) |  Hidden layers: $3$ Nonlinearity: tanh Output nonlinearity: none |  Adam Learning rate: $10^{-3}$ Betas: $0.9$, $0.999$ Epsilon: $10^{-8}$ Batch size: $256$ Batches per epoch: $256$ Epochs: $500$ |  Trajectories (train / val): $100$ / $20$ Total time ($\Omega\,T$): $10$ Time step ($\Omega\,dt$): $0.01$
Table 1: Overview over the architecture, optimizers, training- and system-
parameters for Model I and Model II (top) and the hyper-model (bottom)