# Dynamic programming principle in cost-efficient sequential design:
application to switching measurements

Jeongmin Han
Department of Mathematics and Statistics, University of Jyväskylä
Juha Karvanen
Department of Mathematics and Statistics, University of Jyväskylä
and
Mikko Parviainen
Department of Mathematics and Statistics, University of Jyväskylä  The author
gratefully acknowledges financial support from the National Research
Foundation of Korea (NRF-2021R1A6A3A14045195).

###### Abstract

We study sequential cost-efficient design in a situation where each update of
covariates involves a fixed time cost typically considerable compared to a
single measurement time. The problem arises from parameter estimation in
switching measurements on superconducting Josephson junctions which are
components needed in quantum computers and other superconducting electronics.
In switching measurements, a sequence of current pulses is applied to the
junction and a binary voltage response is observed. The measurement requires a
very low temperature that can be kept stable only for a relatively short time,
and therefore it is essential to use an efficient design. We use the dynamic
programming principle from the mathematical theory of optimal control to solve
the optimal update times. Our simulations demonstrate the cost-efficiency
compared to the previously used methods.

Keywords: optimal design, D-optimality, complementary log-log link, Josephson
junction

##  1 Introduction

We study an optimal design problem motivated by switching measurements in
quantum physics. In these measurements, a Josephson junction (JJ) [8], an
important nonlinear component in superconducting electronics utilized for
example in quantum computers [11], is placed at a temperature close to
absolute zero to obtain superconductivity. When a current pulse is applied to
the junction, a voltage pulse is observed with a probability that depends on
the height of the current pulse. The relation between the height of the
current pulse and the probability of the voltage pulse can be accurately
approximated by a binary response model with a complementary log-log link
[10]. The unknown model parameters $a$ and $b$ determine the slope and the
location of the response curve and depend on the physical dimensions of the
junction. The problem is to choose the covariate values, i.e., the height of
the applied current pulse in our case, in order to estimate the parameters of
interest in the most efficient manner in the sense of D-optimality (see
Section 2). Efficiency is essential because in switching measurements, the
required temperature can be kept stable only for a relatively short time.

Sequential experimentation provides benefits over static designs. In nonlinear
models, the optimal design depends on the parameters. As the true values of
the parameters are unknown, the optimality properties of a design depend on
the accuracy of initial estimates of the parameters. A sequential design
allows us to update the estimates after each new observation and thus improve
the design. In practice, sequential experimentation has its price. In
switching measurements, adjusting the pulse generator to change the covariate
values takes a time that is considerably larger than the time needed for a
single measurement. This leads to cost-efficiency considerations where one
must balance the cost of an update and the information obtained from it. The
sequential design is then replaced by a batch-sequential design, and the key
problem is determining the optimal update times.

The batch-sequential design problem in the context of switching measurements
was formulated by Karvanen, Vartiainen, Timofeev and Pekola [10]. They
implemented an ad hoc method that increased the number of measurements by 10 %
at each stage. Later in [9], an approximate asymptotic model for accumulation
of a quantity approximating the D-criterion was presented. Moreover, an ad hoc
method for solving the optimal update times was implemented. The purpose of
the current paper is in the context of the approximate asymptotic model to
solve the optimal update times in a precise manner using the mathematical
theory of optimal control, and the dynamic programming principle (DPP) in
particular.

The DPP roughly speaking means that instead of considering a process until the
end we look at the problem over one step or stage. A simple example is given
by a random walk where we jump left or right each with probability $1/2$ on
the integers $\\{1,2,\ldots,9\\}$ and stop if we reach $0$ or $10$ where we
get some given payoff. When looking at the expectation in this process when
starting at $x_{0}$ (denoted by $u(x_{0})$), the DPP for the expectation reads
as

$\displaystyle u(x_{0})={\frac{1}{2}}\big{(}u(x_{0}-1)+u(x_{0}+1)\big{)}.$

In other words, this says that the expectation at a point can be considered by
looking at one step and summing up the outcomes with the corresponding
probabilities: we either go left with probability $1/2$ or right with
probability $1/2$ and take the expectations from there. More examples can be
found for example in [6, Section 10.3], [2] or [14, Chapter 1]. In the design
of experiments, dynamic programming has been applied to sequential
experimentation with linear models [3].

The key dynamic programming equations in our case are formulated in Sections
3.1 and 3.2. From these, we then solve the optimal update times using R [15].
In particular, we take into account which target (time or optimal accumulated
$D$) we are optimizing, and present the method to solve both of these
problems. Moreover, the dynamic programming principles automatically take into
account the amount of information we have obtained during the initialization
and the update cost when determining the optimal update times.

In Section 4, we present the results obtained through a simulation testing the
DPP method and previous methods in the literature. The results show that the
suggested method is more cost efficient. Finally, we remark that the use of
optimal control should find a wide variety of applications in cost-efficient
sequential designs with different statistical models. It is by no means
limited to the presented examples.

##  2 Preliminaries: $D$-optimal designs for complementary log-log regression

In this paper, we deal with the DPP in the context of $D$-optimal designs,
where D is the determinant or equivalently the square root of the determinant
of the information matrix. The preliminaries of the DPP are given in the
appendix and the next section: before we can write down the DPP we need to
recall the basic notation for the $D$-optimal design. The process of deriving
the (locally) $D$-optimal design involves maximizing the criterion $D$ under
certain assumptions, in particular assuming knowledge of the true parameter
values. A detailed exposition of this derivation in a generalized form is
presented in [7]. Moreover, the extensively employed logit model has been
explored in various works, including studies by [1], [13], [16], as well as
[12]. These references collectively demonstrate that the $D$-optimal design
represents a two-point design applicable to logit, probit, and cloglog models,
enabling the determination of optimal covariate values. Additionally,
extensions addressing models involving multiple covariates have been
investigated by [17], [5], and [4].

Here we consider an experiment, where our task is to estimate the unknown
parameters $\bm{\theta}=(a,b)$. We assume that the binary random variable is
distributed according to

$\displaystyle P(Y=1|\bm{\theta})=1-e^{-\exp(ax+b)},$ $\displaystyle
P(Y=0|\bm{\theta})=e^{-\exp(ax+b)},$ (2.1)

where $x$ is the covariate.

The task is to select the covariate values $x_{1},\ldots,x_{n}$ in such a way
that the parameters can be estimated from the observed data
$y_{1},\ldots,y_{n}$ as efficiently as possible. Starting from (2) the
likelihood reads for independent observations as

$\displaystyle L(\bm{y}|\bm{\theta})$
$\displaystyle:=\prod_{i=1}^{n}P(Y=y_{i}|\bm{\theta})$
$\displaystyle:=\prod_{i=1}^{n}\Big{\\{}y_{i}(1-e^{-\exp(ax_{i}+b)})+(1-y_{i})e^{-\exp(ax_{i}+b)}\Big{\\}}.$

Thus the log-likelihood is

$\displaystyle l(\bm{\theta}):$ $\displaystyle=\sum_{i=1}^{n}\log
P(Y=y_{i}|\bm{\theta})$
$\displaystyle=\sum_{i=1}^{n}\Big{\\{}y_{i}\log((1-e^{-\exp(ax_{i}+b)}))-(1-y_{i})\exp(ax_{i}+b)\Big{\\}}.$
(2.2)

The Fisher information matrix is defined to be

$\displaystyle J=-{\mathbb{E}\,}\begin{pmatrix}\frac{\partial^{2}l}{\partial
a^{2}}&\frac{\partial^{2}l}{\partial a\partial b}\\\
\frac{\partial^{2}l}{\partial b\partial a}&\frac{\partial^{2}l}{\partial
b^{2}}\end{pmatrix},$

assuming that we know $a$ and $b$ in order to compute the expectation. Observe
that the Fisher information depends on the sample size but not on the
measurements $y_{i}$. Heuristically speaking, the Fisher information measures
the curvature of the graph of log-likelihood and thus tells how stable the
maximum likelihood estimate is.

Computing the Fisher information, we get

$\displaystyle
J=\begin{pmatrix}\sum_{i=1}^{n}g(z_{i})x_{i}^{2}&\sum_{i=1}^{n}g(z_{i})x_{i}\\\
\sum_{i=1}^{n}g(z_{i})x_{i}&\sum_{i=1}^{n}g(z_{i})\end{pmatrix}$

where $z_{i}=ax_{i}+b$ and

$\displaystyle g(z)=\frac{e^{2z}}{e^{\exp(z)}-1}.$

We formulate the $D$-citerion in terms of

$\displaystyle D:=\sqrt{\det(J)}.$

In D-optimal design, our task is to find $x_{1},\ldots,x_{n}$ so that $D$ (or
equivalently $D^{2}$) is maximized. Here we selected the square root instead
of the determinant in order to have correct scaling properties later on. By
re-organizing the terms, we obtain

$\displaystyle D^{2}:=\det(J)$
$\displaystyle=\sum_{i=1}^{n}g(z_{i})x_{i}^{2}\sum_{i=1}^{n}g(z_{i})-\Big{(}\sum_{i=1}^{n}g(z_{i})x_{i}\Big{)}^{2}$
$\displaystyle=\sum_{i=1}^{n-1}\sum_{j=i+1}^{n}g(z_{i})g(z_{j})(x_{i}-x_{j})^{2}.$

The D-optimal design is a two-point design and we can maximize

$\displaystyle D^{2}=\det(J)=\frac{1}{a^{2}}g(z_{1})g(z_{2})(z_{1}-z_{2})^{2}$

to obtain the values

$\displaystyle(z_{1}^{*},z_{2}^{*})=(0.97963269129,-1.337736677)$

as well as $(z_{1}^{*},z_{2}^{*})=(-1.337736677,0.97963269129)$ [7]. The
Fisher information matrix corresponding to these values is denoted by $J^{*}$.
In other words, $J^{*}$ is the information matrix for $D$-optimal covariates
when $n=2$.

We use the approximation from [9] for the expected change in $D$ when two
additional measurements are made. Our intention here is not to study the model
there, but rather show, given such a model, how to solve the optimal update
times using the dynamic programming principle. In any case, we briefly recall
the setting in [9]. Similarly as there, for now, we assume that $(a,b)=(1,0)$
and only comment on the other cases later. We set

$\displaystyle J^{*}$
$\displaystyle=\sum_{i=1}^{2}g(z^{*}_{i})\begin{pmatrix}(x^{*}_{i})^{2}&x^{*}_{i}\\\
x^{*}_{i}&1\end{pmatrix},$

where $z^{*}_{i}=ax^{*}_{i}+b=x^{*}_{i}$, and observe that
$\sqrt{\det(J^{*})}=0.80940268$. Motivated by the asymptotic normality of
maximum likelihood estimates it was assumed that $\hat{a}$ and $\hat{b}$ are
normally distributed with the mean $(a,b)$ and with the covariance matrix
$J^{-1}(D)$, where $J(D)$ represents the accumulated information. Then in the
simulation pairs $(\hat{a},\hat{b})$ were generated from the distribution
$N((1,0),J^{-1}(D))$, and for each simulated $(\hat{a},\hat{b})$, the
following quantity was computed:

$\displaystyle\det\Bigg{(}\sum_{i=1}^{2}g(z_{i})\begin{pmatrix}x_{i}^{2}&x_{i}\\\
x_{i}&1\end{pmatrix}\Bigg{)}^{1/2}$
$\displaystyle=\sqrt{g(z_{1})g(z_{2})}(z_{1}-z_{2}),$

where

$\displaystyle z_{1}$
$\displaystyle=\frac{z_{1}^{*}-\hat{b}}{\hat{a}},\quad\quad
z_{2}=\frac{z_{2}^{*}-\hat{b}}{\hat{a}}.$

The expected change in $D$ was taken to be the average of these simulated
quantities. Since the model is founded on such an accumulation property that
holds asymptotically (similarly to the previous steps) we give a brief
justification for it in Appendix A. The expected change given by the
simulation in [9] was

$\displaystyle h(D)=\frac{0.80940268}{1+e^{\eta+\theta\log(D)}},$ (2.3)

where $\eta=1.88938,\ \theta=-1.51330$. Since $\lim_{D\to 0+}h(D)=0$ we also
need an initialization procedure that gives us the initial $D_{0}>0$ and the
first MLE to get started.

##  3 Optimal update times

Given the model described in the previous section, we now show how to solve
optimal update times using dynamic programming. We consider two problem
settings: given the available time, we want to maximize the accumulated $D$,
or given the desired $D$, we want to minimize the time used to reach this
target. The dynamic programming principles below can easily be solved
numerically by iterating backwards from the final condition.

###  3.1 Given final time, maximize the accumulated $D$

We consider a situation where the termination time $T$ is given, and we want
to maximize the accumulated $D$. For example, it is known that the Josephson
junction remains at a required temperature and thus is stable for a certain
time $T$. Time cost $C_{s}$ for updating the covariate is assumed to be a
nonnegative integer and is here usually large compared to the cost of a single
measurement (the cost of the single measurement is scaled to be 1). For
example, $C_{s}$ could be the time cost for adjusting the current in the pulse
generator. The task is to find the optimal update times so as to accumulate as
large $D$ as possible by the termination time $T$. The problem is to decide
the next update time at the current update time $t$ assuming that the
accumulated $D\geq D_{0}>0$ is given. Let us denote by

_$u(D,t)$ , amount of $D$ we can optimally accumulate during the rest of the
experiment if starting at $(D,t)$ and updating immediately._

This description is enough for now, but the precise definition is given in
Appendix B. Using the terminology from the control theory, the value function
$u$ satisfies the DPP with $C_{s}>0$ and $t<T-C_{s}$

$\displaystyle u(D,t)$
$\displaystyle=\sup_{t_{\text{next}}\in(t+C_{s},T]\cap{\mathbb{Z}}_{+}}\Big{\\{}u(D_{\text{next}},t_{\text{next}})+(D_{\text{next}}-D)\Big{\\}},$
(3.4)

where

$\displaystyle t,C_{s}$ $\displaystyle\in{\mathbb{Z}}_{+},\
{\mathbb{Z}}_{+}=\\{0,1,2,\ldots\\},$ $\displaystyle D_{\text{next}}:$
$\displaystyle=D_{\text{next}}(t_{\text{next}},D):=D+h(D)(t_{\text{next}}-t-C_{s}),$
$\displaystyle u(D,t)$ $\displaystyle=0,\quad
t\in[T-C_{s},T]\cap{\mathbb{Z}}_{+}.$

What the DPP says is that we can compute the optimal value at a point, by
taking one optimal stage of measurements and the optimal value from where we
ended up during this one stage.

Also observe that it does not matter that it was assumed that $a=1$ in $h$
since everything would be just scaled by the true value of $a$, and this does
not affect the update times. However, inside the DPP, the quantity $D$
corresponds to $a=1$ and should not be taken as an approximation of the ’true
$D$’. Nonetheless, the observed $D$ can easily be computed from the sample.

We will solve for $u$ starting from $t=T-C_{s}-1$. When at $(D,t)$, supposing
that we have already solved values $u$ for all $s>t$, we may solve for
$u(D,t)$, using (3.4). In the course of using (3.4), we also solve and save
the stage update times for each $D$. Then we pick the stage update time
$t_{next}$ for $(D,t)=(D_{0},0)$, and by continuing, for the rest of the
stages. Figure 1 gives two examples.

(a) Fixed $D_{0}$ case ($D_{0}=0.5$)

(b) Fixed $C_{s}$ case ($C_{s}=10$)

Figure 1: Update times with $T=1000$ (the point $T=1000$ is not an update time
but shows the final accumulated $D$). The left subfigure indicates the update
times for fixed $D_{0}=0.5$ ($C_{s}=5$ in blue, $C_{s}=30$ in red). The right
one represents those times for fixed $C_{s}=10$ ($D_{0}=0.5$ in cyan,
$D_{0}=5$ in yellow).

The optimization is implemented in the R codes available at optimaltimesCsA.r
and optimaltimesCsB.r. In these R-codes, for example, optimaltimesCs.r, we
have discrete grids for $D_{i},\ i=1,\ldots,N_{D}$ and $t_{j},\
j=1,\ldots,N_{T}$. When at a point $(D_{i_{0}},t_{j_{0}})$, we consider the
time points
$\tilde{t}_{j}\in\\{\tilde{t}_{1},\ldots,\tilde{t}_{N_{t_{j_{0}}}}\\}:=\\{t_{j_{0}}+C_{s},\ldots,t_{N_{T}}=T\\}$
and compute the values

$\tilde{D}_{\tilde{t}_{j}}:=D_{k}+h(D_{k})(\tilde{t}_{j}-t_{j_{0}}-C_{s}).$

On the other hand, we have already computed all the future values for $u$ at
the discretized values of $D_{i},\ i=1,\ldots,N_{D}$. Thus rounding
$\tilde{D}_{t_{j}}$ to the nearest grid value $D_{k}$, we know
$u(D_{k},\tilde{t}_{j})$ and we may form

$\displaystyle
DPP_{i_{0},j}:=u(D_{k},\tilde{t}_{j})+(\tilde{D}_{\tilde{t}_{j}}-D_{i_{0}}),\quad\text{for
all}\quad j=1,\ldots,N_{t_{j_{0}}}.$

Now, according to (3.4), we have solved the values of $u$ at time $t_{j_{0}}$.
Then we continue to time $t_{j_{0}}-1$.

###  3.2 Given desired $D_{\text{final}}$, minimize the time cost

In this section, we consider a different problem where the goal is to reach a
given final value of accuracy $D_{\text{final}}$ in the shortest possible
time. Let us denote by

$v(D)$, _the minimal time to reach $D_{\text{final}}$ if we update
immediately._

Again for a more detailed discussion, see Appendix B. First, as indicated
above, we observe that $v$ does not depend on the time spent so far since we
only look at the future accumulation of $D$. We still denote the change of $D$
in one double measurement by $h(D)$ given in (2.3), and by
$C_{s}\in{\mathbb{Z}}_{+}$ the update cost of the covariates. Now $v$
satisfies the following DPP

$\displaystyle v(D)$ $\displaystyle=\inf_{\Delta
t\in(C_{s},\infty)\cap{\mathbb{Z}}_{+}}\Big{\\{}v\big{(}D+(\Delta
t-C_{s})h(D)\big{)}+\Delta t\Big{\\}},$ (3.5) $\displaystyle v(D)$
$\displaystyle=0\quad\text{ for }\quad D\geq D_{\text{final}}.$

Recall that we only applied $(a,b)=(1,0)$ in the previous subsection, since
the accumulated information in the general case can be obtained by scaling
everything by the factor $a^{-2}$ with the true value of $a$. One should
observe that also in the above DPP $D$ and $D_{\text{final}}$ denote the
rescaled versions with $(a,b)=(1,0)$ and should not be taken as an
approximation of the ’true D’. Also, the updates should not be launched by the
observed $D$ computed from the sample but by the time steps $\Delta t$ that
will be stored while numerically solving for $v$.

To finish the section, observe that the optimal update times depend (as they
should) on whether we are minimizing the time cost or maximizing $D$, what are
the final values of $T$ or $D_{\text{final}}$, the update time cost $C_{s}$,
and $D_{0}$ after the initialization. Thus it is not practical to give update
times in a table, but they can be solved for different values by using the
provided codes.

##  4 Benchmark

In preceding studies [10, 9], several types of batch-sequential $D$-optimal
design methods were presented. They used the measurement data for a JJ circuit
given in https://github.com/JuhaKarvanen/switching_measurements_design. Here
we simulate the experiment from scratch. We also use less intense
initialization with bad initial guesses of $a$ and $b$ to obtain a lower value
for $D_{0}$ and it can be reproduced through the code initialdatagen.r for the
case when the stage size is 100 in the first round under the assumption
$(a,b)=(0.24,-61)$. This is in order to illustrate the difference of the
methods more clearly (the lower $D_{0}$, the more important the updates as
well as cost-efficient designs are). In the code simulation.r, we first solved
the DPP problem (3.4) (given final time, optimize accumulated D) with $T=3500$
and $C_{S}=0.88167/0.00386\approx 228.4$ and then obtained the optimal update
times under our setting. After that, we ran 100 simulations based on those
update times. We also conducted the same number of simulations for the model
suggested in [9] for comparison. The results are illustrated in Figure 2. As a
technical remark, $D$ in the figure is the actual observed $D$ computed from
the sample, not the ‘model $D$’ that appears for example in the DPPs and is
obtained using $h(D)$ and rescaling $(a,b)=(1,0)$. What we observe is that the
DPP method is more cost-efficient than the model in [9].

Figure 2: Simulations for three methods with $D_{0}=0.1408$ and
$(a,b)=(0.24,-61)$. The black line marks the simulated results obtained in the
same way as the ad hoc method used in the original experiment in [10]. The
blue line and cyan line represent the time costs obtained by the approximate
cost-efficient method in [9] and the DPP method of this paper, respectively.
The blue and cyan squares represent the path with the corresponding experiment
design calculated as the median of 100 simulated paths. The dots show the
stages of the 100 simulated paths (red: [9] design, yellow: DPP design). The
dashed line indicates the theoretical limit of the time cost if they are
immediately measured at the optimal points. In the DPP design, we set $T=3500$
and the dotted line represents the total time cost (=17.42268 sec) in this
case.

##  5 Discussion

We considered cost-efficient batch-sequential design in a situation where the
update cost is large compared to the cost of a single measurement. The work
was motivated by switching measurements of Josephson junctions. We made
several improvements to the earlier literature: We formulated and solved a
precise dynamic programming problem for update times if the change in $D$ is
given by the function $h(D)$ obtained from the literature. In particular, the
dynamic programming approach automatically takes into account which target
(target time or target $D$) we are optimizing, and the effect of $D_{0}$,
$C_{S}$ on update times. The simulations illustrate the efficiency of the
proposed approach. Compared to earlier approximations [10, 9], a clear gain in
cost-efficiency was seen.

However, there are still several approximations stemming from the earlier
literature. The particular assumed form of the Fisher information in the
simulation producing $h$, the normality assumption of MLE, and the assumed
accumulation property of $D$, are such approximations and only hold
asymptotically. The reason for these approximations is of course a desire to
formulate the problem in such a way that it can be simulated and solved by
efficient tools and in advance.

The proposed approach is directly applicable to switching measurements. In
addition, a similar approach could be useful in other applications where the
data are collected in batches and a cost is associated with each update. The
functions $h(D)$ have been estimated also for the logit and probit link
functions [9]. One might also want to consider the DPP approach with other
optimality criteria.

##  Appendix A Accumulation of $D$ in the approximative model

To obtain an approximative model for the accumulation of $D$, the expected
change in $D$ in [9] was taken to be the average of the simulated quantities
$\sqrt{\det(J_{i})}$. Since the model is founded on such an accumulation
property that holds asymptotically (similarly to the previous steps in the
model) we give a brief justification for it here. To be more precise, we want
to show in a suitable mode of stochastic convergence that

$\displaystyle\sqrt{\det(J_{1}+\ldots+J_{n})}-\sqrt{\det(J_{1}+\ldots+J_{n-1})}-\sqrt{\det(J_{n})}\to
0$ (1.6)

almost surely as $n\to\infty$, where $J_{i}=J_{i}(\omega)$ denotes the
information obtained in the round $i$ for each sample $\omega$ (we use the
notation $J_{i}$ for the convenience). In other words, the difference between
the first two terms reflects the actual change in $D$ if we obtain the
information $J_{n}$, and the last term the quantity used in the simulation. To
establish the above asymptotic behaviour, let us denote $\Delta
J_{i}:=(J_{i}-J^{*})$ and compute

$\displaystyle\sqrt{\det(J_{1}+\ldots+J_{n})}=\det\Big{(}nJ^{*}+\sum_{i=1}^{n}\Delta
J_{i}\Big{)}^{1/2}$
$\displaystyle=n\det\Big{(}J^{*}+\frac{1}{n}\sum_{i=1}^{n}\Delta
J_{i}\Big{)}^{1/2}$
$\displaystyle=(n-1)\det\Big{(}J^{*}+\frac{1}{n}\sum_{i=1}^{n}\Delta
J_{i}\Big{)}^{1/2}+\det\Big{(}J^{*}+\frac{1}{n}\sum_{i=1}^{n}\Delta
J_{i}\Big{)}^{1/2}.$

To get an approximation of

$\sqrt{\det(J_{1}+\ldots+J_{n})}-\sqrt{\det(J_{1}+\ldots+J_{n-1})},$

in order to prove (1.6) we also observe that

$\displaystyle\sqrt{\det(J_{1}+\ldots+J_{n-1})}=(n-1)\det\Big{(}J^{*}+\frac{1}{n-1}\sum_{i=1}^{n-1}\Delta
J_{i}\Big{)}^{1/2}$
$\displaystyle=(n-1)\det\Big{(}J^{*}+\frac{1}{n}\sum_{i=1}^{n}\Delta J_{i}$
$\displaystyle\hskip 70.0001pt+\frac{1}{n-1}\Big{(}\sum_{i=1}^{n-1}\Delta
J_{i}-\sum_{i=1}^{n}\Delta
J_{i}\Big{)}+\Big{(}\frac{1}{n-1}-\frac{1}{n}\Big{)}\sum_{i=1}^{n}\Delta
J_{i}\Big{)}^{1/2}$
$\displaystyle=(n-1)\det\Big{(}J^{*}+\frac{1}{n}\sum_{i=1}^{n}\Delta
J_{i}-\frac{1}{n-1}\Delta J_{n}+\frac{1}{n(n-1)}\sum_{i=1}^{n}\Delta
J_{i}\Big{)}^{1/2}$ $\displaystyle=:(n-1)\det(A+B)^{1/2},$

where the shorthand notation $A$ contains the first two terms and $B$ the last
two terms. Also, observe that since we are working under the assumption
$J_{i}\to J^{*}$ almost surely as $i\to\infty$, it holds that

$\displaystyle(n-1)\det(B),\ (n-1)\operatorname{tr}(B)\to 0$ (1.7)

almost surely as $n\to\infty$; we will need these facts later. Now we use a
well-known formula for the determinant of a sum of $2\times 2$ matrices, and
the mean value theorem or representation for the remainder term in Taylor’s
theorem as in

$\displaystyle\sqrt{x_{0}+h}=\sqrt{x_{0}}+(2c)^{-1/2}h,\text{ for some
}c\text{ between }\text{$x_{0}$ and $x_{0}+h$}.$

In this way

$\displaystyle\det(A+B)^{1/2}$
$\displaystyle=\big{(}\det(A)+\det(B)+\det(A)\operatorname{tr}(A^{-1}B)\big{)}^{1/2}$
$\displaystyle=\det(A)^{1/2}+(2c)^{-1/2}\big{(}\det(B)+\det(A)\operatorname{tr}(A^{-1}B)\big{)}.$

Using this and (1.7), we get

$\displaystyle(n-1)$
$\displaystyle\bigg{(}\det\Big{(}J^{*}+\frac{1}{n}\sum_{i=1}^{n}\Delta
J_{i}\Big{)}^{1/2}-\det\Big{(}J^{*}+\frac{1}{n-1}\sum_{i=1}^{n-1}\Delta
J_{i}\Big{)}^{1/2}\bigg{)}$
$\displaystyle=(n-1)\bigg{(}\det(A)^{1/2}-\det(A+B)^{1/2}\big{)}$
$\displaystyle=-(n-1)(2c)^{-1/2}\big{(}\det(B)+\det(A)\operatorname{tr}(A^{-1}B)\big{)}\to
0,$

almost surely as $n\to\infty$, since

$\left|\operatorname{tr}(A^{-1}B)\right|\leq\left|\left|A^{-1}\right|\right|_{\infty}\operatorname{tr}(B).$

Combining the above estimates, we deduce

$\displaystyle\sqrt{\det(J_{1}+\ldots+J_{n})}-\sqrt{\det(J_{1}+\ldots+J_{n-1})}-\sqrt{\det(J_{n})}\to
0$

almost surely when $n\to\infty$, as desired.

##  Appendix B Mathematical remarks on dynamic programming

In this section, we review part of the mathematical background for the optimal
control problem in question, see for example [6, 2]. The literature is vast
and it is impossible to give a thorough overview here.

###  B.1 Given final time, maximize the accumulated $D$ problem

Let $C_{s}>0$. The function $u:(0,\infty)\times[0,T]\to{\mathbb{R}}$ denotes
the optimal amount $D$ we can accumulate during the remaining time if updating
at the considered time. To be more precise, the controller chooses the update
times

$\displaystyle{\mathbf{t}}:=(t_{0},\ldots,t_{n+1}),$

where the difference between the update times is larger than $C_{s}$, and
otherwise the number $n$ may vary. Also, observe that $t_{0}$ and $t_{n+1}$
must correspond to the interval we are considering. For example, if we
consider the full interval, then $t_{0}=0$ and $t_{n+1}=T$, but this will be
clear from the context. The definition of $u$ now reads as

$\displaystyle u(D,t)$
$\displaystyle=\begin{cases}\sup_{{\mathbf{t}}}\sum_{i=0}^{n}(t_{i+1}-t_{i}-C_{s})h(D(t_{i})),&t<T-C_{s}\\\
0,&t\in[T-C_{s},T],\end{cases}$ (2.8)

where

$\displaystyle D^{\prime}(s)$ $\displaystyle=\begin{cases}h(D(t_{i})),&\text{
if }s\in[t_{i}+C_{s},t_{i+1}),\\\ 0,&\text{ otherwise}.\end{cases}$ (2.9)

Obviously, $u$ given by (2.8) exists and is unique.

######  Lemma B.1 (DPP)

Let $u$ be as in (2.8). Then for $C_{s}>0$ and $t<T-C_{s}$

$\displaystyle u(D,t)$
$\displaystyle=\sup_{t_{\text{next}}\in(t+C_{s},T]}\Big{\\{}u(D_{\text{next}},t_{\text{next}})+(D_{\text{next}}-D)\Big{\\}},$
(2.10)

where

$\displaystyle D_{\text{next}}:$
$\displaystyle=D+h(D)(t_{\text{next}}-t-C_{s}).$

The validity of the DPP can be readily verified by plugging in the definition
(2.8). One could also study the DPP on its own right without invoking (2.8)
and showing existence and uniqueness using a similar idea as in the numerical
computation: since $C_{s}>0$, iterating backwards from the end $T$, we can
extend the values to each strip of width $C_{s}$.

######  Remark B.2 (Regularity, stability of update times)

Similarly as in [6, Section 10.3], we can also prove that $u$ is bounded and
Lipschitz continuous in $(0,\infty)\times[0,T]$. From the definition of $h$
and (2.8), it is straightforward that

$u(D,t)\leq h^{\ast}(T-C_{s}),$

where $h^{\ast}=\sqrt{\det(J^{*})}=0.80940268$ is the maximum value of $h$.

The basic idea in the proof of Lipschitz continuity is to ‘copy’ the update
times from a nearby point of interest. We will show

$|u(D,t)-u(\tilde{D},\tilde{t})|\leq C(|D-\tilde{D}|+|t-\tilde{t}|)$

for some universal $C>0$ and any
$(D,t),(\tilde{D},\tilde{t})\in(0,\infty)\times[0,T]$. We first estimate

$|u(D,t)-u(\tilde{D},t)|$

and without loss of generality we may assume that $\tilde{D}\geq D$ and
$u(\tilde{D},t)\geq u(D,t)$. To this end, for $t\in[0,T]$, we can choose
update times $\mathbf{t}_{1}=(t_{0},\dots,t_{n})$ with $t=t_{0}$ which
satisfies

$u(\tilde{D},t)\leq\sum_{i=0}^{n}(t_{i+1}-t_{i}-C_{s})h(\tilde{D}(t_{i}))+{\varepsilon}$

for a fixed ${\varepsilon}>0$. Then we have

$\displaystyle u(\tilde{D},t)-u(D,t)$
$\displaystyle\leq\bigg{(}\sum_{i=0}^{n}(t_{i+1}-t_{i}-C_{s})h(\tilde{D}(t_{i}))+{\varepsilon}\bigg{)}-\sum_{i=0}^{n}(t_{i+1}-t_{i}-C_{s})h(D(t_{i}))$
$\displaystyle=\sum_{i=0}^{n}(t_{i+1}-t_{i}-C_{s})\big{(}h(\tilde{D}(t_{i}))-h(D(t_{i}))\big{)}+{\varepsilon}.$

Since

$\displaystyle(\tilde{D}-D)^{\prime}(s)$
$\displaystyle=\begin{cases}h(\tilde{D}(t_{i}))-h(D(t_{i})),&\text{ if
}s\in[t_{i}+C_{s},t_{i+1}),\\\ 0,&\text{ otherwise}\end{cases}$

and $h$ is Lipschitz, we observe that for any $s\in(t,T]$,

$|(\tilde{D}-D)(s)|\leq C|\tilde{D}-D|$

for some $C>0$ independent of ${\varepsilon}$ by Gronwall’s inequality. By
employing this and Lipschitz regularity, we can see that

$u(\tilde{D},t)-u(D,t)\leq C|\tilde{D}-D|$

for some $C>0$.

To get an estimate for $t$-direction, it is enough to show

$u(D,t)-u(D,\tilde{t})\leq C(\tilde{t}-t)$

for $t<\tilde{t}$, since $u$ is decreasing in $t$. We use update times
$\mathbf{t}_{2}=(s_{0},\dots,s_{k}),\ s_{0}=t$ such that

$u(D,t)\leq\sum_{i=0}^{k}(s_{i+1}-s_{i}-C_{s})h(D(s_{i}))+{\varepsilon}$

for some ${\varepsilon}>0$. We set update times
$\mathbf{\tilde{t}}_{2}=(\tilde{s}_{0},\dots,\tilde{s}_{k})$ so that updates
always occur at the same values of $D$ as in the process starting at $t$
before hitting $T$, that is,

$\tilde{s}_{i}=\min\\{\tilde{t}-t+s_{i},T\\}$

for each $i=0,\dots,k$. Now we have

$\displaystyle u(D,t)-u(D,\tilde{t})$
$\displaystyle\leq\bigg{(}\sum_{i=0}^{k}(s_{i+1}-s_{i}-C_{s})h(D(s_{i}))+{\varepsilon}\bigg{)}-\sum_{i=0}^{k}(\tilde{s}_{i+1}-\tilde{s}_{i}-C_{s})h(D(s_{i}))$
$\displaystyle\leq h^{\ast}(\tilde{t}-t)+{\varepsilon}$

and then we can complete the proof by letting ${\varepsilon}\to 0$.

This proof also has an implication that may be of independent interest: Since
the copied strategy provides almost as good (in Lipschitz sense) outcome in
the nearby points, small variations (numerical and other approximation errors)
in the update times are not crucial. In other words, the update times are
stable in this sense and this provides a kind of theoretical basis for the
practical solving of the problem.

It is also rather easy to see that

$\displaystyle h(D)(T-t-C_{s})\leq u_{C_{s}}(D,t)\leq u_{0}(D,t+C_{s}),$
(2.11)

where the subindex in $u_{C_{s}}$ and $u_{0}$ denote the solution to (2.10)
with update cost $C_{s}$ or $0$. Intuitively, we can understand the upper
bound in (2.11) as the case of updating $D$ every time with $C_{s}=0$.
Similarly, the lower bound can be regarded as the case of no update for $D$.

Next, we observe that we cannot derive a partial differential equation from
the DPP by passing to a limit with step size as usual, due to the lower bound
given by the update cost $C_{s}$. Nonetheless, we can derive a partial
differential equation if we consider a continuous time approximation with no
update cost, and this can be solved explicitly. Recall that $u(D,t)$ denotes
the amount of $D$ we can optimally accumulate during the rest of the
experiment if starting at $(D,t)$ and updating immediately. When $C_{s}=0$,
and the time step is 1, it holds that

$\displaystyle u(D,t)=u(D+h(D),t+1)+h(D),$ (2.12)

which is easy to show by induction. This means that updates of $D$ occur after
every measurement. We can instead of (2.12), consider the following
approximation where the time step is ${\varepsilon}>0$

$\displaystyle
u(D,t)=u(D+{\varepsilon}h(D),t+{\varepsilon})+{\varepsilon}h(D),$

and then pass to a continuous time limit. From this we obtain

$\displaystyle\frac{u(D+{\varepsilon}h(D),t+{\varepsilon})-u(D,t)}{{\varepsilon}}=-h(D),$

and if $u\in C^{1}$ passing to ${\varepsilon}\to 0$, we get

$\displaystyle h(D)u_{D}(D,t)+u_{t}(D,t)=-h(D)\qquad\textrm{in}\
(0,\infty)\times(0,T)$ (2.13)

with a boundary condition $u(D,T)=0$, where now the time variables are
continuous. This equation could be solved by the standard method of
characteristics, see for example [6]. The limit is also a solution to the DPP
for all $\delta>0$ such that $0\leq t,t+\delta\leq T$

$\displaystyle\begin{cases}u(D,t)&=\int_{t}^{t+\delta}h(D(s))\,ds+u(D(t+\delta),t+\delta),\\\
u(D,T)&=0,\\\ \end{cases}$

with the dynamics

$\displaystyle D^{\prime}(s)=h(D(s)),\quad D(t)=D.$

In this case, it might be instructive to think $h(D)$ as a density of
accumulation of data so that $h(D(s))\,ds$ gives the accumulated data over a
short time interval.

###  B.2 Given desired $D_{\text{final}}$, minimize the time cost

Here we consider a continuous time version of (3.5) i.e. $v$ now denotes the
minimal time cost to reach the given $D_{\text{final}}$ if updating
immediately. This time cost does not depend on the time spent so far but only
$D$ accumulated so far (this is ’model $D$’ obtained using $h(D)$ with
$(a,b)=(1,0)$). Thus $v$ does not depend on time. We denote

$\displaystyle{\mathbf{D}}=(D_{0},\ldots,D_{n+1})$

where $D_{0}$ and $D_{n+1}$ correspond the interval we are considering. Now
$v:(0,\infty)\to[0,\infty)$ can be defined as

$\displaystyle v(D)$
$\displaystyle=\begin{cases}\inf_{{\mathbf{D}}}\sum_{i=0}^{n}\big{\\{}C_{s}+(D_{i+1}-D_{i})h^{-1}(D_{i})\big{\\}},&D<D_{\text{final}}\\\
0,&t\in[D_{\text{final}},\infty).\end{cases}$ (2.14)

Also observe that we can switch between ${\mathbf{t}}$ and ${\mathbf{D}}$:
moving from ${\mathbf{t}}$ to ${\mathbf{D}}$, recall that $D(s)$ may be solved
from (2.9), and thus we can solve $D_{i}$ corresponding to $t_{i}$. Knowing
$D(t)$, we may also solve the inverse function $t(D)$ by selecting for example
the left continuous version of it at update points where it jumps. Using this
function we may then move from ${\mathbf{D}}$ to ${\mathbf{t}}$. Again $v$
defined by (2.14) exists, is unique, and satisfies the DPP.

######  Lemma B.3 (DPP)

Let $v$ be as in (2.14). Then

$\displaystyle v(D)$ $\displaystyle=\inf_{\Delta
t\in(C_{s},\infty)}\Big{\\{}v\big{(}D+(\Delta t-C_{s})h(D)\big{)}+\Delta
t\Big{\\}},$ (2.15) $\displaystyle v(D)$ $\displaystyle=0\quad\text{ for
}\quad D\geq D_{\text{final}}.$

We see the validity of the DPP directly by plugging in (2.14).

######  Remark B.4 (Regularity, stability of update times)

The solution $v$ is bounded and Lipschitz continuous in $(0,\infty)$. The idea
is again to ’copy’ the update points from a nearby point of interest.
Similarly, the proof also implies that the update points are stable which is
important for approximations.

We also see that the following estimate holds

$v(D)\leq\frac{D_{\text{final}}-D}{h(D)}+C_{s}$

by considering the situation that an update is done only once.

When $C_{s}=0$, we update $D$ every time to minimize the time cost so that the
discrete version (3.5) becomes

$\displaystyle v(D)$ $\displaystyle=v(D+h(D))+1$

since $h$ is increasing. In order to pass to a continuous time, we can
consider

$\displaystyle v(D)$ $\displaystyle=v(D+h(D){\varepsilon})+{\varepsilon}$

for ${\varepsilon}>0$. By passing to a limit ${\varepsilon}\to 0$, we obtain
the equation

$\displaystyle\begin{cases}h(D)v_{D}(D)&=-1\quad\textrm{in}\
(0,D_{\text{final}})\times(0,\infty)\\\ v(D_{\text{final}})&=0.\end{cases}$

The solution also satisfies the DPP

$\displaystyle\begin{cases}v(D)&=t+h+v(D(t+h)),\\\ v(D_{\text{final}})&=0,\\\
\end{cases}$

with the dynamics

$\displaystyle D^{\prime}(s)=h(D(s)),\quad D(t)=D,\quad D<D_{\text{final}}.$

SUPPLEMENTARY MATERIAL

ZIP file supplement:

This supplement contains the R files for reproducing the simulations presented
in the paper. A readme file is provided to connect the files with the
corresponding simulations.

## References

  * [1] K.. Abdelbasit and R.. Plackett  “Experimental design for binary data”  In _Journal of the American Statistical Association_ 78.381 Taylor & Francis, 1983, pp. 90–98
  * [2] M. Bardi and I. Capuzzo-Dolcetta  “Optimal control and viscosity solutions of Hamilton-Jacobi-Bellman equations”  Birkhäuser, 1997
  * [3] I. Ben-Gal and M. Caramanis  “Sequential DOE via dynamic programming”  In _IIE Transactions_ 34.12, 2002, pp. 1087–1100
  * [4] R. Dorta-Guerra, E. González-Dávila and J. Ginebra  “Two-level experiments for binary response data”  In _Computational Statistics & Data Analysis_ 53.1 Elsevier, 2008, pp. 196–208
  * [5] H.. Dror and D.. Steinberg  “Sequential experimental designs for generalized linear models”  In _Journal of the American Statistical Association_ 103.481 Taylor & Francis, 2008, pp. 288–298
  * [6] L.. Evans  “Partial differential equations” 19, Graduate Studies in Mathematics  American Mathematical Society, Providence, RI, 2010
  * [7] I. Ford, B. Torsney and C.F.. Wu  “The use of a canonical form in the construction of locally optimal designs for non-linear problems”  In _Journal of the Royal Statistical Society Series B: Statistical Methodology_ 54.2 Oxford University Press, 1992, pp. 569–583
  * [8] B.. Josephson  “Possible new effects in superconductive tunnelling”  In _Physics letters_ 1.7 Elsevier, 1962, pp. 251–253
  * [9] J. Karvanen  “Approximate cost-efficient sequential designs for binary response models with application to switching measurements”  In _Comput. Statist. Data Anal._ 53.4, 2009, pp. 1167–1176
  * [10] J. Karvanen, J.. Vartiainen, A. Timofeev and J. Pekola  “Experimental designs for binary data in switching measurements on superconducting Josephson junctions”  In _Journal of the Royal Statistical Society: Series C (Applied Statistics)_ 56.2 Wiley Online Library, 2007, pp. 167–181
  * [11] Y. Makhlin, G. Schön and A. Shnirman  “Quantum-state engineering with Josephson-junction devices”  In _Reviews of modern physics_ 73.2 APS, 2001, pp. 357
  * [12] T. Mathew and B.. Sinha  “Optimal designs for binary data under logistic regression”  In _Journal of Statistical Planning and Inference_ 93.1-2 Elsevier, 2001, pp. 295–307
  * [13] S. Minkin  “Optimal designs for binary data”  In _Journal of the American Statistical Association_ 82.400 Taylor & Francis, 1987, pp. 1098–1103
  * [14] M. Parviainen  “Notes on tug-of-war games and the $p$-Laplace equation” To appear in SpringerBriefs on PDEs and Data Science  In _arXiv preprint arXiv:2208.09732_ , 2022
  * [15] R Core Team  “R: A Language and Environment for Statistical Computing”, 2023  R Foundation for Statistical Computing  URL: https://www.R-project.org/
  * [16] R.. Sitter and C… Wu  “Optimal designs for binary response experiments: Fieller, D, and A criteria”  In _Scandinavian Journal of Statistics_ JSTOR, 1993, pp. 329–341
  * [17] D.. Woods, S.. Lewis, J.. Eccleston and K.. Russell  “Designs for generalized linear models with several variables and model uncertainty”  In _Technometrics_ 48.2 Taylor & Francis, 2006, pp. 284–292


# Experimental graybox quantum control

Akram Youssry  Quantum Photonics Laboratory and Centre for Quantum Computation
and Communication Technology, RMIT University, Melbourne, VIC 3000, Australia
Yang Yang  Quantum Photonics Laboratory and Centre for Quantum Computation and
Communication Technology, RMIT University, Melbourne, VIC 3000, Australia
Robert J. Chapman  Quantum Photonics Laboratory and Centre for Quantum
Computation and Communication Technology, RMIT University, Melbourne, VIC
3000, Australia  ETH Zurich, Optical Nanomaterial Group, Institute for Quantum
Electronics, Department of Physics, 8093 Zurich, Switzerland  Ben Haylock
Centre for Quantum Computation and Communication Technology (Australian
Research Council), Centre for Quantum Dynamics, Griffith University, Brisbane,
QLD 4111, Australia  Institute for Photonics and Quantum Sciences, SUPA,
Heriot-Watt University, Edinburgh EH14 4AS, United Kingdom  Francesco Lenzini
Centre for Quantum Computation and Communication Technology (Australian
Research Council), Centre for Quantum Dynamics, Griffith University, Brisbane,
QLD 4111, Australia  Institute of Physics, University of Muenster, 48149
Muenster, Germany  Mirko Lobino  Centre for Quantum Computation and
Communication Technology (Australian Research Council), Centre for Quantum
Dynamics, Griffith University, Brisbane, QLD 4111, Australia  Department of
Industrial Engineering, University of Trento, via Sommarive 9, 38123 Povo,
Trento, Italy  INFN–TIFPA, Via Sommarive 14, I-38123 Povo, Trento, Italy
Alberto Peruzzo  Quantum Photonics Laboratory and Centre for Quantum
Computation and Communication Technology, RMIT University, Melbourne, VIC
3000, Australia

###### Abstract

Understanding and controlling engineered quantum systems is key to developing
practical quantum technology. However, given the current technological
limitations, such as fabrication imperfections and environmental noise, this
is not always possible. To address these issues, a great deal of theoretical
and numerical methods for quantum system identification and control have been
developed. These methods range from traditional curve fittings, which are
limited by the accuracy of the model that describes the system, to machine
learning methods, which provide efficient control solutions but no control
beyond the output of the model, nor insights into the underlying physical
process. Here we experimentally demonstrate a “graybox” approach for
characterizing and controlling quantum devices. We report superior performance
over model fitting, while generating unitaries and Hamiltonians which are
inaccessible with machine learning techniques. Our approach combines physics
principles with high accuracy machine learning and is effective with any
problem where the required controlled quantities are not directly accessible
from the training set. This method naturally extends to time-dependent and
open quantum systems, with applications in quantum noise spectroscopy and
cancellation.

Quantum technology promises to deliver exponentially faster computation,
provably secure communications, and high-precision sensing Leuchs and Bruss
(2019). However, during the fabrication and operation of a quantum device,
there are many factors that can significantly impact its fuctionality,
requiring characterization and control techniques to achieve high-level
performance. Generally, we are interested in uncovering the unknown relation
between the control and the Hamiltonian governing the device, and then utilize
this information to drive the device towards a desired target. Typical targets
include unitary gates, a specific Hamiltonian, or certain output probability
distributions.

The traditional approach to characterizing and controlling physical devices is
based on theoretical models of the underlying processes governing the
relationship between input and output signals. These “whitebox” (WB) models
are based on parameter estimation via curve fitting and can be computationally
expensive, inaccurate, or incomplete. For example, they do not consider
unexpected input parameters or dynamics. Moreover, the models commonly used to
describe the dynamic of an open quantum system (such as Lindblad’s master
equation), is valid only under strict assumptions and approximations of the
noise (such as Markovianity) and the control (such as being ideal impulses).
These assumptions do not hold for many quantum platforms, making the use of WB
models inaccurate.

To increase accuracy and remove some of the limitations of the WB method,
supervised machine learning techniques are proving useful for characterizing
complex physical systems. In particular, techniques such as neural networks
are in many cases superior if not the only viable approach. However, this
approach, referred to as “blackbox” (BB), does not provide any information
about the underlying physics of the system, and cannot be used to study any
quantity beyond the output of the model, posing a major limitation in both
classical and quantum control. Nonetheless, it has been used in many
applications such as quantifying Non-Markovianity of quantum systems Fanchini
_et al._ (2021), characterizing qubits and environments Flurin  _et al._
(2020); Papič and de Vega (2022); Wise  _et al._ (2021); Palmieri  _et al._
(2021), quantum control Ostaszewski  _et al._ (2019); Khait  _et al._ (2022);
Zeng  _et al._ (2020), quantum error correction Bausch and Leditzky (2020);
Baireuther  _et al._ (2018), optimization of experimental quantum measurements
Lennon  _et al._ (2019); Tranter  _et al._ (2018), and calibration of quantum
devices Cimini  _et al._ (2021). Reinforcement learning methods Sivak  _et
al._ (2022); Baum  _et al._ (2021); Niu  _et al._ (2019) have also been
explored for the purposes of identification and control of quantum systems.
They also belong to the category of BB methods.

In situations where the quantities of interest, such as Hamiltonians,
unitaries, and noise operators, cannot be directly accessed experimentally,
the use of a hybrid WB-BB or “graybox” (GB) approach allows access to these
quantities beyond the measurable dataset. The GB approach was proposed for the
characterization and control of quantum devices Youssry  _et al._ (2020a, b,
2021); Youssry and Nudrin (2022), geometric quantum gate synthesis Perrier
_et al._ (2020), and was used to characterize superconducting qubits Genois
_et al._ (2021). However, to date no experimental demonstration of quantum
control based on the GB method has been reported.

A closely-related approach to system identification is the “Hamiltonian
Learning” problem Wiebe  _et al._ (2014); Wang  _et al._ (2017), where a fixed
time-independent Hamiltonian (control is fixed) is inferred from quantum
measurements. Alternatively, there are other approaches that directly aim to
control the quantum device without first identifying it, such as direct
gradient-based optimization. This includes the commonly used GRAPE algorithm
Khaneja  _et al._ (2005) and its variants de Fouquieres  _et al._ (2011);
Ciaramella  _et al._ (2015); Abdelhafez  _et al._ (2019); Leung  _et al._
(2017); Caneva  _et al._ (2011); Haas  _et al._ (2019). These model-free
techniques only work when the dependence of the Hamiltonian on the control is
known, because they are based on optimizing the fidelity to some target with
respect to control.

Here, we experimentally demonstrate how to model and control a quantum device
using the GB architecture when the Hamiltonian dependence on the control is
unknown. We report high-fidelity preparation of arbitrary unitaries and output
probability distributions of a reconfigurable three-mode integrated photonic
device, and uncover the Hamiltonian dependence on the control. Our GB approach
outperforms the traditional model fitting methods and can successfully prepare
unitary operations that are not accessible via BB methods. Our results show a
promising approach to enhance quantum control by understanding the physical
processes, and open the way to improve the engineering of quantum devices.

Figure 1: Physical and machine learning models of the class of quantum devices
considered in this paper. These are described by a time-independent
Hamiltonian in the absence of interaction with the environment. a)
Representation of the class of quantum devices considered in this work. b) A
schematic of an integrated photonic voltage-controlled reconfigurable
waveguide array chip, implementing a noiseless time-independent Hamiltonian.
Photons enter from the input port (on the left), undergo a voltage controlled
propagation along the chip, and are then measured at the output port of the
chip (on the right). c) The structure of the proposed graybox model. The input
to the model is the set of $M$ controls, while the outputs are the quantum
measurements for the set of computational basis as initial states. $P_{a\to
b}$ indicates the transition probability from input port $a$ to output port
$b$. The graybox is a combination of black and white boxes. The blackbox
estimates the real and imaginary components of each matrix element of the
Hamiltonian. The whitebox layers construct the Hamiltonian matrix, and perform
the quantum evolution and measurements. d) A fully whitebox architecture where
a physical model is utilized. The first layer generates predefined Hamiltonian
parameters that follow a known analytical dependence. The remaining layers
perform the quantum evolution and measurements. e) A fully blackbox model
where only a generic neural network is utilized with no physical model.

###  .1 Modeling quantum devices

We consider the class of quantum devices, shown in Fig. 1a, described by a
time-independent Hamiltonian undergoing a closed-system evolution (i.e. in the
absence of quantum noise). Photonic devices are examples of this class, when
the Hamiltonian is not modulated faster than the evolution time of a
propagating photon. We focus on the case where the system is described by a
finite $N$-dimensional Hilbert space (i.e. a qudit). The Hamiltonian governing
the dynamic of the system can be represented in the most general form as an
$N\times N$ complex Hermitian matrix that depends on a set of external
controls. We encode the set of controls in a $M\times 1$-dimensional vector
$\mathbf{V}=[V_{1},V_{2},\cdots V_{M}]^{T}$, where $V_{k}$ is the
$k^{\text{th}}$ control. We assume that during the system evolution the
control vector is fixed. An example of such controls are the voltages applied
to a reconfigurable photonic chip, schematically shown in Fig. 1b. The system
starts in an initial state $\ket{\psi_{0}}$ and evolves to the state
$\ket{\psi_{T}}$ at time $t=T$. The state is then measured in some basis to
obtain a set of probability outcomes. A more general form of time-dependent
evolution in the the presence of unwanted interactions with the environment
have been considered in our previous work Youssry  _et al._ (2020b).

Our aim is to obtain a machine learning (ML) model that describes the
behaviour of the device given a set of controls $\mathbf{V}$, and use it for
controlling the device. The input to the ML model is the $M$-dimensional
control vector $\mathbf{V}$, while the output is the set of measured outcomes
of the state after evolution. Generally, it is required to have
informationally-complete measurements to fully characterize a quantum system.
Here, we restrict the initial states as well as measurement basis to the set
of computational basis states $\\{\ket{0},\ket{1},\cdots\ket{N}\\}$ in order
to be compatible with our experimental setup. For each of these $N$ initial
states we have $N$ possible outcomes with an associated probability
$P_{j\rightarrow k}$ corresponding to the $j^{th}$ input and $k^{th}$ output,
giving a total of $N^{2}$ outputs. The approach, however, is independent of
this choice, and any set of states can be used. In the Supplemental Material,
we provide the results of our method applied to a simulated synthetic dataset
with a more general measurement scheme. It is important to emphasize that the
model input is the set of controls $\mathbf{V}$ applied to the system, and not
the initial state $\ket{\psi_{0}}$.

###  .2 Graybox architecture

The GB structure we propose, shown in Fig. 1c, consists of a BB (in the form
of neural network) followed by a WB part that process the outputs of the BB
into measurable physical quantities. The purpose of the blackbox part is to
map the controls $\mathbf{V}$ (the model inputs) to the Hamiltonian of the
system. The output of the BB then represents the elements of the Hamiltonian
matrix. A general $N\times N$ complex matrix has $2N^{2}$ degree of freedom
($N^{2}$ components with real and imaginary parts). Thus, the output layer of
the BB must consist of $2N^{2}$ neurons. The other BB layers can be designed
arbitrarily, and are custom engineered to provide the best performance for a
given dataset.

The second layer is a Hamiltonian construction layer that arranges the outputs
of the BB into an $N\times N$ matrix. A valid Hamiltonian has to be Hermitian
(i.e. $H^{\dagger}=H$) and this is ensured by calculating the Hermitian part
of the constructed matrix, and discarding the anti-Hermitian part. This can be
done simply by adding the matrix to its Hermitian conjugate. The output of
this layer is a valid control-dependent Hamiltonian. We did not enforce any
structure or constraint on the Hamiltonian, to enable the best fitting allowed
by the rules of quantum mechanics for a given dataset.

Because the Hamiltonian is not directly accessible through quantum
measurements, we cannot build a model with the Hamiltonian as an output. This
is why we need to add subsequent layers to transform the Hamiltonian matrix
into the set of probability outcomes, utilizing the laws of quantum mechanics.
Starting from a valid Hamiltonian, there is no need to further use a BB since
the dynamical equations are known. This saves the algorithm from trying to
learn the rules of quantum mechanics from experimental data, which would
complicate the process and result in a less accurate model. We use WB layers
for the remaining steps. In particular, there is a layer that calculates the
evolution unitary by the matrix exponentiation of the Hamiltonian
$U=e^{-iHT}$, which is the solution of the Schrödinger’s time-independent
equation. The final part of the GB model is a concatenation of $N$-layers
representing the quantum measurement operation for each of the $N$ input
$\ket{j}$. In each layer we calculate $\ket{\psi_{T}}=U\ket{j}$, where
$\ket{j}$ is the initial state of the system. After evolving the state to
$\ket{\psi_{T}}$, the probabilities $P_{j\rightarrow k}$ are calculated by
taking the absolute value squared of each entry of the state, that is applying
the Born’s rule for quantum expectation values.

The added WB layers do not include any trainable parameter; they only exist in
the BB part. Therefore, when we train the model on a dataset, the only updates
occur in the BB, generating a set of outputs that can be interpreted as a
Hamiltonian.

An important aspect of the GB architecture is that it is independent on the
physics of the system, i.e. it provides the most general form of a map between
the control vector and the quantum measurements while keeping the most
important physical quantities accessible via software, namely Hamiltonian,
unitary and evolved state. This is the key aspect needed to perform quantum
control as we are usually interested in implementing a quantum gate
represented by the unitary or Hamiltonian and not by the measured evolved
state. Having access to those quantities, even though the model is trained
with quantum measurements, is the key feature of the GB architecture.

Figure 2: Protocol schematic. The first step is to construct an experimental
dataset by applying controls to the system and measuring the corresponding
outputs. The dataset is then used to train the machine learning models (1).
Next, the trained models are tested (2) for generalization by comparing their
output predictions against a different experimental testing dataset. After
that, the trained models can be used to optimize controls (3) to achieve a
certain target, which could be a Hamiltonian, a unitary gate or an output
probability distribution. Finally, the obtained controls are tested (4)
experimentally and the controlled system output is compared against the
desired target.

###  .3 Whitebox and blackbox architectures

We benchmark the performance of our GB model against the fully WB and fully BB
architectures. The WB approach (Fig. 1d) is equivalent to standard model
(curve) fitting and the details of the architecture depend on the physical
system. The assumption is that all the relations between different dynamical
variables are exactly known except for the parameters we are fitting.
Generally, a WB consists of several layers. The first layer represents the
mathematical relations between the controls and the Hamiltonian, with a set of
unknown parameters. The input of this layer is the control vector, and the
output is a mathematically valid Hamiltonian. The remaining layers are
identical to those of the GB, and represent quantum evolution and quantum
measurements. For some systems we might need more layers (e.g. to model the
fan-in/out in photonic devices Peruzzo  _et al._ (2010)). The WB provides the
same access to hidden quantities as the GB, and even provides more physics as
we know exactly the analytical relations between Hamiltonian and control.
However, if we do not know these relations, or they do not match the physical
reality, the WB will fail and will not be useful for further applications.

The BB architecture (Fig. 1e) is largely different from the WB model since the
relation between the control vector and the quantum measurement is modeled by
a neural network. While any structure can be used, in this paper we consider
fully-connected neural networks with a softmax output layer of $N$ neurons.
This enforces the outputs to form a probability distribution (i.e. positive
numbers in $[0,1]$ whose sum is equal to $1$). This is consistent with what
the model outputs represent, which is probability amplitudes of the evolved
quantum state. Since we are characterising with $N$-initial basis states, we
need $N$ of such layers with the initial state chosen accordingly. This will
give a total of $N^{2}$ outputs. The structure of the other hidden layers can
only be determined and optimized by examining the performance on an actual
dataset. No other physical quantities can be access through this architecture,
but it is very good in fitting a dataset since it gives the maximum freedom in
terms of representation. If we are only interested in controlling the outputs,
a BB would be an efficient solution. But, if the goal is to estimate physical
quantities (like unitary gates), then a BB model is not useful at all.

###  .4 Protocol for training, testing and controlling

Our protocol for training and testing models and controllers is schematically
depicted in Fig. 2. It starts with preparing a dataset that will be used to
train and test the machine learning BB. The dataset consists of examples. Each
example is made of the $M$ control inputs $\mathbf{V}$ and the $N^{2}$ outputs
of the model $P_{j\rightarrow k}$. So we start by generating random values for
our control, let the system evolve, then perform the measurements and obtain
the probabilities $P_{j\rightarrow k}$. We repeat this for the $N$ input
states we consider to obtain all the outputs. The procedure is repeated for
multiple examples. The number of examples of the dataset depends on the
particular structure of the ML model, the noise level in the experiment, and
the acceptable performance level. Generally, the larger dataset is, the better
the ML algorithm will perform.

After the dataset is collected, it is split into the training and testing
subsets, and the ML model is trained. The purpose of training is to minimize a
cost function given by the standard mean square error (MSE) between the
predictions and actual model outputs from the training dataset. Once the model
is trained, its parameters are fixed and do not change in any of the remaining
protocol stages. Next, the model is evaluated using the testing dataset, and
its predictions are compared to the true corresponding values. The testing
examples are not included in the optimization procedure, so they provide an
unbiased evaluation on the performance of the model.

Let’s now consider the trained model for control purposes. In this case, the
model acts as a replacement of the actual experimental setup, and can be
probed via software for any purpose. We consider two types of applications.
The first is for obtaining the values of the control $\mathbf{V}$ to achieve a
target output (i.e. probability amplitudes). In this case we use the MSE as
the loss function, and the optimal control vector $\mathbf{V}^{*}$ can be
expressed as

$\displaystyle\mathbf{V}^{*}=\operatorname*{arg\,min}_{\mathbf{V}\in\mathcal{I}}(\hat{\mathbf{y}}_{P}(\mathbf{V})-\mathbf{\mathbf{y}}_{d})^{T}(\hat{\mathbf{y}}_{P}(\mathbf{V})-\mathbf{y}_{d})$
(1)

where $\hat{\mathbf{y}}_{P}(\cdot)$ is the ML output predictions,
$\mathbf{y}_{d}$ is the desired target, and $\mathcal{I}$ is the control
domain, which reflects the maximum allowed range for the controls. To get
accurate results, the ML model should be trained with dataset examples that
lie in the same control domain as well. Note that the internal ML model
parameters that define $\hat{\mathbf{y}}_{P}(\cdot)$ are not allowed to change
during the optimization as they have been fixed after the training.

The other case is for achieving a target quantum gate (i.e. a target unitary).
For this application, we use the gate fidelity as the loss function for the
controller defined as

$\displaystyle
F(U,W)=\frac{\left|\operatorname{tr}{(U^{\dagger}W)}\right|^{2}}{N^{2}},$ (2)

where $U$ and $W$ are two unitary matrices, and $N$ is their dimension. The
gate fidelity lies in the range $[0,1]$, with $1$ representing maximum overlap
between the two gates. The optimal control voltages can then be represented as

$\displaystyle\mathbf{V}^{*}=\operatorname*{arg\,max}_{\mathbf{V}\in\mathcal{I}}F(\hat{y}_{U}(\mathbf{V}),U_{d})$
(3)

where $\hat{y}_{U}(\cdot)$ is the evolution unitary obtained from the ML
model, and $U_{d}$ is the desired target gate. In this case we can only use
the GB and WB models, since a BB does not provide access to unitaries as
discussed earlier. With this modular approach, the optimization algorithm of
the controller can be chosen arbitrarily. While we choose a gradient-based
method in this paper, other techniques such as genetic algorithms could also
be used.

Finally, the optimal controls for a set of targets are applied experimentally,
and the system is measured to construct the “control” dataset. This dataset is
then assessed and compared against the desired targets. Our main goal is to
control a quantum system, and thus the assessment of any model should not just
rely on its prediction capabilities, but also on how it performs in
conjunction with a controller when tested experimentally.

Figure 3: Experimental performance of the machine learning models. (a) Results
of training the different models on the experimental dataset. The MSE is
plotted versus iteration number (b) The results of evaluating the different
models on the testing set.

###  .5 Experimental results and discussion

We tailor the design of the models described above around the device used for
the experimental verification for our proposal, a voltage controlled quantum
photonic circuit of continuously coupled waveguides based on lithium niobate
technology, schematically shown in Fig. 1b. The details about the chip
physical model as well as fabrication are given in the Supplemental Material.
The chip has 3 waveguides, corresponding to a qutrit system, and is controlled
by 4 electrodes and their respective voltages. We implemented the ML models
using the Tensorflow Python package Abadi  _et al._ (2015); Chollet  _et al._
(2015), and applied the protocol for training the three models, and then
verified the performance of the controllers. The details of the
implementations are also given in the Supplemental Material.

| GB | WB | BB
---|---|---|---
Model |  |  |
average training MSE | $8.3\times 10^{-5}$ | $1.5\times 10^{-2}$ | $9.4\times 10^{-5}$
average testing MSE | $9.1\times 10^{-5}$ | $1.5\times 10^{-2}$ | $1.1\times 10^{-4}$
Output Controller |  |  |
average MSE | $2.6\times 10^{-3}$ | $1.7\times 10^{-2}$ | $2.9\times 10^{-3}$
fidelity |  |  |
average | $99.53\%$ | $97.47\%$ | $99.48\%$
with $>99\%$ | $87.3\%$ | $34\%$ | $86.2\%$
Unitary Controller |  |  |
average MSE | $3.1\times 10^{-3}$ | $1.9\times 10^{-2}$ | –
fidelity |  |  |
average | $99.48\%$ | $97.4\%$ | –
with $>99\%$ | $71.3\%$ | $12\%$ | –
Table 1: Results summary. Comparison of the performance of each model on the
training and testing and control experimental datasets, as well as the
experimental fidelities (average and instances greater than 99%) for 1000
randomly prepared output distributions and unitary gates.

The results of the models training and testing as well as the control
performance are reported in Table 1. The MSE evaluated at each iteration for
training and testing sets is shown in Fig. 3a and Fig. 3b. The plot of the
learning curve in Fig. 3a shows the superior performance of the GB in terms of
accuracy compared to the WB. This is due to the constraints imposed by the
physical model that is used to construct the Hamiltonian in the case of the
WB. Evanescently coupled waveguide arrays, such as the one tested here, are
typically assumed to have only nearest-neighbor coupling. Our results show
that this assumption does not hold for the real device. Furthermore, the
typical assumption of linear dependence on the voltage, does not hold
perfectly either. On the other hand, the GB learns a general mathematically
valid Hamiltonian, and thus is able to better fit the experimental data. The
BB had more similar performance to GB compared to WB, with the main drawback
of completely losing the physical picture. In terms of the testing
performance, Fig. 3b shows that the three models do not overfit, as the final
MSE of the testing set is close to the final MSE of the training set. This
means that the models do not memorize the examples of the training set. In
other words, there’s no significant loss in the prediction accuracy. However,
the GB and BB clearly perform better than the WB.

Figure 4: Experimental quantum control performance. The distribution of the
fidelity between the experimentally measured output power distribution and the
desired target distribution for the three models. The results are for (a) the
output controller, and (b) the unitary controller. The reported values are the
average over the three distributions corresponding to each possible initial
state. (c) Violin plot showing the statistics of the MSE obtained for the
training, testing, and control datasets. The horizontal lines represent from
bottom to top, the minimum, median, and maximum respectively. The plot also
shows an estimated kernel density for the data.

The performance of the models is limited by the size of the dataset. Because
the experimental data always suffer from some level of noise, the minimum MSE
obtainable without overfitting is limited as well. Usually the acceptable
level of MSE depends on the specific application. In this paper, our
application is to control the chip to obtain target power distributions as
well as target unitary operations. The ML models then act as a
replacement/simulator of the actual setup. The experimental assessment of the
optimal control will determine whether the model performance is accepted or
need improvement. In general, the way to improve models is by constructing
very large datasets which is the standard approach in most typical machine
learning applications. For engineering applications, where we characterize and
control a physical device, we are limited by how many measurements we can
obtain. Thus, it becomes a tradeoff between the amount of time and resources
needed to construct the dataset experimentally, and the accuracy of the
trained models, which will also affect the performance of the controller.

The histogram of the controller fidelities between the desired target and the
experimental measurements for 1000 randomly prepared output power
distributions and 1000 randomly prepared unitary gates are shown in Fig. 4a
and Fig. 4b respectively. The plots show that the WB is particularly skewed
towards lower fidelities (minimum is $80.53\%$ compared to a minimum of
approximately $91\%$ for GB and BB for the output controller). Similarly for
the gate controller, the minimum fidelity is $86.6\%$ and $94.95\%$ for WB and
GB. In Fig. 4c, we summarize the statistics of the MSE between the ML model
predictions and actual outputs for the training and testing datasets, in
comparison with the MSE between the experimentally controlled measurements and
the targets for each of the two controllers.

The results of the controller for obtaining a target power distribution, show
once again the superior performance of the GB and BB over the WB in terms of
both MSE and average fidelity. The same controller/optimization algorithm and
cost function are used for the three models. Thus, the lower performance comes
from the lower accuracy of the WB model itself. When considering the
controller for target unitary, it is only possible to use a WB or a GB as they
are the only models that can give access to the overall unitary evolution
matrix. A BB cannot be utilized in this case since it only encodes the
dynamics in an abstract machine-suitable format, and does not provide any
physical picture.

Figure 5: Dependence of the Hamiltonian elements to a subset of input
voltages, as predicted by the graybox model. (a) Real and (b) imaginary parts
of the Hamiltonian matrix elements as function of voltage when all electrodes
are grounded except the first electrode. It should be noted that the imaginary
parts of $H_{11}$, $H_{22}$ and $H_{33}$ are by definition equal to zero.

Finally, the GB provides sufficient physical insights for most purposes, for
example in Fig. 5a and 5b, where the tunability of the Hamiltonian as function
of a single electrode voltage is explored. The figure shows the GB prediction
of the different Hamiltonian elements as a function of the voltage applied to
a single electrode. We can use these predictions more generally when more than
one electrode is tuned (although it would be more difficult to plot in this
situation) for unitary control. We can use the GB to predict the unitary given
the set of voltages even though the dataset originally did not include this
information, but rather the power distribution. It is also important to
realize that this predicted Hamiltonian, besides not being unique
mathematically, represents physically an effective quantity, and so it will
differ in structure (such as the existence of an imaginary part) from the
ideal Hamiltonian that one may expect for a system. Depending on the purpose
of the use of the GB we can control how much we make it “blacker” or “whiter”.
For control applications, the best architecture is to have this effective
Hamiltonian. In another application such as modeling a device for the purpose
of completely understanding the physics, the architecture of the GB might need
to be modified to allow Hamiltonians that are closer to some expected
structure.

###  .6 Conclusion and outlook

We have shown how a GB model can be designed for a general quantum device,
trained on experimental data and verified by generating target unitary
operations and output distributions with high fidelity. The performance was
benchmarked against WB and BB models, showing the superior performance of our
approach. The standard WB model is an order of magnitude less accurate than
the machine leaning methods, and the BB model, as it only provides access to
the model output, cannot generate target unitaries or Hamiltonians. Our
approach is general and can be applied to any quantum system, it can be
extended to time-dependent and open quantum systems Youssry  _et al._ (2020b,
2021); Youssry and Nudrin (2022), for noise spectroscopy and noise
cancellation. There are many possible extensions to this work as well. One
could consider a hybrid approach between Hamiltonian learning using Genetic
Algorithms (such as Schmidt and Lipson (2009); Gentile  _et al._ (2021)) and
our numerical GB for the purpose of obtaining more detailed physical models.
Another aspect related to any ML algorithm is the requirements of the training
dataset size. For complex devices, it can be challenging to collect a large-
sized dataset. However, some emerging techniques can facilitate this process
including incremental learning Losing  _et al._ (2018), transfer learning
Weiss  _et al._ (2016), and adaptive online learning Hoi  _et al._ (2021).
While these techniques are constantly developing in classical machine learning
literature, there is still a gap in porting such methods to physics-based
applications and especially quantum applications.

##  I Acknowledgements

AP acknowledges a RMIT University Vice-Chancellor’s Senior Research Fellowship
and a Google Faculty Research Award. ML was supported by the Australian
Research Council (ARC) Future Fellowship (FT180100055). BH was supported by
the Griffith University Postdoctoral Fellowship. This work was supported by
the Australian Government through the Australian Research Council under the
Centre of Excellence scheme (No: CE170100012), and the Griffith University
Research Infrastructure Program. This work was performed in part at the
Queensland node of the Australian National Fabrication Facility, a company
established under the National Collaborative Research Infrastructure Strategy
to provide nano- and micro-fabrication facilities for Australia’s researchers.

## References

  * Leuchs and Bruss (2019) G. Leuchs and D. Bruss,  _Quantum information: from foundations to quantum technology applications_ (John Wiley & Sons, 2019).
  * Fanchini  _et al._ (2021) F. F. Fanchini, G. Karpat, D. Z. Rossatto, A. Norambuena, and R. Coto, Physical Review A 103, 022425 (2021).
  * Flurin  _et al._ (2020) E. Flurin, L. S. Martin, S. Hacohen-Gourgy, and I. Siddiqi, Physical Review X 10, 011006 (2020).
  * Papič and de Vega (2022) M. Papič and I. de Vega, Physical Review A 105, 022605 (2022).
  * Wise  _et al._ (2021) D. F. Wise, J. J. Morton, and S. Dhomkar, PRX Quantum 2, 010316 (2021).
  * Palmieri  _et al._ (2021) A. M. Palmieri, F. Bianchi, M. G. Paris, and C. Benedetti, Physical Review A 104, 052412 (2021).
  * Ostaszewski  _et al._ (2019) M. Ostaszewski, J. Miszczak, L. Banchi, and P. Sadowski, Quantum Information Processing 18, 1 (2019).
  * Khait  _et al._ (2022) I. Khait, J. Carrasquilla, and D. Segal, Physical Review Research 4, L012029 (2022).
  * Zeng  _et al._ (2020) Y. Zeng, J. Shen, S. Hou, T. Gebremariam, and C. Li, Physics Letters A 384, 126886 (2020).
  * Bausch and Leditzky (2020) J. Bausch and F. Leditzky, New Journal of Physics 22, 023005 (2020).
  * Baireuther  _et al._ (2018) P. Baireuther, T. E. O’Brien, B. Tarasinski, and C. W. Beenakker, Quantum 2, 48 (2018).
  * Lennon  _et al._ (2019) D. T. Lennon, H. Moon, L. C. Camenzind, L. Yu, D. M. Zumbühl, G. A. D. Briggs, M. A. Osborne, E. A. Laird, and N. Ares, npj Quantum Information 5, 1 (2019).
  * Tranter  _et al._ (2018) A. D. Tranter, H. J. Slatyer, M. R. Hush, A. C. Leung, J. L. Everett, K. V. Paul, P. Vernaz-Gris, P. K. Lam, B. C. Buchler, and G. T. Campbell, Nature Communications 9, 4360 (2018).
  * Cimini  _et al._ (2021) V. Cimini, E. Polino, M. Valeri, I. Gianani, N. Spagnolo, G. Corrielli, A. Crespi, R. Osellame, M. Barbieri, and F. Sciarrino, Physical Review Applied 15, 044003 (2021), 2009.07122 .
  * Sivak  _et al._ (2022) V. Sivak, A. Eickbusch, H. Liu, B. Royer, I. Tsioutsios, and M. Devoret, Physical Review X 12, 011059 (2022).
  * Baum  _et al._ (2021) Y. Baum, M. Amico, S. Howell, M. Hush, M. Liuzzi, P. Mundada, T. Merkh, A. R. Carvalho, and M. J. Biercuk, PRX Quantum 2, 040324 (2021).
  * Niu  _et al._ (2019) M. Y. Niu, S. Boixo, V. N. Smelyanskiy, and H. Neven, npj Quantum Inf. 5 (2019), 10.1038/s41534-019-0141-3.
  * Youssry  _et al._ (2020a) A. Youssry, R. J. Chapman, A. Peruzzo, C. Ferrie, and M. Tomamichel, Quantum Science and Technology 5, 025001 (2020a).
  * Youssry  _et al._ (2020b) A. Youssry, G. A. Paz-Silva, and C. Ferrie, npj Quantum Information 6, 1 (2020b).
  * Youssry  _et al._ (2021) A. Youssry, G. A. Paz-Silva, and C. Ferrie, arXiv preprint arXiv:2103.13018 (2021).
  * Youssry and Nudrin (2022) A. Youssry and H. I. Nudrin, arXiv preprint arXiv:2208.03058 (2022).
  * Perrier  _et al._ (2020) E. Perrier, D. Tao, and C. Ferrie, New Journal of Physics 22, 103056 (2020).
  * Genois  _et al._ (2021) É. Genois, J. A. Gross, A. D. Paolo, N. J. Stevenson, G. Koolstra, A. Hashim, I. Siddiqi, and A. Blais, PRX Quantum 2, 040355 (2021), 2106.13126 .
  * Wiebe  _et al._ (2014) N. Wiebe, C. Granade, C. Ferrie, and D. Cory, Phys. Rev. A 89, 042314 (2014).
  * Wang  _et al._ (2017) J. Wang _et al._ , Nat. Phys. 13, 551 (2017).
  * Khaneja  _et al._ (2005) N. Khaneja, T. Reiss, C. Kehlet, T. Schulte-Herbrüggen, and S. J. Glaser, J. Magn. Reson. 172, 296 (2005).
  * de Fouquieres  _et al._ (2011) P. de Fouquieres, S. Schirmer, S. Glaser, and I. Kuprov, J. Magn. Reson. 212, 412 (2011).
  * Ciaramella  _et al._ (2015) G. Ciaramella, A. Borzì, G. Dirr, and D. Wachsmuth, SICOMP 37, A319 (2015).
  * Abdelhafez  _et al._ (2019) M. Abdelhafez, D. I. Schuster, and J. Koch, Phys. Rev. A 99, 052327 (2019).
  * Leung  _et al._ (2017) N. Leung, M. Abdelhafez, J. Koch, and D. Schuster, Phys. Rev. A 95, 042318 (2017).
  * Caneva  _et al._ (2011) T. Caneva, T. Calarco, and S. Montangero, Phys. Rev. A 84, 022326 (2011).
  * Haas  _et al._ (2019) H. Haas, D. Puzzuoli, F. Zhang, and D. G. Cory, New J. Phys. 21, 103011 (2019).
  * Peruzzo  _et al._ (2010) A. Peruzzo, M. Lobino, J. C. F. Matthews, N. Matsuda, A. Politi, K. Poulios, X.-Q. Zhou, Y. Lahini, N. Ismail, K. Wörhoff, Y. Bromberg, Y. Silberberg, M. G. Thompson, and J. L. OBrien, Science 329, 1500 (2010), 1006.4764 .
  * Abadi  _et al._ (2015) M. Abadi, A. Agarwal, P. Barham, E. Brevdo, Z. Chen, C. Citro, G. S. Corrado, A. Davis, J. Dean, M. Devin, S. Ghemawat, I. Goodfellow, A. Harp, G. Irving, M. Isard, Y. Jia, R. Jozefowicz, L. Kaiser, M. Kudlur, J. Levenberg, D. Mané, R. Monga, S. Moore, D. Murray, C. Olah, M. Schuster, J. Shlens, B. Steiner, I. Sutskever, K. Talwar, P. Tucker, V. Vanhoucke, V. Vasudevan, F. Viégas, O. Vinyals, P. Warden, M. Wattenberg, M. Wicke, Y. Yu, and X. Zheng, “TensorFlow: Large-scale machine learning on heterogeneous systems,”  (2015), software available from tensorflow.org.
  * Chollet  _et al._ (2015) F. Chollet _et al._ , “Keras,” https://keras.io (2015).
  * Schmidt and Lipson (2009) M. Schmidt and H. Lipson, Science 324, 81 (2009).
  * Gentile  _et al._ (2021) A. A. Gentile, B. Flynn, S. Knauer, N. Wiebe, S. Paesani, C. E. Granade, J. G. Rarity, R. Santagati, and A. Laing, Nature Physics 17, 837 (2021).
  * Losing  _et al._ (2018) V. Losing, B. Hammer, and H. Wersing, Neurocomputing 275, 1261 (2018).
  * Weiss  _et al._ (2016) K. Weiss, T. M. Khoshgoftaar, and D. Wang, Journal of Big data 3, 1 (2016).
  * Hoi  _et al._ (2021) S. C. Hoi, D. Sahoo, J. Lu, and P. Zhao, Neurocomputing 459, 249 (2021).