# Asymptotics and typicality of sequential generalized measurements

Wen-Long Ma<EMAIL_ADDRESS>State Key Laboratory of Superlattices and
Microstructures, Institute of Semiconductors, Chinese Academy of Sciences,
Beijing 100083, China  College of Materials Science and Opto-Electronic
Technology, University of Chinese Academy of Sciences, Beijing 100049, China
Shu-Shen Li  State Key Laboratory of Superlattices and Microstructures,
Institute of Semiconductors, Chinese Academy of Sciences, Beijing 100083,
China  College of Materials Science and Opto-Electronic Technology, University
of Chinese Academy of Sciences, Beijing 100049, China  Ren-Bao Liu  Department
of Physics, Centre for Quantum Coherence, and The Hong Kong Institute of
Quantum Information Science and Technology, The Chinese University of Hong
Kong, Shatin, New Territories, Hong Kong, China

###### Abstract

The relation between projective measurements and generalized quantum
measurements is a fundamental problem in quantum physics, and clarifying this
issue is also important to quantum technologies. While it has been intuitively
known that projective measurements can be constructed from sequential
generalized or weak measurements, there is still lack of a proof of this
hypothesis in general cases. Here we rigorously prove it from the perspective
of quantum channels. We show that projective measurements naturally arise from
sequential generalized measurements in the asymptotic limit. Specifically, a
selective projective measurement arises from a set of typical sequences of
sequential generalized measurements. We provide an explicit scheme to
construct projective measurements of a quantum system with sequential
generalized quantum measurements. Remarkably, a single ancilla qubit is
sufficient to mediate a sequential weak measurement for constructing arbitrary
projective measurements of a generic system.

Quantum measurements retrieve classical information from quantum states
Wiseman2010 ; Jacobs2014 , and are particularly important to quantum
technologies Nielsen2010 . The traditional description of measurement in
quantum mechanics is through projective measurements (PMs) of observables
represented by Hermitian operators Peres2006 . Measuring an observable
corresponds to statistically projecting the quantum state to one of the
orthogonal eigenspaces of this observable. PMs appear most commonly in quantum
foundation and quantum information theory, and are widely useful for
initialization and readout of quantum systems in quantum technologies
Eizerman2004 ; Vamivakas2010 ; Neumann2010 ; Morello2010 ; Jiang2009 ;
Nakajima2017 ; West2019 .

Recent quantum theories show that there exist more general quantum
measurements, called generalized measurements or positive-operator-valued
measures (POVMs) Kraus1983 ; Andersson2008 ; Chen2018 ; Chen2019 ; Cheong2012
. POVMs can outperform PMs in many tasks in quantum technologies, such as
quantum tomography Renes2004  and quantum state discrimination or estimation
Bergou2010 ; Derka2001 . Moreover, continuous or sequential POVMs can be
exploited for monitoring and maneuvering quantum evolutions Jacobs2006 ;
Gurvitz1997 ; Ashhab2009 ; Korotkov2001 ; Blok2006 ; Jordan2006 ;
Chantasri2013 ; Chantasri2015 ; Presilla1996 ; Diosi2016 . In particular, weak
measurements can extract partial information without projections, and
therefore can help realize optimal qubit tomography Shojaee2018 , reconcile
measurement incompatibility Monroe2021  and extract arbitrary bath
correlations Wang2019 ; Pfender2019 ; Wu2022 .

Substantial efforts have been devoted to illustrating the relation between PMs
and POVMs. A celebrated result is Naimark’s theorem Peres2006 , stating that
any POVM can be implemented as a PM on an enlarged Hilbert space. POVMs can
also be simulated by PMs with classical randomness or postselection
Oszmaniec2017 ; Oszmaniec2019 . In the opposite direction, it has been argued
that sequential weak measurements can generate PMs by analysing the gradual
state collapse Brun2002 ; Lidar2013 ; Oreshkov2005  or the statistics of
measurement results Ma2018 ; Rao2019 ; Liu2017 . However, to our knowledge,
the general relation between PMs and sequential POVMs still remains elusive.

In this paper, we rigorously prove that PMs emerge from sequential POVMs in
the asymptotic limit, when the measurement operators are normal operators and
commuting with each other. The proof is based on the observation that
projective subspaces are fixed points of the quantum channels for such POVMs.
Moreover, from the theory of classical typicality, we find that different
selective PMs arise from different sets of typical sequences of sequential
selective POVMs. These results completely characterize the structures of
sequential commuting POVMs with normal measurement operators. We further
present a general scheme to realize such POVMs with a single qubit ancilla,
and show that sequential POVMs can be used to projectively measure an
arbitrary finite-dimensional quantum system. An example of the scheme is also
provided to measure the modular excitation numbers of an infinite-dimensional
bosonic mode with an ancilla qubit.

Figure 1: (a) Schematic of sequential POVMs in the asymptotic limit. (b)
Emergent PMs as typical sets of the whole set of superoperators for asymptotic
applications of POVMs. (c) The emergent projections in the operator space of
the quantum system.

POVMs and quantum channels. Consider a $d$-level quantum system under non-
selective POVMs with the measurement operators
$\\{M_{\alpha}\\}_{\alpha=1}^{r}$. Its evolution is described by the
completely positive and trace-preserving (CPTP) quantum maps in the Kraus
representation Kraus1983 ,

$\displaystyle\Phi(\rho)=\sum_{\alpha=1}^{r}\mathcal{M}_{\alpha}\rho=\sum_{\alpha=1}^{r}M_{\alpha}\rho
M_{\alpha}^{\dagger},$ (1)

where $\mathcal{M}_{\alpha}=M_{\alpha}(\cdot)M_{\alpha}^{\dagger}$ is a
superoperator in the Hilbert space representing a selective POVM with the
$\alpha$th outcome, and
$\sum_{\alpha=1}^{r}M_{\alpha}^{\dagger}M_{\alpha}=\mathbb{I}$ being the
identity operator ensures the trace-preserving property of such maps. The CPTP
maps are also called quantum operations or quantum channels Caruso2014 ;
Wolf2010 .

Quantum channels have convenient matrix representations in the Hilbert-Schmidt
(HS) space of the quantum system Bengtsson2017 ; SI . While the density
matrices are operators in the Hilbert space with an orthonormal basis
$\\{|i\rangle\\}_{i=1}^{d}$, they are turned into vectors in the HS space,
i.e., $\rho=\sum_{i,j=1}^{d}\rho_{ij}|i\rangle\langle
j|\leftrightarrow|\rho\rangle\rangle=\sum_{i,j=1}^{d}\rho_{ij}|ij\rangle\rangle$,
such that $X\rho Y\leftrightarrow X\otimes Y^{T}|\rho\rangle\rangle$ with $X$,
$Y$ being operators in the Hilbert space and $Y^{T}$ being the transpose of
$Y$. The inner product in the HS space is defined as
$\langle\langle\sigma|\rho\rangle\rangle={\rm Tr}[\sigma^{\dagger}\rho]$. Then
non-selective POVMs act as the following linear operators in the HS space,

$\displaystyle\hat{\Phi}|\rho\rangle\rangle=\sum_{\alpha=1}^{r}\hat{\mathcal{M}}_{\alpha}|\rho\rangle\rangle=\sum_{\alpha=1}^{r}M_{\alpha}\otimes
M_{\alpha}^{*}|\rho\rangle\rangle,$ (2)

where $M_{\alpha}^{*}$ is the complex conjugate of $M_{\alpha}$. Note that we
add hats for operators in the HS space, to distinguish them from the
corresponding superoperators in the Hilbert space. With the HS space, the
probability to get the $\alpha$th outcome is
$\langle\langle\mathbb{I}|\hat{\mathcal{M}}_{\alpha}|\rho\rangle\rangle={\rm
Tr}(M_{\alpha}\rho M_{\alpha}^{\dagger})$.

POVMs with normal and commuting measurement operators. We assume that the set
of measurement operators $\\{M_{\alpha}\\}_{\alpha=1}^{r}$ are all normal
operators and commuting with each other, i.e.,
$[M_{\alpha},M_{\alpha}^{\dagger}]=[M_{\alpha},M_{\beta}]=0$ for all integers
$\alpha,\beta\in[1,r]$, such that $\\{M_{\alpha}\\}_{\alpha=1}^{r}$ can be
simultaneously diagonalized in an orthonormal eigenbasis of the quantum system
Garcia2017 ,

$\displaystyle\left[{\begin{array}[]{*{20}{c}}{{M_{1}}}\\\ \vdots\\\
{{M_{r}}}\end{array}}\right]=\left[{\begin{array}[]{*{20}{c}}{{c_{11}}}&\cdots&{{c_{1d}}}\\\
\vdots&\vdots&\vdots\\\
{{c_{r1}}}&\cdots&{{c_{rd}}}\end{array}}\right]\left[{\begin{array}[]{*{20}{c}}{|1\rangle\langle
1|}\\\ \vdots\\\ {|d\rangle\langle d|}\end{array}}\right].$ (12)

This can be simply denoted as $\mathbf{M}=\mathbf{C}\mathbf{P}$, where
$\mathbf{M}=[M_{1},\cdots,M_{r}]^{T}$, $\mathbf{P}=[|1\rangle\langle
1|,\cdots,|d\rangle\langle d|]^{T}$, and $\mathbf{C}$ is a $r\times d$ complex
matrix ($r$ and $d$ are generally different). We partition $\mathbf{C}$
according to its columns as $[\mathbf{c}_{1},\cdots,\mathbf{c}_{d}]$, then
$\|\mathbf{c}_{j}\|^{2}=\mathbf{c}_{j}^{\dagger}\mathbf{c}_{j}=1$ for any
integer $j\in[1,d]$ due to
$\mathbf{M}^{\dagger}\mathbf{M}=\sum_{\alpha=1}^{r}M_{\alpha}^{\dagger}M_{\alpha}=\mathbb{I}$,
and $\\{\mathbf{c}_{j}\\}_{j=1}^{d}$ is a set of unit vectors in a
$r$-dimensional complex vector space, with $j$ corresponding to the basis
state $|j\rangle$. Note that these unit vectors are not necessarily orthogonal
to each other SI . For a specific POVM, the measurement operators are not
unique, since we can define a new set of measurement operators by
$\mathbf{M}^{\prime}=\mathbf{T}\mathbf{M}$ with $\mathbf{T}$ being a $r\times
r$ unitary matrix, which satisfy
$\mathbf{M}^{\prime\dagger}\mathbf{M}^{\prime}=\mathbb{I}$ and also
characterize the same quantum channel.

The non-selective POVMs above act as diagonal operators in the HS space,

$\displaystyle\hat{\Phi}=\sum_{i,j=1}^{d}\mathbf{c}_{j}^{\dagger}\mathbf{c}_{i}|ij\rangle\rangle\langle\langle
ij|,$ (13)

where $\\{|ij\rangle\rangle\\}_{i,j=1}^{d}$ are the eigenvectors
(eigenmatrices in the Hilbert space) of $\hat{\Phi}$ with the corresponding
eigenvalues $\\{\mathbf{c}_{j}^{\dagger}\mathbf{c}_{i}\\}_{i,j=1}^{d}$. Since
$|\mathbf{c}_{j}^{\dagger}\mathbf{c}_{i}|\leq 1$ (due to the Cauthy-Schwarz
inequality) with equality if and only if
$\mathbf{c}_{i}=e^{i\varphi}\mathbf{c}_{j}$ for some real $\varphi$, all the
eigenvalues of $\hat{\Phi}$ lie within the unit disk of the complex plane. The
eigenvectors with eigenvalue 1 are called fixed points Wolf2010 , and those
with eigenvalues $e^{i\varphi}$ with $\varphi\neq 0$ are rotating points.
Obviously the fixed points are $\\{|jj\rangle\rangle\\}_{j=1}^{d}$, and the
rotating points are $\\{|ij\rangle\rangle|\forall
i,j\in[1,d],\mathbf{c}_{j}^{\dagger}\mathbf{c}_{i}=e^{i\varphi}\neq 1\\}$.

As a simple example, consider $\\{\mathbf{c}_{j}\\}_{j=1}^{d}$ as a set of
orthonormal vectors, then the channel is
$\hat{\Phi}=\sum_{j=1}^{d}|jj\rangle\rangle\langle\langle jj|$, representing
PMs with rank-1 projectors (von Neumann measurements),
${\Phi}(\cdot)=\sum_{j=1}^{d}|j\rangle\langle j|(\cdot)|j\rangle\langle j|$.
This channel has only fixed points but no rotating points. As another example,
consider
$\\{\mathbf{c}_{j}\\}_{j=1}^{d}=\\{\widetilde{\mathbf{c}}e^{i\varphi_{j}}\\}_{j=1}^{d}$,
then
$\hat{\Phi}=\sum_{j=1}^{d}e^{i(\varphi_{i}-\varphi_{j})}|ij\rangle\rangle\langle\langle
ij|$ is a unitary channel ${\Phi}(\cdot)=U(\cdot)U^{\dagger}$ with
$U=\sum_{j=1}^{d}e^{i\varphi_{j}}|j\rangle\langle j|$. For the unitary
channel, $|ij\rangle\rangle$ is a fixed point if $i=j$ or
$\varphi_{i}=\varphi_{j}$, and a rotating point if
$\varphi_{i}\neq\varphi_{j}$.

For general cases, we divide the index set $A=\\{1,\cdots,d\\}$ into $s(\leq
d)$ disjoint subsets $A_{1},\cdots,A_{s}$, with the corresponding
cardinalities (number of elements) being $d_{1},\cdots,d_{s}$, satisfying
$\sum_{i=1}^{s}d_{i}=d$. Then divide the set of unit vectors
$C=\\{\mathbf{c}_{j}\\}_{j=1}^{d}$ into $s$ disjoint subsets
$C_{1},\cdots,C_{s}$ with $C_{k}=\\{\mathbf{c}_{j}|j\in A_{k}\\}$. This
division should ensure that the unit vectors in each subset are the same up to
some phase factors but are different from any other unit vectors in other
subsets, i.e. $C_{k}=\\{\widetilde{\mathbf{c}}_{k}e^{i\varphi_{j}}|j\in
A_{k}\\}$ but
$\widetilde{\mathbf{c}}_{p}\neq\widetilde{\mathbf{c}}_{q}e^{i\varphi}$ for any
$\varphi$ and $p,q\in[1,s]$. This implies that $|ij\rangle\rangle$ with
$i,j\in A_{k}$ is either a fixed point ($\varphi_{i}=\varphi_{j}$) or a
rotating point ($\varphi_{i}\neq\varphi_{j}$).

The division of the index set also partitions the Hilbert space $\mathcal{H}$
of the quantum system into the direct sum of $s$ subspaces,
$\mathcal{H}=\mathcal{H}_{1}\oplus\cdots\oplus\mathcal{H}_{s}$, where
$\mathcal{H}_{k}=\\{{\rm Span}\\{|j\rangle\\}|j\in A_{k}\\}$ with rank-$d_{k}$
projection $P_{k}=\sum_{j\in A_{k}}|j\rangle\langle j|$. Thus the measurement
operators in Eq. (12) can be written in a compact matrix form,
$\mathbf{M}=\widetilde{\mathbf{C}}\widetilde{\mathbf{P}}$, where
$\widetilde{\mathbf{C}}=[\widetilde{\mathbf{c}}_{1},\cdots,\widetilde{\mathbf{c}}_{s}]$
and $\widetilde{\mathbf{P}}=[\widetilde{P}_{1},\cdots,\widetilde{P}_{s}]^{T}$
with $\widetilde{P}_{k}=\sum_{j\in A_{k}}e^{i\varphi_{j}}|j\rangle\langle j|$.
Note that $\widetilde{P}_{k}$ is either a projection operator or a unitary
operator in $\mathcal{H}_{k}$, satisfying
$\widetilde{P}_{k}^{\dagger}\widetilde{P}_{k^{\prime}}=\delta_{kk^{\prime}}P_{k}$
and $\sum_{k=1}^{s}\widetilde{P}_{k}^{\dagger}\widetilde{P}_{k}=\mathbb{I}$.
Such a compact form of $\mathbf{M}$ allows us to extend the above formulation
to infinite-dimensional quantum systems SI , if we divide the identity
operator into a finite set of orthogonal projections.

Asymptotics of sequential POVMs. Sequential POVMs correspond to sequential
applications of the quantum channel $\hat{\Phi}$ [Fig. 1(a)]. Previous works
have studied the asymptotic behaviors of sequential general quantum channels
Albert2019 ; Burgarth2013 ; Novotny2018 ; Blume2010 , mostly trying to find
which information from an initial state can be preserved during the process.

For the channel in Eq. (13), as the number of applications $m$ increases, the
projections to eigenvectors with eigenvalues lying in the interior of the unit
disk ($|\mathbf{c}_{j}^{\dagger}\mathbf{c}_{i}|<1$) gradually vanish, while
the projections to eigenvectors with eigenvalues on the unit circle
($|\mathbf{c}_{j}^{\dagger}\mathbf{c}_{i}|=1$) remain unchanged or change by
some phase factors. So sequential POVMs tend to preserve the quantum coherence
within subspaces $\\{\mathcal{H}_{k}\\}_{k=1}^{s}$ but diminish the coherence
between different subspaces. First assume that the elements in each $C_{k}$
are all the same, i.e., $\varphi_{j}=0$ for all $j\in[1,d]$. In the asymptotic
limit of large $m$,

$\displaystyle\mathop{\lim}\limits_{m\to\infty}\hat{\Phi}^{m}=\sum_{k=1}^{s}\sum_{i,j\in
A_{k}}|ij\rangle\rangle\langle\langle
ij|=\sum_{k=1}^{s}\hat{\mathcal{P}}_{k},$ (14)

corresponding to
$\mathop{\lim}\limits_{m\to\infty}{\Phi}^{m}(\cdot)=\sum_{k=1}^{s}P_{k}(\cdot)P_{k}$
[Figs. 1(b) and (c)], which represents non-selective PMs. Then consider the
possible phase factors in
$C_{k}=\\{\widetilde{\mathbf{c}}_{k}e^{i\varphi_{j}}|j\in A_{k}\\}$, each
application of $\hat{\Phi}$ produces a unitary operation in the Hilbert
subspace $\mathcal{H}_{k}$, i.e., $P_{k}(\cdot)P_{k}$ in the former case
should be replaced by
$\widetilde{P}_{k}^{m}(\cdot)(\widetilde{P}_{k}^{\dagger})^{m}$. For example,
if
$C_{k}=\\{\mathbf{c}_{i},\mathbf{c}_{j}\\}=\\{\widetilde{\mathbf{c}}_{k}e^{i\varphi_{i}},\widetilde{\mathbf{c}}_{k}e^{i\varphi_{j}}\\}$,
then $\widetilde{P}_{k}=e^{i\varphi_{i}}|i\rangle\langle
i|+e^{i\varphi_{j}}|j\rangle\langle j|$.

Typicality of sequential POVMs. Now that sequential non-selective POVMs
produce projections (or unitary operations in the projected subspaces) in the
asymptotic limit, we further ask which sequences of sequential selective POVMs
produce a specific projection. This problem can be perfectly solved by the
theory of classical typicality Wilde2017 ; Cover2006 ; Facchi2015 ;
Goldstein2006 ; Bartsch2009 . Classical typicality mainly cares about the
following problem: if a random variable takes $r$ different values with the
probability distribution $(p_{1},\cdots,p_{r})$, generate $m$ independent
realizations of this variable and find the statistical distributions of the
event sequences with $(m_{1}/m,\cdots,m_{r}/m)$, where $m_{i}$ is the number
of the occurences of the $i$th value. For infinitely large $m$, the event
sequences that are overwhelmingly likely to occur are the set of typical
sequences with $(p_{1},\cdots,p_{r})$. By defining sequences of selective
POVMs as a product of superoperators [Fig. 1(a)], we find that the asymptotic
projections are induced by the sets of typical sequences of selective POVMs.

Since $\hat{\Phi}=\sum_{\alpha=1}^{r}\hat{\mathcal{M}}_{\alpha}$ and
$[\hat{\mathcal{M}}_{\alpha},\hat{\mathcal{M}}_{\beta}]=0$ for
$\alpha,\beta\in[1,r]$, we can expand $\hat{\Phi}^{m}$ according to the
multinomial theorem,

$\displaystyle\hat{\Phi}^{m}=\sum_{\\{F\\}}\frac{m!}{(mf_{1})!\cdots(mf_{r}!)}\hat{\mathcal{M}}_{1}^{mf_{1}}\cdots\hat{\mathcal{M}}_{r}^{mf_{r}},$
(15)

where $F=(f_{1},\cdots,f_{r})$ with $f_{i}\in[0,1]$ (also a rational number
with denominator $m$) satisfying $\sum_{i=1}^{r}f_{i}=1$, and the summation is
over all distributions $\\{F\\}$ in a $(r-1)$-dimensional probability space.
For large $m$, $\hat{\Phi}^{m}$ can be approximated by its projections in the
asymptotic subspaces SI ,

$\displaystyle\hat{\Phi}^{m}\approx\sum_{k=1}^{s}\hat{\mathcal{P}}_{k}\hat{\Phi}^{m}\hat{\mathcal{P}}_{k}\approx\sum_{k=1}^{s}\sum_{\\{F\\}}e^{-mS(F\|F_{k})}\hat{\mathcal{P}}_{k},$
(16)

where
$F_{k}=(f_{k1},\cdots,f_{kr})=(|\widetilde{c}_{1k}|^{2},\cdots,|\widetilde{c}_{rk}|^{2})$
with $\widetilde{c}_{1k},\cdots,\widetilde{c}_{rk}$ being entries of
$\widetilde{\mathbf{c}}_{k}$ satisfying
$\sum_{i=1}^{r}|\widetilde{c}_{ik}|^{2}=1$, and
$S(F\|F_{k})=\sum_{i=1}^{r}f_{i}\ln(f_{i}/f_{ki})$ is the relative entropy
between $F$ and $F_{k}$ (the derivation above uses Stirling’s formula $\ln
m!\approx m\ln m-m$ for large $m$). $S(F\|F_{k})$ takes the minimum when
$F=F_{k}$, so for infinite large $m$, $\\{F_{k}\\}_{k=1}^{s}$ represent $s$
sets of ideal typical sequences of POVMs leading to the projections
$\\{\hat{\mathcal{P}}_{k}\\}_{k=1}^{s}$ correspondingly [Fig. 1(b)].

For large but finite $m$, the distributions of POVM sequences for
$\hat{\mathcal{P}}_{k}$ are concentrated around $F_{k}$, so
$S(F\|F_{k})\approx\sum_{i=1}^{r}(f_{i}-f_{ki})^{2}/(2f_{ki})$. Then Eq. (16)
represents the summation of $s$ Gaussians around $F_{1},\cdots,F_{s}$, with
integration of the $k$th Gaussian over the whole probability space giving rise
to $\hat{\mathcal{P}}_{k}$. For any two Gaussians around $F_{j}$ and $F_{k}$,
they are well separated if the distance between $F_{j}$ and $F_{k}$ is larger
than the sum of the respective Gaussian half widths. This requires
$m>2|\ln\eta|[(\sum_{i=1}^{r}(f_{ji}-f_{ki})^{2}/f_{ji})^{-1/2}+(\sum_{i=1}^{r}(f_{ji}-f_{ki})^{2}/f_{ki})^{-1/2}]^{2}$
SI ; Liu2010 , where $\eta$ is the ratio of the minimum hight to the maximum
hight within the Gaussian width. If all the Gaussians are well separated,
integration of the POVM sequences within a small neighborhood around $F_{k}$
can approximate $\hat{\mathcal{P}}_{k}$ up to arbitrary small error as $m$
increases (see the Supplementary Material SI  for the error rates with finite
$m$).

It may happen that two Gaussians coincide around $F_{j}=F_{k}$ but
$\widetilde{\mathbf{c}}_{j}\neq\widetilde{\mathbf{c}}_{k}$, i.e., partial
elements of $\widetilde{\mathbf{c}}_{j}$ and $\widetilde{\mathbf{c}}_{k}$
differ by some phase factors. Since
$|\widetilde{\mathbf{c}}_{j}^{\dagger}\widetilde{\mathbf{c}}_{k}|<1$, the
coinciding Gaussians actually correspond to different projections, and the
POVMs sequences around $F_{j}$ approximately produce
$\hat{\mathcal{P}}_{j}+\hat{\mathcal{P}}_{k}$. To realize selective
projections, we can get a new set of measurement operators by a unitary
transformation, thus creating different typical sequences of POVMs for
$\hat{\mathcal{P}}_{j}$ and $\hat{\mathcal{P}}_{k}$.

Figure 2: (a) Quantum circuit diagram to realize sequential POVMs on the
target system with PMs of an ancilla qubit. (b) Distributions of eigenvalues
of $U_{\pm}(t)=e^{\pm i\chi a^{\dagger}at/2}$ in the complex unit circle to
detect the $k$ mod $N$ excitation numbers of a bosonic mode for
$t=2\pi/(N\chi)$ and $N$=2, 3.

Physical realization. We present a general physical model to perform PMs on a
$d$-level quantum system (called target system) with sequential POVMs. Without
loss of generality, we assume that the POVMs are realized by PMs of an ancilla
qubit. The coupling Hamiltonian of the composite system (including the ancilla
and target systems) is in the pure-dephasing form Yang2017

$\displaystyle H(t)=\sigma_{z}\otimes B(t),$ (17)

where $\sigma_{i}$ is the Pauli-$i$ operator of the ancilla qubit ($i=x,y,z$),
and $B(t)$ is a time-dependent Hermitian operator of the target system (the
time-dependence of $B(t)$ is due to being in some interaction picture or
external drivings).

The dynamics of the composite system induces a general class of quantum
channels on the target system, which can be written in the Stinespring
representation as Stine1955

$\displaystyle\Phi(\rho)={\rm
Tr}_{a}[U(t)(\rho_{a}\otimes\rho)U^{\dagger}(t)],$ (18)

where
$U(t)=\mathcal{T}e^{-i\sigma_{z}\otimes\int_{0}^{t}B(t^{\prime})dt^{\prime}}$
with $\mathcal{T}$ being the time-ordering operator,
$\rho_{a}=|\psi\rangle_{a}\langle\psi|$ is the initial state of the ancilla,
$\rho$ denotes the density matrix of the target system, and ${\rm Tr}_{a}$
denotes the partial trace over the ancilla. With an orthonormal ancilla basis
$\\{|v_{+}\rangle_{a},|v_{-}\rangle_{a}\\}$, we obtain the Kraus
representation of the quantum channels,
$\Phi(\rho)=\sum_{\alpha\in\\{+,-\\}}M_{\alpha}\rho M_{\alpha}^{\dagger}$ with
$M_{\alpha}=\langle v_{\alpha}|U(t)|\psi\rangle_{a}$ (note that we add
subscripts to the kets only when representing matrix elements or inner
products with respect to the acilla states). With another orthonormal basis
$\\{T|v_{+}\rangle_{a},T|v_{-}\rangle_{a}\\}$ with $T$ being a unitary
operator for the ancilla, the measurement operators become
$\\{M^{\prime}_{\alpha}\\}$ with
$M^{\prime}_{\alpha}=\sum_{\beta\in\\{+,-\\}}T_{\alpha\beta}M_{\beta}$, while
the quantum channels remain unchanged.

We expand $U(t)$ in the ancilla eigenbasis $\\{|+\rangle_{a},|-\rangle_{a}\\}$
of $\sigma_{z}$ as $U(t)=|+\rangle_{a}\langle+|\otimes
U_{+}(t)+|-\rangle_{a}\langle-|\otimes U_{-}(t)$, where
$U_{\pm}(t)=\mathcal{T}e^{\mp i\int_{0}^{t}B(t^{\prime})dt^{\prime}}$. If
$U_{\pm}(t)$ is exactly equal to or well approximated by its first-order
Magnus expansion Ma2016 , i.e., $U_{\pm}(t)=e^{\mp
i\int_{0}^{t}B(t^{\prime})dt^{\prime}}$, then $U_{+}(t)=U_{-}^{\dagger}(t)$
and $[U_{+}(t),U_{-}(t)]=0$, so $U_{+}(t)$ and $U_{-}(t)$ can be
simultaneously diagonalized as $U_{\pm}(t)=\sum_{j=1}^{d}e^{\pm
i\omega_{j}}|j\rangle\langle j|$. So the measurement operators are
$M_{\pm}=\sum_{j=1}^{d}\left(\langle
v_{\pm}|\psi\rangle_{a}\cos\omega_{j}+i\langle
v_{\pm}|\sigma_{z}|\psi\rangle_{a}\sin\omega_{j}\right)|j\rangle\langle j|$.
As a special case, take
$|\psi\rangle_{a}=R_{\phi_{1}}(\frac{\pi}{2})|+\rangle_{a}$ and
$|v_{\pm}\rangle_{a}=R_{\phi_{2}}(-\frac{\pi}{2})|\pm\rangle_{a}$ with
$R_{\phi}(\theta)=e^{-i(\cos\phi\sigma_{x}+\sin\phi\sigma_{y})\theta/2}$, then

$\displaystyle\left[{\begin{array}[]{*{20}{c}}{{{M}_{+}}}\\\
{{{M}_{-}}}\end{array}}\right]=\left[{\begin{array}[]{*{20}{c}}{e^{i\omega_{1}}-e^{i(\Delta\phi-\omega_{1})}}&\cdots&{e^{i\omega_{d}}-e^{i(\Delta\phi-\omega_{d})}}\\\
{e^{i\omega_{1}}+e^{i(\Delta\phi-\omega_{1})}}&\cdots&{e^{i\omega_{d}}+e^{i(\Delta\phi-\omega_{d})}}\end{array}}\right]\frac{\mathbf{P}}{2}$
(23)

where $\Delta\phi=\phi_{1}-\phi_{2}$. Each round of such POVMs corresponds to
a three-step physical process [Fig. 2(a)]: (1) the ancilla starts from
$|+\rangle_{a}$ and is rotated by $R_{\phi_{1}}(\frac{\pi}{2})$; (2) let the
ancilla and target systems evolve under $H(t)$ for time $t$; (3) finally
rotate the ancilla by $R_{\phi_{2}}(\frac{\pi}{2})$ and make a PM of the
ancilla in the basis $\\{|+\rangle_{a},|-\rangle_{a}\\}$. Similar schemes have
been designed to realize single-shot readouts of nuclear spins-1/2 in diamond
Liu2017 , but here we show this scheme can be extended to perform PMs of a
generic system.

Since the POVMs have only two outcomes, the measurement results are solely
determined by the polarization $\Delta f=(m_{-}-m_{+})/m$ Ma2018 , with
$m_{+}/m_{-}$ being the number of outcome $+/-$ in $m$ sequential
measurements. For the spectra $\\{e^{\pm i\omega_{j}}\\}$ of $U_{\pm}(t)$,
calculate $\Delta f_{j}=\cos(2\omega_{j}-\Delta\phi)$ for all $j\in[1,d]$.
Weak measurement corresponds to the regime $|\Delta f_{j}|\ll 1$. If $\Delta
f_{j}\neq\Delta f_{k}$ for any $j,k\in[1,d]$ and $j\neq k$, sequential POVMs
produce von-Neumann measurements of the target system, with the rank-1
projection $P_{j}=|j\rangle\langle j|$ corresponding to typical POVM sequences
with $\Delta f_{j}$. If $\Delta f_{j}=\Delta f_{k}$, then either (I)
$\omega_{j}+\omega_{k}=\Delta\phi+n\pi$ or (II) $\omega_{j}-\omega_{k}=n\pi$
with $n$ being integers. In case-I, the typical POVM sequences for $P_{j}$ and
$P_{k}$ are the same, but selective projections can still be achieved by
choosing a different $\Delta\phi^{\prime}$. In case-II, the typical POVM
sequences with $\Delta f_{j}$ induce the operation $P_{j}+(-1)^{n}P_{k}$.

Example: Modular excitation number measurements of bosonic modes. As an
example, we present a protocol to measure the modular excitation numbers of a
bosonic mode with an ancilla qubit. The ancilla is dispersively coupled to a
bosonic mode with the Hamiltonian $H=-\chi\sigma_{z}a^{\dagger}a/2$, where $a$
($a^{\dagger}$) is the annihilation (creation) operator of the bosonic mode
and $\chi$ is the dispersive coupling strength. The dispersive coupling arises
naturally from the Jaynes-Cumming coupling in cavity quantum electrodynamics
(QED) Raimond2001  and circuit QED Blais2021  when the detuning between the
ancilla and the bosonic mode is much larger than the coupling strength.

We construct the projectors into the sets of bosonic Fock states with modular
excitation number $l$ mod $2N$,
$P_{2N}^{l}=\sum_{j=0}^{\infty}|2jN+l\rangle\langle 2jN+l|$, with
$l\in\\{0,1,\cdots,2N-1\\}$ and $N$ being any positive integer. With the
scheme below Eq. (23) and the evolution time $t=2\pi/(N\chi)$,
$U_{\pm}(t)=e^{\pm i\chi a^{\dagger}at/2}=\sum_{k=0}^{N-1}e^{\pm
ik\pi/N}(P_{2N}^{k}-P_{2N}^{k+N})$, i.e. the eigenvalues of $U_{\pm}(t)$
divides the complex unit circle into $2N$ equal pieces [Fig. 2(b)]. The
measurement operators are $M_{\pm}=\sum_{k=0}^{N-1}(e^{ik\pi/N}\mp
e^{i(\Delta\phi-k\pi/N)})(P_{2N}^{k}-P_{2N}^{k+N})$, and the polarization
$\Delta f_{k}=\cos(2k\pi/N-\Delta\phi)$. We can tune $\Delta\phi$ so that
$\Delta f_{k}$ is maximally distinguishable for different $k\in[0,N-1]$. For
$N=1$, $\Delta\phi=0$ is optimal as $\Delta f_{0}=-\Delta f_{1}=1$; while for
$N\geq 2$, we can choose $\Delta\phi=\pi/(2N)$ so that $\Delta
f_{k}=\cos[(2k-1/2)\pi/N]$. Then for a large and even $m$, sequential POVMs
induce the $k$ mod $N$ excitation number measurement of the bosonic mode. The
modular excitation numbers are the error syndromes of rotation-symmetric error
correction codes of bosonic modes Grimsmo2020 , such as cat codes Leghtas2013
; Mirrahimi2014 ; Li2017 ; Bergmann2016  and binomial codes Michael2016 . So
this protocol is useful for quantum non-demolition measurements in bosonic
quantum information processing Blais2020 ; Cai2021 ; Joshi2021 ; Ma2021 ,
especially for tracking the error syndromes of high-order bosonic error
correction codes Sun2014 ; Ofek2016 ; Hu2019 .

Summary. We have revealed the elegant structures of sequential POVMs by
studying their asymptotic behaviors and typical sequences. We rigorously prove
that PMs emerge from sequential POVMs when the measurement operators are
normal and commuting with each other. Each selective projective measurement
comes from a set of typical sequences of sequential POVMs, which is determined
solely by the structures of the measurement operators. While the POVMs here
are restricted to have normal and commuting measurement operators, they
describes a large class of quantum channels on a quantum system induced by a
pure-dephasing coupling between this system and an ancilla system. For future
works, it will be interesting to relax this restriction, and study the
asymptotics and typicality of sequential POVMs with general measurement
operators.

###### Acknowledgements.

We acknowledge support from the Startup Foundation of Institute of
Semiconductors, Chinese Academy of Sciences (No. E0SEBB11), and National
Natural Science Foundation of China (No. 12174379). R.B.L. was supported by
the Hong Kong Research Grants Council - General Research Fund Project
14300119.

## References

  * (1) H. M. Wiseman and G. J. Milburn, Quantum Measurement and Control (Cambridge University Press, Cambridge, England, 2010).
  * (2) K. Jacobs, Quantum Measurement Theory and Its Applications (Cambridge University Press, Cambridge, England, 2014).
  * (3) M. A. Nielsen and I. L. Chuang, Quantum Computation and Quantum Information (Cambridge University Press, Cambridge, England, 2010).
  * (4) A. Peres, Quantum Theory: Concepts and Methods (Springer Science & Business Media, New York, 2006), Vol. 57.
  * (5) J. M. Eizerman, R. Hanson, L. H. W. Van Beveren, B. Witkamp, L. M. K. Vandersypen, and L. P. Kouwenhoven, Single-Shot Read-out of an Individual Electron Spin in a Quantum Dot, Nature 430, 431 (2004).
  * (6) A. N. Vamivakas, C.-Y. Lu, C. Matthiesen, Y. Zhao, S. Fält, A. Badolato, and M. Atat re, Observation of spin-dependent quantum jumps via quantum dot resonance fluorescence, Nature 467, 297 (2010).
  * (7) P. Neumann, J. Beck, M. Steiner, F. Rempp, H. Fedder, P. R. Hemmer, J. Wrachtrup, and F. Jelezko, Single-shot readout of a single nuclear spin, Science 329, 542 (2010).
  * (8) A. Morello, J. J. Pla, F. A. Zwanenburg, K. W. Chan, K. Y. Tan, H. Huebl, M. Möttönen, C. D. Nugroho, C. Yang, J. A. van Donkelaar, A. D. C. Alves, D. N. Jamieson, C. C. Escott, L. C. L. Hollenberg, R. G. Clark, and A. S. Dzurak, Single-shot readout of an electron spin in silicon, Nature 467, 687 (2010).
  * (9) L. Jiang, J. S. Hodges, J. R. Maze, P. Maurer, J. M. Taylor, D. G. Cory, P. R. Hemmer, R. L. Walsworth, A. Yacoby, A. S. Zibrov, and M. D. Lukin, Repetitive readout of a single electronic spin via quantum logic with nuclear spin ancillae, Science 326, 267 (2009).
  * (10) T. Nakajima et al., Robust Single-Shot Spin Measurement with 99.5% Fidelity in a Quantum Dot Array, Phys. Rev. Lett. 119, 017701 (2017).
  * (11) A. West et al., Gate-Based Single-Shot Readout of Spins in Silicon, Nat. Nanotechnol. 14, 437 (2019).
  * (12) K. Kraus, States, Effects, and Operations: Fundamental Notions of Quantum Theory, Lecture Notes in Physics Vol. 190 (Springer-Verlag, Berlin, 1983).
  * (13) E. Andersson and D. K. L. Oi, Binary search trees for generalized measurements, Phys. Rev. A 77, 052108 (2008).
  * (14) Y. -H. Chen and T. A. Brun, Decomposing qubit positive-operator-valued measurements into continuous destructive weak measurements, Phys. Rev. A 98, 062113 (2018).
  * (15) Y. -H. Chen and T. A. Brun, Qubit positive-operator-valued measurements by destructive weak measurements, Phys. Rev. A 99, 062121 (2019).
  * (16) Y. W. Cheong and S. -W. Lee, Balance between Information Gain and Reversibility in Weak Measurement, Phys. Rev. Lett. 109, 150402 (2012).
  * (17) J. M. Renes, R. Blume-Kohout, A. J. Scott, and C. M. Caves, Symmetric informationally complete quantum measurements, J. Math. Phys. 45, 2171 (2004).
  * (18) J. A. Bergou, Discrimination of quantum states, J. Mod. Opt. 57, 160 (2010).
  * (19) R. Derka, V. Buz̆ek, and E. A. Kert, Universal Algorithm for Optimal Estimation of Quantum States from Finite Ensembles via Realizable Generalized Measurement, Phys. Rev. Lett. 80, 1571 (2001).
  * (20) K. Jacobs and D. A. Steck, A straightforward introduction to continuous quantum measurement, Contemp. Phys. 47, 279 (2006).
  * (21) S. A. Gurvitz, Measurements with a noninvasive detector and dephasing mechanism, Phys. Rev. B 56, 15215 (1997).
  * (22) S. Ashhab, J. Q. You, and F. Nori, The information about the state of a qubit gained by a weakly coupled detector, New J. Phys. 11, 083017 (2009).
  * (23) A. N. Korotkov, Selective quantum evolution of a qubit state due to continuous measurement, Phys. Rev. B 63, 115403 (2001).
  * (24) M. S. Blok, C. Bonato, M. L. Markham, D. J. Twitchen, V. V. Dobrovitski, and R. Hanson, Manipulating a qubit through the backaction of sequential partial measurements and real-time feedback, Nat. Phys. 10, 189 (2014)
  * (25) A. N. Jordan and A. N. Korotkov, Qubit feedback and control with kicked quantum nondemolition measurements: A quantum Bayesian analysis, Phys. Rev. B 74, 085307 (2006).
  * (26) A. Chantasri, J. Dressel, and A. N. Jordan, Action principle for continuous quantum measurement, Phys. Rev. A 88, 042110 (2013).
  * (27) A. Chantasri and A. N. Jordan, Stochastic path-integral formalism for continuous quantum measurement, Phys. Rev. A 92, 032125 (2015)
  * (28) C. Presilla, R. Onofrio, and U. Tambini, Measurement quantum mechanics and experiments on quantum Zeno effect, Ann. Phys. 248, 95 (1996).
  * (29) L. Diósi, Structural features of sequential weak measurements, Phys. Rev. A 94, 010103 (2016).
  * (30) E. Shojaee, C. S. Jackson, C. A. Riofrío, A. Kalev, and I. H. Deutsch, Optimal Pure-State Qubit Tomography via Sequential Weak Measurements, Phys. Rev. Lett. 121, 130404 (2018).
  * (31) J. T. Monroe, N. Yunger Halpern, T. Lee, and K. W. Murch, Weak Measurement of a Superconducting Qubit Reconciles Incompatible Operators, Phys. Rev. Lett. 126, 100403 (2021).
  * (32) P. Wang, C. Chen, X. Peng, J. Wrachtrup, and R. B. Liu, Characterization of Arbitrary-Order Correlations in Quantum Baths by Weak Measurement, Phys. Rev. Lett. 123, 50603 (2019).
  * (33) M. Pfender et al., High-resolution spectroscopy of single nuclear spins via sequential weak measurements, Nat. Commun. 10, 594 (2019).
  * (34) Z. Wu, P. Wang, T. Wang, Y. Li, R. Liu, Y. Chen, X. Peng, R.-B. Liu, and J. Du, Detection of Arbitrary Quantum Correlations via Synthesized Quantum Channels, arXiv: 2206. 05883 (2022).
  * (35) M. Oszmaniec, L. Guerini, P. Wittek, and A. Acín, Simulating Positive-Operator-Valued Measures with Projective Measurements, Phys. Rev. Lett. 119, 190501 (2017).
  * (36) M. Oszmaniec, F. B. Maciejewski, and Z. Puchała, Simulating all quantum measurements using only projective measurements and postselection, Phys. Rev. A 100, 01235 (2019).
  * (37) T. A. Brun, A simple model of quantum trajectories, Am. J. Phys. 70, 719 (2002).
  * (38) D. Lidar and T. D. Brun, Quantum Error Correction (Cambridge University Press, Cambridge, England, 2013).
  * (39) O. Oreshkov and T. A. Brun, Weak Measurements Are Universal, Phys. Rev. Lett. 95, 110409 (2005).
  * (40) W. -L. Ma, P. Wang, W. -H. Leong, and R. -B. Liu, Phase transitions in sequential weak measurements, Phys. Rev. A 98, 012117 (2018).
  * (41) G. -Q. Liu, J. Xing, W. -L. Ma, P. Wang, C. -H. Li, H. C. Po, Y. -R. Zhang, H. Fan, R. -B. Liu, and X. -Y. Pan, Single-Shot Readout of a Nuclear Spin Weakly Coupled to a Nitrogen-Vacancy Center at Room Temperature, Phys. Rev. Lett. 118, 150504 (2017).
  * (42) D. D. Bhaktavatsala Rao, S. Yang, S. Jesenski, E. Tekin, F. Kaiser, and J. Wrachtrup, Observation of nonclassical measurement statistics induced by a coherent spin environment, Phys. Rev. A 100, 22307 (2019).
  * (43) F. Caruso, V. Giovannetti, C. Lupo, and S. Mancini, Quantum channels and memory effects, Rev. Mod. Phys. 86, 1203 (2014).
  * (44) M. M. Wolf, Quantum Channels & Operations Guided Tour, 2010.
  * (45) I. Bengtsson and K. Życzkowski, Geometry of Quantum States: An Introduction to Quantum Entanglement (Cambridge University Press, Cambridge, England, 2017).
  * (46) See Supplementary Material for details about the HS space and POVM, structral propeties of POVMs with normal and commuting measurement operators, matrix representation of quantum channels for POVMs, derviations in typicality of sequential POVMs including the conditions and error rates for approximating sequential POVMs with PMs.
  * (47) S. R. Garcia and R. A. Horn, A Second Course in Linear Algebra (Cambridge University Press, Cambridge, England, 2017).
  * (48) V. V. Albert, Asymptotics of quantum channels: conserved Quantities, an adiabatic Limit, and matrix product states, Quantum 3, 151 (2019).
  * (49) D. Burgarth, G. Chiribella, V. Giovannetti, P. Perinotti, and K. Yuasa, Ergodic and mixing quantum channels in finite dimensions, New J. Phys. 15, 073045 (2013).
  * (50) J. Novotný, J. Maryška, and I. Jex, Quantum Markov processes: From attractor structure to explicit forms of asymptotic states: Asymptotic dynamics of quantum Markov processes, Eur. Phys. J. Plus 133, 310 (2018).
  * (51) R. Blume-Kohout, H. K. Ng, D. Poulin, and L. Viola, Information-Preserving Structures: A General Framework for Quantum Zero-Error Information, Phys. Rev. A 82, 062306 (2010).
  * (52) M. W. Wilde, Quantum Information Theory (Cambridge University Press, Cambridge, England, 2017).
  * (53) T. M. Cover and J. A. Thomas, Elements of Information Theory (Wiley-Interscience, New York, New Jersey, 2006).
  * (54) P. Facchi, S. Pascazio, and F. V. Pepe, Quantum typicality and initial conditions, Phys. Scr. 90, 074057 (2015).
  * (55) S. Goldstein, J. L. Lebowitz, R. Tumulka, and N. Zanghì, Canonical Typicality, Phys. Rev. Lett. 96, 050403 (2006).
  * (56) C. Bartsch and J. Gemmer, Dynamical Typicality of Quantum Expectation Values, Phys. Rev. Lett. 102, 110403 (2009).
  * (57) R. B. Liu, W. Yao, and L. J. Sham, Quantum Computing by Optical Control of Electron Spins, Adv. Phys. 59, 703 (2010).
  * (58) W. Yang, W. -L. Ma, and R. -B. Liu, Quantum Many-Body Theory for Electron Spin Decoherence in Nanoscale Nuclear Spin Baths, Rep. Prog. Phys. 80, 16001 (2017).
  * (59) W. F. Stinespring, Positive functions on $C^{*}$-algebras, Proc. Am. Math. Soc. 6, 211 (1955).
  * (60) W. -L. Ma and R. -B. Liu, Angstrom-Resolution Magnetic Resonance Imaging of Single Molecules via Wave-Function Fingerprints of Nuclear Spins, Phys. Rev. Appl. 6, 024019 (2016).
  * (61) J. M. Raimond, M. Brune, and S. Haroche, Colloquium: Manipulating Quantum Entanglement with Atoms and Photons in a Cavity, Rev. Mod. Phys. 73, 565 (2001).
  * (62) A. Blais, A. L. Grimsmo, S. M. Girvin, and A. Wallraff, Circuit Quantum Electrodynamics, Rev. Mod. Phys. 93, 025005 (2021).
  * (63) A. L. Grimsmo, J. Combes, and B. Q. Baragiola, Quantum Computing with Rotation-Symmetric Bosonic Codes, Phys. Rev. X 10, 11058 (2020).
  * (64) Z. Leghtas, G. Kirchmair, B. Vlastakis, R. J. Schoelkopf, M. H. Devoret, and M. Mirrahimi, Hardware-Efficient Autonomous Quantum Memory Protection, Phys. Rev. Lett. 111, 120501 (2013).
  * (65) M. Mirrahimi, Z. Leghtas, V. V. Albert, S. Touzard, R. J. Schoelkopf, L. Jiang, and M. H. Devoret, Dynamically Protected Cat-Qubits: A New Paradigm for Universal Quantum Computation, New J. Phys. 16, (2014).
  * (66) L. Li, C. L. Zou, V. V. Albert, S. Muralidharan, S. M. Girvin, and L. Jiang, Cat Codes with Optimal Decoherence Suppression for a Lossy Bosonic Channel, Phys. Rev. Lett. 119, 030502 (2017).
  * (67) M. Bergmann and P. Van Loock, Quantum Error Correction against Photon Loss Using Multicomponent Cat States, Phys. Rev. A 94, 042332 (2016).
  * (68) M. H. Michael, M. Silveri, R. T. Brierley, V. V. Albert, J. Salmilehto, L. Jiang, and S. M. Girvin, New Class of Quantum Error-Correcting Codes for a Bosonic Mode, Phys. Rev. X 6, 031006 (2016).
  * (69) A. Blais, S. M. Girvin, and W. D. Oliver, Quantum Information Processing and Quantum Optics with Circuit Quantum Electrodynamics, Nat. Phys. 16, 247 (2020).
  * (70) W. Cai, Y. Ma, W. Wang, C.-L. Zou, and L. Sun, Bosonic Quantum Error Correction Codes in Superconducting Quantum Circuits, Fundam. Res. 1, 50 (2021).
  * (71) A. Joshi, K. Noh, and Y. Y. Gao, Quantum Information Processing with Bosonic Qubits in Circuit QED, Quantum Sci. Technol. 6, 033001 (2021).
  * (72) W. -L. Ma, S. Puri, R. J. Schoelkopf, M. H. Devoret, S. M. Girvin, and L. Jiang, Quantum Control of Bosonic Modes with Superconducting Circuits, Sci. Bull. 66, 1789 (2021).
  * (73) L. Sun, A. Petrenko, Z. Leghtas, B. Vlastakis, G. Kirchmair, K. M. Sliwa, A. Narla, M. Hatridge, S. Shankar, J. Blumoff, L. Frunzio, M. Mirrahimi, M. H. Devoret, and R. J. Schoelkopf, Tracking Photon Jumps with Repeated Quantum Non-Demolition Parity Measurements, Nature 511, 444 (2014).
  * (74) N. Ofek, A. Petrenko, R. Heeres, P. Reinhold, Z. Leghtas, B. Vlastakis, Y. Liu, L. Frunzio, S. M. Girvin, L. Jiang, M. Mirrahimi, M. H. Devoret, and R. J. Schoelkopf, Extending the Lifetime of a Quantum Bit with Error Correction in Superconducting Circuits, Nature 536, 441 (2016).
  * (75) L. Hu, Y. Ma, W. Cai, X. Mu, Y. Xu, W. Wang, Y. Wu, H. Wang, Y. P. Song, C. L. Zou, S. M. Girvin, L. M. Duan, and L. Sun, Quantum Error Correction and Universal Gate Set Operation on a Binomial Bosonic Logical Qubit, Nat. Phys. 15, 503 (2019).


# Artificial neural network as a universal model of nonlinear dynamical
systems

Pavel V. Kuptsov1,∗, Anna V. Kuptsova2, Nataliya V. Stankevich1

( 1 Laboratory of topological methods in dynamics, National Research
University Higher School of Economics, Nizhny Novgorod, 25/12 Bolshay
Pecherskaya str., Nizhny Novgorod 603155, Russia
2 Institute of electronics and mechanical engineering, Yuri Gagarin State
Technical University of Saratov, Politekhnicheskaya 77, Saratov 410054, Russia
∗ Send correspondence to P.V.K. E-mail<EMAIL_ADDRESS>
August 27, 2024)

###### Abstract

We suggest a universal map capable to recover a behavior of a wide range of
dynamical systems given by ODEs. The map is built as an artificial neural
network whose weights encode a modeled system. We assume that ODEs are known
and prepare training datasets using the equations directly without computing
numerical time series. Parameter variations are taken into account in the
course of training so that the network model captures bifurcation scenarios of
the modeled system. Theoretical benefit from this approach is that the
universal model admits using common mathematical methods without needing to
develop a unique theory for each particular dynamical equations. Form the
practical point of view the developed method can be considered as an
alternative numerical method for solving dynamical ODEs suitable for running
on contemporary neural network specific hardware. We consider the Lorenz
system, the Röessler system and also Hindmarch–Rose neuron. For these three
examples the network model is created and its dynamics is compared with
ordinary numerical solutions. High similarity is observed for visual images of
attractors, power spectra, bifurcation diagrams and Lyapunov exponents.

Keywords: neural network, dynamical system, numerical solution, universal
approximation theorem, Lyapunov exponents

##  1 Introduction

In this paper we suggest a map, i.e., a discrete time dynamical system,
capable to recover dynamics of systems given by ODEs. The map is built as an
artificial neural network whose encode the modeled system. Using neural
networks for dynamical systems reconstruction is a long standing problem. But
typically networks are used to predict dynamics when governing equations are
unknown and only time series are available [1, 2, 3, 4]. We assume that ODEs
are known and create their neural network model. The structure of the network
is the same for all cases while network weights are trained to fit the modeled
dynamical system. To prepare training datasets we do not use system time
series. Instead we feed the network and the modeled system by random time
series sampled form normal distributions and update network weights comparing
the outputs. Parameter variations are taken into account so that the network
captures bifurcation scenarios of the modeled system.

The motivations of this study are the following. We want to develop a method
of training an artificial neural network that can operates as a discrete time
system and can reproduce behavior of a wide variety of dynamical systems. We
consider a perceptron with one hidden layer and sigmoidal activation. Also
such network is said to consists of two dense layers. Many contemporary
investigations deal with deep networks whose number of layers is much more
then two and whose neuron interconnections is much more complicated. We prefer
a classical architecture because of a solid mathematical background behind it,
that is the universal approximation theorem. According to this theorem,
considered network is the simplest universal approximator, i.e., is able to
reproduce any function of multiple variables on a compact set. Such simple
universal model can be interesting for theoretical studies. Theoretical
analysis of a dynamical system often requires developing highly specialized
mathematical approaches unique for the system. Considering the universal model
that covers a wide range of systems one can extend the theoretical results for
these range of systems without needing always recreate a special mathematical
approaches.

From the practical point of view the development of methods of creation of the
universal neural network that can models dynamics can be considered an
alternative to the existing numerical methods for solving dynamical ODEs.
Although a large variety of well established and effective numerical methods
is available for computer simulation of dynamics, these methods are basically
developed for single-thread computation. But contemporary trend in
computational hardware development is in increasing of a number of computation
cores instead of increasing of single-core speed. In particular many hardware
are known today specialized for implementing artificial neural networks. In
this situation it seems to be very important to develop new numerical
approaches good fitted to a powerful contemporary hardware. Our model operates
as a neural network that can be run either using various available today
network software like TensorFlow [5] and PyTorch [6] or it can be downloaded
to a dedicated computer chip called AI accelerator (AI stands for artificial
intelligence) [7, 8].

##  2 Mathematical background: the universal approximation theorem

The problem of a universal construct for approximation of functions with many
variables has a long story. First can be mentioned the Weierstrass theorem [9]
that states that any continuous function over a closed interval on the real
axis can be expressed in that interval as an absolutely and uniformly
convergent series of polynomials. David Hilbert in the International Congress
of Mathematicians in Paris in the year 1900 outlined 23 major mathematical
problems for in the coming new century. His 13th problem is whether solutions
to 7th degree polynomial equation can be written as the composition of
finitely many two-variable functions. Hilbert believed they could not be. In
1956-57 years, Kolmogorov and Arnold proved that each continuous function of
$N$ variables — including the case in which $N=7$ — can be written as a
composition of continuous functions of two variables [10, 11, 12]. This is
called Kolmogorov–Arnold representation theorem.

Research interest in the virtues of multilayer perceptrons as devices for the
representation of arbitrary continuous functions was perhaps first put into
focus by Hecht-Nielsen [13]. In the context of traditional multilayer
perceptrons, it was Cybenko who demonstrated rigorously for the first time
that a single hidden layer is sufficient to uniformly approximate any con-
tinuous function with support in a unit hypercube [14]. In 1989, two other
papers were published independently on multilayer perceptrons as universal
approximators [15, 16]. For subsequent contributions to the approximation
problem, see [17]. Review on this topic can also be found in [1]

To sum up, universal approximation theorem states that a feed-forward network
with a single hidden layer containing a finite number of neurons can
approximate continuous functions on compact subsets of $\mathbb{R}^{N}$, under
mild assumptions on the activation function. The theorem thus states that
simple neural networks can represent a wide variety of interesting functions
when given appropriate parameters.

##  3 The network and training details

Assume that we have ODE

$\dot{u}=f(u,p),$ (1)

where $u\in\mathbb{R}^{N_{u}}$ is a vector of $N_{u}$ dynamical variables, and
$p\in\mathbb{R}^{N_{p}}$ is a vector of $N_{p}$ parameters.

We consider a perceptron with one hidden layer, or using more contemporary
terms, a network with two dense layers. Formally the network can be
represented as a function that maps vectors $u\in\mathbb{R}^{N_{u}}$ to
vectors $d\in\mathbb{R}^{N_{u}}$,

$d=F(u,p,w),$ (2)

where $w$ is a vector of network weights. Our purpose is to tune $w$ in such a
way that

$u(t+\Delta t)=u(t)+d(t)$ (3)

where $u(t)$ is a solution to ODE (1) and $\Delta t$ is a time step. The size
of the time step is defined before training the network. We take $\Delta
t=0.01$.

Consider a semi implicit numerical scheme of ODE solution:

$u(t+\Delta t)=u(t)+\frac{\Delta t}{2}\\{f[u(t)]+f[u(t+\Delta t)]\\}$ (4)

Compute the difference between (3) and (4):

$e=d(t)-\frac{\Delta t}{2}\\{f[u(t)]+f[u(t+\Delta t)]\\},$ (5)

where $e$ is the approximation error. Substituting $u(t+\Delta t)$ as
$u(t)+d(t)$ from Eq. (3) and omitting $t$ we obtain

$e=d-\frac{\Delta t}{2}[f(u)+f(u+d)]$ (6)

The network approximation (2), (3) works well if the approximation error tends
to zero $e\to 0$ for any $u$ and $p$ from the domain of interest.

Before training we need to define a localization areas for $u$ and $p$. This
is done empirically via testing various numerical solutions of Eq. (1). We
define in this way a mean value $\mu_{u}$ and a standard deviation $s_{u}$ of
$u$ and the corresponding $\mu_{p}$ and $s_{p}$ varying parameters $p$. The
training occurs on a random $u$ and $p$ sampled from normal distribution
defined by given mean values $\mu_{u}$, $\mu_{p}$ and standard deviations
$s_{u}$ and $s_{p}$. Since we use a random number generator to produce dataset
its size is limited only by a period of random number generator that is very
large.

Let us now discuss the structure of the network denoted above as $F(u,p,w)$.
The network includes liner and nonlinear data transformations. The linear one
is done via multiplication of data vectors by a matrix of neuron weights. For
neural networks the usual order of vector-matrix manipulation is the reversed:
Typically we multiply a matrix by a vector-column and in the neural network
context a vector-row is multiplied by a matrix. This is done because in the
course of training a batch of vectors is processed in parallel. A rectangular
data matrix with the vectors stowed in rows is multiplied by a matrix of
weights. Thus we assume that $u$ and $p$ are vector-rows of dimension $N_{u}$
and $N_{p}$ respectively.

The training data vectors $u$ and $p$ are sampled from a normal distribution
and elements of $u$ and $p$ can have different scales. Thus the first
transformation of the network inputs $u$ and $p$ is a non-trainable
normalization layer that rescales inputs to a standard normal distribution.

$\displaystyle\mathop{\mathrm{Norm}}(x)=(x-\mu_{x})/s_{x},$ (7)
$\displaystyle\mathop{\mathrm{Denorm}}(x)=xs_{x}+\mu_{x}.$ (8)

Here $x$, $\mu_{x}$ and $s_{x}$ are vectors-rows and operations are performed
element-wise. Also we define here the layer performing backward operation
$\mathop{\mathrm{Denorm}}()$. It will be done at the very end of the network
to fit the values to an appropriate range. It might be seem that the layers
(7) and (8) are superfluous — one can expect that the network is able to fit
these scales itself in the course of training. But in fact this is not the
case. All network training methods are developed in the assumption that both
inputs and outputs do not deviate much from a standard range. So the training
is efficient if we know in advance what are the ranges of the inputs and the
outputs and rescale them appropriately.

After normalization we concatenate two resulting row vector into the one
vector:

$\mathop{\mathrm{Concat}}(x,y)=(x,y)$ (9)

Here $x$ and $y$ are vectors of $N_{x}$ and $N_{y}$ elements, respectively,
and $(x,y$) is a row vector of $N_{x}+N_{y}$ elements.

The next step is a dense layer. This is mere a affine transformation:

$\mathop{\mathrm{Dense}}(x,N)=xW_{x,N}+b_{x,N}$ (10)

Here $W_{x,N}$ is a rectangular matrix whose number of rows equals to the
number of columns of $x$ and the number of columns of $W_{x,N}$ is $N$,
$b_{x,N}$ is a vector-row with $N$ elements.

After that a nonlinear transformation is applied that is called activation:

$\mathop{\mathrm{Activ}}(x)=\sigma(x)$ (11)

Here $\sigma()$ is a scalar function of a scalar argument and if a vector is
passed to it the element-wise operation is assumed.

Subsequent transformations are done using already defined operators so that
the whole network $d=F(u,p,w)$ can be described as follows:

$\displaystyle
z=\mathop{\mathrm{Concat}}(\mathop{\mathrm{Norm}}(u),\mathop{\mathrm{Norm}}(p))$
(12) $\displaystyle
h=\mathop{\mathrm{Activ}}(\mathop{\mathrm{Dense}}(z,N_{h}))$ (13)
$\displaystyle g=\mathop{\mathrm{Dense}}(h,N_{u})$ (14) $\displaystyle
d=\Delta t\mathop{\mathrm{Denorm}}(g)$ (15)

Variable $w$ in $F(u,p,w)$ represents a set of trainable parameters of the
networks. As follows from the equations above

$w=\left\\{W_{z,N_{h}},b_{z,N_{h}},W_{h,N_{u}},b_{h,N_{u}}\right\\}$ (16)

In the very beginning the network parameters $w$ are initialized at random.
Then the training process is performed as follows. We generate an input batch
$\\{U,P\\}$ of $N_{\text{batch}}$ random $u$ and $p$ sampled form a normal
distribution. Here $U$ and $P$ matrices with $N_{\text{batch}}$ rows and their
number of columns are $N_{u}$ and $N_{p}$, respectively. This batch is feed to
the network (12)-(15) and the matrix $D$ with $N_{\text{batch}}$ rows and
$N_{u}$ columns is obtained. Then the input matrix $U$ and the output one $D$
is substituted into (6) to compute an error matrix $E$ of $N_{\text{batch}}$
rows and $N_{u}$ columns. Finally a mean squared error (MSE) is computed for
the elements of $E$ as:

$\ell=\frac{1}{N_{\text{batch}}N_{u}}\sum_{i=1}^{N_{\text{batch}}}\sum_{j=1}^{N_{u}}e_{ij}^{2}$
(17)

This $\ell$ is the loss function for our training. To update the network
parameters a gradient of $\ell$ is computed with respect of each of the
network parameter gathered in $w$, see (16), and then it is used in a gradient
descent step that computes corrections to the network parameters with respect
to the minimization of $\ell$. The simplest version of the gradient descent
step reads

$w\leftarrow w-\gamma\nabla_{w}\ell$ (18)

where the step size scale $\gamma$ is a small parameter controlling the
convergence.

The iteration that starts from a random batch generation and ends after
updating the network parameters via the gradient descent is repeated
$t_{\text{epoch}}$ times. This is considered as an epoch. Notice that usually
the epoch has a bit different meaning. A neural network is trained on a large
unaltered dataset it cannot not be passed to the network at once due to the
lack of a computer memory. In this case the whole dataset is split into
batches (they are also called mini-batches) and they are passed one by one.
The parameter updates are computed for each batch. The optimization method
applied not to the whole dataset at once but to its batches is called
stochastic gradient descent and the epoch ends when all the bathes have shown
to the network. In our case the batches are always generated at random so that
dividing training process into epochs is required only to interrupt the
training and to compute metrics to see the progress of the network
performance.

We use two metrics: the loss function (17) and the mean relative norm error
(MRNE) that is defined as follows:

$m=\frac{1}{N_{\text{batch}}}\sum_{i=1}^{N_{\text{batch}}}\left(\sum_{j=1}^{N_{u}}|e_{ij}|/\sum_{j=1}^{N_{u}}|u_{ij}|\right),$
(19)

where $e_{ij}$ as above are the elements of the error matrix $E$ and $u_{ij}$
are the elements of the network input batch $U$. To estimate the network
performance, after each epoch we perform $t_{\text{valid}}$ validation steps:
Generate a new random batch $\\{U,P\\}$, feed the network (12)-(15), obtain
$E$ and compute $\ell$ (17) and $m$ (19) without network parameters updating;
finally the computed metrics are averaged over the validation steps
$t_{\text{valid}}$. The dependence of the average metrics on the number of
epochs passed is called learning curves.

For actual computations instead of the simplest one (18) we the use more
sophisticated version of the gradient descent method called Adam. The
difference is that the step size scale $\gamma$ is not a constant, but is
tuned according to the accumulated gradients on the previous steps [18]. This
method has a meta-parameter learning rate $\alpha$ that control the overall
scale of the computed step size. We decrease it in the course of the
computations according to the inverse time decay rule:

$\alpha=\frac{0.1}{1+0.96t/(30t_{\text{epoch}})}$ (20)

where $t$ is the gradient descent step, and $t_{\text{epoch}}$ is a number of
steps comprising one epoch. The particular numerical values of the
coefficients in this formula are chosen empirically to provide the fastest
convergence.

At the activation layer $\mathop{\mathrm{Activ}}()$ in Eq. (13) we apply the
sigmoid function

$\sigma(x)=\frac{1}{1+\mathrm{e}^{-x}}$ (21)

We will train the neural network models to achieve the mean relative error
MRNE at level $10^{-5}$.

The transformation that is done by the network under consideration (2),
(12)-(15) can be represented as a map. Normalization operator in Eq. (12) can
be taken into account inside the dense layer in (13) by an appropriate
rescaling and shift of the elements of $W_{z,N_{h}}$ and $b_{x,N_{h}}$.
Similarly, denormalization operator in Eq. (15) can be merged with the dense
layer in (14). Also instead of concatenating the normalized vectors $u$ and
$p$ we split the matrix $W_{z,N_{h}}$ into two blocks corresponding to $u$ and
$p$ respectively. As a result we obtain the following map that models
solutions to Eq. (1):

$u_{n+1}=u_{n}+\sigma(u_{n}A_{0}+pB_{0}+a_{0})A_{1}+a_{1}$ (22)

where $A_{0}$ is a matrix with $N_{u}$ rows and $N_{h}$ columns, $B_{0}$ has
$N_{p}$ rows and $N_{h}$ columns, $A_{1}$ is a matrix with $N_{h}$ rows and
$N_{u}$ columns. Vector-row $a_{0}$ has $N_{h}$ elements and $a_{1}$ has
$N_{u}$ elements.

Equation (22) is a universal model of a solution to ODE (1). Particular system
is selected by choosing an appropriate size $N_{h}$ of the hidden layer and by
numerical values of the elements of matrices $A_{0}$, $B_{0}$, $A_{1}$ and
vectors $a_{0}$ and $a_{1}$.

For Eq. (22) we can find the variational equation suitable for applying to
this system the Lyapunov analysis, in particular, for computing Lyapunov
exponents. Differentiating the elements of $u_{n+1}$ by the elements $u_{n}$
one obtains the Jacobian matrix:

$J_{n}=I+A_{0}H_{n}A_{1}$ (23)

where $I$ is the identity matrix, and $H_{n}$ is a diagonal square matrix
$N_{h}$ by $N_{h}$:

$H_{n}=\mathop{\mathrm{diag}}(h_{n}(1-h_{n}))$ (24)

and $h_{n}$ is a row-vector computed as

$h_{n}=\sigma(u_{n}A_{0}+pB_{0}+a_{0})$ (25)

In the other words it is computted according to Eqs. (12), (13) when $u_{n}$
and $p$ corresponding to the current trajectory point are substituted there.

Thus the variational equation for the system (22) reads:

$\delta u_{n+1}=(I+A_{0}H_{n}A_{1})\delta u_{n}$ (26)

This variational equation can be used to compute Lyapunov exponents. For this
purpose we apply the standard algorithm [19, 20]: Iterate the main system (22)
simultaneously with the required number of copies of the variational Eq. (26)
with periodic orthogonalization and normalization of the set of vectors
$\delta u_{n}$. Accumulated and averaged in time logarithms of the norms of
variational vectors converge to the Lyapunov exponents.

Since the training and running of neural network is usually done in a
multithread computation environment the preferable way of finding of the
exponents is to iterated vary many trajectories simultaneously for not very
large time cuts and then average the resulting exponents over the
trajectories.

All the computations including training and running are preformed using
TensorFlow [5] and CUDA [21] software.

##  4 Models

###  4.1 Lorenz system

First we consider the Lorenz system [22, 23, 24]:

$\begin{gathered}\dot{x}=\sigma(y-x),\\\ \dot{y}=x(r-z)-y,\\\ \dot{z}=xy-
bz\end{gathered}$ (27)

To train this model we choose the vectors of mean $\mu_{u}$ and standard
deviation $s_{u}$ as follows:

$\begin{gathered}\mu_{u}=(0,0,0),\;s_{u}=(10,10,20)\\\
\mu_{p}=(0,0,0),\;s_{p}=(5,20,2)\\\
\mu_{g}=(0,0,0),\;s_{g}=(70,280,110)\end{gathered}$ (28)

These vectors are used in Eq. (12), see also Eq. (7).

The vectors $\mu_{g}$ and $s_{g}$ are computed as mean and standard deviation
of elements of $f(u,p)$ that is the right hand side of Eq. (27) when $u$ and
$p$ are sampled from normal distribution with mean and standard deviations
$\mu_{u}$, $s_{u}$, $\mu_{p}$ and $s_{p}$. In this case the network output
$d$, see Eq. (15) will approximately have the range of $\frac{\Delta
t}{2}[f(u)+f(u+d)]$, see Eqs. (4) and (5).

One more parameter that we need to define is $N_{h}$, the size of the hidden
layer, see (13). We consider different values to check which one is preferred.
Figure 1 shows the learning curves for the Lorenz system. Its panels (a) and
(b) show MSE (17) and MRNE (19), respectively, computed for validation data.
Three cases are shown corresponding to $N_{h}=50$, $100$ and $200$. We see
that all there curves decay that means that the performance of the network
improves. The fastest decay is observed for $N_{h}=200$. In what follows we
will consider the network with $N_{h}=200$. In the course of the training
after each epoch we compare the attained MRNE level with the previously
smallest. And if the new one is smaller, we save the corresponding network
parameters $w$, see Eq. (16). Up to 20000 epoch we were able to find a network
whose MRNE is approximately $3\times 10^{-5}$. We use this metric as a
criterion of the performance because it is normalized by the dynamical
variables scale so that we can compare the performance of different systems.
Since MRNE is already sufficiently small we did not considered larger values
of $N_{h}$.

Figure 1: Learning curves for the Lorenz network model corresponding to ODEs
(27). The curves are computed at validation steps, i.e., for inputs that did
not used for updating the network parameters. Panel (a) represents MSE, Eq.
(17). Panel (b) shows MRNE, Eq. (19). Different curves correspond to different
$N_{h}$.

The training result is shown in Figs. 2 and 3. Figure 2(a, b) demonstrates the
Lorenz attractor computed for the standard set of parameters $\sigma=10$,
$r=28$, $b=8/3$ using fourth order Runge-Kutta method (a) and the network
model (22) (b). Observe very high coincidence of two plots.

Neural networks architecture is very good suited for parallel computations. So
doing computations with the network model we employ it considering multiple
trajectories at once: to plot Fig. 2(a) via the Runge-Kutta method we compute
$10000$ steps with the time interval $\Delta t=0.01$, while in Fig. 2(b) we
compute $100$ trajectories at once, each of the length $100$ steps $\Delta
t=0.01$.

a) b)

Figure 2: Lorenz attractor computed (a) as a numerical solution of Eqs. (27)
using the forth order Runge-Kutta method and (b) as iteration of the network
model (22). Parameters are $\sigma=10$, $r=28$, $b=8/3$.

Figure 3 shows Fourier spectra computed for $x$ at the parameters $\sigma=10$,
$r=28$, $b=8/3$, panels (a, b), and $\sigma=16$, $r=45.92$, $b=4$, panels (c,
d). Left panels (a) and (c) are computed for the Runge-Kutta data and the
right ones are obtained for the network model. The spectra coincide very well
that indicates that the obtained network (22) models the Lorenz dynamics very
well.

a) b)
c) d)

Figure 3: Fourier spectra of the Lorenz attractor. Data series for panels (a)
and (c) are computed numerically using the forth order Runge-Kutta method, and
data for the panels (b) and (d) are obtained after iterations of the network
model (22). Parameter for the panels (a) and (b) are $\sigma=10$, $r=28$,
$b=8/3$, and panels (c) and (d) are obtained with $\sigma=16$, $r=45.92$,
$b=4$

Now compute Lyapunov exponents using the standard algorithm[19, 20]. Using the
Runge-Kutta method and at $\sigma=10$, $r=28$, $b=8/3$ we obtain the values of
$\lambda_{i}$ in Eq. (29). Lyapunov exponents $\tilde{\lambda}_{i}$ computed
for the network model (22) and corresponding variational equation (26) are in
Eq. (30). Observe the very good coincidence. Notice that $\lambda_{2}$ is
expected to be zero since describes marginally stable perturbations along
trajectories. However actual values in computations are never exact zero.
Their closeness to zero indicates the quality of the computation. In our case
both $\lambda_{2}$ and $\tilde{\lambda}_{2}$ are very small.

$\displaystyle\lambda_{1}$ $\displaystyle=0.906$ $\displaystyle\lambda_{2}$
$\displaystyle=8.26\times 10^{-6}$ $\displaystyle\lambda_{3}$
$\displaystyle=-14.6$ (29) $\displaystyle\tilde{\lambda}_{1}$
$\displaystyle=0.905$ $\displaystyle\tilde{\lambda}_{2}$
$\displaystyle=1.26\times 10^{-5}$ $\displaystyle\tilde{\lambda}_{3}$
$\displaystyle=-14.6$ (30)

Similarly the Lyapunov exponents are computed for the parameters $\sigma=16$,
$r=45.92$, $b=4$. Observe again the very high similarity of $\lambda_{i}$ with
network model exponents $\tilde{\lambda}_{2}$:

$\displaystyle\lambda_{1}$ $\displaystyle=1.50$ $\displaystyle\lambda_{2}$
$\displaystyle=-1.89\times 10^{-5}$ $\displaystyle\lambda_{3}$
$\displaystyle=-22.5$ (31) $\displaystyle\tilde{\lambda}_{1}$
$\displaystyle=1.49$ $\displaystyle\tilde{\lambda}_{2}$
$\displaystyle=4.34\times 10^{-5}$ $\displaystyle\tilde{\lambda}_{3}$
$\displaystyle=-22.7$ (32)

###  4.2 Röessler system

Another system that we consider is the Röessler system [25, 24, 26]:

$\begin{gathered}\dot{x}=-y-z,\\\ \dot{y}=x+ay,\\\
\dot{z}=b+z(x-c),\end{gathered}$ (33)

For this system we choose the following $\mu_{u,p}$ and $s_{u,p}$ and compute
the corresponding $\mu_{g}$ and $s_{g}$:

$\begin{gathered}\mu_{u}=(0,0,0),\;s_{u}=(10,10,10),\\\
\mu_{p}=(0,0,0),\;s_{p}=(10,10,10),\\\
\mu_{g}=(0,0,0),\;s_{g}=(14,101,142).\end{gathered}$ (34)

Figure 4: Learning curves for the Röessler network model corresponding to ODEs
(33).

Figure 4 demonstrates the learning curves for the Röessler system. We observe
that the training now goes much faster then for the Lorenz system, see Fig. 1.
Inspecting the learning curves we can conclude that the network models with
$N_{h}=100$ and $200$ do not differ much. So, unlike the Lorenz system we will
consider the network model for the Röessler system with $N_{h}=100$. After
20000 epochs of the training we save the model with the smallest MRNE equals
to approximately $1.0\times 10^{-5}$.

Figure 5 demonstrates chaotic Röessler attractor and the corresponding Fourier
spectra computed for ODEs (33) (left column) and for the trained network model
(right column). Numerical solutions of ODEs here and below are obtained using
the forth order Runge-Kutta method. A very high similarity of the graphs
indicates the high quality of approximation of the network model. Another
example of dynamics is in Fig. 6. Parameters here correspond to the period 2
oscillations. Limit cycles in Figs. 6(a) and (c) looks almost identical. The
spectrum for the network model in Figs. 6(d) also repeats the spectrum in Fig.
6(b) in location and relative heights of harmonics. The difference between
these two spectra is in small fluctuations. Since the regime of the considered
system is periodic the fluctuations are mere artifacts related in particular
with the computation method. The methods of computations are different and so
the fluctuations are.

a) b) |  c) d)
---|---
Figure 5: Röessler chaotic attractor at $a=0.15$, $b=0.2$, and $c=10$, panels (a) and (c) and the corresponding Fourier spectra, panels (b) and (d). Panels (a) and (b): numerical solution of ODEs (33). Panels (c) and (d): iterations of a trained network model. a) b) |  c) d)
---|---
Figure 6: Period 2 oscillations of the Röessler system at $a=0.1$, $b=0.1$,
and $c=6$. Limit cycle (a) and Fourier spectrum (b) are computed for ODEs and
the corresponding panels (c) and (d) are obtained for the network model.

To demonstrate that the trained network model reproduces the dynamics of the
modeled ODEs in a wide range of parameter in Fig. 7 we show a bifurcation
diagram for the Röessler system. Parameters $a$ and $b$ are fixed and $c$ is
varying. For each $c$ we compute a trajectory then find its Poincaré section
at $x=0$. Absent values of variables between the time discretization points
are obtained via linear interpolation. The diagram obtained for the numerical
solution of (33), see Fig. 7(a) is reproduced very well by the network model,
see Fig. 7(b). Notice however that the bifurcation points for the network
model are little bit shifted to the right. Nevertheless the overall
correspondence is very high.

a)
b)

Figure 7:  Bifurcation diagrams for the Röessler system at $a=0.1$ and
$b=0.1$. Panel (a) corresponds to a numerical solution of ODEs (33) and panel
(b) is computed for the network model. Bifurcation diagrams are obtained as
$y$ values at Poincaré sections at $x=0$. The sections are computed for
linearly interpolated time series.

Let us now compare Lyapunov exponents applying the standard algorithm for ODEs
(33) and for the corresponding network model. We demonstrate two cases. For
parameters $a=0.15$, $b=0.2$, $c=10$ the Lyapunov exponents $\lambda_{i}$ for
ODEs are shown in Eq. (35). For comparison the exponents $\tilde{\lambda}_{i}$
for the corresponding network model at the same parameters are shown in Eq.
(36). The values coincide very well. Because the considered system is
autonomous the value of $\lambda_{2}$ must be zero. Actually computed values
are indeed very close to zero.

$\displaystyle\lambda_{1}$ $\displaystyle=0.0886$ $\displaystyle\lambda_{2}$
$\displaystyle=-8.66\times 10^{-7}$ $\displaystyle\lambda_{3}$
$\displaystyle=-9.80$ (35) $\displaystyle\tilde{\lambda}_{1}$
$\displaystyle=0.0839$ $\displaystyle\tilde{\lambda}_{2}$
$\displaystyle=2.70\times 10^{-5}$ $\displaystyle\tilde{\lambda}_{3}$
$\displaystyle=-9.64$ (36)

One more example is considered at $a=0.1$, $b=0.1$, $c=13$ for that the
Röessler systems also has a chaotic attractor. From Eqs. (37) and (38) we
again observe that the exponents for the network model $\tilde{\lambda}_{i}$
are close to those obtained for the numerical solution of ODEs $\lambda_{i}$.

$\displaystyle\lambda_{1}$ $\displaystyle=0.0116$ $\displaystyle\lambda_{2}$
$\displaystyle=1.87\times 10^{-5}$ $\displaystyle\lambda_{3}$
$\displaystyle=-12.8$ (37) $\displaystyle\tilde{\lambda}_{1}$
$\displaystyle=0.0189$ $\displaystyle\tilde{\lambda}_{2}$
$\displaystyle=8.53\times 10^{-5}$ $\displaystyle\tilde{\lambda}_{3}$
$\displaystyle=-12.8$ (38)

However we must notice that the correspondence of the Lyapunov exponents for
the Röessler system is not so good as for the Lorenz system, see. Eqs. (29)-
(32). We address it to the parameter mismatch observed in the bifurcation
diagrams.

###  4.3 Hindmarch–Rose neuron

Now we consider the Hindmarsh–Rose model of neuronal activity [27, 28]:

$\begin{gathered}\dot{x}=y-ax^{3}+bx^{2}-z+I,\\\ \dot{y}=c-dx^{2}-y,\\\
\dot{z}=r(s(x-\alpha)-z)\end{gathered}$ (39)

Totally this system has eight parameters. However often the system is
considered when six of them have standard values: $a=1.0$, $b=3.0$, $c=1.0$,
$d=5.0$, $s=4.0$, $\alpha=-1.6$. Parameters $I$ and $r$ are varied.

The Hindmarch–Rose model (39) is a simplified model for biological neurons
presenting bursting oscillations. In this regime, bursts of fast spikes are
followed by quiescent periods. Typical values of parameters where the bursts
are observed are $I=2.7$ and $r=0.003$. Thus we select the normalization to be
close to these values:

$\begin{gathered}\mu_{u}=(0,-5,2.5),\;s_{u}=(0.8,2.5,0.5),\\\
\mu_{p}=(0.012,2.7),\;s_{p}=(0.024,0.3),\\\
\mu_{g}=(-2.9,2.8,0.047),\;s_{g}=(4.2,5.2,0.13)\end{gathered}$ (40)

Learning curves for the network model of the system (39) are shown in Fig. 8.
Unlike two previous cases the training occurs much faster: it takes 7000
epochs for the loss function at $N_{h}=100$ to reach values about $10^{-10}$.
The model at $N_{h}=50$ also demonstrates a very good convergence, and the
model at $N_{h}$ behave much worse. Thus we will consider a model with
$N_{h}=50$.

Figure 8: Learning curves for the Hindmarch-Rose network model corresponding
to ODEs (39).

Figure 9(a, b) demonstrates typical solutions of the Hindmarch-Rose model
(39): Panel (a) demonstrates periodical the bursts of spikes and in the panel
(b) we see chaotic spikes. The system contains fast and slow variables, i.e.,
is stiff. To solve it numerically the method LSODA is used [29].

Figures 9(c, d) shows the corresponding time series obtained from the neural
network model. The behavior of the network model is very similar, but the
close inspection revels that in the panel (c) there are seven spikes in each
burst, while the “original” curve contains only six of them. It means that
although the model demonstrates a neural dynamics as well as in original ODEs,
its parameters do not coincide exactly. The chaotic regimes in the panels (b)
and (d) obviously represent the same regime.

a) b)

c) d)

Figure 9: Time series of $x(t)$ obtained as numerical solution of ODEs (39),
panels (a) and (b), and corresponding iterations of the neural network model,
panels (c) and (d). Parameters are $r=0.003$, $I=2.7$ for panels (a) and (d),
and $r=0.013$, $I=2.9$ in panels (b) and (d).

In Fig. 10(a) and (b) bifurcation diagrams provides a more detailed comparison
of the neural network model with ODEs. In both panels the diagrams are
computed for Poincaré sections at $x=0$ computed for linearly interpolated
times series. The diagrams have similar global structure. One can see areas of
bursts in their left parts and chaotic areas to the right. However the
detailed arrangement is different. The diagram for the neural network model
looks less regular along parameter axis. Often changes of the regimes are
observed.

a)
b)

Figure 10:  Bifurcation diagrams for the Röessler system at $a=0.1$ and
$b=0.1$. Panel (a) corresponds to a numerical solution of ODEs (33) and panel
(b) is computed for the network model. Bifurcation diagrams are obtained as
$y$ values at Poincaré sections at $x=0$. The sections are computed for
linearly interpolated time series.

Finally we compare Lyapunov exponents computed for numerical solutions of ODEs
(39) and for iterations of the neural network model. For chaotic spikes at
$r=0.013$, $I=2.9$, see Fig. 9(b), the Lyapunov exponents $\lambda_{i}$ are
given by Eq. (41). The corresponding exponents for the network model
$\tilde{\lambda}_{i}$ are gathered in Eq. (42). The exponents are pairwise
close but do not coincide.

$\displaystyle\lambda_{1}$ $\displaystyle=8.39\times 10^{-3}$
$\displaystyle\lambda_{2}$ $\displaystyle=-1.32\times 10^{-5}$
$\displaystyle\lambda_{3}$ $\displaystyle=-9.55$ (41)
$\displaystyle\tilde{\lambda}_{1}$ $\displaystyle=7.95\times 10^{-3}$
$\displaystyle\tilde{\lambda}_{2}$ $\displaystyle=9.77\times 10^{-6}$
$\displaystyle\tilde{\lambda}_{3}$ $\displaystyle=-9.63$ (42)

Similar situations is obtained for another parameter values $r=0.012$,
$I=2.7$: the exponents $\lambda_{i}$ and $\tilde{\lambda}_{i}$ are similar,
but the difference is notable.

$\displaystyle\lambda_{1}$ $\displaystyle=5.41\times 10^{-3}$
$\displaystyle\lambda_{2}$ $\displaystyle=7.59\times 10^{-6}$
$\displaystyle\lambda_{3}$ $\displaystyle=-10.2$ (43)
$\displaystyle\tilde{\lambda}_{1}$ $\displaystyle=3.20\times 10^{-3}$
$\displaystyle\tilde{\lambda}_{2}$ $\displaystyle=3.04\times 10^{-5}$
$\displaystyle\tilde{\lambda}_{3}$ $\displaystyle=-10.3$ (44)

We address this not very good coincides of the Lyapunov exponents to the very
weak chaos. The first exponents are very small by magnitude so that the
numerical routine converges poorly and is strongly affected by numerical
errors. The parameter mismatch observed above for the Röessler system also
includes in the not very good coincides of the Lyapunov exponents.

Thus, we observe that the discussed neural network model for the Hindmarch-
Rose system provides good qualitative approximation of this system, however
the quantitative correspondence is not high.

##  5 Conclusion

We discussed the universal neural network, a perceptron with one hidden level,
that can be trained to model behavior of various dynamical systems given by
ODEs. Mathematically the universal neural network model is a discrete time
system, see (22). We aware of contemporary success in using of so called deep
networks. Our network on contrary is not deep. We have preferred it because
there is a rigorous mathematical evidence, the Universal Approximation
Theorem, that the network with such architecture is able to approximate
various dependencies. Another reason to apply a classical perceptron is its
simple structure. We believe that it will help to trigger new theoretical
studies of dynamical systems. From the practical point of view this simple
network can be effectively simulated using so called AI accelerators, a
hardware dedicated to deal with artificial neural networks. The approach
developed in this paper can be considered as an alternative numerical method
of modeling dynamical system that is able to utilize contemporary parallel
hard- and software.

The universal network model was trained to reproduce the dynamics of the three
systems: Lorenz and the Röessler systems and Hindmarch-Rose model. It was very
successful for the Lorenz system. This is confirmed by visual inspection of
attractors, and by coincidence of Fourier spectra and Lyapunov exponents. For
the Röessler system the correspondence is also high. However a certain
mismatch of the bifurcation points is observed on bifurcation diagrams
computed for the numerical solution of Röessler ODEs and for the network
model.

For the Hindmarch-Rose system good qualitative correspondence is achieved
however quantitative characteristics are sometimes differ. This system allows
to reveal the limitations of the suggested approach. Hindmarch-Rose system is
stiff and also its regimes changes fast withing a narrow range parameters.
Probably for such cases like this system a more subtle approach is required.

Work of PVK on theoretical formulation and numerical computations and work of
NVS on results analysis was supported by grant of Russian Science Foundation
No 20-71-10048.

## References

  * [1] Simon Haykin.  Neural networks and learning machines.  Pearson Prentice Hall, third edition, 2009.
  * [2] G. Peter Zhang.  Neural networks for time-deries forecasting, pages 461–477.  Springer Berlin Heidelberg, Berlin, Heidelberg, 2012.
  * [3] Nigel Da Costa Lewis.  Deep time series forecasting with Python: An intuitive introduction to deep learning for applied time series modeling.  ND Lewis, 2016.
  * [4] Jason Brownlee.  Deep learning for time series forecasting: predict the future with MLPs, CNNs and LSTMs in Python.  Machine Learning Mastery, 2018.
  * [5] Martín Abadi, Ashish Agarwal, Paul Barham, Eugene Brevdo, Zhifeng Chen, Craig Citro, Greg S. Corrado, Andy Davis, Jeffrey Dean, Matthieu Devin, Sanjay Ghemawat, Ian Goodfellow, Andrew Harp, Geoffrey Irving, Michael Isard, Yangqing Jia, Rafal Jozefowicz, Lukasz Kaiser, Manjunath Kudlur, Josh Levenberg, Dandelion Mané, Rajat Monga, Sherry Moore, Derek Murray, Chris Olah, Mike Schuster, Jonathon Shlens, Benoit Steiner, Ilya Sutskever, Kunal Talwar, Paul Tucker, Vincent Vanhoucke, Vijay Vasudevan, Fernanda Viégas, Oriol Vinyals, Pete Warden, Martin Wattenberg, Martin Wicke, Yuan Yu, and Xiaoqiang Zheng.  TensorFlow: Large-scale machine learning on heterogeneous systems, 2015.  Software available from tensorflow.org.
  * [6] Adam Paszke, Sam Gross, Soumith Chintala, Gregory Chanan, Edward Yang, Zachary DeVito, Zeming Lin, Alban Desmaison, Luca Antiga, and Adam Lerer.  Automatic differentiation in PyTorch.  In NIPS-W, 2017.
  * [7] Y. Wei, J. Zhou, Y. Wang, Y. Liu, Q. Liu, J. Luo, C. Wang, F. Ren, and L. Huang.  A review of algorithm & hardware design for AI-based biomedical applications.  IEEE Transactions on Biomedical Circuits and Systems, 14(2):145–163, 2020.
  * [8] Manar Abu Talib, Sohaib Majzoub, Qassim Nasir, and Dina Jamal.  A systematic literature review on hardware implementation of artificial intelligence algorithms.  The Journal of Supercomputing, 77:1897–1938, 2021.
  * [9] Morris Kline.  Mathematical thought from ancient to modern times: Volume 2, volume 2.  Oxford university press, 1990.
  * [10] A. N. Kolmogorov.  On the representation of continuous functions of several variables by superpositions of continuous functions of a smaller number of variables.  Doklady Akademii Nauk SSSR, 108:179–182, 1956.  English translation: Amer. Math. Soc. Transl., 17 (1961), pp. 369-373.
  * [11] A. N. Kolmogorov.  On the representation of continuous functions of many variables by superposition of continuous functions of one variable and addition.  Doklady Akademii Nauk SSSR, 114:953–956, 1957.  English translation: Amer. Math. Soc. Transl., 28 (1963), pp. 55–59.
  * [12] V. I. Arnold.  On functions of three variables.  Doklady Akademii Nauk SSSR, 114:679–681, 1957.  English translation: Amer. Math. Soc. Transl., 28 (1963), pp. 51–54.
  * [13] R. Hecht-Nielson.  Kolmogorov’s mapping neural network eistence theorem.  In First IEEE international conference on neural networks, volume III, pages 11–14. San Diego, CA, 1987.
  * [14] G. Cybenko.  Approximation by superpositions of a sigmoidal function.  Mathematics of Control, Signals and Systems, 2(4):303–314, December 1989.
  * [15] K. Funahashi.  On the approximate realization of continues mappings by neural networks.  Neural Networks, 2:183–192, 1989.
  * [16] Kurt Hornik, Maxwell Stinchcombe, and Halbert White.  Multilayer feedforward networks are universal approximators.  Neural Networks, 2(5):359–366, 1989.
  * [17] W. Light.  Ridge functions, sigmoidal functions and neural networks, volume VII of Approximation Theory, pages 163–206.  Academic Press, Boston, 1992.
  * [18] Diederik P. Kingma and Jimmy Ba.  Adam: A method for stochastic optimization.  arXiv preprint arXiv:1412.6980, 2014.  Published as a conference paper at International Conference on Learning Representations (ICLR) 2015.
  * [19] Giancarlo Benettin, Luigi Galgani, Antonio Giorgilli, and Jean-Marie Strelcyn.  Lyapunov characteristic exponents for smooth dynamical systems and for Hamiltonian systems: a method for computing all of them. Part 1: Theory.  Meccanica, 15(1):9–20, 1980.
  * [20] Ippei Shimada and Tomomasa Nagashima.  A numerical approach to ergodic problem of dissipative dynamical systems.  Prog. Theor. Phys., 61(6):1605–1616, 1979.
  * [21] John Nickolls, Ian Buck, Michael Garland, and Kevin Skadron.  Scalable parallel programming with CUDA.  ACM Queue, 6(2):40–53, 2008.
  * [22] E. N. Lorenz.  Deterministic nonperiodic flow.  Journal of the atmospheric sciences, 20(2):130–141, 1963.
  * [23] C. Sparrow.  The Lorenz equations: bifurcations, chaos, and strange attractors.  Springer-Verlag, NY, Heidelberg, Berlin, 1982.
  * [24] H. G. Schuster and W. Just.  Deterministic chaos: an introduction.  Wiley-VCH, 2005.
  * [25] O. E. Rössler.  An equation for continuous chaos.  Physics Letters A, 57(5):397–398, 1976.
  * [26] S. P. Kuznetsov.  Dynamical chaos.  Moscow: Fizmatlit, 2006.
  * [27] J. L. Hindmarsh and R. M. Rose.  A model of neuronal bursting using three coupled first order differential equations.  Proc. R. Soc. Lond. B, 221:87–102, 1984.
  * [28] X.-J. Wang.  Genesis of bursting oscillations in the Hindmarsh-Rose model and homoclinicity to a chaotic saddle.  Physica D: Nonlinear Phenomena, 62(1):263–274, 1993.
  * [29] L. Petzold.  Automatic selection of methods for solving stiff and nonstiff systems of ordinary differential equations.  SIAM Journal on Scientific and Statistical Computing, 4(1):136–148, 1983.