# On the stability of second order parametric ordinary differential equations
and applications

Z. Mazgouri1,∗, A. El Ayoubi2

1Applied Mathematics Engineering Department, National School of Applied
Sciences, Sidi Mohamed Ben Abdellah University, Laboratory LAMA, 30000 Fez,
Morocco
2Department of Mathematics, Faculty of Sciences, Chouaib Doukkali University,
El Jadida Morocco

Abstract. This work deals with Lipschitz stability for a parametric version of
the general second order Ordinary Differential Equation (ODE) initial-value
Cauchy problem. We first establish a Lipschitz stability result for this
problem under a partial variation of the data. Then, we apply our abstract
result to second order differential equations governed by cocoercive
operators. Furthermore, we discuss more concrete applications of the stability
for two specific applied mathematical models inherent in electricity and
control theory. Finally, we provide numerical tests based on the software
source Scilab, which are done with respect to parametric linear control
systems, illustrating henceforth the validity of our abstract theoretical
result.

Keywords. Cauchy problem; Differential equation initial-value problem;
Parametric perturbation; Lipschitz stability; Perov inequality.

2020 Mathematics Subject Classification. 47D09, 34A12, 34D10, 34D20.

††footnotetext:  ∗Corresponding author. E-mail addresses:
<EMAIL_ADDRESS>(Z. Mazgouri<EMAIL_ADDRESS>(A. El
Ayoubi). Received ; Accepted

##  1\. Introduction

Throughout this paper, unless otherwise is specified, the Euclidian space
$\mathbb{R}^{n}$ is equipped with the supremum norm $\|.\|.$ We denote by
$B(x_{0},r)$ the open ball of radius $r>0$ centered at
$x_{0}\in\mathbb{R}^{n}$. Let $T>0$ and $\dot{x}_{0}\in\mathbb{R}^{n}$ and let
$f:[0,T]\times B(x_{0},r)\times B(\dot{x}_{0},r)\longrightarrow\mathbb{R}^{n}$
be a given function. We consider the following second order ordinary
differential equation (ODE) initial-value Cauchy problem:

$S(f,x_{0},\dot{x}_{0})\begin{cases}\ddot{x}(t)=f(t,x(t),\dot{x}(t)),\,\,\mbox{for}\,\,a.e\,\,t\in[0,T]\\\
x(0)=x_{0}\\\ \dot{x}(0)=\dot{x}_{0},\end{cases}$ (1.1)

where $T$ stands for the final time of the interval of interest and
$(x_{0},\dot{x}_{0})\in\mathbb{R}^{2n}$ is the initial-value condition.

It is well known that the mathematical formulation of second order ODE plays a
central role in modelling a variety of interesting problems inherent in
optimal control systems, mechanics and quantum mechanics, dynamical systems
subject to classical mechanical laws, damped simple harmonic motion, the
motion of a particle in one dimensional box besides standard concrete
applications to electric circuits, the charge of capacitors in related RLC
series and many other applied models of a particular relevance well described
in the literature as for example in [7, 17, 21] and the references quoted
there.

For these reasons, the classic ODE problems always attract in the widest sense
the attention of large community of scientists not only from mathematics and
physics or else economic disciplines but also in all scientific fields as well
as engineers and highly qualified technicians. In fact, the literature on the
topic available today presents a plethora of contributions which fundamentally
regard the existence theory, asymptotic analysis in time and numerical
approximation procedures by the use of several methods such as fixed point
approaches, discretization, regularization procedures, equilibrium techniques
and others, see [14, 18, 20] and the references therein for a detailed
presentation on these elements.

A further interesting research direction for ODE problems regards the study of
the impact of parametric perturbations on the corresponding solutions. Such
parameters could be of an internal type steaming from endogen properties of
the considered systems or of an external/exogen one coming from an independent
source, which always exposes trajectories of the systems under consideration
to be a subject to variations. In this context, the challenge is to find
conditions trough which those perturbations will not affect the objective of
the systems under consideration. For such purposes, we refer to the very nice
contributions of [1, 3, 6] on the mathematical management of different
perturbations within the dynamic feature of the considered differential
equation in the same illuminating stability techniques of [2, 4, 5, 8, 9, 10,
11].

Let us underline the interest of the quantitative stability analysis for
general dynamic problems. It, actually, allows to compute the distance between
the solution of the initial system and the one of its perturbed counterpart,
which is doubly meaningful since on the one hand ensures the qualitative
convergence of solutions with respect to parameters at any time, and on the
other hand provides Lipschitz estimates of trajectories highly required in
auxiliary steps arising in the corresponding algorithms for such evolution
problems. Motivated by these numerous factors, the authors of [7], being
inspired by [8, 9], presented quantitative stability results for first order
ODE parametric Cauchy problem, which have been proved to be extendable to
second order framework as confirmed in the more recent paper [3], wherein the
authors considered second order differential equations governed by maximally
monotone operators by using two different approaches. In our opinion, such
important stability topics still deserve further attention due to the
increasing involvement of different old and/or new perturbations in the
systems modelled by ODE problems.

Accordingly, this paper could be viewed as a continuation of [3] and
constitutes by the meantime a natural extension of the results in [7]. Thus,
being very encouraged by the aforementioned motivations and practical
applications, we first investigate parametric stability of the system (1.1)
under perturbation of the function of the second member of the problem and the
initial conditions as well.

The perturbation of our purpose involves an external parameter $\lambda$ that
belongs to another space different from the one of the state variable of the
system. Precisely, $\lambda$ is localized in some open subset $U$ of a some
normed space $(\Lambda,\|\,.\,\|)$. In this way, we denote by
$\overline{\lambda}\in\Lambda$, the reference value of the parameter
$\lambda$, that is
$f_{\overline{\lambda}}=f,\,\,x_{\overline{\lambda}}=x,\,\,x_{\overline{\lambda}}(0)=x_{0}\,\,\,\mbox{and}\,\,\,\dot{x}_{\overline{\lambda}}(0)=\dot{x}_{0}$,
and consider a neighborhood $\mathcal{V}(\overline{\lambda})\subset\Lambda$ of
$\overline{\lambda}$. Thus, for each
$\lambda\in\mathcal{V}(\overline{\lambda})$, the second order parametric
Cauchy problem we consider is stated as follows:

$S(f_{\lambda},x_{0,\lambda},\dot{x}_{0,\lambda})\begin{cases}\ddot{x}_{\lambda}(t)=f_{\lambda}(t,x_{\lambda}(t),\dot{x}_{\lambda}(t)),\,\,\mbox{for}\,\,a.e\,\,t\in[0,T]\\\
x_{\lambda}(0)=x_{0,\lambda}\\\
\dot{x}_{\lambda}(0)=\dot{x}_{0,\lambda},\end{cases}$ (1.2)

where $f_{\lambda}:[0,T]\times B(x_{0},r)\times
B(\dot{x}_{0},r)\longrightarrow\mathbb{R}^{n}$ is the perturbed function.

In this setting, we seek the conditions under which, for each $\lambda$ near
$\overline{\lambda},$ the trajectory $x_{\lambda}:=x(\lambda)$ is Lipschitz
continuous, that is the following Lipschitz type estimate holds: for all
$\lambda\in\mathcal{V}(\overline{\lambda})$ and for a.e $t\in[0,T]$:

$\displaystyle\|x_{\lambda}(t)-x(t)\|\leq
c_{1}(t)\|x_{0,\lambda}-x_{0}\|+c_{2}(t)\|\dot{x}_{0,\lambda}-\dot{x}_{0}\|+c_{3}(t)\|\lambda-\overline{\lambda}\|,$
(1.3)

where $c_{1},c_{2},c_{3}$ are positive functions to be determined.

Our second concern in this paper is the application of our quantitative
stability on the abstract formulation of (ODE) in (1.2) to concrete examples
quoted above. In this sense, our attention is restricted to control systems
and electric circuits though our contribution demonstrates indeed the
applicability to a large class of similar or other problems that fit in the
second order ODE framework.

In what follows we give an outline of the paper. In Section 2, we present the
ingredients and tools we need for establishing our main results. In Section 3,
based on direct computations with the use of a Perov inequality, we formulate
and prove our principal abstract result on Lipschitz stability for the
parametric system (1.2) (see, Theorem 3.1). In Subsection 4.1, we first apply
this result to second order dynamical systems governed by cocoercive operators
(see, Theorem 4.1). Afterwards, in subsection 4.2, a further interesting
applied model is discussed with respect to differential equations modelling
RLC circuit’s current (see, Theorem 4.2). Furthermore, in subsection 4.3,
based on our main result of Theorem 3.1, we investigate sensitivity analysis
for second order optimal control systems (see, Theorem 4.3). Finally, in
Section 5, we report a numerical example to support and illustrate the
theoretical result given in Theorem 4.3.

##  2\. Preliminaries

Let $T>0$ and $x_{0}\in\mathbb{R}^{n}$. For a real-valued function
$g:[0,T]\times B(x_{0},r)\longrightarrow\mathbb{R}^{n}$, we first consider the
associated Cauchy problem of initial-value first order ODE which is expressed
as follows:

$S(g,x_{0})\begin{cases}\dot{x}(t)=g(t,x(t)),\,\,\mbox{for}\,\,a.e\,\,t\in[0,T]\\\
x(0)=x_{0}.\end{cases}$

Let $U$ be an open subset of $\mathbb{R}^{n}$. We recall that a function
$g:[0,T]\times U\longrightarrow\mathbb{R}^{n}$ is said to be
$L^{1}$-Carathéodory if the following conditions are satisfied:

  1. (i)

for each $x\in U,$ the map $t\longmapsto g(t,x)$ is measurable;

  2. (ii)

for a.e $t\in[0,T]$, the map $x\longmapsto g(t,x)$ is continuous;

  3. (iii)

there exists a function $m\in L^{1}([0,T],\mathbb{R}^{+})$ such that for all
$x\in U$ and for a.e $t\in[0,T]$

$\|g(t,x)\|\leq m(t)$.

The following two results give existence and uniqueness of the solution for a
Cauchy problem of the above form in the Carathéodory sense.

######  Theorem 2.1 ([15, Theorem 1]).

Let $T>0,$ $x_{0}\in\mathbb{R}^{n}$ and let $g:[0,T]\times
B(x_{0},r)\longrightarrow\mathbb{R}^{n}$ be a $L^{1}$-Carathéodory function.
Then, for any real number $d$ such that $0<d\leq T$ and
$\int_{0}^{d}m(s)ds\leq r,$ the Cauchy problem $S(g,x_{0})$ admits a unique
solution on $[0,d].$

######  Theorem 2.2 ([15, Theorem 2]).

Assume that there exists an integrable function $l$ such that

$\|g(t,x)-g(t,y)\|\leq l(t)\|x-y\|\quad\mbox{for
all}\,\,(t,x),(t,y)\in[0,T]\times B(x_{0},r).$

Then, the Cauchy problem $S(g,x_{0})$ admits at most one solution on
$[0,T]\times B(x_{0},r).$

Let us now return to the problem (1.1) subject to our treatment in the present
contribution. We make use of the change of variable $u:=(x,\dot{x})$ and
define the corresponding function $\tilde{f}(t,x):=(\dot{x},f(t,x,\dot{x}))$.
With this change of variable, by augmenting the dimension of the state
variable, the second order ODE (1.1) becomes a first order ODE as follows:

$S(\tilde{f},u_{0})\;\begin{cases}\dot{u}(t)=\tilde{f}(t,u(t)),\,\,\mbox{for}\,\,a.e\,\,t\in[0,T],\,\,u\in\mathbb{R}^{2n}\\\
u(0)=u_{0}.\end{cases}$ (2.1)

Consequently, the next result is straight from Theorem 2.1.

######  Corollary 2.1.

Let $f:[0,T]\times B(x_{0},r)\times
B(\dot{x}_{0},r)\longrightarrow\mathbb{R}^{n}$ be a function that satisfies
the assumptions below:

  1. (1)

for every $(x,y)\in\mathbb{R}^{2n}$, $f(.,x,y)$ is measurable;

  2. (2)

for a.e $t\in[0,T]$, $f(t,.,.)$ is continuous;

  3. (3)

for every $r>0$ there exists $g_{r}\in L^{1}([0,T],\mathbb{R}^{+})$ such that
whenever $\|x-x_{0}\|\leq r$ and $\|y-\dot{x}_{0}\|\leq r$, it holds

$\|f(t,x,y)\|\leq g_{r}(t)$, for a.e $t\in[0,T]$.

Then, for any real number $d$ such that $0<d\leq T$ and
$\int_{0}^{d}g_{r}(s)ds\leq r,$ the Cauchy problem (1.1) admits a unique
solution on $[0,d].$

We also need the first part of Theorem 1 on page 360 of the book [22], which
is the statement of Perov’s inequality.

######  Theorem 2.3 ([22, Theorem 1]).

Let $u(.)$ be a nonnegative function that satisfies the integral inequality

$u(t)\leq c+\int_{t_{0}}^{t}(a(s)u(s)+b(s)u^{\alpha}(s))ds,\quad c,\alpha\geq
0,$ (2.2)

where $a(t)$ and $b(t)$ are continuous nonnegative functions for $t\geq
t_{0}$.
Then, for $0\leq\alpha<1$ we have

$u(t)\leq\left\\{c^{1-\alpha}e^{(1-\alpha)\int_{t_{0}}^{t}a(s)ds}+(1-\alpha)\int_{t_{0}}^{t}b(s)e^{(1-\alpha)\int_{s}^{t}a(r)dr}ds\right\\}^{\frac{1}{1-\alpha}}.$
(2.3)

We close this section with the following lemma that will be needed to prove
our main result of the next section.

######  Lemma 2.1.

Let $f\in C^{2}([0,T],\mathbb{R}^{n})$ such that $f(0)=0$. Then, for each
$t\in[0,T]$, we have

$\int_{0}^{t}\|f(s)\|\|\dot{f}(s)\|ds\leq\frac{t}{2}\int_{0}^{t}\|\dot{f}(s)\|^{2}ds.$

###### Proof.

Consider the function $g$ defined on $[0,T]$ by
$g(t)=\int_{0}^{t}\|\dot{f}(u)\|du$.
We have for every $t\in[0,T]$, $\dot{g}(t)=\|\dot{f}(t)\|$. Fix $t\in[0,T]$,
thus for $s\in[0,t]$, we have

$\|f(s)\|\|\dot{f}(s)\|\leq g(s)\dot{g}(s).$

Integrating this inequality in $[0,t]$ and using the Cauchy Schwarz
inequality, immediately follows

$\int_{0}^{t}\|f(s)\|\|\dot{f}(s)\|ds\leq\int_{0}^{t}g(s)\dot{g}(s)ds=\frac{1}{2}g^{2}(t)\leq\frac{t}{2}\int_{0}^{t}\|\dot{f}(s)\|^{2}ds.$

This completes the proof. ∎

##  3\. Main result

In this section, we consider the perturbed format of the system
$S(f,x_{0},\dot{x}_{0})$ given by (1.2). Based on a Perov’s argument, we
provide conditions under which the Lipschitz estimate in (1.3) is satisfied.

######  Theorem 3.1.

Assume that for some $L>0,\,\,L^{\prime}>0$ the following conditions hold:

  1. (1)

the function $f_{\lambda}$ is continuous;

  2. (2)

the function $f_{\lambda}$ is $L$-Lipschitz in $(x,y)$ uniformly in $t$ and
$\lambda$, that is

$\|f_{\lambda}(t,x,y)-f_{\lambda}(t,x^{\prime},y^{\prime})\|\leq
L(\|x-x^{\prime}\|+\|y-y^{\prime}\|);$

  3. (3)

the function $f_{\lambda}$ is $L^{\prime}$-Lipschitz in $\lambda$ uniformly in
$t$ and $(x,y)$, that is

$\|f_{\lambda}(t,x,y)-f_{\lambda^{\prime}}(t,x,y)\|\leq
L^{\prime}\|\lambda-\lambda^{\prime}\|.$

Let $x_{\lambda}$ $($resp., $x)$ be the solution of the problem
$S(f_{\lambda},x_{0,\lambda},\dot{x}_{0,\lambda})$ $($resp.,
$S(f,x_{0},\dot{x}_{0}))$. Then, for each
$\lambda\in\mathcal{V}(\overline{\lambda})$ and for a.e $t\in[0,T]$, we have

$\begin{split}&\|x_{\lambda}(t)-x(t)\|\leq\left(1+L\left(\frac{2}{2+LT}\right)^{2}\left(e^{\frac{(2+LT)t}{2}}-1-t\right)\right)\|x_{0,\lambda}-x_{0}\|\\\
&\quad\quad+\frac{2}{2+LT}\left(e^{\frac{(2+LT)t}{2}}-1\right)\|\dot{x}_{0,\lambda}-\dot{x}_{0}\|+L^{\prime}\left(\frac{2}{2+LT}\right)^{2}\left(e^{\frac{(2+LT)t}{2}}-1-t\right)\|\lambda-\overline{\lambda}\|.\end{split}$
(3.1)

###### Proof.

Consider the function
$\phi_{\lambda}(t)=\|\dot{x}_{\lambda}(t)-\dot{x}(t)\|^{2}$. We have

$\dot{\phi}_{\lambda}(t)=2\langle\dot{x}_{\lambda}(t)-\dot{x}(t),\ddot{x}_{\lambda}(t)-\ddot{x}(t)\rangle.$

Using Cauchy Schwarz and the triangular inequality we obtain

$\dot{\phi}_{\lambda}(t)\leq
2\|\dot{x}_{\lambda}(t)-\dot{x}(t)\|\left(\|f(t,x(t),\dot{x}(t))-f(t,x_{\lambda}(t),\dot{x}_{\lambda}(t))\|+\|f(t,x_{\lambda}(t),\dot{x}_{\lambda}(t))-f_{\lambda}(t,x_{\lambda}(t),\dot{x}_{\lambda}(t))\|\right).$

Since $f_{\lambda}$ is $L$-Lipschitz in $(x,y)$ uniformly in $t$ and $\lambda$
and $L^{\prime}$-Lipschitz in $\lambda$ uniformly in $t$ and $(x,y)$, we get

$\dot{\phi}_{\lambda}(t)\leq
2L\phi_{\lambda}(t)+2L\|\dot{x}_{\lambda}(t)-\dot{x}(t)\|\|x_{\lambda}(t)-x(t)\|+2L^{\prime}\|\lambda-\overline{\lambda}\|\|\dot{x}_{\lambda}(t)-\dot{x}(t)\|.$

On the other hand, we have

$\|x_{\lambda}(t)-x(t)\|\leq\int_{0}^{t}\|\dot{x}_{\lambda}(s)-\dot{x}(s)\|ds+\|x_{\lambda}(0)-x(0)\|.$
(3.2)

By combining the last two inequalities we obtain

$\dot{\phi}_{\lambda}(t)\leq
2\phi_{\lambda}(t)+2L(\phi_{\lambda}(t))^{\frac{1}{2}}\int_{0}^{t}(\phi_{\lambda}(s))^{\frac{1}{2}}ds+\left(2L\|x_{\lambda}(0)-x(0)\|+2L^{\prime}\|\lambda-\overline{\lambda}\|\right)(\phi_{\lambda}(t))^{\frac{1}{2}}.$
(3.3)

Applying Lemma 2.1, we obtain

$\int_{0}^{t}\left((\phi_{\lambda}(s))^{\frac{1}{2}}\int_{0}^{s}(\phi_{\lambda}(u))^{\frac{1}{2}}du\right)ds\leq\frac{t}{2}\int_{0}^{t}\phi_{\lambda}(s)ds.$

Integrating the above inequality from $0$ to $t$ and using the inequality
(3.3) yields

$\phi_{\lambda}(t)\leq\phi_{\lambda}(0)+2\int_{0}^{t}\phi_{\lambda}(s)ds+Lt\int_{0}^{t}\phi_{\lambda}(s)ds+\left(2L\|x_{\lambda}(0)-x(0)\|+2L^{\prime}|\lambda-\overline{\lambda}\|\right)\int_{0}^{t}(\phi_{\lambda}(s))^{\frac{1}{2}}ds.$

Thus, for all $t\in[0,T]$

$\phi_{\lambda}(t)\leq\phi_{\lambda}(0)+(2+LT)\int_{0}^{t}\phi_{\lambda}(s)ds+\left(2L\|x_{\lambda}(0)-x(0)\|+2L^{\prime}\|\lambda-\overline{\lambda}\|\right)\int_{0}^{t}(\phi_{\lambda}(s))^{\frac{1}{2}}ds.$

Applying Theorem 2.3 with

$t_{0}=0,\alpha=\frac{1}{2},c=\phi_{\lambda}(0),a(.)=2+LT,b(.)=2L\|x_{\lambda}(0)-x(0)\|+2L^{\prime}\|\lambda-\overline{\lambda}\|$
and $u(.)=\phi_{\lambda}(.)$,

we see that

$\phi_{\lambda}(t)\leq\left\\{\phi_{\lambda}(0)^{\frac{1}{2}}e^{\frac{1}{2}\int_{0}^{t}(2+LT)ds}+\left(L\|x_{\lambda}(0)-x(0)\|+L^{\prime}\|\lambda-\overline{\lambda}\|\right)\int_{0}^{t}e^{\frac{1}{2}\int_{s}^{t}(2+LT)dr}ds\right\\}^{2}.$

Thus,

$\|\dot{x}_{\lambda}(t)-\dot{x}(t)\|\leq\|\dot{x}_{\lambda}(0)-\dot{x}(0)\|e^{\frac{(2+LT)t}{2}}+2\frac{L\|x_{\lambda}(0)-x(0)\|+L^{\prime}\|\lambda-\overline{\lambda}\|}{2+LT}(e^{\frac{(2+LT)t}{2}}-1).$
(3.4)

Integrating this inequality from $0$ to $t$ and using the inequality (3.2), we
obtain

$\begin{split}\|x_{\lambda}(t)-x(t)\|&\leq\left(1+L\left(\frac{2}{2+LT}\right)^{2}\left(e^{\frac{(2+LT)t}{2}}-1-t\right)\right)\|x_{\lambda}(0)-x(0)\|\\\
&\quad+\frac{2}{2+LT}\left(e^{\frac{(2+LT)t}{2}}-1\right)\|\dot{x}_{\lambda}(0)-\dot{x}(0)\|+L^{\prime}\left(\frac{2}{2+LT}\right)^{2}\left(e^{\frac{(2+LT)t}{2}}-1-t\right)\|\lambda-\overline{\lambda}\|.\end{split}$

The proof is hence complete. ∎

##  4\. Applications

###  4.1. Second order dynamical systems governed by cocoercive operators

The object of this subsection is to study the sensitivity analysis for a
parametric second-order dynamical system governed by a given maximal monotone
operator $A:H\rightarrow H$ defined on a Hilbert space $H$. For $T>0$, the
problem under consideration is stated as follows:

$\begin{cases}\;\ddot{x}(t)+\gamma\dot{x}(t)+Ax(t)=0,\;\;\mbox{for}\,\,a.e\,\,t\in[0,T]\\\
\;x(0)=u_{0}\hbox{ and }\;\dot{x}(0)=v_{0},\end{cases}$ $S(A,u_{0},v_{0})$

where $\gamma>0$ is a damping coefficient and $u_{0},v_{0}\in H$ are initial
conditions.
We consider here the notion of strong solution of $S(A,u_{0},v_{0})$, that is
$x:[0,T]\rightarrow H$ such that $x$ and $\dot{x}$ are absolutely continuous
on the interval $[0,T]$, $x(0)=u_{0}$, $\dot{x}(0)=v_{0}$ and
$\ddot{x}(t)+\gamma\dot{x}(t)+Ax(t)=0$ for almost every $t\in[0,T]$. For
interesting particular cases which motivate the above dynamical system, we
refer to [3, 12, 13] and the references therein. In our concern, under
perturbation of both the operator $A$ and the initial values, the parametric
form of the above system is stated as follows:

$\quad\begin{cases}\ddot{x}_{\lambda}(t)+\gamma\dot{x}_{\lambda}(t)+A_{\lambda}x_{\lambda}(t)=0,\;\;\mbox{for}\,\,a.e\,\,t\in[0,T]\\\
x_{\lambda}(0)=u_{0,\lambda}\,\,\,\mbox{and}\,\,\,\dot{x}_{\lambda}(0)=v_{0,\lambda}.\end{cases}$
$S(A_{\lambda},u_{0,\lambda},v_{0,\lambda})$

In our sensitivity analysis of $S(A_{\lambda},u_{0,\lambda},v_{0,\lambda})$,
the following two hypotheses will be considered:

  * $(C_{1})$

Cocoercivity assumption: for each $\lambda\in\mathcal{V}(\overline{\lambda}),$
there exists $\alpha_{\lambda}>0$ such that $\forall u,v\in H$

$\langle
A_{\lambda}(u)-A_{\lambda}(v),u-v\rangle\geq\alpha_{\lambda}\|A_{\lambda}(u)-A_{\lambda}(v)\|^{2};$

  * $(C_{2})$

Lipschitz property with respect to parameter $\lambda$: there exists
$L^{\prime}>0$ such that for all $u\in H$ and all
$\lambda,\lambda^{\prime}\in\mathcal{V}(\overline{\lambda}),$

$\|A_{\lambda}(u)-A_{\lambda^{\prime}}(u)\|\leq
L^{\prime}\|\lambda-\lambda^{\prime}\|.$

######  Theorem 4.1.

Assume that conditions $(C_{1})-(C_{2})$ are satisfied. Then, for each
$\lambda\in\mathcal{V}(\overline{\lambda})$, the problem
$S(A_{\lambda},u_{0,\lambda},v_{0,\lambda})$ admits a unique solution
$x_{\lambda}.$ Furthermore, if $\frac{1}{\alpha_{\lambda}}<\gamma$ for all
$\lambda\in\mathcal{V}(\overline{\lambda})$, then for a.e $t\in[0,T]$, the
following Lipschitz estimate holds:

$\begin{split}&\|x_{\lambda}(t)-x(t)\|\leq\left(1+\gamma\left(\frac{2}{2+\gamma
T}\right)^{2}\left(e^{\frac{(2+\gamma
T)t}{2}}-1-t\right)\right)\|u_{0,\lambda}-u_{0}\|\\\ &\quad+\frac{2}{2+\gamma
T}\left(e^{\frac{(2+\gamma
T)t}{2}}-1\right)\|v_{0,\lambda}-v_{0}\|+L^{\prime}\left(\frac{2}{2+\gamma
T}\right)^{2}\left(e^{\frac{(2+\gamma
T)t}{2}}-1-t\right)\|\lambda-\overline{\lambda}\|.\end{split}$ (4.1)

###### Proof.

Clearly, the formulation in the problem
$S(A_{\lambda},u_{0,\lambda},v_{0,\lambda})$ may be expressed as (1.1) by
considering the function $f_{\lambda}$ defined, for all
$\lambda\in\mathcal{V}(\overline{\lambda})$ and all $x,y\in H$, by
$f_{\lambda}(t,x,y)=-A_{\lambda}x-\gamma y$. Owing to the fact that the
$\alpha_{\lambda}$-cocoercivity of $A_{\lambda}$ in $(C_{1})$ implies that
$A_{\lambda}$ is $\frac{1}{\alpha_{\lambda}}$-Lipschitz continuous, the
existence and uniqueness of a strong solution to
$S(A_{\lambda},u_{0,\lambda},v_{0,\lambda})$ can be shown via the classical
Cauchy-Lipschitz-Picard theorem (see, [16]) by rewriting this system as a
first order dynamical system in the product space $H\times H$.

Now, for proving the above estimate, we will make use of Theorem 3.1.
Combining the $\frac{1}{\alpha_{\lambda}}$-Lipschitz continuity of
$A_{\lambda}$ with the hypothesis $\frac{1}{\alpha_{\lambda}}<\gamma$ for all
$\lambda\in\mathcal{V}(\overline{\lambda})$, one can conclude that the
function $f_{\lambda}$ is $\gamma$-Lipschitz in $(x,y)$ uniformly in $t$ and
$\lambda$. On the other hand, it is easy to check that the condition $(C_{2})$
implies that $f_{\lambda}$ is $L^{\prime}$-Lipschitz in $\lambda$ uniformly in
$t$ and $(x,y).$ Consequently, assumptions of Theorem 3.1 are all satisfied.
Applying this theorem, we immediately obtain the required estimate. ∎

######  Remark 4.1.

  1. (1)

For some instances for which the assumptions $(C_{1})$ and $(C_{2})$ are
verified, we refer to [3, Example 3.3].

  2. (2)

If the initial conditions are not subject to perturbation, (i.e., if
$(u_{0,\lambda},v_{0,\lambda})=(u_{0},v_{0})$ for all
$\lambda\in\mathcal{V}(\overline{\lambda})$), then the estimation of Theorem
4.1 reduces to: for all $\lambda\in\mathcal{V}(\overline{\lambda})$ and for
a.e. $t\in[0,T]$,

$\|x_{\lambda}(t)-x(t)\|\leq L^{\prime}\left(\frac{2}{2+\gamma
T}\right)^{2}\left(e^{\frac{(2+\gamma
T)t}{2}}-1-t\right)\|\lambda-\overline{\lambda}\|.$

  3. (3)

In [3, Theorem 3.1], based on a Gronwall-Bellman argument, the authors
obtained an estimate similar to the one of our Theorem 4.1 with same
hypotheses by reformulating the second order dynamical system
$S(A_{\lambda},u_{0,\lambda},v_{0,\lambda})$ to a first order one. While in
our case here, thanks to the Perov’s inequality, the estimate (4.1) is rather
obtained by a direct computation without recourse to the first order case.

###  4.2. Differential equation of RLC circuit’s current

In this paragraph, we apply our main result on Lipschitz stability of the
previous section to a concrete problem mathematically modelled as an initial-
value problem for second order ODE, which interests in particular the audience
from physics and engineering. In this respect, we consider the following
equation modelling the electric current in an RLC parallel circuit, also known
as a tuning circuit (see [17]):

$\begin{cases}\frac{d^{2}x}{dt^{2}}+\frac{R}{L}\frac{dx}{dt}=g(t,x(t)),\,\,\,t\in[0,1]\\\
x(0)=\dot{x}(0)=0,\end{cases}$ (4.2)

where $g:[0,1]\times\mathbb{R}^{+}\longrightarrow\mathbb{R}$ is a continuous
function, $R$ is the resistance of the resistor and $L$ is the inductance of
the inductor.

Figure 1. RLC parallel circuit

Clearly, the problem (4.2) is equivalent to the following integral equation

$x(t)=\int_{0}^{1}G(t,s)g(s,x(s))ds,\,\,t\in[0,1],$ (4.3)

where $G$ is the Green’s function defined by

$G(t,s)=\begin{cases}\frac{1}{\tau}(1-e^{\tau(s-t)})&if\quad 0\leq s\leq t\leq
1\\\ 0\qquad\qquad&if\quad 0\leq t\leq s\leq 1,\end{cases}$ (4.4)

and $\tau=\frac{R}{L}$. Following the same approach used in the proof of [19,
Theorem 4.1], the problem (4.3) admits a unique solution under the following
hypotheses:

  * •

there exists a continuous function $w:[0,1]\longrightarrow(0,\infty)$ such
that

$|g(s,x(s))|\leq\frac{\tau}{2}w(s)|x(s)|,$

for all $s\in[0,1]$;

  * •

$\max_{s\in[0,1]}w(s)=e^{-\alpha}$, where $\alpha>e$.

In our concern, we consider that the function $g$ is also subject to
perturbation, the perturbed circuit is modelled by the following equation:

$\begin{cases}\frac{d^{2}x_{\lambda}}{dt^{2}}+\frac{R}{L}\frac{dx_{\lambda}}{dt}=g_{\lambda}(t,x_{\lambda}(t)),\,\,\mbox{for}\,\,a.e\,\,t\in[0,1]\\\
x_{\lambda}(0)=x_{0,\lambda}\\\
\dot{x}_{\lambda}(0)=\dot{x}_{0,\lambda}.\end{cases}$ (4.5)

The initial value of the parameter $\lambda$ is denoted by
$\overline{\lambda}$, that is,
$g_{\overline{\lambda}}=g,x_{\overline{\lambda}}=x,x_{\overline{\lambda}}(0)=x_{0}\,\mbox{and}\,\dot{x}_{\overline{\lambda}}(0)=\dot{x}_{0}.$

As a direct application of Theorem 3.1, we obtain the following parametric
stability result for the system (4.5).

######  Theorem 4.2.

Suppose that the following assumptions are satisfied:

  1. (i)

the function $g_{\lambda}$ is continuous;

  2. (ii)

for some $\beta>0$, the function $g_{\lambda}$ is $\beta$-Lipschitz in $x$
uniformly in $t$ and $\lambda$;

  3. (iii)

for some $L^{\prime}>0$, the function $g_{\lambda}$ is $L^{\prime}$-Lipschitz
in $\lambda$ uniformly in $t$ and $x$;

  4. (iv)

there exists a continuous function $w:[0,1]\longrightarrow(0,\infty)$ such
that

$|g_{\lambda}(s,x(s))|\leq\frac{\tau}{2}w(s)|x(s)|,$

for all $s\in[0,1]$ and all $\lambda\in\mathcal{V}(\overline{\lambda})$;

  5. (v)

$\max_{s\in[0,1]}w(s)=e^{-\alpha}$, where $\alpha>e$.

Then, for each $\lambda\in\mathcal{V}(\overline{\lambda})$, the problem (4.5)
admits at least a solution and the following Lipschitz estimate holds: for a.e
$t\in[0,1]$

$\begin{split}\|x_{\lambda}(t)-x(t)\|\leq&\left(1+L\left(\frac{2}{2+L}\right)^{2}\left(e^{\frac{(2+L)t}{2}}-1-t\right)\right)\|x_{0,\lambda}\|\\\
&+\frac{2}{2+L}\left(e^{\frac{(2+L)t}{2}}-1\right)\|\dot{x}_{0,\lambda}\|+L^{\prime}\left(\frac{2}{2+L}\right)^{2}\left(e^{\frac{(2+L)t}{2}}-1-t\right)\|\lambda-\overline{\lambda}\|,\end{split}$
(4.6)

where $L=\sup(\beta,\tau)$.

###### Proof.

Under assumptions $(iv)$ and $(v)$, the existence of a solution to the problem
(4.5) follows from [19, Theorem 4.1]. Now, to prove the Lipschitz estimate
(4.6), it suffices to apply Theorem 3.1. To this end, let us consider the
function
$f_{\lambda}:[0,1]\times\mathbb{R}^{+}\times\mathbb{R}^{+}\longrightarrow\mathbb{R}$
defined, for each $\lambda\in\mathcal{V}(\overline{\lambda})$, by

$f_{\lambda}(t,x,y)=g_{\lambda}(t,x)-\tau y.$

Using hypothesis $(ii)$, we see that the function $f_{\lambda}$ is
$L=\sup(\beta,\tau)$-Lipschitz in $(x,y)$ uniformly in $t$. On the other hand,
assumption $(iii)$ implies that $f_{\lambda}$ is $L^{\prime}$-Lipschitz in
$\lambda$ uniformly in $(x,y)$. Consequently, the desired estimate (4.6)
immediately follows from inequality (3.1) with $T=1$ and
$x_{0}=\dot{x}_{0}=0$. ∎

######  Remark 4.2.

Our approach is of course also applicable to an RLC series circuit of the
following type:

Figure 2. A series RLC circuit

Indeed, since the charge on the capacitor in such circuit can be modelled by
the following second order ODE:

$\frac{d^{2}q(t)}{dt^{2}}+\frac{R}{L}\frac{dq(t)}{dt}+\frac{1}{LC}q(t)=V(t),$

where, $L$ is the inductance, $R$ is the resistance, $C$ is the capacitance
and $V(t)$ is the voltage source.
We easily see that the above equation can be equivalently written as

$\ddot{q}(t)=f(t,q(t),\dot{q}(t)),$

where $f(t,x,y)=g(t,x)-\frac{R}{L}y$ and
$g(t,x)=\frac{1}{L}V(t)-\frac{1}{LC}x$.
Hence, with initial conditions $q(0)=\dot{q}(0)=0$, when the function $V$ is
subject to perturbation, our result in Theorem 4.2 can be applied whenever the
perturbed form of the function $g$ satisfies all the assumptions of that
Theorem.

###  4.3. Linear control systems

In this subsection, we give a further interesting application of our main
result in Theorem 3.1. In this respect, for given 2$\times$2 matrices
$A,B,C,D$, a nonnegative real number $\gamma$ and an initial condition
$(x_{0},\dot{x}_{0})$, we aim to tackle linear control systems of second order
ODE of the form:

$\;\begin{cases}\ddot{x}(t)=Ax(t)+\gamma\dot{x}(t)+Bu(t),\\\
z(t)=C(x(t)+\dot{x}(t))+Du(t),\\\
x(0)=x_{0},\,\,\,\dot{x}(0)=\dot{x}_{0},\,\,\,z(0)=C(x_{0}+\dot{x}_{0})+Du(0),\\\
\end{cases}$ ($LCS$)

where $x(.),z(.)$ and $u(.)$ are respectively the state variable, the
observability variable and the control of the system. In this example, the
function $f$ is nothing else but

$f(t,x(t),y(t))=Ax(t)+\gamma y(t)+Bu(t).$

We consider here a linear perturbation of the matrix $A$, i.e.,
$A=(a_{ij})_{i,j=1,2}$ and

$A_{\lambda}=\begin{pmatrix}a_{11}+\lambda-\overline{\lambda}&a_{12}\\\
a_{21}&a_{21}+\lambda-\overline{\lambda}\end{pmatrix},$

$\overline{\lambda}$ being the initial value of $\lambda$. Therefore, the
perturbed linear control system under consideration is stated as follows:

$\;\begin{cases}\ddot{x}_{\lambda}(t)=A_{\lambda}x(t)+\gamma\dot{x}_{\lambda}(t)+B_{\lambda}u(t),\\\
z_{\lambda}(t)=C(x_{\lambda}(t)+\dot{x}_{\lambda}(t))+D_{\lambda}u(t),\\\
x_{\lambda}(0)=x_{0,\lambda},\,\,\,\dot{x}_{\lambda}(0)=\dot{x}_{0,\lambda},\,\,\,z_{\lambda}(0)=C(x_{0,\lambda}+\dot{x}_{0,\lambda})+D_{\lambda}u(0),\\\
\end{cases}$ ($LCS_{\lambda}$)

and the associated function $f_{\lambda}$ is given by

$f_{\lambda}(t,x(t),y(t))=A_{\lambda}x(t)+\gamma y(t)+B_{\lambda}u(t).$

In the next, we will make use of the following assumptions:

  * $(H_{1})$

there exist $r>0$ and a neighborhood $\mathcal{V}(\overline{\lambda})$ of
$\overline{\lambda}$ such that for all
$\lambda\in\mathcal{V}(\overline{\lambda})$, $x,\,x_{\lambda}\in
B(x_{0},r)\cap B(x_{0,\lambda},r)$ and $x^{\prime},\,x_{\lambda}^{\prime}\in
B(\dot{x}_{0},r)\cap B(\dot{x}_{0,\lambda},r)$;

  * $(H_{2})$

it is further assumed that on the same neighborhood
$\mathcal{V}(\overline{\lambda})$ the norm $\|A_{\lambda}\|$ is independently
bounded in $\lambda$, that is there exists $\alpha>0$ such that
$\|A_{\lambda}\|\leq\alpha$;

  * $(H_{3})$

there exists a constant $\beta>0$ such that for all
$\lambda,\lambda^{\prime}\in\mathcal{V}(\overline{\lambda})$ we have,

$\|B_{\lambda}-B_{\lambda^{\prime}}\|\leq\beta\|\lambda-\lambda^{\prime}\|\quad\text{and}\quad\|D_{\lambda}-D_{\lambda^{\prime}}\|\leq\beta\|\lambda-\lambda^{\prime}\|;$

  * $(H_{4})$

the control $u$ is such that $\|u\|_{\infty}\leq 1$.

In the Banach space $C([0,T];H)$, we consider the norm $\|\,.\,\|_{\infty}$
given for a function $x$ by
$\|x\|_{\infty}=\displaystyle\sup_{t\in[0,T]}\|x(t)\|$ and state the following
stability result for ($LCS_{\lambda}$):

######  Theorem 4.3.

Under the assumptions above, the following estimation holds:

$\begin{split}\|z-z_{\lambda}\|_{\infty}\leq
c_{1}\|x_{0,\lambda}-x_{0}\|+c_{2}\|\dot{x}_{0,\lambda}-\dot{x}_{0}\|+c_{3}\|\lambda-\overline{\lambda}\|,\end{split}$
(4.7)

where,

$c_{1}=\|C\|\left(\frac{2+LT-2L}{2+LT}+\frac{2L(4+LT)}{(2+LT)^{2}}e^{\frac{(2+LT)T}{2}}\right),\;\;c_{2}=\frac{\|C\|(4+LT)}{2+LT}e^{\frac{(2+LT)T}{2}}$

and

$c_{3}=\beta+\frac{4\|C\|L^{\prime}(1-L^{\prime})}{(2+LT)^{2}}e^{\frac{(2+LT)T}{2}},\;\;\text{with}\;\;L=\sup(\alpha,\gamma)\;\;\text{and}\;\;L^{\prime}=r+\|x_{0}\|+\beta.$

###### Proof.

We have for all $\lambda\in\mathcal{V}(\overline{\lambda})$,

$\|z-z_{\lambda}\|_{\infty}\leq\|C\|(\|x-x_{\lambda}\|_{\infty}+\|\dot{x}-\dot{x}_{\lambda}\|_{\infty})+\|D-D_{\lambda}\|\|u\|_{\infty}.$

Then, by assumptions $(H_{3})$ and $(H_{4})$, we get

$\|z-z_{\lambda}\|_{\infty}\leq\|C\|(\|x-x_{\lambda}\|_{\infty}+\|\dot{x}-\dot{x}_{\lambda}\|_{\infty})+\beta\|\lambda-\overline{\lambda}\|.$
(4.8)

Clearly, assumption $(H_{1})$ implies that

$\|A_{\lambda}x-A_{\lambda^{\prime}}x\|\leq(r+\|x_{0}\|)\|\lambda-\lambda^{\prime}\|.$

Combining this inequality with hypotheses $(H_{3})$ and $(H_{4})$, we conclude
that $f_{\lambda}$ is $L^{\prime}=(r+\|x_{0}\|+\beta)$-Lipschitz in $\lambda$
uniformly in $t$ and $(x,y)$.
On the other hand, under assumption $(H_{2})$, the function $f_{\lambda}$ is
$L$-Lipschitz in $(x,y)$ uniformly in $t$ and $\lambda$ with
$L=\sup(\alpha,\gamma)$.
Consequently, assumptions of Theorem 3.1 are all satisfied. Applying this
theorem, immediately follows from the estimate (3.1):

$\begin{split}\|x-x_{\lambda}\|_{\infty}\leq&\left(1+L\left(\frac{2}{2+LT}\right)^{2}e^{\frac{(2+LT)T}{2}}\right)\|x_{0,\lambda}-x_{0}\|\\\
&+\frac{2}{2+LT}e^{\frac{(2+LT)T}{2}}\|\dot{x}_{0,\lambda}-\dot{x}_{0}\|+L^{\prime}\left(\frac{2}{2+LT}\right)^{2}e^{\frac{(2+LT)T}{2}}\|\lambda-\overline{\lambda}\|.\end{split}$

Using thus inequality (3.4), we also have

$\|\dot{x}_{\lambda}-\dot{x}\|_{\infty}\leq\|\dot{x}_{0,\lambda}-\dot{x}_{0}\|e^{\frac{(2+LT)T}{2}}+2\frac{L\|x_{0,\lambda}-x_{0}\|+L^{\prime}\|\lambda-\overline{\lambda}\|}{2+LT}(e^{\frac{(2+LT)T}{2}}-1).$

Combining the last two inequalities with (4.8) yields the desired estimate
(4.7). ∎

##  5\. Numerical experiment

In this section, we present numerical tests of our stability based on the
Scientific Laboratory, Scilab, programming language, so we take back to the
linear control system ($LCS$) considered in the previous section. The
objective here is to justify the validity of Theorem 4.3 when the parametric
functions satisfy the conditions of that theorem. In this regard, for the sake
of simplicity, we make the following standard choices:

$A=\begin{pmatrix}0&-3\\\ 1&-4\end{pmatrix}$ and $C=\begin{pmatrix}1&1\\\
1&1\end{pmatrix}$.

We take the matrix $B=D=0$, the control $u(t)=\left(\begin{array}[]{c}1\\\
1\\\ \end{array}\right)$ doesn’t dependent on time and start from the
following initial conditions $x_{0}=\left(\begin{array}[]{c}1\\\ 1\\\
\end{array}\right)$ and $\dot{x}_{0}=\left(\begin{array}[]{c}0\\\ 1\\\
\end{array}\right).$
For the associated parametric form ($LCS_{\lambda}$) of our example, we
consider here a linear perturbation of the matrices $A,B$ and $D$ as follows:

$A_{\lambda}=\begin{pmatrix}\lambda&-3\\\
1&-4+\lambda\end{pmatrix}\quad\text{and}\quad
B_{\lambda}=D_{\lambda}=\begin{pmatrix}\lambda&0\\\ 0&\lambda\end{pmatrix},$

$\overline{\lambda}=0$ being the initial value of $\lambda$, i.e.,
$x_{\overline{\lambda}}=x,\,A_{\overline{\lambda}}=A,\,B_{\overline{\lambda}}=B=0,\,C_{\overline{\lambda}}=C$
and $D_{\overline{\lambda}}=D=0.$

Figure 3 below describes the behavior of the norm of the difference between
the observability solutions: $z_{\lambda}$ of the parametric system
($LCS_{\lambda}$) and $z$ of the initial one ($LCS$), i.e., the solution for
$\overline{\lambda}=0.$

Figure 3. The uniform convergence of $z_{\lambda}$ to $z$.

The numerical results in the left hand side of Figure 3 illustrate the rate of
convergence of $\|z-z_{\lambda}\|_{\infty}$ for various choices of $\lambda$,
when the initial condition is not subject to perturbation, while the second
panel displays the convergence rate of $\|z-z_{\lambda}\|_{\infty}$ for the
perturbed system starting from perturbed initial condition
$(x_{0,\lambda},\dot{x}_{\lambda,0})$ with:
$x_{0,\lambda}=\left(\begin{array}[]{c}1+\lambda\\\ 1+\lambda\\\
\end{array}\right)$ and $\dot{x}_{\lambda,0}=\dot{x}_{0}$ is keeping constant,
by using the same values of the perturbation parameter $\lambda$. Clearly,
both plots indicate that the norm $\|z-z_{\lambda}\|_{\infty}$ decreases
linearly with respect to the parameter $\lambda$ which confirm the validity of
the inequality (4.7).

Acknowledgements. The authors would like to thank Prof. M. Ait Mansour for the
nice discussion on the topic and the valuable comments on the content of the
paper. The second author would like to thank his supervisor, Prof. J.
Lahrache, for the continuous support during the preparation of his thesis. The
authors wish equally to warmly thank the editor in chief of ASVAO Journal
Prof. Christiane Tammer, and as well by the meantime, the guest editors Profs:
Valentina Sessa, Elisabeth Köbis, Jaafar Lahrache and Mohamed Ait Mansour for
their very nice cooperation during the workshops ASEM 23-El Jadida and ASSVAS
23-Safi to which is dedicated this issue.

## References

  * [1] M. Ait Mansour, On the perturbation of the observability equation in linear control systems, J. Global Optim., 106 (2008), 169-174.
  * [2] M. Ait Mansour, Two remarks on the stability of generalized hemivariational inequalities, J. Ineq. Pure Appl. Math., 3 (2002), 1-15.
  * [3] M. Ait Mansour, H. Attouch and Z. Mazgouri, Second order dynamical systems governed by maximally monotone operators: Quantitative stability, Evol. Eq. Control Th., to appear, DOI: 10.3934/eect.2023056.
  * [4] M. Ait Mansour, M. A. Bahraoui and A. El Bekkali, Sharp estimates for approximate exact solutions to quasi-optimization problems, Optimization, 71 (2022), 1331-1355.
  * [5] M. Ait Mansour, R-.A. Elakri and M. Laghdir, Quantitative stability for equilibrium problems under relaxed upper-sign properties: Application to quasiconvex programming, J. Appl. Anal. Optim., 3 (2019), 307-331.
  * [6] M. Ait Mansour, A. El Ayoubi and J. Lahrache, Sensitivity analysis of an optimal control problem under Lipschitzian perturbations, J. Appl. Numer. Optim., to appear.
  * [7] M. Ait Mansour, J. Lahrache and A. El Ayoubi, The stability of the parametric Cauchy problem of initial-value ordinary differential equations revisited, J. Appl. Numer. Optim., 5 (2023), 111-124.
  * [8] M. Ait Mansour, Z. Mazgouri and H. Riahi, A dynamical approach for the quantitative stability of bilevel parametric equilibrium problems and applications, Optimization, 71 (2021), 1389-1408.
  * [9] M. Ait Mansour, Z. Mazgouri and H. Riahi, Stability of parametric time dependent and dynamical equilibrium problems: application to hemivariational inequalities, J. Nonlinear Convex Anal., 21 (2020), 605-628.
  * [10] M. Ait Mansour and H. Riahi, On the eigenvalues problem for hemivariational inequalities: Existence and Stability, J. Convex Anal., 9 (2002), 1-17.
  * [11] M. Ait Mansour and H. Riahi, Sensitivity analysis for abstract equilibrium problems, J. Math. Anal. Appl., 306 (2005), 684-691.
  * [12] F. Alvarez and H. Attouch, An inertial proximal method for maximal monotone operators via discretization of a nonlinear oscillator with damping, Set-Valued Anal., 9 (2001), 3-11.
  * [13] F. Alvarez and H. Attouch, The heavy ball with friction dynamical system for convex constrained minimization problems. In: Optimization (Namur, 1998). Vol. 481, Lecture notes in economics and mathematical systems. Berlin: Springer; (2000). 25–35.
  * [14] M. Braun, Differential Equations and Their Applications. Springer, (1978).
  * [15] A. F. Filippov, Differential equations with discontinuous right-hand side, Matematicheskii sbornik 93 (1960), 99-128.
  * [16] A. Haraux, Systèmes Dynamiques Dissipatifs et Applications, Recherches en Mathématiques Appliquées 17, Masson, Paris, (1991).
  * [17] K. Harwood, Modeling a RLC circuit’s current with differential equations (2011).
  * [18] V. Henner, A. Nepomnyashchy, T. Belozerova and M. Khenner, Ordinary Differential Equations: Analytical Methods and Applications. Springer Nature, (2023).
  * [19] I. Iqbal, M. Rizwan, Existence of the solution to second order differential equation through fixed point results for nonlinear $F$-contractions involving $w_{0}$-distance, Filomat, 34 (2020), 4079-4094.
  * [20] J.D. Lambert, Numerical Methods for Ordinary Differential Systems. Wiley, Chichester (1991)
  * [21] A. I. Lurie, Analytical mechanics. Springer Science and Business Media, (2002).
  * [22] D. S. Mitrinović, J. E. Pec̆arić and A. M. Fink, Inequalities involving functions and their integrals and derivatives, Vol. 53. Springer Science and Business Media, (1991).


# Deep Bayesian Experimental Design for Quantum Many-Body Systems

Leopoldo Sarra1,2, and Florian Marquardt1,2 0000-0001-7504-8656
0000-0003-4566-1753 1 Max Planck Institute for the Science of Light,
Staudtstraße 2, 91058 Erlangen, Germany  2 Department of Physics, Friedrich-
Alexander Universität Erlangen-Nürnberg, Staudtstraße 5, 91058 Erlangen,
Germany<EMAIL_ADDRESS>
###### Abstract

Bayesian experimental design is a technique that allows to efficiently select
measurements to characterize a physical system by maximizing the expected
information gain. Recent developments in deep neural networks and normalizing
flows allow for a more efficient approximation of the posterior and thus the
extension of this technique to complex high-dimensional situations. In this
paper, we show how this approach holds promise for adaptive measurement
strategies to characterize present-day quantum technology platforms. In
particular, we focus on arrays of coupled cavities and qubit arrays. Both
represent model systems of high relevance for modern applications, like
quantum simulations and computing, and both have been realized in platforms
where measurement and control can be exploited to characterize and counteract
unavoidable disorder. Thus, they represent ideal targets for applications of
Bayesian experimental design.

Keywords: active learning, bayesian optimal experimental design, quantum many-
body systems, quantum technologies

##  1 Introduction

We are currently witnessing rapid scaling in the number of components for
quantum technology platforms. Fulfilling the promise of fault-tolerant quantum
computation will eventually require millions of qubits, and while current
implementations still fall far short of that goal, the specific challenges of
scaling are already apparent for the present-day devices including on the
order of hundred qubits [1]. Other areas like integrated photonic circuits [2]
are following a similar trajectory, for applications such as neuromorphic
computing [3] or sensing and more fundamental studies in areas like
topological transport [4]. There, large networks of beamsplitters, waveguides
and resonators, again with component counts on the order of dozens or
hundreds, are being fabricated and deployed. Similar large-scale networks have
now been fabricated and investigated for coupled mechanical resonators,
producing phononic circuits with local access to vibrational modes.

Characterizing any of these devices is a very important but nontrivial task,
especially if it is to be done in an efficient manner, within a limited time
budget [5]. Active learning [6], in the form of optimal experimental design
can help, provided that one can employ techniques able to deal with the large
number of parameters that are to be determined. These parameters comprise
resonance frequencies of optical or microwave modes, couplings between
components (e.g. between resonators and qubits or waveguides, or beamsplitter
transparencies), nonlinearities, propagation phase shifts, matrix elements for
the effect of external drives, and decay rates. Many different measurement
approaches can be drawn upon, each of them coming with its respective
parameters that can be adjusted prior to each new measurement. In linear
devices, specific components of the scattering matrix can be measured by
injecting waves in some port and performing a homodyne measurement on another
port. Here, the frequency would be a continuous parameter, while, depending on
the setup, the choice of ports could be another, discrete parameter. In
nonlinear systems, such as circuits comprised of qubits, manipulation via
pulses in Ramsey-type schemes is the natural choice, with a final projective
measurement of one or several qubits. The drive amplitude and pulse times
would then represent the measurement parameters to be optimized. Systems that
couple qubits and resonators or waveguides can also be characterized via
nonlinear transmission, with the amplitude and drive frequency being treated
as adjustable.

In all these cases, there is unavoidable fundamental noise in the measurement
outcomes, namely shot noise in the case of wave transmission measurements and
quantum projection noise in the case of qubit measurements [7]. This noise can
naturally be reduced by extending the measurement duration, increasing the
wave amplitude, or repeating multiple times the qubit pulse sequence together
with its final projective qubit measurement (multi-shot measurement). As a
consequence, the information gain per such extended measurement increases.
However, there are limits to the usefulness of this naive optimization: wave
amplitudes can be increased only so far before entering nonlinear regimes or
heating up the device. In addition, one cannot keep measuring at only one
single choice of the measurement parameter, since that will eventually only
pinpoint a particular function of the many underlying setup parameters and not
allow to resolve all the parameters individually. That is where active
learning, i.e. choosing next measurement settings efficiently, can offer true
benefits.

In this paper, we investigate recent and promising techniques from machine
learning [8, 9] to efficiently and accurately make the best parameter
prediction given past measurement, and propose the next best measurement to
perform. We show applications of these techniques to the above-mentioned
quantum devices. We focus especially on quantum systems, such as chains of
coupled cavities and arrays of qubits. The presented framework is nonetheless
completely general and appliable in a similar way to other settings.

[width=0.9]figures/algorithm_sketch

Figure 1: Sketch of the model of a system’s response and summary of the
parameter estimation procedure. At the top: the observed system produces a
noisy observation given the measurement settings and the system’s parameter
$\lambda$. At the bottom: starting from the prior we have about the system, we
can perform multiple steps. In each step, given the prior and the system’s
likelihood, we approximate the posterior distribution by minimizing the
Barber-Agakov bound (3). At the same time, the minimization also gives us the
optimal measurement setting to perform the measurement at. After the
measurement has been performed, the approximated posterior, conditioned on the
measured value and outcome, becomes the new prior, and a new iteration can be
performed.

##  2 Related Work

The idea of finding the best experiment to perform is known as Active Learning
in the machine learning literature [10, 6, 11] and as Bayesian optimal
experimental design [12] in more specific parameter estimation applications
[13, 14]. It is very common in science to have a class of possible models,
dependent on a set of parameters, and to perform experiments to find those
that better describe the true system. The outcome of each experiment does not
give direct knowledge on the parameters themselves because of noise and
measurement errors, but it provides partial information. Therefore, it is
generally not enough to perform as many measurements as the number of unknown
parameters of the system, and more experiments are required. When enough
information is collected through measurements, we can identify the true
parameters with a certain confidence.

When experiments are expensive, either in terms of cost, time and effort, it
can be important to be efficient in the number of experiments required to
characterize the given physical system. In those cases, it is desirable to
exploit the information obtained from each experiment as much as possible, and
to choose the sequence of experiments in such a way that the smallest number
of experiments is required.

Information Theory provides satisfactory answers to this problem, at least
from the theoretical point of view. In particular, the Bayes theorem shows the
proper way to include a new measurement into our knowledge and update our
predictions of the parameters of the system, and expected information gain
[15] is the quantity to look at to find the best possible experiments to
perform.

There have already been many remarkable applications in science [16, 17] and
in physics in particular, ranging from devising quantum metrological
procedures [18] to characterizing quantum dots [19, 20], superconducting
quantum processors [21], or sensors [22], multi-phase estimation [23, 24], and
their use is becoming more and more common. Parameter estimation can then be
further generalized to the complete design of the experiment. For example,
[25, 26, 27] completely automate the search for the experiment that solves a
given task, in the case of quantum optics.

Nevertheless, optimal experimental design techniques have not been used in
their full power up to recently, even though these techniques have been known
for decades. The main reason is that the exact evaluation of those statistical
quantities can be very expensive. Therefore, very rough approximations were
usually made, including the use of empirical heuristics, the use of Gaussian
processes [28] and Gaussian posterior approximation [29, 30] and the
optimization through maximum likelihood estimation [31]. The approximation
usually depends on where the bottleneck is in practice, according to the
specific problem: high-dimensional parameter space, large number of different
possible experiments to perform, large number of measured quantities of each
experiment, experiment execution time and cost, and total allowed number of
experiments.

Modern neural networks [32] that employ normalizing flows [33, 34] for the
approximation of the posterior distribution, developed in the past few years,
recently allowed for quite more efficient estimations [8, 9] that do not
require such rough (and often unjustified) assumptions. The price to pay for
the increased precision of the approximation is clearly a larger computational
effort.

##  3 Methods

Bayesian experimental design can be implemented as follows. The settings of
the experiment, which could be, for example, a measurement at a particular
frequency, are described by a vector $x$. The physical system is identified by
hidden parameters $\lambda$, which we would like to determine. Because of
measurement noise, the outcome of the experiment, $y$, is not deterministic,
but distributed according to a distribution $P(y|\lambda,x)$, called
likelihood of an observation. We can imagine, for example,
$y=f(x,\lambda)+\epsilon$, where $f$ is a deterministic function and
$\epsilon$ is a random variable, e.g. Gaussian distributed. We can update our
knowledge of the parameters of the system with the Bayes rule [35]: given our
inferred distribution $P_{n}(\lambda|\mathcal{M}_{n})$ after $n$ measurements
$\mathcal{M}_{n}=\\{(x_{1},y_{1}),\ldots,(x_{n},y_{n})\\}$ (to which we will
refer to as prior at step $n+1$), and the subsequent measurement
$(x_{n+1},y_{n+1})$, the updated distribution (posterior at step $n+1$) is

$\displaystyle
P_{n+1}(\lambda|\mathcal{M}_{n+1})=\frac{P(y_{n+1}|\lambda,x_{n+1})P_{n}(\lambda|\mathcal{M}_{n})}{P_{n}(y_{n+1}|x_{n+1})}.$
(1)

The initial prior distribution $P_{0}(\lambda)$ is chosen arbitrarily, and it
reflects our initial assumptions on the parameters of the system.

To choose the next measurement to perform, it is possible to define a query
function, which assigns to each possible experiment $x$ its expected value:
the $x$ for which the query function is maximized is the one that is expected
to be the most useful to measure [6]. A possible choice of the query function
is the expected information gain when measuring at $x$:

$\displaystyle I_{n}[\lambda,y](x)=$ $\displaystyle\;\;=\int dyd\lambda
P_{n}(\lambda)P(y|\lambda,x)\log\frac{P_{n+1}(\lambda|y,x)}{P_{n}(\lambda)}.$
(2)

This can be interpreted as the mutual information between $y$ and $\lambda$
given a measurement at $x$, or in other words, as the entropy reduction after
measuring $(x,y)$.

However, Eq. (3) is hard to estimate. First, it requires to be able to sample
from the prior distribution $P_{0}(\lambda)$, and second, it requires the
value itself of the posterior probability distribution
$P_{1}(\lambda|y_{0},x_{0})$, since it appears inside the logarithm. In
principle, we could obtain the posterior distribution with the Bayes rule, Eq.
(1), but it is not efficient to calculate its normalization factor
$P_{0}(y|x)=\int d\lambda P_{0}(\lambda)\log P_{0}(\lambda)$. Overall, to get
Eq. (3) we would require two nested Monte Carlo estimates, one for the
evidence and one for the mutual information estimation. Then, we should
optimize over $x$ to get the best measurement to perform.

Modern neural network techniques allow implementing variational bounds and
perform a much more efficient estimate [36]. For example, we can avoid
evaluating the posterior [37] explicitly by introducing a new function
$Q(\lambda|y,x)$ and calculating the quantity

$I(x)=\int dyd\lambda
P_{0}(y|x)P_{1}(\lambda|y,x)\log\frac{Q(\lambda|y,x)}{P_{0}(\lambda)}.$ (3)

If $Q$ is a probability distribution, we know that
$KL(P_{1}(\lambda|y,x)||Q(\lambda|y,x))\geq 0$ thus

$\displaystyle I(x)\leq I_{0}[\lambda,y](x).$ (4)

This is the so-called Barber-Agakov bound of mutual information [38]. We can
parametrize $Q(\lambda|y,x)$ with a neural network, so that we do not need to
find a different $Q$ function for every different $(x,y)$, but we can
interpolate, amortizing the costs of the evaluation.

In particular, we use a conditional normalizing flow [39]: we start from a
normal Gaussian distribution $\mathcal{G}(z)$ and perform a series of
invertible transformations (parametrized by neural networks) having a Jacobian
that is simple to estimate. Optimizing over $Q(\lambda|y,x)$ means to optimize
over the parameters of the neural networks that implement the transformations.
Normalizing flows can be easily implemented in common deep learning frameworks
like TensorFlow [40]. In particular, we use TensorFlow Probability [41]. The
procedure is summarized in figure 1. Please refer to A for more details.

There are multiple ways to estimate the performance of the proposed result
(and compare with alternative techniques). The sum of the information gained
at step $n+1$ after measuring $(x_{n+1},y_{n+1})$,

$\int d\lambda
P_{n}(\lambda|y_{n},x_{n})\log\frac{Q(\lambda|y_{n},x_{n})}{P_{n}(\lambda)},$
(5)

can be a measure for how much we understood about the system. However, in
contrast to (3), which represents the expected information gain averaged over
all the possible measurement outcomes, this quantity can also be negative.
Indeed, it is possible that an unexpected outcome of a measurement increases
the global uncertainty on the parameters. Properties of the final posterior
distribution such as its variance, or the comparison of its mean (or most
likely value) with the true parameters are also a possible metric to consider.
On the other hand, the best prediction for an observation is actually the
average over all the allowed parameters

$P(y|x)=\int d\lambda P_{n}(\lambda|\mathcal{M}_{n})P(y|\lambda,x).$ (6)

It is possible to compare this prediction with the likelihood calculated at
the true parameters $\lambda^{*}$, for example using the Kullback-Leibler
divergence. This gives an idea of how much the predictions of our model can
differ from the observations. It is important to emphasize that if the
physical model is correct, we would expect only a single parameter set to be
relevant, corresponding to the true parameters of the system. However, if the
model does not represent well the true system, the use of multiple values may
be helpful to produce a better approximation.

##  4 Applications and discussion

We have already mentioned in the introduction the wide variety of quantum
technology platforms, system parameters, and potential measurement approaches
that could benefit in principle from active learning approaches. In the
following, we have chosen two illustrative examples that are sufficiently
distinct in their characteristics. As a first example, we consider a linear
device, where wave transmission is measured to extract the resonance
frequencies of a coupled-cavity array. Such cavity arrays [42, 43, 4, 44] have
been considered in studies of transport as coupled-resonator optical
waveguides in optical setups or quantum many-body physics of photons, when
combined with nonlinear elements like qubits. In our second example, we turn
to an array of qubits, which is a widespread platform used for quantum
simulations and quantum computing [43, 45, 46, 47, 44, 48]. There, we consider
a typical pulsed scheme followed by a projective measurement to extract qubit
parameters. B provides additional technical details about the implementations
of the examples.

###  4.1 Coupled cavities

[width=]figures/measurements

Figure 2: Deep Bayesian optimal experimental design applied to estimating the
frequencies of an array of $6$ coupled cavities. (a) Sketch of the system. (b)
Measurements at each step. The red symbols show the chosen frequency $x$ to
measure at each step. In the background, the value of the expected information
gain for each possible $x$ (brighter is higher), normalized at each step, i.e.
$IG(x)/\max_{x}IG(x)$. (c) Expected information gain values for each possible
measurement $x$. A peak represents an optimal value to measure. Different
lines represent different measurement steps. (d) Inferred parameter
distribution after $15$ measurements. The diagonal shows the marginal
distribution $P(\lambda_{i})$, the off-diagonal plots the two-variable slice
$P(\lambda_{i},\lambda_{j})$. (e) Evolution of the marginalized posterior
distribution of the last cavity frequency, until converging to a sharp peak
around the true value.

We consider a cavity array with the following Hamiltonian:

${\hat{H}}=\sum_{j}\Omega_{j}{\hat{a}}_{j}^{\dagger}{\hat{a}}_{j}+\sum_{j}J_{j+1}({\hat{a}}_{j}^{\dagger}{\hat{a}}_{j+1}+{\rm
h.c.}),$ (7)

where the sum runs over all but the last of the cavity modes in this chain
with open boundary conditions. Here and in what follows, we consider
$\hbar=1$.

Since the setup we are dealing with is linear, it is sufficient to solve the
classical equations of motion for the coupled modes, driven by a wave entering
from a waveguide coupled to the first cavity. To this end, we consider the
classical coherent-state amplitudes $a_{j}$ corresponding to the quantum
operators ${\hat{a}}_{j}$, including drive and decay as prescribed by input-
output theory [measurement]. For brevity of our notation, we collect all these
amplitudes in a vector and all frequencies in a matrix, where first-neighbour
interactions are $J$ and proper frequencies $\Omega_{i}$. Let our system also
have some internal decay $\kappa_{\text{int}}$ and external decay
$\kappa_{\text{ext}}$, which only applies to the cavities coupled to the
environment, for example the first and last.

Given the entering fluctuating field $a_{\text{in}}$, which can also contain a
laser drive, the output field $a_{\text{out}}$ is given by [49]

$a_{\text{out}}=\mathcal{S}a_{\text{in}}=\left(\mathbf{1}-\frac{\sqrt{\kappa_{\text{ext}}}}{-i(\omega-\Omega)+\frac{\kappa}{2}}\right)a_{\text{in}}.$
(8)

In our example, we imagine driving the first cavity and looking at the
response at the last cavity. Our goal is to extract the scattering matrix
element for transmission from the first to the last cavity, given by the
$\mathcal{S}_{0N}$ in (8). Being it a complex number, its knowledge
corresponds to measuring both the amplitude and phase of the emitted light. By
performing the smallest possible number of measurements at a frequency
$\omega$ close to the resonance frequency, we want to discover the frequencies
of all the cavities. For example, they may all be close to a common resonance
frequency, but slightly detuned. Each measurement is affected by noise
$\epsilon$, which can be assumed to be Gaussian [7]. This noise, being
observed in a measurement of the field amplitude, ultimately arises from
quantum vacuum noise being injected into the input ports of the device (and
through the dissipation channels) and propagating linearly through the setup.

In figure 2, we show our results. We consider $6$ coupled cavities, with each
random coupling uniformly chosen in $J_{j}/\kappa\in[1,3]$. We assume a
setting where the couplings $J$ and the dissipation $\kappa$ have already been
calibrated and are therefore known, and our goal is to find the detunings. In
the general notation of section 3, the measurement setting $x$ corresponds in
this case to the choice of the drive frequency $\omega$, the response function
$f(x,\lambda)$ corresponds to the scattering matrix element
$\mathcal{S}_{0N}$, and the unknown parameters of the system $\lambda$ are the
frequencies of the cavities $\Omega_{i}$.

As more measurements are performed, we improve our estimate of the cavity
frequencies. We see that the information provided by new measurements
decreases as we make more observations. At the same time, the posterior
distribution converges to a sharp peak around the true values. For example,
figure 2e shows the evolution of the last cavity frequency distribution. It is
important to emphasize that even if the posterior eventually converges to a
Gaussian, it is not Gaussian at the beginning, as shown in figure 2d. This is
why it is necessary to approximate the posterior using a neural network
approximation.

By looking at the measurement frequencies that the active approach selects, we
can identify a pattern, or measurement strategy: it measures the frequencies
where the slope of the response is larger, alternating between different
points. These points give indeed the largest information on the intrinsic
frequencies of the cavities.

[width=0.95]figures/strategies

Figure 3: Comparison of different measurement strategies for the coupled
cavities system. (a) Final response functions after inferring the parameters
with the active, random, or fixed strategy. The orange curves represent the
response function induced by various samples of parameters $\lambda$ from the
final posterior distribution. In blue, the true response. Red symbols mark the
performed measurements and their outcome. (b) Evolution of the KL-divergence
between the inferred response function distribution $P(y|x)$ and the system
likelihood during the various measurements. Brighter means higher. (c)
Cumulative sum of information gained after each measurement, i.e.
$\sum_{n=1}^{N}\log{P_{n}(\lambda|y,x)}-\log{P_{n-1}(\lambda)}$. Higher values
represent learning more about the system. Results are averaged over $3$ runs.

We compare the results of the inference with an active choice of the
measurement with two other strategies: fixed and random. The fixed strategy is
the simplest: all the measurements are performed at the same $x$ value, and it
is clearly often not possible to learn all the parameters from that. There
will usually be a region of the response function, close to the measurement
region, that is very well learned, while other regions will not be very
accurate. A better naive measurement selection approach is the random one. In
this case, measurements are chosen uniformly within the specified measurement
range (here, $x=\omega/\kappa\in[-12,12]$). We can see that in figure 3: the
active strategy, i.e. the one that chooses the next measurement greedily by
maximizing the expected information gain, learns the parameters of the system
faster than the random one. Asymptotically, both the random and active
strategy will converge to the right values. However, active selection allows
to save a lot of measurements.

###  4.2 Array of qubits

We consider an array of $N$ qubits, whose many-body Hamiltonian reads like

$H=\sum_{i=0}^{N-1}\omega_{i}\sigma_{z}^{(i)}+J\sum_{i=0}^{N-2}\sigma_{x}^{(i)}\sigma_{x}^{(i+1)}.$
(9)

At the beginning, the system starts in its ground state $|\psi_{0}\rangle$. We
assume that qubits can be individually addressed, and can be subjected to
microwave pulses and projective measurements. In particular, we excite the
first qubit of the array by applying a microwave pulse with pulse area $\pi$.
This corresponds to implementing the unitary

$U=e^{-i\frac{\pi}{2}\sigma_{x}^{(0)}}.$ (10)

We let the system evolve freely for some time $t$, and then we perform a
projective measurement on the first qubit. The evolution allows for the
propagation of the external excitation. The outcome of the measurement is
binary, according to probability

$p_{\uparrow}^{(0)}=\langle\psi(t)|\hat{\mathcal{P}}_{\uparrow}^{(0)}|\psi(t)\rangle,$
(11)

where
$\hat{\mathcal{P}}_{\uparrow}^{(0)}=|\uparrow\rangle_{0}\langle\uparrow|$ is
the projection operator on the ”up”-state of the first qubit (qubit number
$0$). In each experiment, we can choose the duration of the free evolution
after the initial pulse, which we will denote as the ”measurement pulse
duration” $x\in[0,5]$. By performing the smallest possible number of
experiments, we want to discover the frequencies $\omega_{i}$ of all the
qubits.

We considered a $4$-qubit system and performed $100$ measurements in series,
choosing the pulse duration for each of them. In this case, each of those
measurements is in itself a multishot measurement, being repeated $100$ times
to give better statistics (otherwise the outcome would only be either $0$ or
$1$). Again, as shown in figure 4, we compare the active strategy with a
random and a fixed one, and show the efficiency advantages of the active one
in terms of gained information. In this case, we also show the performance of
a uniform strategy, which takes sequential equally-spaced measurements, i.e.
sweeping the pulse duration from $0$ to the maximum allowed value. It is
interesting to observe that this strategy eventually learns the parameters of
the system, but it requires many more measurements than the random one. Also,
in figure 4d we see that the active strategy usually requires, on average,
sensibly fewer measurements than the random strategy to provide the same total
information gained. For example, to reach a total information gained of $10$,
we require about $10$ measurements with the active approach or about $20$
measurements with the random one.

The fixed strategy always measures at $x=0$, i.e. does not allow any
relaxation after the microwave pulse on the first qubit. As a consequence,
this measurement effectively only provides information about the interacting
ground state of the system, which already contains some limited information
about the qubit frequencies. We emphasize, indeed, that since the interacting
quantum many-body ground state is not an eigenstate of $\sigma_{z}$ acting on
the first qubit, the probability of measuring $\uparrow$ at $t=0$ is not $1$.

Furthermore, it is possible to spot a sort of strategy in the pattern of the
measurements performed by the active selection: a few measurement regions are
alternated cyclically. Indeed, we can imagine that some regions give the
largest information about some parameters. After we measure in one region, we
increase the accuracy on some parameters, and another region becomes more
relevant, until it is useful to measure again at the initial spot. This cycle
is repeated until we reach the desired accuracy on the parameters.

[width=]figures/strategies-qubits

Figure 4: Deep Bayesian experimental design applied to the estimation of the
frequencies of an array of $4$ qubits. (a) Sketch of the system. (b) Chosen
pulse duration $x$ at each step. The red symbols show the chosen $x$ to
measure at each step. In the background, the value of the expected information
gain for each possible $x$ (brighter is higher), normalized at each step
maximum. (c) Response functions after $10$ measurements, comparing the active,
random, fixed and uniformly spaced approaches. The orange curves represent the
response function induced by various samples of parameters $\lambda$ from the
final posterior distribution. In blue, the true response. Red symbols mark the
performed measurements and their outcome. (d) Sum of information gained after
each measurement, i.e.
$\sum_{n=1}^{N}\log{P_{n}(\lambda|y,x)}-\log{P_{n-1}(\lambda)}$. Higher values
imply learning more about the system. Results are averaged over $3$ runs.

##  5 Outlook

In this paper, we have introduced deep optimal Bayesian experimental design
for modern quantum technologies. This approach approximates the posterior
update with a variational bound and a deep neural network and allows
extracting the optimal measurement to perform at each step. We have shown the
application of this technique to two promising quantum platforms, cavity
arrays and qubit chains. In both cases, the active measurement selection
technique allows learning the parameters of the system with fewer measurements
than other strategies.

Depending on the application, different assumptions should be made, leading to
different tradeoffs between the approximation precision and the time
efficiency. The approach we propose is especially advantageous in the setting
where each measurement is very expensive and it is worth spending a long
calculation time to really find the optimal next measurement setting.

The main challenge at present is the time it takes to update the neural
network representing the Bayes posterior distribution as well as optimizing
over possible measurement settings. This time is still too large in order for
this technique to be deployed economically “as is” in the given scenarios,
where individual measurements can happen on microsecond time scales, and even
extended sequences used to reduce shot noise will probably not last longer
than a millisecond for one measurement setting, except when extreme accuracy
is called for. In the present (not yet very optimized) technique, an
optimization run still takes on the order of hundreds of seconds per
measurement step, which is currently a drawback of Bayesian Optimal
Experimental Design based on deep neural networks in general. However, a
remedy would consist in performing the optimization on a number of different
example scenarios, in simulations, and then train a neural network in a
supervised fashion to learn the suitable choice of the next measurement
setting based on all the previous settings and outcomes. In this way, we could
really develop a sort of policy that can be deployed to characterize unknown
systems of the same class.

A specific cost function may also be added to take into account that some
measurement values may require more resources and thus be more expensive than
others [6].

A further generalization can be to drop the likelihood assumption and try to
develop a likelihood-free active inference technique. This would be similar to
[24]. However, we do not want to discretize the input space $x$, since it
would not scale up to large dimensions. Finally, we notice that we are
employing a greedy strategy which always chooses the single next best
measurement to perform. Alternative approaches may include providing an
overall measurement budget (e.g. total number of measurements) or set a target
total information gain to achieve. It would be interesting to study how these
sequential strategies [50], perhaps approximated with a reinforcement learning
agent [51, 22], may suggest better measurement strategies than the greedy one
we applied in this paper.

From the broader perspective, we can imagine in the future to generalize the
measurement setting $x$ to become an entire experimental setup, and to perform
the search on a broader class of experiments. This would be a very useful
aspect of a future “artificial scientist” that tries to explore the world and
learn a model by performing experiments autonomously.

##  Appendix A Bayesian Experimental Design

In this Appendix, we describe in more detail the standard framework of
Bayesian experimental design.

###  A.1 The Bayesian parameter update

Starting from some a priori distribution $P_{0}(\lambda)$, that reflects all
our expectations on the parameters of the system (e.g. parameter range,
symmetries, constraints), we can update it according to the Bayes rule [35]

$P_{1}(\lambda|x_{0},y_{0})=\frac{P(y_{0}|\lambda,x_{0})P_{0}(\lambda)}{P_{0}(y_{0}|x_{0})},$
(12)

where $(x_{0},y_{0})$ is the first measurement and its outcome and

$P_{0}(y_{0}|x_{0})=\int d\lambda P(y|\lambda,x)P_{0}(\lambda).$ (13)

The distribution $P_{1}(\lambda)$, which depends on the performed experiment
and on its outcome $(x_{0},y_{0})$, is called posterior distribution. It
represents the updated knowledge on $\lambda$ after taking into account the
result of the experiment. It is easy to show that in case of multiple
(independent) experiments, $\left\\{x_{i},y_{i}\right\\}^{n}$, this equation
can be easily generalized to a recursive form:

$\displaystyle
P_{n}\left(\lambda|\\{x_{i},y_{i}\\}^{n}\right)=\frac{\prod^{n}P(y_{i}|\lambda,x_{i})P_{0}(\lambda)}{\int
d\lambda^{\prime}\prod^{n}P(y_{i}|\lambda^{\prime},x_{i})P_{0}(\lambda^{\prime})}=$
(14)
$\displaystyle\;\;=\frac{P(y_{n}|\lambda,x_{n})P_{n-1}\left(\lambda|\\{x_{i},y_{i}\\}^{n-1}\right)}{P_{n-1}(y_{n}|x_{n})}$
(15)

with

$P_{n}(y_{n}|x_{n})=\int d\lambda
P(y_{n}|\lambda,x_{n})P_{n-1}\left(\lambda|\\{x_{i},y_{i}\\}^{n-1}\right).$
(16)

As one would expect, we get exactly the same updated prior if we perform many
experiments in parallel and then update our knowledge or if we perform them
one after the other and update our prior after each single experiment.

###  A.2 The information gain query function

We define a query function to choose which experiment to perform. This
function assigns to each possible experiment $x$ its expected usefulness: the
$x$ for which the query function is maximized is the one expected to be most
useful to measure [6]. The query function we use in this paper is the expected
Kullback-Leibler divergence between the updated parameter distribution and the
prior. We start from

$KL(P_{1}(\lambda|y,x)||P_{0}(\lambda))=\int d\lambda
P_{1}(\lambda|y,x)\log\frac{P_{1}(\lambda|y,x)}{P_{0}(\lambda)},$ (17)

which estimates how much the prior differs from the updated distribution as a
function of the experiment $x$ and its outcome $y$. Then, we calculate its
average over the possible $y$, which gives the expected information gain when
measuring at $x$ i

$\displaystyle IG(x)=\int
dyP_{0}(y|x)KL(P_{1}(\lambda|y,x)||P_{0}(\lambda))=\noindent$ (18)
$\displaystyle\;\;=\int dyd\lambda
P_{0}(y|x)P_{1}(\lambda|y,x)\log\frac{P_{1}(\lambda|y,x)}{P_{0}(\lambda)}.$
(19)

This equation can be rewritten as

$\displaystyle\int dyd\lambda
P_{0}(\lambda)P(y|\lambda,x)\log\frac{P_{1}(\lambda|y,x)}{P_{0}(\lambda)}=\noindent$
(20) $\displaystyle\;\;=H_{0}(\lambda)-H_{1}(\lambda|y)(x)=\noindent$ (21)
$\displaystyle\;\;=I_{0}(\lambda,y)(x),$ (22)

with $H_{0}(\lambda)=-\int d\lambda P_{0}(\lambda)\log P_{0}(\lambda)$ and
$H_{1}(\lambda|y)(x)=-\int P_{0}(\lambda)P(y|\lambda,x)\log
P_{1}(\lambda|y,x)$. This can be interpreted as the mutual information between
$y$ and $\lambda$, or, in other words, as the entropy reduction after
measuring $(x,y)$. The optimal $x$ to measure next, if we consider a greedy
strategy, is the one that maximizes this quantity.

##  Appendix B Applications

###  B.1 Implementation

In our experiments, we employ a $4$-layer normalizing flow, each employing a
$2$-layer MADE network (Masked Autoregressive flow for Density Estimation
[52]) with $64$ neurons each. After each layer of the normalizing flow, the
input vector is shifted by one with a permutation layer, so that the
subsequent layer applies a transformation to a different subset of the input.
We optimize the network with an Adam optimizer with learning rate
$\eta=10^{-3}$.

###  B.2 Coupled-cavity arrays

We consider a cavity array with the following Hamiltonian:

${\hat{H}}=\sum_{j}\Omega_{j}{\hat{a}}_{j}^{\dagger}{\hat{a}}_{j}+\sum_{j}J_{j+1}({\hat{a}}_{j}^{\dagger}{\hat{a}}_{j+1}+{\rm
h.c.}),$ (23)

where the sum runs over all but the last of the cavity modes in this chain
with open boundary conditions. We consider the classical coherent-state
amplitudes $a_{j}$ corresponding to the quantum operators ${\hat{a}}_{j}$,
including drive and decay as prescribed by input-output theory [53]. Given
first-neighbour interactions $J_{i}$ and proper frequencies $\Omega_{i}$, we
can define the matrix

$\Omega=\begin{pmatrix}\Omega_{0}&J_{1}&0&\ldots&\\\
J_{1}&\Omega_{1}&J_{2}&0&\ldots\\\ 0&J_{2}&\Omega_{2}&J_{3}&\ldots\\\
\vdots&&&\ddots&\\\ &\ldots&0&J_{N}&\Omega_{N}\end{pmatrix}.$ (24)

We also include some internal decay $\kappa_{\text{int}}$ and external decay
$\kappa_{\text{ext}}$, which only applies to the cavities coupled to the
environment, for example the first and last. The overall decay can be
represented by a vector

$\kappa=\begin{pmatrix}\kappa_{\text{int}}+\kappa_{\text{ext}}\\\
\kappa_{\text{int}}\\\ \vdots\\\ \kappa_{\text{int}}\\\
\kappa_{\text{int}}+\kappa_{\text{ext}}\end{pmatrix}.$ (25)

Let $a_{\text{in}}$ be the entering fluctuating field, which can also contain
a laser drive, and $a_{\text{out}}$ the output field. The behavior of the
system is described by the input-output relations [49]

$\begin{cases}\dot{a}&=\left(-i\Omega-\frac{\kappa}{2}\right)a+\sqrt{\kappa_{\text{ext}}}a_{\text{in}}\\\
a_{\text{out}}&=a_{\text{in}}-\sqrt{\kappa_{\text{ext}}}a\end{cases}$ (26)

From the first equation, we can write the Green function (i.e. consider a
drive $a_{\text{in}}\propto e^{-i\omega t}$ and assume also $a\propto
e^{-i\omega t}$)

$a=\frac{\sqrt{\kappa_{\text{ext}}}}{-i(\omega-\Omega)+\frac{\kappa}{2}}a_{\text{in}}.$
(27)

Therefore, the final input-output relation becomes

$a_{\text{out}}=\left(\mathbf{1}-\frac{\sqrt{\kappa_{\text{ext}}}}{-i(\omega-\Omega)+\frac{\kappa}{2}}\right)a_{\text{in}}.$
(28)

For the case of spatially constant cavity mode frequencies, the $N$
eigenfrequencies of the open-boundary array span a bandwidth $4J$. All the
frequencies in the bulk of the spectrum are doubly degenerate, and the spacing
becomes smaller near the boundaries of the spectrum. This picture changes, and
the degeneracies are broken, when we introduce variation in the onsite
frequencies $\Omega_{j}$, which will be the generic situation we want to
explore.

The system we discuss in the main text has couplings $J$ uniformly sampled in
$[1,3]$ for each cavity, and proper frequencies sampled from a Gaussian
distribution with zero mean and unit variance. The values for the specific
example in figure 2 and figure 3 are shown in table 1. We optimize the
posterior neural network for $8000$ steps after each measurement, using a
batch size of $1500$ samples in (3).

Parameters | Value
---|---
(freq.) $\Omega_{0}$ | $1.040$
$\Omega_{1}$ | $0.326$
$\Omega_{2}$ | $0.520$
$\Omega_{3}$ | $0.900$
$\Omega_{4}$ | $-0.466$
$\Omega_{5}$ | $0.004$
Known par. | Value
---|---
(coupling) $J_{1}$ | $2.733$
$J_{2}$ | $2.615$
$J_{3}$ | $1.956$
$J_{4}$ | $1.568$
$J_{5}$ | $2.620$
(int. diss.) $\kappa_{int}$ | $0.5$
(ext. diss.) $\kappa_{ext}$ | $0.5$
(msmt noise) $\epsilon$ | $0.05$

Table 1:  Parameters for the coupled cavity example discussed in the main
text. On the left, the true frequencies we want to discover; on the right, the
known parameters of the system.

###  B.3 Qubit chain

Table 2 shows the parameters employed for the qubit chain example. We optimize
the posterior neural network for at least $2500$ steps after each measurement,
using a batch size of $700$ samples in (3). When also optimizing over $x$
(i.e. optimal measurement setting), we keep optimizing the loss until the
change in $500$ steps is smaller than $0.05$. This choice allows decreasing
the training time for this more computationally-demanding system, since we
train for a longer time only when required.

Parameters | Value
---|---
(freq.) $\omega_{0}$ | $1.163$
$\omega_{1}$ | $1.003$
$\omega_{2}$ | $1.045$
$\omega_{3}$ | $0.910$
Known par. | Value
---|---
(coupling) $J$ | $1.7$
(multishot) $n_{\text{msmt}}$ | $100$

Table 2:  Parameters for the qubit chain example discussed in the main text.
On the left, the true frequencies we want to discover; on the right, the known
parameters of the system.

It is also interesting to look at how the prediction of the response function
of a qubit system as in figure 4 improves after the measurements suggested by
an active strategy. As shown in figure 5, starting from a very uncertain
prediction, the distribution of possible response functions converges to the
true curve.

[width=]figures/app-response-evolution

Figure 5: Prediction of the response function after a different number of
steps. In orange, curves given by sampling parameters from the posterior; in
blue, the true response function. The initial best guess is given by sampling
from the prior distribution, before any measurement is performed.

## References

## References

  * [1] Kim Y, Eddins A, Anand S, et al. 2023 Nature 618 500–505 ISSN 0028-0836, 1476-4687 URL https://www.nature.com/articles/s41586-023-06096-3
  * [2] Wang J, Sciarrino F, Laing A, and Thompson M G 2020 Nature Photonics 14 273–284 ISSN 1749-4885, 1749-4893 URL https://www.nature.com/articles/s41566-019-0532-1
  * [3] Wetzstein G, Ozcan A, Gigan S, et al. 2020 Nature 588 39–47 ISSN 0028-0836, 1476-4687 URL https://www.nature.com/articles/s41586-020-2973-6
  * [4] Ningyuan J, Owens C, Sommer A, et al. 2015 Physical Review X 5 021031 ISSN 2160-3308 URL https://link.aps.org/doi/10.1103/PhysRevX.5.021031
  * [5] Gebhart V, Santagati R, Gentile A A, et al. 2023 Nature Reviews Physics 5 141–156 ISSN 2522-5820 URL https://www.nature.com/articles/s42254-022-00552-1
  * [6] Settles B 2012 Synthesis Lectures on Artificial Intelligence and Machine Learning 6 1–114 ISSN 1939-4608 URL https://www.morganclaypool.com/doi/abs/10.2200/S00429ED1V01Y201207AIM018
  * [7] Clerk A A, Devoret M H, Girvin S M, et al. 2010 Rev. Mod. Phys. 82 1155–1208 publisher: American Physical Society URL https://link.aps.org/doi/10.1103/RevModPhys.82.1155
  * [8] Foster A, Jankowiak M, Bingham E, et al. 2019 Variational Bayesian Optimal Experimental Design Advances in Neural Information Processing Systems vol 32 ed Wallach H, Larochelle H, Beygelzimer A, et al. (Curran Associates, Inc.) URL https://proceedings.neurips.cc/paper/2019/file/d55cbf210f175f4a37916eafe6c04f0d-Paper.pdf
  * [9] Kleinegesse S and Gutmann M U 2020 arXiv:2002.08129 [cs, stat] ArXiv: 2002.08129 URL http://arxiv.org/abs/2002.08129
  * [10] Cohn D, Atlas L, and Ladner R 1994 Machine Learning 15 201–221 ISSN 1573-0565 URL https://link.springer.com/content/pdf/10.1007/BF00993277.pdf
  * [11] Tong S 2001 ACTIVE LEARNING: THEORY AND APPLICATIONS (Stanford University) URL http://www.robotics.stanford.edu/ stong/papers/tong_thesis.pdf
  * [12] Chaloner K and Verdinelli I 1995 Statistical Science 10 ISSN 0883-4237 URL https://projecteuclid.org/journals/statistical-science/volume-10/issue-3/Bayesian-Experimental-Design-A-Review/10.1214/ss/1177009939.full
  * [13] Fredlund R, Everson R M, and Fieldsend J E 2010 A Bayesian framework for active learning The 2010 International Joint Conference on Neural Networks (IJCNN) (Barcelona, Spain: IEEE) pp 1–8 ISBN 978-1-4244-6916-1 URL http://ieeexplore.ieee.org/document/5596917/
  * [14] Ryan E G, Drovandi C C, McGree J M, and Pettitt A N 2016 International Statistical Review 84 128–154 ISSN 03067734 URL http://doi.wiley.com/10.1111/insr.12107
  * [15] Lindley D V 1956 The Annals of Mathematical Statistics 27 986–1005 ISSN 0003-4851 URL http://projecteuclid.org/euclid.aoms/1177728069
  * [16] Gal Y, Islam R, and Ghahramani Z 2017 arXiv:1703.02910 [cs, stat] ArXiv: 1703.02910 URL http://arxiv.org/abs/1703.02910
  * [17] Sanchez‐Lengeling B, Roch L M, Perea J D, et al. 2019 Advanced Theory and Simulations 2 1800069 ISSN 2513-0390, 2513-0390 URL https://onlinelibrary.wiley.com/doi/10.1002/adts.201800069
  * [18] Hentschel A and Sanders B C 2011 Physical Review Letters 107 233601 ISSN 0031-9007, 1079-7114 URL https://link.aps.org/doi/10.1103/PhysRevLett.107.233601
  * [19] Lennon D T, Moon H, Camenzind L C, et al. 2019 npj Quantum Information 5 79 ISSN 2056-6387 URL http://www.nature.com/articles/s41534-019-0193-4
  * [20] Nguyen V, Orbell S B, Lennon D T, et al. 2021 npj Quantum Information 7 100 ISSN 2056-6387 URL http://www.nature.com/articles/s41534-021-00434-x
  * [21] Hangleiter D, Roth I, Eisert J, and Roushan P 2021 Precise Hamiltonian identification of a superconducting quantum processor arXiv:2108.08319 [cond-mat, physics:physics, physics:quant-ph] URL http://arxiv.org/abs/2108.08319
  * [22] Fiderer L J, Schuff J, and Braun D 2021 PRX Quantum 2 020303 ISSN 2691-3399 URL https://link.aps.org/doi/10.1103/PRXQuantum.2.020303
  * [23] Valeri M, Polino E, Poderini D, et al. 2020 npj Quantum Information 6 92 ISSN 2056-6387 URL http://www.nature.com/articles/s41534-020-00326-6
  * [24] Cimini V, Valeri M, Polino E, et al. 2023 Advanced Photonics 5 ISSN 2577-5421
  * [25] Krenn M, Malik M, Fickler R, et al. 2016 Physical Review Letters 116 090405 ISSN 0031-9007, 1079-7114 URL https://link.aps.org/doi/10.1103/PhysRevLett.116.090405
  * [26] Melnikov A A, Nautrup H P, Krenn M, et al. 2018 Proceedings of the National Academy of Sciences 115 1221–1226 ISSN 0027-8424, 1091-6490 arXiv: 1706.00868 URL http://arxiv.org/abs/1706.00868
  * [27] Krenn M, Erhard M, and Zeilinger A 2020 Nature Reviews Physics 2 649–661 ISSN 2522-5820 URL https://www.nature.com/articles/s42254-020-0230-4
  * [28] Liang Q, Gongora A E, Ren Z, et al. 2021 npj Computational Materials 7 188 ISSN 2057-3960 URL https://www.nature.com/articles/s41524-021-00656-9
  * [29] Vargas-Hernández R A, Sous J, Berciu M, and Krems R V 2018 Physical Review Letters 121 255702 ISSN 0031-9007, 1079-7114 URL https://link.aps.org/doi/10.1103/PhysRevLett.121.255702
  * [30] Duris J, Kennedy D, Hanuka A, et al. 2020 Physical Review Letters 124 124801 ISSN 0031-9007, 1079-7114 URL https://link.aps.org/doi/10.1103/PhysRevLett.124.124801
  * [31] Leclercq F 2018 Physical Review D 98 063511 ISSN 2470-0010, 2470-0029 URL https://link.aps.org/doi/10.1103/PhysRevD.98.063511
  * [32] Goodfellow I, Bengio Y, and Courville A 2016 Deep Learning (MIT Press) ISBN 978-0-262-03561-3 URL 10.5555/3086952
  * [33] Dinh L, Krueger D, and Bengio Y 2015 arXiv:1410.8516 [cs] ArXiv: 1410.8516 URL http://arxiv.org/abs/1410.8516
  * [34] Kobyzev I, Prince S J D, and Brubaker M A 2021 IEEE Transactions on Pattern Analysis and Machine Intelligence 43 3964–3979 ISSN 0162-8828, 2160-9292, 1939-3539 arXiv: 1908.09257 URL http://arxiv.org/abs/1908.09257
  * [35] Papoulis A and Pillai S U 2009 Probability, random variables, and stochastic processes 4th ed (Boston, Mass.: McGraw-Hill) ISBN 978-0-07-122661-5 978-0-07-112256-6 oCLC: 255904469 URL https://doi.org/10.1115/1.3269815
  * [36] Kleinegesse S and Gutmann M U 2021 arXiv:2105.04379 [cs, stat] ArXiv: 2105.04379 URL http://arxiv.org/abs/2105.04379
  * [37] Foster A, Jankowiak M, O’Meara M, et al. 2020 arXiv:1911.00294
  * [38] Barber D and Agakov F V 2003 The IM algorithm: a variational approach to Information Maximization NIPS 2003
  * [39] Winkler C, Worrall D, Hoogeboom E, and Welling M 2019 arXiv:1912.00042 [cs, stat] ArXiv: 1912.00042 URL http://arxiv.org/abs/1912.00042
  * [40] Martín Abadi, Paul Barham, Eugene Brevdo, et al. 2015 TensorFlow: Large-Scale Machine Learning on Heterogeneous Systems (Google) URL https://www.tensorflow.org/
  * [41] Dillon J V, Langmore I, Tran D, et al. 2017 arXiv:1711.10604 [cs, stat] ArXiv: 1711.10604 URL http://arxiv.org/abs/1711.10604
  * [42] Hartmann M, Brandão F, and Plenio M 2008 Laser & Photonics Reviews 2 527–556 ISSN 18638880 URL https://onlinelibrary.wiley.com/doi/10.1002/lpor.200810046
  * [43] Houck A A, Türeci H E, and Koch J 2012 Nature Physics 8 292–299 ISSN 1745-2473, 1745-2481 URL https://www.nature.com/articles/nphys2251
  * [44] Altman E, Brown K R, Carleo G, et al. 2021 PRX Quantum 2 017003 ISSN 2691-3399 URL https://link.aps.org/doi/10.1103/PRXQuantum.2.017003
  * [45] Bernien H, Schwartz S, Keesling A, et al. 2017 Nature 551 579–584 ISSN 0028-0836, 1476-4687 URL https://www.nature.com/articles/nature24622
  * [46] Zhang J, Pagano G, Hess P W, et al. 2017 Nature 551 601–604 ISSN 0028-0836, 1476-4687 URL https://www.nature.com/articles/nature24654
  * [47] Ma R, Saxberg B, Owens C, et al. 2019 Nature 566 51–57 ISSN 0028-0836, 1476-4687 URL https://www.nature.com/articles/s41586-019-0897-9
  * [48] Zhang X, Kim E, Mark D K, et al. 2023 Science 379 278–283 ISSN 0036-8075, 1095-9203 URL https://www.science.org/doi/10.1126/science.ade7651
  * [49] Cohadon P F, Harris J, Marquardt F, and Cugliandolo L F 2020 Quantum optomechanics and nanomechanics: Lecture Notes of the Les Houches Summer School 1st ed Lnlh:c (New York: Oxford University Press) ISBN 978-0-19-882814-3
  * [50] Foster A, Ivanova D R, Malik I, and Rainforth T 2021 arXiv:2103.02438 [cs, stat] ArXiv: 2103.02438 URL http://arxiv.org/abs/2103.02438
  * [51] Sutton R S and Barto A G 2018 Reinforcement learning: an introduction second edition ed Adaptive computation and machine learning series (Cambridge, Massachusetts: The MIT Press) ISBN 978-0-262-03924-6 URL 10.5555/3312046
  * [52] Papamakarios G, Pavlakou T, and Murray I 2017 Masked Autoregressive Flow for Density Estimation Proceedings of the 31st International Conference on Neural Information Processing Systems NIPS’17 (Red Hook, NY, USA: Curran Associates Inc.) pp 2335–2344 ISBN 978-1-5108-6096-4 event-place: Long Beach, California, USA
  * [53] Aspelmeyer M, Kippenberg T J, and Marquardt F (eds) 2014 Cavity Optomechanics: Nano- and Micromechanical Resonators Interacting with Light 1st ed Quantum Science and Technology (Berlin, Heidelberg: Springer Berlin Heidelberg : Imprint: Springer) ISBN 978-3-642-55312-7