Learning dynamical systems: an example from open
quantum system dynamics.

Pietro Novelli
Istituto Italiano di Tecnologia
pietro.novelli@iit.it

Abstract
Machine learning algorithms designed to learn dynamical systems from data canbe used to forecast, control and interpret the observed dynamics. In this work weexemplify the use of one of such algorithms, namely Koopman operator learning,in the context of open quantum system dynamics. We will study the dynamics of asmall spin chain coupled with dephasing gates and show how Koopman operator
learning is an approach to efficiently learn not only the evolution of the density
matrix, but also of every physical observable associated to the system. Finally,leveraging the spectral decomposition of the learned Koopman operator, we show
how symmetries obeyed by the underlying dynamics can be inferred directly from
data.

1

A primer on Koopman operator learning

The Koopman operator is a mathematical object concerned with the evolution of a dynamical system.
Formally it is defined as the linear operator evolving any observable1 in a pre-specified set. It owes
its name to the works of Bernard Koopman on the evolution of Hamiltonian systems [Koo31], butcurrent applications of the Koopman operator embrace the whole field of dynamical systems [BMM12,
MMS20, BBKK22]. We now provide a minimalistic introduction to the Koopman operator, referring
to [MMS20, KNM+ 22] for a mathematically rigorous presentation.
Given a discrete dynamical system over a state space X , the Koopman operator U is defined as(U f ) (xt ) := f (xt+1 ) for all t ∈ N,
for every observable f : X → R in a suitable set. The Koopman operator, therefore, acts by evolving
the observables of the system by one step of the dynamics. The above definition is appropriate
for deterministic dynamical systems, but can be consistently extended to stochastic systems via
conditional expectation operators [MT93]. In this work we hinge on the formulation of Koopman
operator learning in [KNM+ 22], in which the space of observables is restricted to a reproducing
kernel Hilbert space H [Aro50, SC08]. Within this setting, an estimator Û : H → H of the Koopman
operator can be learned from a single dataset of observations (xi , yi )ni=1 consisting of sampled statesxi ∈ X along with their one-step-forward evolutions yi ∈ X .

An obvious use of the learned operator Û is forecasting both future states and every observable in H
(see Remark 1). The linearity of Û , however, also allows to compute its spectral decomposition, andif (λi , ξi , ψi )ri=1 are eigen-triplets2 of Û , the evolution of an observable f : X → R can be dissected
as a sum of modesr
X(Û f )(xt ) =λti γif ψi (x0 ),
(1)
i=1
1
2

That is, a scalar function of the state of the dynamical system.
Namely: eigenvalues, left and right eigenfunctions.

Machine Learning and the Physical Sciences workshop, NeurIPS 2022.

where γif is a linear functional of the observable f linked to the i-th left eigenfunction ξi . The
Koopman mode decomposition [BMM12, AM17] can be used to infer properies of the dynamicalsystem such as constant of motions, correlation and dephasing between different observables, mixing
times, normal modes, metastable states and many others. These quantities are of paramount importance in the physical sciences and machine-learning-enabled tools to compute the Koopman mode
decomposition unlock the possibility to study complex physical systems too hard to be approachedwith analytical theories alone.
Remark 1 (Forecasting the observables). Suppose we have observed T states (xt )t≤T of a (possibly
noisy) dynamical system and we are interested in learning the evolution of its observables f : X → R.
A straightforward approach is to first learn an map Ŝ between the states such that xt+1 ≈ Ŝ(xt ) andthen predict any observable by simply composing it with the learned map Ŝ, i.e. f (xt+1 ) ≈ f (Ŝ(xt )).
Here, Ŝ : X → X can be anything: a neural network, a vector-valued kernel etc. The framework
developed in [KNM+ 22] for the Koopman operator learning, on the contrary, actively capitalizeon the structure of the space of observables H to directly provide the best forecast as f (xt+1 ) ≈(Û f )(xt ) for any observable f ∈ H, without the middle step of learning the dynamical map Ŝ.
The two approaches coincide if the underlying dynamics is deterministic, whereas for stochastic
dynamical systems the first approach is sub-optimal (i.e. leading to larger errors).

2

Quantum dynamics and the Lindblad equation

In the context of open quantum system dynamics, the state of a quantum mechanical system is
described by a density matrix (or operator) ρ : H → H, where H is a Hilbert space characterizing
the quantum system [NC12, Chapter 2]. Physical requirements prescribe that the density matrix
is a positive definite, hermitian operator satisfying tr(ρ) = 1. Any physical observable A is anhermitian operator acting on H, and its quantum mechanical expectation value for a given state ρ ishAiρ := tr(Aρ). Noticing that tr(Aρ) is the scalar product between A and ρ in the Hilbert-Schmidt
sense, we conclude that the observables of physical interest are linear functionals of ρ ∈ HS(H),where HS(H) is the space of Hilbert-Schmidt operators on H.
For quantum systems with discrete energy levels, the Hilbert space H is isomorphic to Cd for someN(possibly large) d. In the case of N two-level systems such as spins, one has H ∼ C2 .

If the physical system is closed, i.e. isolated from any interaction with the environment, its evolution isunitary and governed by the Schrödinger equation. Conversely, when the system is open its equationof motion is also affected by the environment. Under proper assumptions and approximations3 , thedynamics of an open quantum system is described by the Lindblad equation [Lin76, WM09, Chapter
3]"#X
ρ(t)L†k Lk + L†k Lk ρ(t)di†
ρ(t) = − [H, ρ(t)] +Lk ρ(t)Lk −
,
dt
~
2
k

where H is the Hamiltonian of the system and Lk : H → H are collapse operators describing
disspiative processes. The Lindblad equation is a linear ordinary differential equation for the density
matrix ρ. Its form guarantees that for all t, ρ(t) is a density matrix (i.e. satisfy the aforementioned
requirements of positivity, hermiticity and unit trace).
Remark 2 (Relation with Heisenberg and Schrödinger pictures). It is well known that the evolutionof quantum systems can be described by two equivalent pictures. In the Schrödinger picture, which
we employed above, the states evolve in time, while the operators are fixed. In the Heisenberg picture
the opposite is true: states are fixed and operators evolve. Clearly, because of theirs equivalence,the choice of the picture do not affect the evolution of the expected value hAiρ (t) = tr(AH (t)ρ) =tr(AρS (t)). The Koopman operator acts directly on the expected value hAiρ (t) and is therefore welldefined in both pictures. Koopman operator learning, however, is based on a dataset of evolved statesand is arguably more natural to adopt the Schrödinger picture to study its properties.
3

Namely: separability, Markovianity, Born and secular approximations.

2

3

Koopman operator learning on a spin chain

We now illustrate Koopman operator learning on a simple spin chain model by generating the training
dataset from a simulated Lindblad evolution. This synthetic setting is mostly useful as a benchmark,as we can compare the learned dynamical system with the ground truth. A more sensible use of
Koopman operator learning that we envision, however, is to make use of experimental measurementson complex quantum systems (e.g. quantum computers) to distill an accurate and interpretable modelof theirs dynamical evolution.
We consider a spin chain with N spins described by the Hamiltonian
N −1

1 X
y 
x
z
H := −
Jk σix σi+1
+ σiy σi+1
+ J⊥ σiz σi+1
,
2 i=1

(2)

where σix,y,z are Pauli matrices corresponding to the i-th spin. The model (2) is also known as
quantum Heisenberg model, and under proper assumptions on the coupling constants Jk , J⊥ can beexactly solved using the Bethe ansatz [KMGT97, KHM98].
We let the spin chain (2) interact with the environment via dephasing channels [MP08]. p
To this end,z
we choosethecollapse
operatorsoftheLindbladequationasN
puredephasingsL=σΓ/2 withiippzrate Γ/2. As explained in [MP08], a pure dephasing σ Γ/2 for a two level system is responsible
for the decay of the off-diagonal elements of the density matrix at a rate ∝ e−Γt .
Experimental details. The data are generated by simulating a chain of N = 5 spins with the
QuTiP framework [JNN12]. We construct a dataset of quantum states by saving the density matrixof the simulated system ρ(t) at times t = k∆t with k = 1, . . . , 200 and ∆t := 0.5. The first half ofthe data, corresponding to the initial 100 timesteps, is used to train the Koopman operator estimatorÛ , while the second half is used for testing. The initial state ρ(0) for the integration of the Lindblad
equation is chosen as the pure state | ↓↑↑↑↑i, corresponding to the first spin down and the rest up. The
coupling constants and dephasing rate are set to Jk := 0.1π, J⊥ := 0.2π, and Γ = 0.2, respectively.
The estimator Û of the Koopman operator is computed using the reduced rank regression algorithm
presented in [KNM+ 22] with rank r = 19, Tikhonov regularization λ = 10−6 and a linear kernel.
The reduced rank regression algorithm with linear kernel on the dataset (xi , yi )ni=1 returns theminimizern1Xkyi − U ∗ xi k2 + λkU k2HS ,Û := arg min
U ∈HSr (X ) n i=1
where HSr (X ) is the set of Hilbert-Schmidt operators on X with rank smaller or equal than r andU ∗ is the hermitian conjugate of U .
The code to reproduce the experiments is available at the link https://github.com/CSML-IITUCL/kooplearn. The experiments have been conducted on a workstation equipped with an Intel(R)CoreTM i9-9900X CPU @ 3.50GHz, 48GB of RAM and a NVIDIA GeForce RTX 2080 Ti GPU.
3.1

Forecasting the physical observables

We now compare the dynamics of selected physical observables as forecasted by the Koopman
operator estimator Û against the ground truth obtained by direct integration of the Lindblad equation. We focus on the spin polarization and on the spin current. The spin polarization (along z) of
the i-th spin is the observable associated with the spin operator Siz := ~σiz /2. A non-vanishingspin polarization tr(Siz ρ) across many sites signals the onset of a magnetic order for the state ρ.
The Heisenberg model (2), indeed, is a rudimentary model of a magnetic material and depending on the coupling constants Jk and J⊥ might exhibit either ferromagnetic or antiferromagnetic
order [KMGT97, KHM98].
The spin current at site i can be derived from the continuity equation as detailed in [SKK04] andreduces to the operator
jiz :=

iJk  + −++−(σi σi+1 − σi− σi+1
) − (σi−1
σi− − σi−1
σi+ ) ,
4
3

Average spin current [Jk ]Average spin polarization [h̄]

0.4
0.2
0.0
hS1z i

−0.2
−0.4
0

25

50
Time

75

0.1
0.0

−0.1

hS5z iTraining

0.2

−0.2

100

h j3z iTraining
0

25

50
Time

75

100

Figure 1: Left panel: forecast of the average spin polarization of the first (1) and last (5) spin of thechain in units of ~. Right panel: forecast of the average spin current of the central (3) spin of thechain in units of Jk . In both panels the gray traces correspond to the ground truth and the first half ofthe data were used to train the Koopman operator estimator.

where σi± := σix ± iσiy are the usual ladder operators for the i-th spin. The first term in roundparenthesis in the definition of jiz accounts for outgoing currents while the second term for incoming
currents at site i. Manipulation and control of spin currents is of great interest in the field of
spintronics, dealing with the use electrons’ spin degrees of freedom for information processing. In
spintronic devices, data storage and transmission can be achieved at higher efficiency compared tostandard electronic devices [HYN+ 20].
In Figure 1 we show the forecast of, respectively, the average spin polarization and the average spincurrent of selected sites along the chain. In both panels, the left portion of the plot corresponds tothe data used to learn the Koopman operator estimator Û , while in the right portion of the plot arecompared the forecasted dynamics on unseen data against the true dynamics. The forecast is accurate
both for spin polarization and current at short times (i.e. for 50 ≤ t / 70). At later times, however,the forecasted spin polarization diverge from the ground truth, while the forecasted spin current stays
remarkably close to the true spin current.
3.2

Decay rates, normal modes, constants of motion and symmetries

We now highlight how the Koopman mode decomposition (1) might be used to infer quantities ofphysical interest such as decay rates, frequency of the normal modes and constants of motion of the
dynamics.
Decay rates and frequency of the normal modes are easily obtained noticing that the time-dependencein (1) is expressed only by the eigenvalues λi of Û . The theory of dynamical systems dictatesthat |λi | ≤ 1 for stable systems. The powered eigenvalues λti appearing in (1) can be conveniently expressed in polar form, yielding λti := e−t/τi ei2πωi t , with τi = −t/ log(|λi |) andωi = arg(λi )/(2πt).
Computing the decay rates of the learned Û we get that 16 out of 19 eigenvalues have decay rates
in the interval [0.018, 0.021], which is quite close to the true rate Γ = 0.02 selected for the puredephasing gates. The remaining three decay rates are, respectively 0.014, 0.026 and 0.003, the last
one of which corresponds to the steady-state solution. Similarly, the computed frequencies ωi lie in
the range [0.05, 0.2], of the same order of the (properly normalized) coupling constants Jk /2π = 0.05
and J z /2π = 0.1.
The last point we would like to comment concerns the mode corresponding to the steady-stateof the system, i.e. the one with eigenvalue one. This mode, in practice, is characterized by an
estimated Koopman eigenvalue very close to 1 and in our calculation is given by λ1 = 0.9985. By
construction, the right eigenfunction associated to the steady-state mode is a constant of motion.
Indeed, ψ1 (ρt ) = U tψ1 (ρ0 ) = λt1 ψ1 (ρ0 ) = ψ1 (ρ0 ) for all t.
4

In the present case we are dealing with linear observables4 , and every eigenfunction ψi (ρ) can be
expressed, by construction, as ψi (ρ) = tr(Ψi ρ) for some operator Ψi . For the steady-state mode,the operator Ψ1 is a symmetry transformation for which the dynamics is invariant. We verified thatthe estimated Ψ1 commutes with the total spin operator along z, as expected for the pure dephasing
channel [MP08].

Broader impact
The abstract focus on two areas of active research: machine learning and quantum technologies.
The techniques discussed in this abstract do not pose any direct ethical concern. Machine learning,
here, is used only as a tool to aid and foster scientific research in the field of quantum technologies.
The results presented might impact (hopefully in a positive way) the cohort of highly specialized
researchers working on quantum technologies, and are very unlikely to impact the society at large.

Checklist
1. (a) Do the main claims made in the abstract and introduction accurately reflect the paper’s
contributions and scope? [Yes](b) Did you describe the limitations of your work? [No](c) Did you discuss any potential negative societal impacts of your work? [N/A](d) Have you read the ethics review guidelines and ensured that your paper conforms tothem? [Yes]2. (a) Did you state the full set of assumptions of all theoretical results? [Yes](b) Did you include complete proofs of all theoretical results? [No]3. (a) Did you include the code, data, and instructions needed to reproduce the main experimental results (either in the supplemental material or as a URL)? [Yes](b) Did you specify all the training details (e.g., data splits, hyperparameters, how they
were chosen)? [Yes](c) Did you report error bars (e.g., with respect to the random seed after running experiments multiple times)? [N/A](d) Did you include the total amount of compute and the type of resources used (e.g., typeof GPUs, internal cluster, or cloud provider)? [Yes]4. (a) If your work uses existing assets, did you cite the creators? [N/A](b) Did you mention the license of the assets? [N/A](c) Did you include any new assets either in the supplemental material or as a URL? [N/A](d) Did you discuss whether and how consent was obtained from people whose data you’reusing/curating? [N/A](e) Did you discuss whether the data you are using/curating contains personally identifiable
information or offensive content? [N/A]
5. (a) Did you include the full text of instructions given to participants and screenshots, if
applicable? [N/A](b) Did you describe any potential participant risks, with links to Institutional Review
Board (IRB) approvals, if applicable? [N/A](c) Did you include the estimated hourly wage paid to participants and the total amountspent on participant compensation? [N/A]

References[AM17] Hassan Arbabi and Igor Mezić. Ergodic theory, dynamic mode decomposition, andcomputation of spectral properties of the koopman operator. SIAM Journal on Applied
Dynamical Systems, 16(4):2096–2126, 2017.
4That is, linear functionals of the density matrix.

5

[Aro50] Nachman Aronszajn. Theory of reproducing kernels. Transactions of the Americanmathematical society, 68(3):337–404, 1950.
[BBKK22] Steven L. Brunton, Marko Budišić, Eurika Kaiser, and J. Nathan Kutz. Modern Koopman
Theory for Dynamical Systems. SIAM Review, 64(2):229–340, May 2022.
[BMM12] Marko Budišić, Ryan Mohr, and Igor Mezić. Applied Koopmanism. Chaos: An
Interdisciplinary Journal of Nonlinear Science, 22(4):047510, December 2012.
[HYN+ 20] Atsufumi Hirohata, Keisuke Yamada, Yoshinobu Nakatani, Ioan-Lucian Prejbeanu,
Bernard Diény, Philipp Pirro, and Burkard Hillebrands. Review on spintronics:
Principles and device applications. Journal of Magnetism and Magnetic Materials,
509:166711, September 2020.
[JNN12] J.R. Johansson, P.D. Nation, and Franco Nori. QuTiP: An open-source python framework for the dynamics of open quantum systems. Computer Physics Communications,
183(8):1760–1772, August 2012.
[KHM98] Michael Karbach, Kun Hu, and Gerhard Muüller. Introduction to the bethe ansatz II.
Computers in Physics, 12(6):565, 1998.
[KMGT97] Michael Karbach, Gerhard Muüller, Harvey Gould, and Jan Tobochnik. Introduction to
the bethe ansatz i. Computers in Physics, 11(1):36, 1997.
[KNM+ 22] Vladimir Kostic, Pietro Novelli, Andreas Maurer, Carlo Ciliberto, Lorenzo Rosasco,and Massimiliano Pontil. Learning dynamical systems via koopman operator regression
in reproducing kernel hilbert spaces, 2022.
[Koo31] B. O. Koopman. Hamiltonian systems and transformation in hilbert space. Proceedingsof the National Academy of Sciences, 17(5):315–318, May 1931.
[Lin76] G. Lindblad. On the generators of quantum dynamical semigroups. Communications inMathematical Physics, 48(2):119–130, June 1976.
[MMS20] Alexandre Mauroy, Igor Mezić, and Yoshihiko Susuki, editors. The Koopman Operator
in Systems and Control. Springer International Publishing, 2020.
[MP08] Florian Marquardt and Annett Püttmann. Introduction to dissipation and decoherence inquantum systems, 2008.
[MT93] Sean P. Meyn and Richard L. Tweedie. Markov Chains and Stochastic Stability. Communications and Control Engineering. Springer London, 1993.
[NC12] Michael A. Nielsen and Isaac L. Chuang. Quantum Computation and Quantum Information. Cambridge University Press, June 2012.
[SC08] Ingo Steinwart and Andreas Christmann. Support Vector Machines. Springer New York,2008.
[SKK04] F. Schütz, P. Kopietz, and M. Kollar. What are spin currents in heisenberg magnets?
The European Physical Journal B, 41(4):557–560, October 2004.
[WM09] Howard M. Wiseman and Gerard J. Milburn. Quantum Measurement and Control.
Cambridge University Press, November 2009.

6