# Loophole-free Bell tests with randomly chosen subsets of measurement
settings

Jaskaran Singh<EMAIL_ADDRESS>Departamento de Física Aplicada II, Universidad
de Sevilla, E-41012 Sevilla, Spain  Adán Cabello<EMAIL_ADDRESS>Departamento de
Física Aplicada II, Universidad de Sevilla, E-41012 Sevilla, Spain  Instituto
Carlos I de Física Teórica y Computacional, Universidad de Sevilla, E-41012
Sevilla, Spain

###### Abstract

There are bipartite quantum nonlocal correlations requiring very low detection
efficiency to reach the loophole-free regime but that need too many
measurement settings to be practical for actual experiments. This leads to the
general problem of what can be concluded about loophole-free Bell nonlocality
if only a random subset of these settings is tested. Here we develop a method
to address this problem. We show that, in some cases, it is possible to detect
loophole-free Bell nonlocality testing only a small random fraction of the
settings. The prize to pay is a higher detection efficiency. The method allows
for a novel approach to the design of loophole-free Bell tests in which, given
the dimension of the local system, the visibility, and the detection
efficiency available, one can calculate the fraction of the contexts needed to
reach the detection-loophole-free regime. The results also enforce a different
way of thinking about the costs of classically simulating quantum nonlocality,
as it shows that the amount of resources that are needed can be made
arbitrarily large simply by considering more contexts.

##  I Introduction

###  I.1 Motivation

Quantum advantage in applications based on the violation of Bell inequalities
[1, 2, 3, 4, 5, 6, 7, 8, 9] requires assuring that the observed correlations
cannot be simulated with local realistic models [10, 11]. This demands a
detection efficiency above a certain threshold $\eta_{\rm crit}$ [12] that
depends on the nonlocal correlations. Recently [13, 14], several bipartite
quantum correlations between high-dimensional quantum systems have been
identified requiring $\eta_{\rm crit}$ smaller than the ones needed for qubits
[15, 16] and ququarts [17], which are the systems used in the detection-
loophole-free Bell tests performed so far [18, 19, 20, 21, 22, 23, 24, 25].
However, for all these cases, achieving a small detection efficiency requires
a large number of settings. For example, to achieve $\eta_{\rm crit}=0.510$ in
Ref. [13], each party has to measure $2^{8}$ settings. Or, for example, to
achieve $\eta_{\rm crit}=0.324$ in Ref. [14], each party has to measure
$2^{30}$ settings. These correlations need too many settings to be useful for
actual experiments.

This begs the question of whether it is possible to detect loophole-free Bell
nonlocality when the parties only choose, randomly, a fraction of these
settings.

###  I.2 Aim

In a bipartite Bell experiment, a context is one of the possible pairs of
local settings used to evaluate the Bell inequality. Here, our aim is to
investigate whether the parties can detect loophole-free Bell nonlocality with
a certain confidence level using only a randomly chosen strict subset of the
contexts. If they can, we want to know how large $\eta$ must be to reach the
loophole-free regime. Reciprocally, for a given detection efficiency, we want
to know what is the minimum fraction of contexts needed to detect loophole-
free Bell nonlocality. In particular, we want to elucidate whether this
approach is useful when applied to the correlations in [13, 14] and similar
examples that may appear in the future.

###  I.3 Structure

The paper is organized as follows. In Sec. II we develop a general framework
for deciding whether some given correlations can exhibit loophole-free Bell
violation by evaluating only a fraction of the contexts. There, we analyze the
case of a general bipartite Bell inequality and construct an estimator of the
value of the Bell parameter as a function of the fraction of randomly chosen
contexts. Then, by using Chebyshev’s inequality [26], we bound the minimum
fraction of contexts required to detect loophole-free Bell nonlocality for a
given confidence value.

In Secs. III and IV, we apply the framework to the quantum correlations
maximally violating the so-called penalized $N$-product (PNP) Clauser-Horne-
Shimony-Holt (CHSH) Bell inequalities [13] and to some graph-theoretic Bell
inequalities [14]. We focus on these correlations because they require very
low detection efficiency. For these cases, we will obtain how the detection
efficiency depends on the fraction of contexts.

Finally, in Sec. V, we summarize the pros and cons of the approach, explain
why it offers a new way to design loophole-free Bell tests, and discuss the
implications of our results for the attempts to simulate quantum nonlocality
with classical resources.

##  II Loophole-free Bell violation using a subset of the measurement settings

###  II.1 Method

Consider two spatially separated parties, Alice and Bob, each of them
performing measurements on a subsystem of a composite system. Let us denote
Alice’s measurement setting by $x$ and Bob’s measurement setting by $y$. Let
us denote the outcome of Alice’s measurement by $a$ and the outcome of Bob’s
measurement by $b$. A Bell inequality is a linear combination of joint
conditional probabilities $p(a,b|x,y)$ for the outcomes of different
combinations of measurement settings of Alice and Bob. Specifically, a Bell
inequality is an expression of the form

$\beta=\sum_{a,b,x,y}c_{a,b}^{x,y}p(a,b|x,y)=\sum_{j=1}^{M}\beta_{j}\leq C,$
(1)

where $j=\left(x,y\right)$ is the context corresponding to the settings $x$
and $y$, thus $\beta_{j}=\sum_{a,b}c_{a,b}^{x,y}p(a,b|j)$, $M$ is the total
number of contexts, and $C$ is the maximum achievable value by local realistic
models. For example, in the CHSH Bell inequality [27],
$c_{a,b}^{x,y}=\delta_{a\oplus b,xy}$, $M=4$, and $C=3$.

Alice and Bob randomly select a subset of contexts of cardinality $L$ and
evaluate the corresponding joint probability distributions. In principle, we
assume that $L<M$. However, as we will see, not all nonlocal correlations
allow us to certify nonlocality using only a subset of settings. Now, Alice
and Bob want to estimate the value of $\beta$. That is, they want to obtain
the most likely value of $\beta$ (whose detailed calculation requires
measuring all contexts) from the probabilities of contexts that have been
chosen.

For each of the chosen contexts, Alice and Bob evaluate the corresponding
$\beta_{j}$. We assume that, for each of the contexts chosen, Alice and Bob
conduct $K$ rounds of the experiment. This allows them to determine
$\beta_{j}$ up to some finite precision $\epsilon^{\prime}$ and probability of
failure $\delta^{\prime}$. Once all the terms $\beta_{j}$ are evaluated, Alice
and Bob can estimate $\beta$ with some finite precision $\epsilon$ and
probability of failure $\delta$.

To estimate $\beta$, Alice and Bob can proceed as follows. They select
$j\in\\{1,\ldots,M\\}$ at random with probability $p(j)=\frac{1}{M}$. Here,
for simplicity, we will assume a uniform probability distribution for $p(j)$.
However, $p(j)$ has to be tailored according to $\beta$ and may be not
uniform. For simplicity, we will also assume that it is possible to evaluate
$\beta_{j}$ with infinite precision (i.e., that $K\rightarrow\infty$). The
case of finite precision will be discussed later. Then, an estimator of
$\beta$ is

$X=M\beta_{j}.$ (2)

By construction, $\langle X\rangle=\beta$, where the mean value is averaged
over $j$.

Then, the parties choose $L$ contexts. If the first context is context $p$,
then we define $X_{1}=M\beta_{p}$, if the second context is context $q\neq p$,
then $X_{2}=M\beta_{q}$, etc. Each $X_{l}$ is an estimator of the value of
$\beta$. Let $Y$ be the average value of these estimators. That is,

$\displaystyle Y$ $\displaystyle=\frac{1}{L}\sum_{l=1}^{L}X_{l}$ (3)
$\displaystyle=\frac{M}{L}\sum_{i}\beta_{i},$

where the second sum also has $L$ terms. The variance of each $X_{l}$ is
Var$(X_{l})\leq M\beta$, which implies that Var$(Y)\leq\frac{M\beta}{L}$.
Therefore, using Chebyshev’s inequality [26],

$p\left(|Y-\beta|\geq\lambda\sqrt{\frac{M\beta}{L}}\right)\leq\frac{1}{\lambda^{2}},$
(4)

where $\lambda>0$ is a real number.

Since we want $Y$ with probability of failure $\delta$ and error $\epsilon$,
we take

$\lambda=\frac{1}{\sqrt{\delta}}$ (5)

and

$L=\frac{M\beta}{\epsilon^{2}\delta}$ (6)

to obtain

$p\left(|Y-\beta|\geq\epsilon\right)\leq\delta.$ (7)

Then, the fraction of contexts needed is

$\nu=\frac{L}{M}=\frac{\beta}{\epsilon^{2}\delta}.$ (8)

We now deal with the fact that $K$, the number of rounds used to evaluate each
of the terms $\beta_{j}$, must be finite. In this case, the average value of
an estimator of $\beta_{j}$ is

$B_{j}=\frac{1}{K}\sum_{k=1}^{K}\beta_{j}^{(k)},$ (9)

where $\beta_{j}^{(k)}$ is the value of $\beta_{j}$ obtained in round $k$. Let
$\delta^{\prime}$ be the probability of failure and $\epsilon^{\prime}$ the
error in evaluating $\beta_{j}$. Then, applying Hoeffding’s inequality [28],
we have

$p\left(B_{j}-\beta_{j}\geq\epsilon^{\prime}\right)\leq\delta^{\prime},$ (10)

where

$\displaystyle\delta^{\prime}$
$\displaystyle=\exp\left\\{-\frac{2\epsilon^{\prime
2}}{\sum_{k=1}^{K}\left[\sum_{a,b}c_{a,b}^{x,y}p(a,b|x,y)\right]^{2}}\right\\}$
(11) $\displaystyle=\exp\left\\{-\frac{2\epsilon^{\prime
2}}{K\left[\beta_{j}^{(k)}\right]^{2}}\right\\}.$

This yields

$K\geq\frac{-\ln(\delta^{\prime})\beta_{j}^{(k)^{2}}}{2\epsilon^{\prime 2}}.$
(12)

That is, the minimum number of rounds needed to estimate $\beta_{j}$.

###  II.2 The method does not always work

In many cases, all the contexts will be needed. For example, consider the CHSH
Bell inequality written so it has local value $C=3$ and maximum quantum
violation $2+\sqrt{2}\approx 3.414$. Let us suppose that the visibility is
$V=0.9$. This yields $\beta=3.272$ (thus $\epsilon=0.272$). Let us assume that
$\delta=0.00003$. Then, $L\geq 5,896,771$. However, the total number of
contexts in the CHSH Bell inequality is $M=4$. Therefore, $L>M$ indicates that
it is not possible to consider a strict subset of contexts and obtain a
loophole-free Bell violation.

##  III The PNP Bell inequalities using a subset of settings

Here, we apply the tools described in Sec. II to the case of the PNP Bell
inequalities of Ref. [13].

###  III.1 The PNP Bell inequalities

Given a Bell inequality with local realistic bound $C$, its PNP version is the
product of $n$ copies of that Bell inequality with an extra penalization term
chosen to guarantee that the local realistic bound of the PNP version is
$C^{n}$. That is, the PNP Bell inequality can be written as

$\beta_{PNP}=\sum_{\textbf{a},\textbf{b},\textbf{x},\textbf{y}}p(\textbf{a},\textbf{b}|\textbf{x},\textbf{y})\prod^{n}_{i=1}c_{a_{i},b_{i}}^{x_{i},y_{i}}-\kappa(\text{A}+\text{B})\leq
C^{n},$ (13)

where $\textbf{x}=(x_{1},\ldots,x_{n})$ is Alice’s measurement setting (which
can be seen as $n$ measurement settings, one for each copy of the original
Bell inequality), $\textbf{y}=(y_{1},\ldots,y_{n})$ is Bob’s measurement
setting, $\textbf{a}=(a_{1},\ldots,a_{n})$ is Alice’s outcome (which can be
seen as $n$ outcome, one for each copy of the original Bell inequality), and
$\textbf{b}=(b_{1},\ldots,b_{n})$ is Bob’s outcome.
$\kappa=2^{n-1}(\Sigma_{n}-C^{n})$, where $C$ and $\Sigma_{n}$ are the maximum
local bound and the algebraic bound of the original Bell inequality,
respectively. In addition,

A
$\displaystyle=\sum_{i=1}^{n}\sum_{\textbf{x}}\sum_{\textbf{x}\neq\textbf{x}^{\prime}}^{x_{i}=x^{\prime}_{i}}\sum_{a_{i}=0}|p(a_{i}|\textbf{x})-p(a_{i}|\textbf{x}^{\prime})|,$
(14a) B
$\displaystyle=\sum_{i=1}^{n}\sum_{\textbf{y}}\sum_{\textbf{y}\neq\textbf{y}^{\prime}}^{y_{i}=y^{\prime}_{i}}\sum_{b_{i}=0}|p(b_{i}|\textbf{y})-p(b_{i}|\textbf{y}^{\prime})|,$
(14b)

where the sum over $\textbf{x}^{\prime}$ is taken such that
$\textbf{x}^{\prime}\neq\textbf{x}$, except for the $i$th element. Similarly
for the sum over $\textbf{y}^{\prime}$. Let $Q$ be the maximum quantum value
of the original Bell inequality. Then, the maximum quantum value of the
corresponding PNP Bell inequality is simply $Q^{n}$.

Figure 1: Fraction of the settings required for a loophole-free test of the
PNP CHSH Bell inequalities, as a function of $n$, the number of CHSH Bell
inequalities conducted in parallel.

###  III.2 The PNP CHSH Bell inequalities with a subset of settings

We focus on the PNP Bell inequalities constructed by taking $n$ copies of the
CHSH Bell inequality. In this case,
$c_{a_{i},b_{i}}^{x_{i},y_{i}}=\delta_{a_{i}\oplus b_{i},x_{i}y_{i}}$ $\forall
a_{i},b_{i},x_{i},y_{i}$, $C=3$, and $Q=2+\sqrt{2}$. There are $2^{n}$
measurement settings per party and the total number of contexts is $M=4^{n}$.
For the time being, we will assume that $\text{A}+\text{B}=0$. Later on, we
will consider the case of $\text{A}+\text{B}\neq 0$.

The parties randomly select $L$ contexts. Each of them is tested $K$ times.
Applying the tools in Sec. II, we obtain the fraction of local settings that
each party must choose randomly to obtain a loophole-free Bell violation of
the PNP CHSH Bell inequality such that the probability of failure is
$\delta=0.00003$ and under the assumption that the visibility of each qubit-
qubit subspace is $V=0.9$. The results are shown in Fig. 1. There, we observe
that the fraction is smaller than $1$ only after $n=9$.

Next, we calculate the minimum detection efficiency, $\eta_{\nu}$, for a
loophole-free Bell violation as a function of the fraction of contexts chosen,
$\nu=\frac{L}{M}$. In order to do so, we model detection inefficiencies as
follows. We bin the no clicks to one of the outcomes (always the same one).
Then,

$\beta_{PNP}=\eta^{2}Q^{n}+\eta(1-\eta)(A^{n}+B^{n})+(1-\eta)^{2}C^{n},$ (15)

with

$\displaystyle Q$
$\displaystyle=\sum_{a,b,x,y}c_{a,b}^{x,y}\text{tr}\left(\mathcal{A}_{a}^{x}\otimes\mathcal{B}_{b}^{y}\rho_{AB}\right),$
(16a) $\displaystyle A$
$\displaystyle=\sum_{a,b,x,y}c_{a,b}^{x,y}\text{tr}\left(\mathcal{A}_{a}^{x}\rho_{A}\right),$
(16b) $\displaystyle B$
$\displaystyle=\sum_{a,b,x,y}c_{a,b}^{x,y}\text{tr}\left(\mathcal{B}_{b}^{y}\rho_{B}\right),$
(16c)

where the state shared by Alice and Bob is of the form
$\bigotimes_{i=1}^{n}\rho_{AB}$, and
$\bigotimes_{i=1}^{n}\mathcal{A}_{a_{i}}^{x_{i}}$ and
$\bigotimes_{i=1}^{n}\mathcal{B}_{b_{i}}^{y_{i}}$ are the elements of the
positive operator valued measures. We assume that
$\mathcal{A}_{a_{i}}^{x_{i}}=\mathcal{A}_{a}^{x}$ and
$\mathcal{B}_{b_{i}}^{y_{i}}=\mathcal{B}_{b}^{y}$, $~{}\forall i$.

Fig. 2 shows the minimum detection efficiency $\eta_{\nu}$ needed for a
loophole-free Bell violation, when using a fraction $\nu$ of all the contexts
in the PNP CHSH Bell inequality for $n=10,11,12,13$, and $14$, for
$\delta=0.00003$ and $V=1$. Notice that, in all cases, $\nu=1$ occurs before
the respective $\eta_{\rm crit}$ is reached ($\eta_{\rm crit}$ for
$n=10,11,12,13$, and $14$ are $0.43,0.38,0.34,0.31$, and $0.28$,
respectively). In addition, Table 1 presents the fraction of contexts needed
for a given detection efficiency $\eta$, for different values of $n$.

Figure 2: Minimum detection efficiency $\eta_{\nu}$ needed to reach the
loophole-free regime as a function of the fraction $\nu$ of the contexts for
the PNP CHSH Bell inequality with different values of $n$.

Finally, we analyze the case when $\text{A}+\text{B}\neq 0$. The value of this
sum cannot be very large, otherwise a violation will not be observed. This is
because $\kappa$ increases exponentially with $n$. As an example, consider
$\text{A}+\text{B}=10^{-6}$, for which the parties can observe a violation up
to $n=14$. Then, the maximum quantum violation for the PNP Bell inequality is
$\beta^{n}-10^{-6}\kappa$. As it can be seen from Fig. 1, the number of random
local settings chosen by either of the parties do not increase significantly
from the case when $\text{A}+\text{B}=0$. Therefore, it is possible for Alice
and Bob to detect nonlocality even when $\text{A}+\text{B}\neq 0$.

Bell inequality |  $d$ | $M$ |  $\eta$ | $\nu$
---|---|---|---|---
PNP CHSH with $n=14$ |  $2^{14}$ | $4^{14}$ |  $0.28$ | $1$
$0.40$ | $0.081$
$0.60$ | $0.008$
$0.80$ | $0.003$
$0.95$ | $0.001$
PNP CHSH with $n=13$ |  $2^{13}$ | $4^{13}$ |  $0.31$ | $1$
$0.40$ | $0.528$
$0.60$ | $0.036$
$0.80$ | $0.011$
$0.95$ | $0.006$
PNP CHSH with $n=12$ |  $2^{12}$ | $4^{12}$ |  $0.34$ | $1$
$0.60$ | $0.157$
$0.80$ | $0.045$
$0.95$ | $0.025$
PNP CHSH with $n=11$ |  $2^{11}$ | $4^{11}$ |  $0.38$ | $1$
$0.60$ | $0.728$
$0.80$ | $0.176$
$0.95$ | $0.093$
PNP CHSH with $n=10$ |  $2^{10}$ | $4^{10}$ |  $0.43$ | $1$
$0.80$ | $0.710$
$0.95$ | $0.353$
Table 1: Fraction $\nu$ of the contexts needed to detect loophole-free Bell
nonlocality for different values of the detection efficiency $\eta$, for the
correlations violating some of the PNP CHSH Bell inequalities. “Bell
inequality” indicates the Bell inequality considered. $d$ is the dimension of
each of the local quantum systems. $M$ is the total number of contexts. For
each Bell inequality, the smallest value of $\eta$ corresponds to the critical
detection efficiency when all the contexts are evaluated (i.e., when $\nu=1$).
$\nu$ is obtained assuming $V=1$, $\text{A}+\text{B}=0$, and $\delta=0.00003$.

##  IV Graph-theoretic Bell inequalities using a subset of settings

Here, we apply the tools described in Sec. II to the case of the graph-
theoretic Bell inequalities of Ref. [14].

###  IV.1 The graph-theoretic Bell inequalities

Consider two separated parties, Alice and Bob, each of them having access to a
set of measurement settings corresponding in quantum mechanics to projectors
$\\{\Pi_{i}\\}$, having each of them two outcomes, and such that the graph $G$
describes the relations of orthogonality between the members of
$\\{\Pi_{i}\\}$: each element of $\\{\Pi_{i}\\}$ is represented by a vertex of
$G$ and orthogonal projectors correspond to adjacent vertices. Then, the
corresponding graph-theoretic Bell inequality can be written as

$\displaystyle\beta_{G}$
$\displaystyle=\sum_{i\in\mathcal{V}}p(\Pi_{i}^{A}=\Pi_{i}^{B}=1)$ (17)
$\displaystyle-\sum_{(i,j)\in\mathcal{E}}\frac{1}{2\Xi}\left[p(\Pi_{i}^{A}=\Pi_{j}^{B}=1)+p(\Pi_{j}^{A}=\Pi_{i}^{B}=1)\right]\leq
C,$

where $\mathcal{V}$ is the vertex set of $G$, $\mathcal{E}$ is the edge set of
$G$, and $\Xi$ is the xi number of $G$ (which, for simplicity, can be taken to
be $1$), $p(\Pi_{i}^{A}=\Pi_{j}^{B}=1)$ is the probability for Alice and Bob
to obtain the outcome $1$ when they measure $\Pi_{i}^{A}$ and $\Pi_{j}^{B}$,
respectively. Interestingly, both the maximum local realistic value $C$ and
the maximum quantum value $Q$ can be related to properties of $G$.
Specifically, $C$ is the independence number of $G$ and $Q$ is
$\frac{|\mathcal{V}|}{\xi}$, where $\xi$ is the orthogonal rank of the graph
$G$. See Ref. [14] for details.

###  IV.2 The graph-theoretic Bell inequalities with a subset of settings

For any graph $G$, the number of settings per party in $\beta_{G}$ in (17) is
$|\mathcal{V}|$ and the total number of contexts is
$M=|\mathcal{V}|+2|\mathcal{E}|$. To estimate $\beta_{G}$, Alice and Bob
choose settings $i$ and $j$ from the set $\mathcal{V}$ with the following
nonuniform probability distribution

$p(i,j)=\begin{cases}\frac{|\mathcal{V}|}{|\mathcal{V}|+2|\mathcal{E}|}\quad
i=j\\\ \frac{2|\mathcal{E}|}{|\mathcal{V}|+2|\cal{E}|}\quad(i,j)\in\cal{E}\\\
0\quad\quad\quad\quad\text{otherwise}.\end{cases}$ (18)

Then, the estimator is

$X_{i,j}=\begin{cases}\frac{1}{p(i,j)}\beta_{i,j}\quad i=j\\\
\frac{-1}{p(i,j)}\beta_{i,j}\quad(i,j)\in\mathcal{E}\\\
0\quad\quad\quad\quad\text{otherwise},\end{cases}$ (19)

where $\beta_{i,j}=p(\Pi_{i}^{A}=\Pi_{j}^{B}=1)$. It can be seen that the mean
value of the estimator $X_{i,j}$ when averaged over the variables $i$ and $j$
is simply $\langle X_{i,j}\rangle=\beta_{G}$.

Then, Alice and Bob choose $L$ contexts according to the distribution $p(i,j)$
to evaluate the corresponding $L$ estimators $X_{i,j}$. In this case, we apply
Hoeffding’s inequality instead of Chebyshev’s inequality to evaluate the
number of contexts required. The reason is that Chebyshev’s inequality
requires evaluating the variance of each estimator in terms of the Bell value,
which is not possible in this case because of the inherent asymmetry in the
Bell inequality. Instead, we can use the fact that each estimator $X_{i,j}$ is
a function of a probability distribution and can thus be bounded.

Following Eq. (3), let $Y$ be the average value of the estimators
corresponding to the $L$ contexts chosen by the parties. Since $Y$ is a sum of
independent random variables $\frac{X_{l}}{L}$, we can apply Hoeffdings’s
inequality to bound its value. In order to do so, we can bound each of the
estimators by noting that their achievable maximum and minimum value is
$\frac{(|\mathcal{V}|+2|\mathcal{E}|)}{|\mathcal{V}|}$ and
$\frac{-(|\mathcal{V}|+2|\mathcal{E}|)}{2|\mathcal{E}|}$, respectively. Using
Hoeffding’s inequality, we have

$p(Y-\beta_{G}\geq\epsilon)\leq\delta.$ (20)

For a fixed precision $\epsilon$ and probability of failure $\delta$, the
number of contexts that Alice and Bob would need to evaluate is

$L=\frac{-\ln(\delta)\left(|\mathcal{V}|+2|\mathcal{E}|\right)^{4}}{8\epsilon^{2}|\mathcal{E}|^{2}|\mathcal{V}|^{2}}.$
(21)

Now we calculate the minimum detection efficiency $\eta_{\nu}$ needed for
loophole-free Bell violation as a function of the fraction of the contexts
$\nu$. In order to model the detection inefficiency, we bin the no clicks to
the outcome $0$. This outcome is chosen because no terms corresponding to it
appear in the Bell inequality (17). This simplifies the evaluation of the
critical detection efficiency, which, in this case, is now only a function of
the maximum local realistic and quantum values. Then, the minimum detection
efficiency is

$\eta_{\nu}=\left[\frac{1}{Q}\left(\frac{-\ln\delta\left(|\mathcal{V}|+2|\mathcal{E}|\right)^{3}}{2\nu|\mathcal{V}|^{2}|\mathcal{E}|^{2}}\right)^{\frac{1}{2}}+\frac{C}{Q}\right]^{\frac{1}{2}}.$
(22)

Next, we consider $8$ different correlations associated to the maximum quantum
violation of a different Bell inequality, each of them corresponding to a
different graph. The details of these graphs can be found in Ref. [14] but are
not relevant for our purpose. For each case, we calculate the fraction of
contexts needed for a given detection efficiency $\eta$. The results are
presented in Table 2. As it can be seen in Table 2, in all the cases, by
randomly selecting a small strict subset of the contexts, the parties can
claim loophole-free Bell nonlocality with a very small probability of failure.

Bell inequality |  $d$ | $M$ |  $\eta$ | $\nu$
---|---|---|---|---
$Y_{44}$ |  $44$ | $4.62\times 10^{24}$ |  $0.163$ | $1$
$0.200$ | $7.01\times 10^{-19}$
$0.400$ | $7.01\times 10^{-21}$
$0.600$ | $1.12\times 10^{-21}$
$0.800$ | $3.31\times 10^{-22}$
$0.950$ | $1.62\times 10^{-22}$
$Y_{36}$ |  $36$ | $7.79\times 10^{19}$ |  $0.260$ | $1$
$0.400$ | $7.07\times 10^{-16}$
$0.600$ | $7.06\times 10^{-17}$
$0.800$ | $1.84\times 10^{-17}$
$0.950$ | $8.66\times 10^{-18}$
$Y_{32}$ |  $32$ | $3.22\times 10^{17}$ |  $0.326$ | $1$
$0.400$ | $4.51\times 10^{-13}$
$0.600$ | $2.03\times 10^{-14}$
$0.800$ | $4.54\times 10^{-15}$
$0.950$ | $2.02\times 10^{-15}$
$Y_{28}$ |  $28$ | $1.34\times 10^{15}$ |  $0.407$ | $1$
$0.600$ | $7.19\times 10^{-12}$
$0.800$ | $1.20\times 10^{-12}$
$0.950$ | $4.99\times 10^{-13}$
$\mathcal{P}_{4}(\mathds{R})$ |  $16$ | $8752320$ |  $0.516$ | $1$
$0.600$ | $3.84\times 10^{-4}$
$0.800$ | $2.40\times 10^{-4}$
$0.950$ | $8.29\times 10^{-5}$
$\mathcal{P}_{3}(\mathds{C})$ |  $8$ | $341280$ |  $0.730$ | $1$
$0.750$ | $0.098$
$0.850$ | $0.002$
$0.950$ | $6.17\times 10^{-4}$
$\mathcal{P}_{3}(\mathds{R})$ |  $8$ | $25440$ |  $0.730$ | $1$
$0.850$ | $0.072$
$0.950$ | $0.019$
$\mathcal{P}_{2}(\mathds{C})$ |  $4$ | $960$ |  $0.894$ | $1$
$0.950$ | $0.668$
$\mathcal{P}_{2}(\mathds{R})$ |  $4$ | $240$ |  $0.912$ | $1$
Table 2: Fraction $\nu$ of the contexts needed to detect loophole-free Bell
nonlocality for different values of the detection efficiency $\eta$, for the
correlations violating some of the graph-theoretic Bell inequalities. “Bell
inequality” indicates the Bell inequality considered. $d$ is the dimension of
each of the local quantum systems. $M$ is the total number of contexts. Here,
we assume $V=1$ and $\delta=0.00003$.

##  V Discussion

###  V.1 Pros and cons of the approach

So far, we have shown that, for some bipartite correlations, a small fraction
of the contexts is enough to detect loophole-free Bell nonlocality with a
given confidence. For example, in the PNP CHSH Bell inequality with $n=13$, we
have shown that the parties need only a fraction of $0.528$ of the contexts to
detect loophole-free nonlocality (see Table 1). However, then the detection
efficiency must be $\eta\geq 0.4$, which is higher than the critical detection
efficiency had the parties measured all contexts (which is $\eta_{\rm
crit}=0.313$). Moreover, the number of settings ($0.528\times 10^{13})$ that
the parties have to measure is still too large to be practical.

Therefore, measuring only a fraction of the contexts is sometimes enough to
detect loophole-free Bell nonlocality. However, this comes at the cost of a
higher detection efficiency and, at least in the example, the reduction in the
number of settings is not enough to be practical for standard experiments.

Then the question is: What is this approach useful for? In the following, we
will argue that: (i) It provides a new way to design loophole-free Bell tests.
(ii) It enforces a different way to look at the classical simulation of
quantum correlations, as it shows that the amount of classical resources may
depend on how large the fraction of random contexts is, which is a choice that
can be modified during the experiment. (iii) It stimulates the search for
correlations with many settings and low detection efficiency, to which no
attention has been paid so far.

###  V.2 Designing loophole-free Bell tests

The tools introduced here offer a new approach to think of loophole-free Bell
tests. Suppose that we want to achieve loophole-free nonlocality with quantum
systems of dimension, e.g., $32$. So far, the only possibility was conducting
experiments targeting correlations (and Bell inequalities) having a critical
detection efficiency smaller than the detection efficiency that we have in our
laboratory. The problem is that we do not have practical (i.e., with a
sufficiently small number of settings) examples to base our experiment on.

However, there are quantum correlations in $d=32$ which, if we would have
infinite time, would only require $\eta_{\rm crit}=0.326$ (see Table 2). We
may not have infinite time but we may be curious about what would be the
result of the Bell test if we would have infinite time. In that respect, it is
useful to notice that, in all the cases explored here, all measurements are
equally difficult from an experimental point of view. For example, all the
projectors to be measured in the graph-theoretic Bell inequalities associated
to the graphs $Y_{d}$ only differ in phases (see [14] for details). Therefore,
by measuring a random subset of the contexts, we will get a reliable
estimation of what we would get if we measured all the contexts.

Not only that. Now, if we know the visibility and the detection efficiency
that we can achieve in the laboratory, we can answer the question of what is
the minimum fraction of contexts needed to detect loophole-free nonlocality.
For example, consider the Bell inequality $Y_{32}$. As said, it requires a
minimum detection efficiency $\eta_{\rm crit}=0.326$ when all contexts are
tested. Suppose that we have, e.g., $\eta\geq 0.5$, and we want to actually
conduct a loophole-free Bell test. Then, the tools introduced here show that
we can do it using a fraction of contexts $\nu\geq 2.5\times 10^{-13}$.
Similarly, if we have a higher detection efficiency, e.g., $\eta\geq 0.6$,
then the fraction of contexts further reduces to $\nu\geq 2\times 10^{-14}$,
which can be further reduced to $\nu\geq 4.5\times 10^{-15}$ if the detection
efficiency is $\eta\geq 0.8$, etc. This way, knowing the visibility and
detection efficiency, we can choose the corresponding fraction of contexts
$\nu$ needed for a loophole-free Bell test.

###  V.3 Making classical simulations asymptotically impossible

The tools introduced here also enforce a new way to look at the cost of
simulating Bell nonlocality. In a local realistic model that simulates certain
nonlocal correlations, the detection efficiency is also a target of the
simulation (see, e.g., [29, 30, 31]). That is, the local realistic model
should not only reproduce the quantum statistics but should also simulate a
fixed detection efficiency $\eta<\eta_{\rm crit}$. However, there are
correlations that can be classically simulated for a fixed $\eta$ that cannot
be simulated when more contexts are added. The problem is identifying them.

For example, consider a local realistic model which simulates $\eta=0.9$ and
the correlations maximally violating the Bell inequality in
$\mathcal{P}_{2}(\mathds{R})$ for the $240$ contexts in Table 2. Such a local
realistic model is possible because the simulated $\eta$ is smaller than the
corresponding critical detection efficiency $\eta_{\rm crit}=0.912$, otherwise
the simulation would be impossible. Now notice that the correlations maximally
violating the $\mathcal{P}_{2}(\mathds{R})$ Bell inequality are a subset of
the correlations maximally violating $\mathcal{P}_{2}(\mathds{C})$ (see [14]
for details). If the parties decide to test more contexts until they cover all
the $960$ contexts needed for the correlations maximally violating the
$\mathcal{P}_{2}(\mathds{C})$ Bell inequality, then no local realistic model
simulating $\eta=0.9$ can also simulate the correlations. The reason is that
simultaneously simulating the correlations for $\mathcal{P}_{2}(\mathds{C})$
and $\eta=0.9$ is impossible since, in this case, $\eta_{\rm crit}=0.894$ (see
Table 2).

However, it is not true that this is always the case. For example, the
correlations maximally violating the $\mathcal{P}_{3}(\mathds{R})$ Bell
inequality are a subset of the correlations maximally violating
$\mathcal{P}_{3}(\mathds{C})$. However, the critical detection efficiency is
the same in both cases (see Table 2).

Moreover, both examples above refer to nonlocal correlations in which we can
identify a subset of them which is also nonlocal and we can compute the
corresponding $\eta_{\rm crit}$. The problem is that this is something that is
very difficult in general. Nevertheless, in this work we have found
correlations (all those in Tables 1 and 2) for which a smaller randomly chosen
fraction always requires a larger $\eta_{\rm crit}$ than a larger randomly
chosen fraction. As mentioned, this is not true in general, even in the case
when the parties are allowed to choose specific subsets. However, the power of
our approach is that now the parties can calculate $\eta_{\rm crit}$ for any
fraction of contexts. In all the cases in Tables 1 and 2 (except for the last
correlation in Table 2; this is why there is only a single row there), by
simply increasing the fraction of contexts, the parties can make it impossible
for a local realistic model to simulate the correlations.

This brings us to a final thought: How does the classical simulation know in
advance what fraction of contexts to simulate? For not failing, the simulation
should work for all contexts. There are probably cases in which the resources
needed for the simulation rapidly tend to infinity as the number of contexts
increases. This question and this conjecture indicate that it would be
interesting to investigate how fast these resources can grow with the number
of contexts.

###  V.4 Further research

While our tools allow for the detection of loophole-free nonlocality by using
only fractions of contexts (at the cost of a higher detection efficiency),
when we apply these tools to the correlations in Tables 1 and 2, then the
number of measurements required is either still too large for practical Bell
tests or not too large but then they offer critical detection efficiencies
that are comparable with the ones needed in existing detection-loophole-free
experiments with smaller dimension and number of settings. That is, the
examples used in this work do not yet provide practical targets with
sufficiently low critical detection efficiency and number of settings.

However, it is important to emphasize that the examples we have used to
illustrate our tools are based on correlations obtained by sophisticated but
arguably sub optimal methods with the purpose of showing that the critical
detection efficiency can be arbitrary low without the need of quantum systems
of dimension impossible to achieve in the laboratory (as it was thought).
Presumably, there are many quantum correlations with sufficiently low critical
detection and number of settings waiting to be discovered. Our hope is that
the tools introduced here stimulate the search for such correlations.

## Acknowledgments

This work is supported by the QuantERA project SECRET by MCINN/AEI (Project
No. PCI2019-111885-2). A.C. is also supported by the EU-funded project FoQaCiA
and MCINN/AEI (Project No. PID2020-113738GB-I00).

## References

  * Nadlinger  _et al._ [2022] D. P. Nadlinger, P. Drmota, B. C. Nichol, G. Araneda, D. Main, R. Srinivas, D. M. Lucas, C. J. Ballance, K. Ivanov, E. Y. Z. Tan, P. Sekatski, R. L. Urbanke, R. Renner, N. Sangouard, and J.-D. Bancal, Experimental quantum key distribution certified by Bell’s theorem, Nature (London) 607, 682 (2022).
  * Zhang  _et al._ [2022] W. Zhang, T. van Leent, K. Redeker, R. Garthoff, R. Schwonnek, F. Fertig, S. Eppelt, W. Rosenfeld, V. Scarani, C. C. W. Lim, and H. Weinfurter, A device-independent quantum key distribution system for distant users, Nature (London) 607, 687 (2022).
  * Liu  _et al._ [2022] W.-Z. Liu, Y.-Z. Zhang, Y.-Z. Zhen, M.-H. Li, Y. Liu, J. Fan, F. Xu, Q. Zhang, and J.-W. Pan, Toward a photonic demonstration of device-independent quantum key distribution, Phys. Rev. Lett. 129, 050502 (2022).
  * Primaatmaja  _et al._ [2023] I. W. Primaatmaja, K. T. Goh, E. Y.-Z. Tan, J. T.-F. Khoo, S. Ghorai, and C. C.-W. Lim, Security of device-independent quantum key distribution protocols: a review, Quantum 7, 932 (2023).
  * van Leent  _et al._ [2023] T. van Leent, R. Arnon-Friedman, W.-Z. Liu, Q. Zhang, H. Weinfurter, and M. Curty, Advances in device-independent quantum key distribution, NPJ Quantum Inf. 9, 10 (2023).
  * Pironio  _et al._ [2010] S. Pironio, A. Acín, S. Massar, A. B. de la Giroday, D. N. Matsukevich, P. Maunz, S. Olmschenk, D. Hayes, L. Luo, T. A. Manning, and C. Monroe, Random numbers certified by Bell’s theorem, Nature (London) 464, 1021 (2010).
  * Šupić and Bowles [2020] I. Šupić and J. Bowles, Self-testing of quantum systems: a review, Quantum 4, 337 (2020).
  * Acín and Masanes [2016] A. Acín and L. Masanes, Certified randomness in quantum physics, Nature (London) 540, 213 (2016).
  * Shalm  _et al._ [2021] L. K. Shalm, Y. Zhang, J. C. Bienfang, C. Schlager, M. J. Stevens, M. D. Mazurek, C. Abellán, W. Amaya, M. W. Mitchell, M. A. Alhejji, H. Fu, J. Ornstein, R. P. Mirin, S. W. Nam, and E. Knill, Device-independent randomness expansion with entangled photons, Nat. Phys. 17, 452 (2021).
  * Freedman and Clauser [1972] S. J. Freedman and J. F. Clauser, Experimental test of local hidden-variable theories, Phys. Rev. Lett. 28, 938 (1972).
  * Weihs  _et al._ [1998] G. Weihs, T. Jennewein, C. Simon, H. Weinfurter, and A. Zeilinger, Violation of Bell’s inequality under strict Einstein locality conditions, Phys. Rev. Lett. 81, 5039 (1998).
  * Pearle [1970] P. M. Pearle, Hidden-variable example based upon data rejection, Phys. Rev. D 2, 1418 (1970).
  * Miklin  _et al._ [2022] N. Miklin, A. Chaturvedi, M. Bourennane, M. Pawłowski, and A. Cabello, Exponentially decreasing critical detection efficiency for any Bell inequality, Phys. Rev. Lett. 129, 230403 (2022).
  * Xu  _et al._ [2023] Z.-P. Xu, J. Steinberg, J. Singh, A. J. López-Tarrida, J. R. Portillo, and A. Cabello, Graph-theoretic approach to Bell experiments with low detection efficiency, Quantum 7, 922 (2023).
  * Garg and Mermin [1987] A. Garg and N. D. Mermin, Detector inefficiencies in the Einstein-Podolsky-Rosen experiment, Phys. Rev. D 35, 3831 (1987).
  * Eberhard [1993] P. H. Eberhard, Background level and counter efficiencies required for a loophole-free Einstein-Podolsky-Rosen experiment, Phys. Rev. A 47, R747 (1993).
  * Vértesi  _et al._ [2010] T. Vértesi, S. Pironio, and N. Brunner, Closing the detection loophole in Bell experiments using qudits, Phys. Rev. Lett. 104, 060401 (2010).
  * Rowe  _et al._ [2001] M. A. Rowe, D. Kielpinski, V. Meyer, C. A. Sackett, W. M. Itano, C. Monroe, and D. J. Wineland, Experimental violation of a Bell’s inequality with efficient detection, Nature (London) 409, 791 (2001).
  * Giustina  _et al._ [2013] M. Giustina, A. Mech, S. Ramelow, B. Wittmann, J. Kofler, J. Beyer, A. Lita, B. Calkins, T. Gerrits, S. W. Nam, R. Ursin, and A. Zeilinger, Bell violation using entangled photons without the fair-sampling assumption, Nature (London) 497, 227 (2013).
  * Christensen  _et al._ [2013] B. G. Christensen, K. T. McCusker, J. B. Altepeter, B. Calkins, T. Gerrits, A. E. Lita, A. Miller, L. K. Shalm, Y. Zhang, S. W. Nam, N. Brunner, C. C. W. Lim, N. Gisin, and P. G. Kwiat, Detection-loophole-free test of quantum nonlocality, and applications, Phys. Rev. Lett. 111, 130406 (2013).
  * Shalm  _et al._ [2015] L. K. Shalm, E. Meyer-Scott, B. G. Christensen, P. Bierhorst, M. A. Wayne, M. J. Stevens, T. Gerrits, S. Glancy, D. R. Hamel, M. S. Allman, K. J. Coakley, S. D. Dyer, C. Hodge, A. E. Lita, V. B. Verma, C. Lambrocco, E. Tortorici, A. L. Migdall, Y. Zhang, D. R. Kumor, W. H. Farr, F. Marsili, M. D. Shaw, J. A. Stern, C. Abellán, W. Amaya, V. Pruneri, T. Jennewein, M. W. Mitchell, P. G. Kwiat, J. C. Bienfang, R. P. Mirin, E. Knill, and S. W. Nam, Strong loophole-free test of local realism, Phys. Rev. Lett. 115, 250402 (2015).
  * Giustina  _et al._ [2015] M. Giustina, M. A. M. Versteegh, S. Wengerowsky, J. Handsteiner, A. Hochrainer, K. Phelan, F. Steinlechner, J. Kofler, J.-Å. Larsson, C. Abellán, W. Amaya, V. Pruneri, M. W. Mitchell, J. Beyer, T. Gerrits, A. E. Lita, L. K. Shalm, S. W. Nam, T. Scheidl, R. Ursin, B. Wittmann, and A. Zeilinger, Significant-loophole-free test of Bell’s theorem with entangled photons, Phys. Rev. Lett. 115, 250401 (2015).
  * Hensen  _et al._ [2015] B. Hensen, H. Bernien, A. E. Dréau, A. Reiserer, N. Kalb, M. S. Blok, J. Ruitenberg, R. F. L. Vermeulen, R. N. Schouten, C. Abellán, W. Amaya, V. Pruneri, M. W. Mitchell, M. Markham, D. J. Twitchen, D. Elkouss, S. Wehner, T. H. Taminiau, and R. Hanson, Loophole-free Bell inequality violation using electron spins separated by 1.3 kilometres, Nature (London) 526, 682 (2015).
  * Rosenfeld  _et al._ [2017] W. Rosenfeld, D. Burchardt, R. Garthoff, K. Redeker, N. Ortegel, M. Rau, and H. Weinfurter, Event-ready Bell test using entangled atoms simultaneously closing detection and locality loopholes, Phys. Rev. Lett. 119, 010402 (2017).
  * Hu  _et al._ [2022] X.-M. Hu, C. Zhang, B.-H. Liu, Y. Guo, W.-B. Xing, C.-X. Huang, Y.-F. Huang, C.-F. Li, and G.-C. Guo, High-dimensional Bell test without detection loophole, Phys. Rev. Lett. 129, 060402 (2022).
  * Tchébychef [1867] P.-L. Tchébychef, Des valeurs moyennes, J. Math. Pures Appl. 2e série 12, 177 (1867).
  * Clauser  _et al._ [1969] J. F. Clauser, M. A. Horne, A. Shimony, and R. A. Holt, Proposed experiment to test local hidden-variable theories, Phys. Rev. Lett. 23, 880 (1969).
  * Hoeffding [1963] W. Hoeffding, Probability inequalities for sums of bounded random variables, J. Am. Stat. Assoc. 58, 13 (1963).
  * Larsson [1999] J.-Å. Larsson, Modeling the singlet state with local variables, Phys. Lett. A 256, 245 (1999).
  * Aerts  _et al._ [1999] S. Aerts, P. Kwiat, J.-Å. Larsson, and M. Żukowski, Two-photon Franson-type experiments and local realism, Phys. Rev. Lett. 83, 2872 (1999).
  * Cabello  _et al._ [2009] A. Cabello, A. Rossi, G. Vallone, F. De Martini, and P. Mataloni, Proposed Bell experiment with genuine energy-time entanglement, Phys. Rev. Lett. 102, 040401 (2009).


# Control-Oriented Identification for the Linear Quadratic Regulator:
Technical Report

Sean Anderson, João P. Hespanha This material is based upon work supported by
the National Science Foundation Graduate Research Fellowship under Grant No.
2139319 and by the U.S. Office of Naval Research MURI grant No.
N00014-23-1-2708. Any opinion, findings, and conclusions or recommendations
expressed in this material are those of the authors(s) and do not necessarily
reflect the views of the National Science Foundation.The authors are with
Department of Electrical and Computer Engineering, University of California
Santa Barbara, Santa Barbara, CA, 93106 USA<EMAIL_ADDRESS><EMAIL_ADDRESS>
###### Abstract

Data-driven control benefits from rich datasets, but constructing such
datasets becomes challenging when gathering data is limited. We consider an
offline experiment design approach to gathering data where we design a control
input to collect data that will most improve the performance of a feedback
controller. We show how such a control-oriented approach can be used in a
setting with linear dynamics and quadratic objective and, through design of a
gradient estimator, solve the problem via stochastic gradient descent. We
contrast our method with a classical A-optimal experiment design approach and
numerically demonstrate that our method outperforms A-optimal design in terms
of improving control performance.

##  I Introduction

Model-based control methods benefit from useful models of the controlled
system. Consider a setting in which there is uncertainty in the model
parameters and there is an opportunity to collect experimental data to learn
more about the system. The data collection involves selecting a control input,
and this paper focuses on the optimal selection of such an input with the goal
of eventually generating a “high-performance” controller for the uncertain
process. In this context, “high-performance” is defined in terms of a pre-
specified criterion that we use to evaluate the performance of a given
controller. This motivates the following control-oriented experiment design
problem: select a control input for a data-collection experiment so that the
feedback controller designed using the data acquired will lead to the best
possible performance. Our approach is general in terms of the control design
procedure used to generate the controller from the experimental data
collected. However, in this paper we focus our attention on controllers
generated through certainty equivalence, which in this context means
constructing an _a-posteriori_ estimate for the process and designing a
controller for this estimate.

This paper includes two key contributions: First, we show how such a control-
oriented approach to experiment design can be carried out for the control of a
linear system with uncertain matrix dynamics and a quadratic objective
function. While this problem does not have a closed-form solution, we show
that it can be efficiently solved by stochastic gradient descent. The second
contribution lies in the observation that our formulation of a control-
oriented approach to experiment design generally leads to closed-loop
controllers that exhibit higher performance than what would be achieved by
classical forms of experiment design, which are typically aimed at minimizing
_a-posteriori_ estimation error (or other generic forms of post-experiment
uncertainty), rather than focusing directly on improving closed-loop
performance.

In this work, we focus on the offline experiment design setting for data-
driven control. In Section II we present a general formulation for experiment
design that takes into account a-priori parameter uncertainty in discrete-time
dynamics. We consider an experiment that generates a dataset consisting of
state and input measurements, where the dataset is used in the construction of
a controller. The resulting optimal experiment design problem aims to improve
the expected performance of the controller with respect to possible parameter
values as well process noise during the experiment. This formulation is
contrasted with classical experiment design formulations that aim to minimize
a measure of the parameter uncertainty [1].

We then turn to the linear quadratic regulator (LQR) setting in Section III,
where we define the linear dynamics and quadratic control objective. We
address the system identification step and how to handle exploding
trajectories during (simulated) experiments. For the data-driven controller,
we use an LQR controller with certainty equivalence in the parameters [2]. The
solutions available in LQR are amenable to fast computations, which is
conducive to the numerical approach we take later on.

We go on to present the solution method to the experiment design approach in
Section IV. We use a first-order approach by designing a pathwise gradient
estimator for the purposes of stochastic gradient descent. By leveraging the
structure of the problem, we find a gradient estimate that proves to be
effective numerically as illustrated in Section V. In order to validate our
method, we compare our approach against a classical experiment design that
aims to reduce a measure of the uncertainty in the parameters. We first
illustrate our method’s advantages when compared against the classical
experiment design and then go on to show how the problem can be scaled to
systems of a moderate order.

Related work: The issue of how to gather data through well-planned experiments
has traditionally been addressed through the framework of optimal experiment
design. Modern optimal experiment design is often attributed to Gustav
Elfving, who designed experiments to minimize measures of parameter error
covariance [3]. Later on, researchers worked on aligning experiment design
with particular criteria, including control objectives [4]. See [5] for a
thorough review of early work in this area. Recent relevant work in this area
includes [6], which proposes a stochastic gradient descent approach to
designing experiments that minimizes a post-experiment optimal control
objective. Work in experiment design for control in the statistical learning
community includes [7, 8, 9], all of which emphasize theoretical aspects of
learning linear systems. As an alternative paradigm, online learning or
adaptive control allow for improvement during the experiment trial. Adaptive
control has been well-covered in [10], and online linear quadratic regulation
has been explored from learning theory in [11].

As we focus on the linear quadratic regulator setting, the most relevant
papers are from [7], in which the authors derive fundamental limits for
learning the linear quadratic regulator via offline experiments, and [8], in
which the authors propose an experiment design method for learning linear
systems for particular tasks. The former focuses on the statistical learning
guarantees associated with the LQR problem but do not broach how to construct
an optimal experiment. In the latter, the authors propose a weighted-trace
optimal experiment design that bears similarity to $L$-optimal experiment
design for control-oriented identification as noted in [5]; the authors
consider convergence of sequential experiment designs to the underlying
system.

Gradient estimation has seen particular attention in the machine learning
community [12], and two dominant methods are via the pathwise gradient and the
score function (or log score) method [13]. Score function estimators benefit
from only taking the gradient of the density function, but the structure of
our problem is not amenable: as such, we focus on the pathwise gradient
estimate.

##  II Experiment Design for Data-Driven Control

We consider a discrete-time system with dynamics of the form

$\displaystyle x_{t+1}=f(x_{t},u_{t},w_{t};\theta),$ (1)

with state $x_{t}\in\mathbb{R}^{n_{x}}$, control input
$u_{t}\in\mathbb{R}^{n_{u}}$, and an unmeasured stochastic disturbance
$w_{t}\in\mathbb{R}^{n_{w}}$ independent and identically distributed across
time. The dynamics depend on stochastic parameters $\theta$ that are unknown,
but for which we have an _a-priori_ distribution.

An experiment will be performed to provide additional information about the
parameter $\theta$. State and input measurements are collected throughout the
experiment providing a sequence of $M$ triples
$\mathcal{D}:=\\{(x_{i}^{+},x_{i},u_{i}):i=1,...,M\\}$ that satisfy the basic
model of the dynamical system:

$\displaystyle x_{i}^{+}=f(x_{i},u_{i},w_{i};\theta),\quad\forall
i\in\\{1,\dots,M\\}$ (2)

with the $w_{i}$ independent across indices $i$ and with the same distribution
as the disturbance. If the experiment consists of a single run of (1) over a
time horizon $t=1$ through $t=T$, then the index $i$ is simply time and $M=T$.
However, in general, “an experiment” may include multiple runs of (1) over
different time horizons, in which case (2) include all the data collected. To
simplify notation we collect all the columns vectors
$x_{i}^{+},x_{i}\in\mathbb{R}^{n_{x}},u_{i}\in\mathbb{R}^{n_{u}}$ into
matrices with $M$ columns that we denote by $X^{+}$,
$X\in\mathbb{R}^{n_{x}\times M},U\in\mathbb{R}^{n_{u}\times M}$, respectively.

Our goal is to design a controller $\pi$ that optimizes a given cost function
$J(\pi;\theta)$ that depends both on the controller and on the actual (but
unknown) value of the parameter $\theta$. We also take as given a control
design procedure that maps the experiment design data $\mathcal{D}$ to a
specific controller $\pi$, with the goal of minimizing the cost
$J(\pi;\theta)$. A reasonable procedure would be to select a controller that
minimizes the expected value of the cost $J(\pi;\theta)$, given the data
$\mathcal{D}$ collected during the experiment, as in:

$\displaystyle
K(\mathcal{D})=\operatorname*{arg\,min}_{\pi}\operatorname*{\mathbb{E}}[J(\pi;\theta)\,|\,\mathcal{D}].$
(3)

However, and because this optimization is often intractable, our presentation
considers a general control design method $K:\mathcal{D}\mapsto\pi$, which may
or may not be optimal.

The experiment design problem arises from the observation that the data
$\mathcal{D}$ collected depends on the control inputs $U$ used during the
experiment as well as on the actual realizations of the the random
disturbances (also during the experiment) that we denote by
$W\in\mathcal{R}^{n_{w}\times M}$. To emphasize this dependence, we use the
notation $\mathcal{D}_{U,W}$ to express the dependence of the dataset on these
variables. The optimal experiment design problem can then be formulated as

$\displaystyle\min_{U\in\mathcal{U}}\operatorname*{\mathbb{E}}\big{[}J\big{(}\pi;\theta\big{)}\big{]},\quad\pi:=K\big{(}\mathcal{D}_{U,W}\big{)},$
(4)

where the expectation refers to an integration over (i) the _a-priori_
distribution of the parameter $\theta$, and (ii) the realization of the
disturbance $W$ during the actual experiment. The minimization is performed
over a set of admissible controls that we denote generically by $\mathcal{U}$.
The experiment design criterion (4) should be contrasted with more classical
formulations that set the goal of the experiment to minimize some measure of
uncertainty about the unknown parameter $\theta$. For example, an $A$-optimal
experiment design essentially tries to minimize

$\displaystyle\min_{U\in\mathcal{U}}\operatorname*{\mathbb{E}}\Big{[}\big{\|}\hat{\theta}\big{(}\mathcal{D}_{U,W}\big{)}-\theta\big{\|}^{2}\Big{]},\quad\hat{\theta}\big{(}\mathcal{D}_{U,W}\big{)}:=\operatorname*{\mathbb{E}}\big{[}\theta\,|\,\mathcal{D}_{U,W}\big{]}.$
(5)

The key distinction is that (5) ignores the impact of uncertainty on the
control objective $J(\pi;\theta)$ and therefore does not take advantage of the
fact that reducing uncertainty on some parameters may be much more important
than on others, _for our ultimate control objective of minimizing
$J(\pi;\theta)$._

##  III Experiment Design for the Linear Quadratic Regulator

We now specialize the general setup described above to the finite-horizon
linear quadratic regulator setup. Specifically, we consider the process

$\displaystyle x_{t+1}=Ax_{t}+Bu_{t}+w_{t},$ (6)

such that $\theta$ contains the elements of $A$ and $B$, and $w_{t}$ is
Gaussian noise identically distributed across time with zero mean and
covariance $\Sigma_{w}$.

We consider a quadratic optimization criterion of the form:

$\displaystyle
J(\pi;\theta):=\mathbb{E}\bigg{[}x_{N}^{T}Q_{N}x_{N}+\sum_{t=0}^{N-1}x_{t}^{T}Qx_{t}+u_{t}^{T}Ru_{t}\,|\,\theta\bigg{]},$
(7)

where the expectation refers to an integration over the disturbances
encountered by the controller $\pi$; $Q_{N},Q$ are positive semidefinite
matrices; and $R$ a positive definite matrix.

We consider a common option for control design generally known as _certainty
equivalence (CE)_ : certainty equivalence design $K_{CE}(\mathcal{D})$
computes the _a-posteriori_ expected value of the unknown parameters
$\hat{\theta}:=\operatorname*{\mathbb{E}}[\theta\,|\,\mathcal{D}]$ and
computes the linear optimal controller $u_{t}=K_{t}x_{t}$ that minimizes (7),
_assuming that the estimate $\hat{\theta}$ is correct_.

####  III-1 System identification

In order to generate the _a-posteriori_ estimate of the system $\hat{\theta}$
for $K_{CE}(\mathcal{D})$, we employ linear regression on a dataset
$\mathcal{D}$, which in our case will be the dataset generated under the
experiment decision variable $U$. For identification, we express (6) as:

$\displaystyle X^{+}=\Theta Z+W$ (8)

with $Z=[X;U]\in\mathbb{R}^{(n_{x}+n_{u})\times M}$, and
$\Theta:=[A,B]\in\mathbb{R}^{n_{x}\times(n_{x}+n_{u})}$. For ease of notation,
we use $\theta\in\mathbb{R}^{n_{x}(n_{x}+n_{u})}$ to denote the vectorized
version of $\Theta$ via stacking its columns.

Corollary 1: [14] Consider a Gaussian prior on the parameters with mean
$\Theta_{0}\in\mathbb{R}^{n_{x}\times(n_{x}+n_{u})}$ and covariance of the
$(i,j)$th element with the $(k,l)$th element of $\Theta$ given by
$\mathbb{E}[(\hat{\Theta}-\Theta)_{ij}(\hat{\Theta}-\Theta)_{kl}]=(\Sigma_{w})_{ki}(\Lambda_{0}^{-1})_{jl}$,
where $\Sigma_{w}$ is the known noise covariance and
$\Lambda_{0}^{-1}\in\mathbb{R}^{(n_{x}+n_{u})\times(n_{x}+n_{u})}$ is a prior
on the parameter covariance. The weighted Bayesian estimator for $\Theta$ is

$\displaystyle\hat{\Theta}=(\Theta_{0}\Lambda_{0}+X^{+}SZ^{T})\Lambda_{n}^{-1},$
(9a) and the error covariance of the estimate $\hat{\Theta}$ is
$\displaystyle\mathbb{E}[(\hat{\Theta}$
$\displaystyle-\Theta)_{ij}(\hat{\Theta}-\Theta)_{kl}]=(\Sigma_{w})_{ki}(\Lambda_{n}^{-1})_{jl},$
(9b) where $\Lambda_{n}:=\Lambda_{0}+ZSZ^{T}$, and $S\in\mathbb{R}^{M\times
M}$ is the weight matrix.

###### Proof:

See [14]. ∎

The weight matrix $S$ improves the numerics of the regression problem,
particularly since simulating unstable systems can lead to exponential growth
in the state that, due to large numbers, lead to deleterious performance in
the inversion of $\Lambda_{n}$. In particular, let

$\displaystyle S(X):=\text{diag}([s(x_{0}),...,s(x_{N})])$ (10)

where $s(x)\in[0,1]$ ensures that the weight matrix assigns zero weight to
points on trajectories that are numerically too large. For this work, we
choose $s(x):=\arctan(\norm{x_{t}-\alpha_{1}}\alpha_{2})/\pi+0.5$, and
$\alpha_{1},\alpha_{2}$ are design parameters.

####  III-2 Certainty equivalent control

Given an estimate of the parameters $\theta$ from (9) with means $\hat{A}$ and
$\hat{B}$, respectively, we construct our controller $K_{CE}(\mathcal{D})$ by
recursively solving the Riccati difference equations given by

$\displaystyle K_{t}$
$\displaystyle=-(R+\hat{B}^{T}P_{t+1}\hat{B})^{-1}\hat{B}^{T}P_{t+1}\hat{A},$
(11a) $\displaystyle P_{t}$
$\displaystyle=Q+K_{t}^{T}RK_{t}-(\hat{A}+\hat{B}K_{t})^{T}P_{t+1}(\hat{A}+\hat{B}K_{t}),$
(11b)

with $P_{N}=Q_{N}$; $Q,Q_{N}$ are positive semidefinite matrices and $R$ a
positive definite matrix.

Corollary 2: [15] For a sequence of linear feedback gains,
$\pi:=\\{K_{0},...,K_{N-1}\\}$ from $K_{CE}(\mathcal{D})$, we can express the
finite-horizon LQR cost (7) for the system in (6) parameterized by $\theta$ as

$\displaystyle J(\pi;\theta)$
$\displaystyle=x_{0}^{T}P_{0}x_{0}+\sum_{t=0}^{N-1}\tr(P_{t+1}\Sigma_{w}),$
(12a) where $\displaystyle P_{t}$
$\displaystyle=Q+K_{t}^{T}RK_{t}-(A+BK_{t})^{T}P_{t+1}(A+BK_{t}),$ (12b) with
boundary condition $P_{N}=Q_{N}$.

###### Proof:

See, for example, [15]. ∎

##  IV Experiment design via gradient descent

In order to solve the experiment design optimization (4), we take a gradient-
descent approach:

$\displaystyle
U_{i+1}=\text{Proj}_{\mathcal{U}}(U_{i}-\eta_{i}\hat{\nabla}_{U}),$ (13)

where $\text{Proj}_{\mathcal{U}}(\cdot)$ projects $U$ onto the set of
admissible inputs $\mathcal{U}$, $\eta_{i}$ is step size, and
$\hat{\nabla}_{U}$ is an approximation of the true gradient $\nabla_{U}$ given
by differentiating the experiment criteria (4):

$\displaystyle\begin{split}\nabla_{U}\operatorname*{\mathbb{E}}\Big{[}&J\Big{(}K\big{(}\mathcal{D}\big{)};\theta\Big{)}\Big{]}=\\\
&\nabla_{U}\int
J\big{(}K\big{(}\mathcal{D}\big{)};\theta\big{)}p(X\,|\,U,x_{0},\theta)p(\theta)d{\theta}dX,\end{split}$
(14)

where $x_{0}$ is the initial state for the experiment. The gradient of the
integral is analytically intractable, motivating the use of a Monte Carlo
gradient estimate. A gradient estimate generally requires an exchange of the
gradient and integral

$\displaystyle\nabla_{U}\int
F(X,U)p(X|U)dX=\int\nabla_{U}\big{(}F(X,U)p(X|U)\big{)}dX,$ (15)

for an arbitrary $F(\cdot)$ and density $p(\cdot)$, where the exchange is
valid when (i) the magnitude of the gradient of the integrand is bounded by a
function (ii) that is integrable with respect to the random variable. See [16]
for further treatment of the exchange.

By introducing a change of variable $g(\cdot)$ and using the Law of the
Unconscious Statistician [17], we remove the need to differentiate the density
such that

$\displaystyle\begin{split}&\nabla_{U}\int F(X,U)p(X\,|\,U)dX=\\\ &\int
p(\epsilon)\nabla_{U}F(X,U)\rvert_{X=g(\epsilon;U,x_{0},\theta)}d\epsilon,\leavevmode\nobreak\
\quad\epsilon\sim p(\epsilon),\end{split}$ (16)

and only $F(\cdot)$ needs to differentiated. This leads to the Monte Carlo
pathwise gradient estimator [12]:

$\displaystyle\frac{1}{L}\sum_{l=1}^{L}\nabla_{U}F(X,U)\rvert_{X=g(\epsilon^{(l)};U)},\quad\epsilon^{(l)}\sim
p(\epsilon).$ (17)

For our problem in (14), we can find a suitable change of variable using the
dynamics.

Theorem 1: The change of variable $g:\mathbb{R}^{n_{x}\times
M}\times\mathbb{R}^{n_{u}\times
M}\times\mathbb{R}^{n_{x}}\times\mathbb{R}^{n_{\theta}}\to\mathbb{R}^{n_{x}\times
M}$ for (14) given by the recursion of the dynamics (1):

$\displaystyle x_{1}$ $\displaystyle=f(x_{0},u_{0},w_{0};\theta),$ (18a)
$\displaystyle x_{2}$
$\displaystyle=f(f(u_{0},x_{0},w_{0};\theta),u_{1},w_{1};\theta),$ (18b)
$\displaystyle\vdots$ $\displaystyle x_{M}$
$\displaystyle=f(f(x_{M-2},u_{M-2},w_{M-2};\theta),u_{M-1},w_{M-1};\theta),$
(18c) where $x_{M}$ is the final state, recursively dependent on all preceding
states, satisfies (16) such that $X=g(W;U,x_{0},\theta)$ with $W$ distributed
according to $p(W)$, the process noise.

###### Proof:

Under the process (1), use the Markov property to write the density of $X$ as
$\displaystyle
p(X\,|\,U,x_{0},\theta)=\Pi_{t=0}^{T-1}p(x_{t+1}\,|\,x_{t},u_{t},\theta).$
(19a) For any time $t$, $\displaystyle
p(x_{t+1}\,|\,x_{t},u_{t},\theta)=p(x_{t+1}\,|\,x_{t},u_{t},w_{t},\theta)p(w_{t}),$
(19b) where $p(x_{t+1}\,|\,x_{t},u_{t},w_{t},\theta)$ is deterministic given
the dynamics, $f(\cdot)$, in (1) such that sampling from
$p(x_{t+1}\,|\,x_{t},u_{t},\theta)$ is equivalent to sampling from
$f(x_{t},u_{t},w_{t};\theta)$, $w_{t}\sim p(w)$. Extend this recursively by
expressing $\displaystyle
p(x_{t+1}\,|\,x_{t},u_{t},w_{t},\theta)=\Pi_{i=0}^{t}p(x_{i+1}\,|\,x_{i},u_{i},\theta)p(w_{i})$
(19c) such that $x_{t+1}$ can be expressed as
$f\big{(}...f(x_{0},u_{0},w_{0};\theta),u_{t},w_{t};\theta\big{)}$ where the
first argument is taken to mean that the state is recursively defined given an
initial state $x_{0}$, input sequence, and $w_{0},...,w_{t}\sim p(w)$. Define
this recursion as in (18) by $X=g(W;U,x_{0},\theta)$ where $W\sim p(W)$ to get
the desired result.

∎

Lemma 1: The linear system (6) has a change of variable $G(W;U,x_{0},\theta)$,
such that the $t$th column of $X$ in (18) is given by

$\displaystyle x_{t}=A^{t}x_{0}+\sum_{l=0}^{t-1}A^{t-1-l}(Bu_{l}+w_{l}),$
(20a) $w_{l}\sim p(W)$, the process noise distribution.

###### Proof:

See Appendix VII-A. ∎

Lemma 2: The change of variable (20a) and LQR experiment criteria (7), with
$K_{CE}(\mathcal{D})$, yields a valid form of (17):

$\displaystyle\hat{\nabla}_{U}=\frac{1}{L}\sum_{l=1}^{L}\nabla_{U}J\Big{(}K\big{(}\mathcal{D}\rvert_{X=G(W^{(l)};U,x_{0},\theta^{(l)})}\big{)};\theta^{(l)}\Big{)},$
(20b)

with $W^{(l)}\sim p(W).$

###### Proof:

See Appendix VII-B for satisfying the differentiability of the value function.
∎

###  IV-A Algorithm for Experiment Design Problem

In the LQR setting we derived a pathwise estimator (20). For more general
problems as in (4), if we assume the exchange of integral and gradient is
valid, we can express the pathwise gradient estimate as

$\displaystyle\hat{\nabla}_{U}=\frac{1}{L}\sum_{l=1}^{L}\nabla_{U}J\Big{(}K\big{(}\mathcal{D}\rvert_{X=G(W^{(l)};U,x_{0},\theta^{(l)})}\big{)};\theta^{(l)}\Big{)}.$
(21)

For each sample, $l$, we obtain a single experiment realization under sampled
noise $W$ for a sampled system, $\theta$. For this realization, we compute an
_a-posteriori_ system estimate, $\hat{\theta}$, and compute the control,
$\pi$, as in (3). Computing the gradient $\nabla_{U}J$ can be done
analytically in some cases. Given the structure of $J$, it may be easier to
use automatic differentiation as we go on to do in Section V.

Algorithm 1 Control-Oriented Experiment Design

Input $p_{\theta}\text{ (prior on }\theta\text{)},U_{0}\leavevmode\nobreak\
(\text{initialization}),L\leavevmode\nobreak\ \text{(batch size)},$

$\mathcal{U}\leavevmode\nobreak\ \text{(feasible
set)},x_{0}\leavevmode\nobreak\ \text{(initial
condition)},p_{W}\leavevmode\nobreak\ \text{(noise dist.)}$

Output $U^{*}$

function ConstructControl($\mathcal{D},p_{\theta}$)

$\hat{\theta}\leftarrow$ estimate system given $\mathcal{D},p_{\theta}$

$K\leftarrow$ Solve control problem given $\hat{\theta}$

end function

function SampleGrad($U,\theta,W,p_{\theta},x_{0}$)

$X\leftarrow G(W;U,x_{0},\theta)$

$\mathcal{D}\leftarrow X,U$

$K\leftarrow\textsc{ConstructControl}(\mathcal{D},p_{\theta})$

$\nabla_{U}J\leftarrow$ Compute gradient of $J(K;\theta)$

end function

while not converged do

for i=1 to L do

$\theta\sim p_{\theta},W\sim p_{W}$

$\nabla_{U}J_{i}\leftarrow\textsc{SampleGrad}(U_{j},\theta,W,p_{\theta},x_{0})$

end for

$U_{j+1}\leftarrow\text{Proj}_{\mathcal{U}}\big{(}U_{j}-\eta_{j}\frac{1}{L}\sum_{l=1}^{L}\nabla_{U}J_{l}\big{)}$

end while

##  V Numerical Experiments

###  V-A Car String

Figure 1: We focus on a scenario in which cars need to regulate to a specified
gap $L$ at a desired reference velocity $v$. Deviations from the desired
position are $\Delta w^{(n)}$ and deviations from the reference velocity are
$\Delta v^{(n)}$.

We consider the problem of maintaining a fixed distance, $\bar{L}$, between
$n$ cars at a desired velocity $v$ as depicted in Figure 1. We adapt the
continuous-time dynamics for _relative_ position as given in [18] to discrete-
time dynamics with sampling time $T_{s}$:

$\displaystyle\Delta v_{t+1}^{(n)}$
$\displaystyle=\bigg{(}-\frac{\alpha^{(n)}T_{s}}{m^{(n)}}+1\bigg{)}\Delta
v_{t}^{(n)}+\frac{T_{s}}{m^{(n)}}\Delta u_{t}^{(n)},$ (22)
$\displaystyle\Delta w_{t+1}^{(n)}$ $\displaystyle=T_{s}(\Delta
v_{t}^{(n)}-\Delta v_{t}^{(n+1)})+\Delta w_{k}^{(n)},$ (23)

where $\Delta v^{(n)}$ is the deviation from the reference velocity at car $n$
and $\Delta w^{(n)}$ is the deviation of the gap between cars $n+1$ and $n$
from the desired gap $L$. $\Delta u$ is a change in force input for each car.
This leads to an $n$ car state-vector $x_{t+1}:=[\Delta v_{t+1}^{(1)},\Delta
w_{t+1}^{(1)},\Delta v_{t+1}^{(2)},...,\Delta v_{t+1}^{(n)}]^{T}$. As such $A$
and $B$ are given by:

$\displaystyle A$
$\displaystyle=\begin{bmatrix}-\frac{\alpha^{(1)}T_{s}}{m^{(1)}}+1&0&0&0&\ldots\\\
T_{s}&1&-T_{s}&0&\ldots\\\
0&0&-\frac{\alpha^{(2)}T_{s}}{m^{(2)}}+1&0&\ldots\\\ 0&0&T_{s}&1&\ldots\\\
\vdots&\vdots&\vdots&\vdots&\ddots\end{bmatrix}$ (24) $\displaystyle B$
$\displaystyle=\begin{bmatrix}\frac{T_{s}}{m^{(1)}}&0&\ldots\\\ 0&0&\ldots\\\
0&\frac{T_{s}}{m^{(2)}}&\ldots\\\ 0&0&\ldots\\\
\vdots&\vdots&\ddots\end{bmatrix}.$ (25)

While there is a specific structure to $(A,B)$ here, we assume we do not know
the structure and estimate all $(2n-1)(3n-1)$ entries. We specify the noise
covariance in the dynamics (6) as $\Sigma_{w}=1\mathrm{e}{-2}\times I_{5}$.
The prior on the parameters (9) is $\Theta_{0}=[A,B]$ with
$m^{(1)}=m^{(2)}=m^{(3)}=1$, $\alpha^{(1)}=\alpha^{(2)}=\alpha^{(3)}=1$,
$T_{s}=0.1$; $\lambda_{0}^{-1}=1\mathrm{e}{-2}\times I_{8}$. Horizon $N=30$.

###  V-B Experiment Design Setup

In the results that follow, we use an experiment horizon of $M=20$ time steps,
and batch size $L=1000$ in Algorithm 1. As in [18], for the criteria in (7)
$Q$ includes penalties of magnitude $10$ on the positions $\Delta w$ and zero
on the velocity $\Delta v$. $R$ is the identity matrix. The weight matrix $S$
has parameters $\alpha_{1}=10^{3},\alpha_{2}=10^{6}$. $U$ is initialized with
$u_{t}\sim U[10^{-}3,10^{-}2]$ and is fixed across experiments.

To solve the A-optimal experiment design in (5), we use a pathwise gradient
estimator of the same form as (17) to perform stochastic gradient descent,
where the objective is to minimize the trace of the posterior error covariance
(9).

###  V-C Results and Discussion

We start by comparing the performance of our method against A-optimal design
(5) in terms of post-experiment LQR control performance (7). For the
experiment design (4), we consider a feasible input set $\mathcal{U}$ where
the Frobenius norm of the input sequence, $U$, is bounded by a design
parameter $\beta$ such that
$\mathcal{U}=\left\\{U\,|\,\|U\|_{F}\leq\beta\right\\}$. We vary the allowed
magnitude, $\beta$, in Figure 2 and observe our method outperforms the
A-optimal design uniformly. For any experiment design, we have a lower and
upper bound on the performance. If we knew the values of $(A,B)$, we would
achieve the lowest possible control cost such that this is a lower bound on
achievable performance. The upper bound comes from the expected control
performance associated with using a controller that uses the _a-priori_ system
estimate instead of the _a-posteriori_ as in (3). Intuitively, as the budget
$\beta$ for an experiment increases, so should the experiment performance as
the experiment can “explore more”. For small $\beta$ the methods are nearly
comparable as the additional information in minimal, but as $\beta$ increases
the performance of our method approaches that of the perfect knowledge case.
The A-optimal design only slowly decreases even though in the limit it should
reach the lower bound.

Figure 2: We compare the performance of our control-oriented system
identification against A-optimal experiment design for a system with five
states and three inputs, and known initial condition
$x_{0}=[0.,-4.3,0.,2.1,2.5]^{T}$ as in [18]. Any experiment design performance
is lower-bounded by the optimal control given knowledge of $(A,B)$ and upper
bounded by the performance of the controller given the _a-priori_ system
estimate. We include 95$\%$ confidence intervals using $10^{5}$ samples.

In Figure 3, we observe the optimal experiment inputs ($\Delta u$) for a
three-car system. The experiment inputs for our control-oriented approach
exhibit more excitation over the time horizon than the A-optimal design, which
has a relatively smooth experiment input sequence, suggesting our method would
perform better, which is verified by Figure 2. We also consider what the input
sequence would look like if the controller given the _a-priori_ system
estimate is used during an experiment trial. Since the controller is closed-
loop (11), we show the average input sequence. Using the optimal control may
seem like a natural way to conduct an experiment, but the norm bound $\beta$
on the input is not active in this case, indicating that simply performing the
optimal control leaves experiment budget unused.

Figure 4 shows how the problem scales with the system dimension. In the first
subplot we see the convergence of the experiment criteria in (4) as a function
of iterations. The criteria is normalized by the lower bound (given by the
performance if we knew $A,B$) for illustrative purposes. The number of
iterations until the criteria stabilizes is roughly constant across problem
dimension suggesting that the number of iterations required is independent of
the system size in this case. In the second subplot, the time to compute each
gradient sample is shown as a function of the state dimension. The compute
time is dominated by “control solution”–the time to solve the LQR problem–thus
suggesting that the time-complexity is dominated by the control problem.
“Overhead” refers to the remainder of the compute tasks such as automatic
differentiation and initialization of objects in the python library JAX [19].

Figure 3: The experiment input sequences for each car in a 3 car system are
compared under the proposed method, A-optimal experiment design, and by
conducting the experiment with a feedback controller given the prior. Figure
4: In the upper subplot, the number of iterations to converge for the car
string problem is essentially the same regardless of system size suggesting
good scaling properties of our method. In the lower subplot, we observe the
average time to compute a single gradient sample. The time is dominated by
solving the control problem and “overhead” refers to tasks such as automatic
differentiation, initial compile time, etc.

##  VI Conclusion

We proposed a control-oriented identification approach that in expectation
improves any data-driven controller by construction. Our solution method via
stochastic gradient descent is shown in the LQR setting to provide solutions
that outperform a typical experiment design objective in the sense that the
post-experiment control performance is better with our method.

Our experiment design approach extends beyond the LQR setting, and it would be
interesting to apply this to more general parametric problems. Establishing
analytical results on the convergence rate and sample complexity of the
stochastic gradient descent is an important direction for future research.

## References

  * [1] H. Chernoff, “Locally Optimal Designs for Estimating Parameters,” _The Annals of Mathematical Statistics_ , vol. 24, no. 4, pp. 586–602, Dec. 1953, publisher: Institute of Mathematical Statistics.
  * [2] H. A. Simon, “Dynamic Programming Under Uncertainty with a Quadratic Criterion Function,” _Econometrica_ , vol. 24, no. 1, pp. 74–81, 1956, publisher: [Wiley, Econometric Society].
  * [3] G. Elfving, “Optimum Allocation in Linear Regression Theory,” _The Annals of Mathematical Statistics_ , vol. 23, no. 2, pp. 255–262, Jun. 1952, publisher: Institute of Mathematical Statistics.
  * [4] K. Lindqvist and H. Hjalmarsson, “Identification for control: adaptive input design using convex optimization,” in _Proceedings of the 40th IEEE Conference on Decision and Control (Cat. No.01CH37228)_ , vol. 5. Orlando, FL, USA: IEEE, 2001, pp. 4326–4331.
  * [5] M. Gevers, “Identification for Control: From the Early Achievements to the Revival of Experiment Design,” in _Proceedings of the 44th IEEE Conference on Decision and Control_ , Dec. 2005, pp. 12–12, iSSN: 0191-2216.
  * [6] S. Anderson, K. Byl, and J. P. Hespanha, “Experiment design with gaussian process regression with applications to chance-constrained control,” in _2023 62nd IEEE Conference on Decision and Control (CDC)_ , 2023, pp. 3931–3938.
  * [7] B. D. Lee, I. Ziemann, A. Tsiamis, H. Sandberg, and N. Matni, “The Fundamental Limitations of Learning Linear-Quadratic Regulators,” Mar. 2023, arXiv:2303.15637 [cs, eess].
  * [8] A. J. Wagenmaker, M. Simchowitz, and K. Jamieson, “Task-Optimal Exploration in Linear Dynamical Systems,” in _Proceedings of the 38th International Conference on Machine Learning_. PMLR, Jul. 2021, pp. 10 641–10 652, iSSN: 2640-3498.
  * [9] S. Dean, H. Mania, N. Matni, B. Recht, and S. Tu, “On the Sample Complexity of the Linear Quadratic Regulator,” _Foundations of Computational Mathematics_ , vol. 20, no. 4, pp. 633–679, Aug. 2020.
  * [10] K. J. Åström and B. Wittenmark, _Adaptive control_. Courier Corporation, 2008.
  * [11] M. Simchowitz and D. Foster, “Naive Exploration is Optimal for Online LQR,” in _Proceedings of the 37th International Conference on Machine Learning_. PMLR, Nov. 2020, pp. 8937–8948, iSSN: 2640-3498.
  * [12] S. Mohamed, M. Rosca, M. Figurnov, and A. Mnih, “Monte Carlo Gradient Estimation in Machine Learning,” Sep. 2020.
  * [13] J. Peters and S. Schaal, “Reinforcement learning of motor skills with policy gradients,” _Neural Networks_ , vol. 21, no. 4, pp. 682–697, May 2008.
  * [14] P. Rossi, G. Allenby, and R. McCulloch, _Bayesian Statistics and Marketing_. Wiley, Oct. 2006.
  * [15] D. Bertsekas, _Dynamic programming and optimal control: Volume I_. Athena scientific, 2012, vol. 4.
  * [16] P. Glasserman, _Monte Carlo methods in financial engineering_. Springer, 2004, vol. 53.
  * [17] G. Grimmett and D. Stirzaker, _Probability and random processes_. Oxford university press, 2020.
  * [18] W. Levine and M. Athans, “On the optimal error regulation of a string of moving vehicles,” _IEEE Transactions on Automatic Control_ , vol. 11, no. 3, pp. 355–361, Jul. 1966, conference Name: IEEE Transactions on Automatic Control.
  * [19] “JAX: Autograd and XLA,” Feb. 2023, original-date: 2018-10-25T21:25:02Z.

##  VII Appendix

###  VII-A Change of variable for linear system

We want to find a change of variable for the linear system (6). We start by
showing the case for $t=2$ such that

$\displaystyle x_{1}$ $\displaystyle=Ax_{0}+Bu_{0}+w_{0}$ (26a) $\displaystyle
x_{2}$ $\displaystyle=Ax_{1}+Bu_{1}+w_{1}$ (26b) and expressing $x_{1}$ in
terms of $x_{0}$ $\displaystyle=A^{2}x_{0}+ABu_{0}+Aw_{0}+Bu_{1}+w_{1}$ (26c)
Then assuming this holds for time $t$: $\displaystyle
x_{t}=A^{t}x_{0}+\sum_{l=0}^{t-1}A^{t-1-l}(Bu_{l}+w_{l}),$ (26d) at time $t+1$
$\displaystyle x_{t+1}$ $\displaystyle=Ax_{t}+Bu_{t}+w_{t}$ (26e) substituting
the recursion (26d) $\displaystyle x_{t+1}$
$\displaystyle=A(A^{t}x_{0}+\sum_{l=0}^{t-1}A^{t-1-l}(Bu_{l}+w_{l}))+B_{t}+w_{t},$
(26f) $\displaystyle x_{t+1}$
$\displaystyle=A^{t+1}x_{0}+\sum_{l=0}^{t}A^{t-l}(Bu_{l}+w_{l}),$ (26g)

the desired result, where $w_{t}$ is distributed according to the process
noise and is independent of $u_{t}$.

###  VII-B Differentiability of the value function

Here we show the value function is differentiable with respect to the decision
variable $U$. In particular, with $\pi:=\\{K_{0},...,K_{N-1}\\}$ we write the
gradient as

$\displaystyle\partialderivative{J}{u_{ij}}=\sum\partialderivative{J}{K_{t,k,l}}\partialderivative{K_{t,k,l}}{\hat{{\Theta}}_{m,n}}\partialderivative{\hat{\Theta}_{m,n}}{u_{ij}}$
(27)

####  VII-B1 Gradient of $J$ with respect to $K_{t}$

First, we observe that the gradient of (12a) at $t=0$ is given by

$\displaystyle\partialderivative{J}{K_{0}}$
$\displaystyle=2\big{(}(R+B^{T}P_{1}B)K_{0}+B^{T}P_{1}A\big{)}x_{0}x_{0}^{T}$
(28a) and for $t>0$ as
$\displaystyle\begin{split}&\partialderivative{J}{K_{t}}=2\big{[}(R+B^{T}P_{t+1}B)K_{t}+B^{T}P_{t+1}A\big{]}\\\
&\times\bigg{(}\Pi_{i=0}^{t-1}(A+BK_{i})x_{0}x_{0}^{T}\Pi_{i=t-1}^{0}(A+BK_{i})^{T}+\Sigma_{w}\\\
&+\sum_{j=1,t>1}^{t-1}\Pi_{i=j}^{t-1}(A+BK_{i})\Sigma_{w}\Pi_{i=t-1}^{j}(A+BK_{i})^{T}\bigg{)}\end{split}$
(28b)

For a finite horizon, the entries of $P_{t}$ are finite even if the cost grows
exponentially in time such that the gradient itself will be finite.
Furthermore, the the gradient is polynomial in the Gaussian random variables
$(A,B)$ such that there exists a polynomial function of the random variables,
which is integrable.

####  VII-B2 Derivation: Gradient of $J$ with respect to $K_{t}$

We want to take the gradient of the value function

$\displaystyle
J(\Theta,\pi):=x_{0}^{T}P_{0}x_{0}+\sum_{t=0}^{N-1}\tr(P_{t+1}\Sigma_{W})$
(29a) with respect to $K_{t}$. For $K_{0}$ we expand $P_{0}$ to see the
dependence
$\displaystyle\begin{split}&\partialderivative{J(\Theta,\pi)}{K_{0}}=\partialderivative{K_{0}}\bigg{(}x_{0}^{T}(Q+K_{0}^{T}RK_{0}+\\\
&(A+BK_{0})^{T}P_{1}(A+BK_{0})x_{0}+\sum_{t=0}^{N-1}\tr(P_{t+1}\Sigma_{W})\bigg{)}.\end{split}$
(29b) Evaluting this, we get
$\displaystyle\partialderivative{J(\Theta,\pi)}{K_{0}}=(2RK_{0}+2B^{T}P_{1}BK_{0})x_{0}x_{0}^{T}+2B^{T}P_{1}Ax_{0}x_{0}^{T},$
(29c) which can be rearranged to give the desired result. For $t>0$, there is
dependence in both the initial condition and the process noise term. For the
initial condition term, recursively expand $P_{i}$ until $i=t$, and then take
the gradient as for $K_{0}$. If we define the current state as
$x_{t}:=\Pi_{i=0}^{t-1}(A+BK_{i})x_{0}$, then we can express this relationship
as
$\displaystyle\partialderivative{(x_{0}^{T}P_{0}x_{0})}{K_{t}}=2\big{[}(R+B^{T}P_{t+1}B)K_{t}+B^{T}P_{t+1}A\big{]}x_{t}x_{t}^{T}.$
(29d) This gives us the first part of the gradient. The second part is due to
the process noise and follows a similar pattern. Start by expanding $P_{t+1}$
to get terms of $K_{t+1}$:
$\displaystyle\begin{split}&\tr(P_{t}\Sigma_{w})=\\\
&\tr\big((Q+K_{t}^{T}RK_{t}+(A+BK_{t})^{T}P_{t+1}(A+BK_{t}))\Sigma_{w}\big{missing})\end{split}$
(29e) Expanding $P_{t+1}$, we need to take gradients of the following terms
(here given at $t$): $\displaystyle\partialderivative{K_{t}}\tr(K_{t}\Sigma
K_{t}^{T}R)=2RK_{t}\Sigma_{w}.$ (29f)
$\displaystyle\partialderivative{K_{t}}\tr(K_{t}\Sigma_{w}K_{t}^{T}B^{T}P_{t+1}B)=2B^{T}P_{t+1}BK_{t}\Sigma_{w}.$
(29g)
$\displaystyle\partialderivative{K_{t}}2\tr(\Sigma_{w}A^{T}P_{t+1}BK_{t})=2B^{T}P_{t+1}A\Sigma_{w}.$
(29h)

Using these gradients and algebraic manipulations, we get the desired result
for one step for the process noise term. This can be repeated for all time
steps to get the overall result. Combining the initial condition terms with
the noise terms gives us the gradient for $t>0$.

####  VII-B3 Gradient of $K_{t}$ with respect to estimate $\Theta$

Next, we examine the gradient of the data-driven control with respect to the
estimated system, denoted above as $(\hat{A},\hat{B})$ as we use certainty
equivalence in the dynamics parameters. The gradient of

$\displaystyle
K_{t}=-(R+\hat{B}^{T}P_{t+1}\hat{B})^{-1}\hat{B}^{T}P_{t+1}\hat{A}$ (30a) with
respect to $\hat{\Theta}$ is most easily written in terms of the elements of
$\hat{A},\hat{B}$.

The gradient is recursively computed as

$\displaystyle\begin{split}\partialderivative{K_{t}}{A_{ij}}&=(R+B^{T}P_{t+1}B)^{-1}B^{T}\partialderivative{P_{t+1}}{A_{ij}}B(R+B^{T}P_{t+1}B)^{-1}\\\
&-(R+B^{T}P_{t+1}B)^{-1}\bigg{(}B^{T}\partialderivative{P_{t+1}}{A_{ij}}A+B^{T}P_{t+1}e_{ij}\bigg{)}\end{split}$
(30b)
$\displaystyle\begin{split}\partialderivative{K_{t}}{B_{ij}}&=(R+B^{T}P_{t+1}B)^{-1}\bigg{(}2e_{ij}^{T}P_{t+1}B+B^{T}\partialderivative{P_{t+1}}{B_{ij}}B\bigg{)}\\\
&\times(R+B^{T}P_{t+1}B)^{-1}\\\
&-(R+B^{T}P_{t+1}B)^{-1}\bigg{(}B^{T}\partialderivative{P_{t+1}}{B_{ij}}B+2e_{ij}^{T}P_{t+1}B\bigg{)}\end{split}$
(30c)

with

$\displaystyle\begin{split}\partialderivative{P_{t}}{A_{ij}}&=2\partialderivative{K_{t}}{A_{ij}}^{T}RK_{t}+A^{T}\partialderivative{P_{t+1}}{A_{ij}}A\\\
&+2e_{ij}^{T}P_{t+1}A+2e_{ij}^{T}P_{t+1}BK_{t}\\\
&+2A^{T}\partialderivative{P_{t+1}}{A_{ij}}BK_{t}+2A^{T}P_{t+1}B\partialderivative{K_{t}}{A_{ij}}\\\
&+2\partialderivative{K_{t}}{A_{ij}}^{T}B^{T}P_{t+1}BK_{t}+K_{t}^{T}B^{T}\partialderivative{P_{t+1}}{A_{ij}}BK_{t}\end{split}$
(30d)
$\displaystyle\begin{split}\partialderivative{P_{t}}{B_{ij}}&=2\partialderivative{K_{t}}{B_{ij}}^{T}RK_{t}+A^{T}\partialderivative{P_{t+1}}{B_{ij}}A+2A^{T}P_{t+1}e_{ij}K_{t}\\\
&+2A^{T}\partialderivative{P_{t+1}}{B_{ij}}BK_{t}+2A^{T}P_{t+1}B\partialderivative{K_{t}}{B_{ij}}\\\
&+2\partialderivative{K_{t}}{B_{ij}}^{T}B^{T}P_{t+1}BK_{t}+K_{t}^{T}B^{T}\partialderivative{P_{t+1}}{B_{ij}}BK_{t}\\\
&+2K_{t}^{T}e_{ij}^{T}P_{t+1}BK_{t}.\end{split}$ (30e)

with $P_{N}=Q_{N}$. Again, we appeal to the fact that the elements of $P_{t}$
as governed by (11) will be finite for a finite horizon, leading to finite
values for the gradients. Furthermore, the gradient is again polynomial in the
parameters such that there exists a polynomial function that upper bounds the
gradient and is integrable.

####  VII-B4 Derivation: Gradient of $K_{t}$ with respect to estimate $\Theta$

For a posterior distribution with mean $\hat{\Theta}=[\hat{A},\hat{B}]$, and
the controller defined by the Riccati difference equations:

$\displaystyle K_{t}$
$\displaystyle=-(R+\hat{B}^{T}P_{t+1}\hat{B})^{-1}\hat{B}^{T}P_{t+1}\hat{A},$
(31a) $\displaystyle P_{t}$
$\displaystyle=Q+K_{t}^{T}RK_{t}-(\hat{A}+\hat{B}K_{t})^{T}P_{t+1}(\hat{A}+\hat{B}K_{t}),$
(31b)

we want to find the gradient with respect to elements of $\hat{A}$ and
$\hat{B}$. Starting with $K_{t}$:

$\displaystyle\partialderivative{K_{t}}{\hat{A}_{ij}}$
$\displaystyle=\partialderivative{\hat{A}_{ij}}\bigg{(}-(R+\hat{B}^{T}P_{t+1}\hat{B})^{-1}\hat{B}^{T}P_{t+1}\hat{A}\bigg{)}$
(32a)
$\displaystyle\begin{split}&=\partialderivative{\hat{A}_{ij}}\big{(}-(R+\hat{B}^{T}P_{t+1}\hat{B})^{-1}\big{)}\hat{B}^{T}P_{t+1}\hat{A}\\\
&+\big{(}-(R+\hat{B}^{T}P_{t+1}\hat{B})^{-1}\big{)}\partialderivative{\hat{A}_{ij}}\big{(}\hat{B}^{T}P_{t+1}\hat{A}\big{)}\end{split}$
(32b) For the gradient of the first component:
$\displaystyle\partialderivative{\hat{A}_{ij}}\big{(}-(R+\hat{B}^{T}P_{t+1}\hat{B})^{-1}\big{)}=$
(32c)
$\displaystyle(R+\hat{B}^{T}P_{t+1}\hat{B})^{-1}\partialderivative{(R+\hat{B}^{T}P_{t+1}\hat{B})}{\hat{A}_{ij}}\big{(}(R+\hat{B}^{T}P_{t+1}\hat{B})^{-1}\big{)}$
(32d)
$\displaystyle=(R+\hat{B}^{T}P_{t+1}\hat{B})^{-1}\hat{B}^{T}\partialderivative{P_{t+1}}{\hat{A}_{ij}}\hat{B}\big{(}(R+\hat{B}^{T}P_{t+1}\hat{B})^{-1}\big{)}$
(32e) and the second
$\displaystyle\partialderivative{\hat{A}_{ij}}\big{(}\hat{B}^{T}P_{t+1}\hat{A}\big{)}=\hat{B}^{T}\bigg{(}\partialderivative{P_{t+1}}{\hat{A}_{ij}}\hat{A}+P_{t+1}e_{ij}\bigg{)}.$
(32f) The results for the partial with respect to $B_{ij}$ follows similarly.
In each case, we need to compute the partial of $P_{t}$:
$\displaystyle\begin{split}&\partialderivative{P_{t}}{\hat{A}_{ij}}=\\\
&\partialderivative{\hat{A}_{ij}}\big{(}Q+K_{t}^{T}RK_{t}-(\hat{A}+\hat{B}K_{t})^{T}P_{t+1}(\hat{A}+\hat{B}K_{t})\big{)}.\end{split}$
(32g) $Q$ is independent of $\hat{A}$ (and $\hat{B}$). The rest of terms are:
$\displaystyle\partialderivative{\hat{A}_{ij}}\big{(}K_{t}^{T}RK_{t}\big{)}=2\partialderivative{K_{t}}{\hat{A}_{ij}}^{T}RK_{t},$
(32h)
$\displaystyle\partialderivative{\hat{A}_{ij}}\big{(}(\hat{A}+\hat{B}K_{t})^{T}P_{t+1}(\hat{A}+\hat{B}K_{t})\big{)}=$
(32i) $\displaystyle
2e_{ij}^{T}P_{t+1}\hat{A}+\hat{A}^{T}\partialderivative{P_{t+1}}{\hat{A}_{ij}}\hat{A}+\partialderivative{\hat{A}_{ij}}\big{(}\hat{A}^{T}P_{t+1}\hat{B}K_{t}\big{)},$
(32j) and
$\displaystyle\begin{split}&\partialderivative{\hat{A}_{ij}}\big{(}\hat{A}^{T}P_{t+1}\hat{B}K_{t}\big{)}=2e_{ij}^{T}P_{T+1}\hat{B}K_{t}\\\
&+2\hat{A}^{T}\partialderivative{P_{t+1}}{\hat{A}_{ij}}\hat{B}K_{t}+2\hat{A}^{T}P_{t+1}\hat{B}\partialderivative{K_{t}}{\hat{A}_{ij}},\end{split}$
(32k) A similar derivation follows with respect to $\hat{B}_{ij}$.

####  VII-B5 Gradient of estimate $\Theta$ with respect to $U$ with derivation

Finally, to address the gradient of the estimated value with respect to the
design variable $U\in\mathbb{R}^{n_{u}\times T}$ with entries $u_{ij}$, we
first rewrite the estimator in (9) using sums as

$\displaystyle\hat{\Theta}$
$\displaystyle=(\Theta_{0}\Lambda_{0}+\sum_{t=0}^{T-1}y_{t}s_{t}z_{t}^{T})(\Lambda_{0}+\sum_{t=0}^{T-1}z_{t}s_{t}z_{t}^{T})^{-1},$
(33a) $\displaystyle=:\Psi\Delta,$ (33b) where $\displaystyle y_{t}$
$\displaystyle=x_{t+1}=Ax_{t}+B\gamma_{t},$ (33c) $\displaystyle x_{t}$
$\displaystyle=A^{t}x_{0}+\sum_{l=0}^{t-1}A^{t-1-l}(Bu_{l}+w_{l}),$ (33d)
$\displaystyle z_{t}$ $\displaystyle=[x_{t};u_{t}].$ (33e)

$y_{t},s_{t},z_{t}$ all depend on $U$. We use the $\mathrm{vec}\left(\right)$
operator, which stacks the columns on tops of each other, to simplify the
derivation. As such,

$\displaystyle\begin{split}\mathrm{vec}\left(\partialderivative{\hat{\Theta}}{u_{ij}}\right)&=(\Delta^{T}\otimes
I)\mathrm{vec}\left(\partialderivative{\Psi}{u_{ij}}\right)+(I\otimes\Psi)\mathrm{vec}\left(\partialderivative{\Delta}{u_{ij}}\right),\end{split}$
(33f)

where

$\displaystyle\begin{split}&\mathrm{vec}\left(\partialderivative{\Psi}{u_{ij}}\right)=\\\
&\mathrm{vec}\left(\partialderivative{\sum_{t=0}^{T-1}y_{t}s_{t}z_{t}^{T}}{u_{ij}}\right)=\sum_{t=0}^{T-1}\mathrm{vec}\left(\partialderivative{y_{t}s_{t}z_{t}^{T}}{u_{ij}}\right),\\\
&=(z_{t}s_{t}^{T}\otimes
I)\mathrm{vec}\left(\partialderivative{y_{t}}{u_{ij}}\right)+(z_{t}\otimes
y_{t})\mathrm{vec}\left(\partialderivative{s_{t}}{u_{ij}}\right)\\\
&+(I\otimes
y_{t}s_{t})\mathrm{vec}\left(\partialderivative{z_{t}^{T}}{u_{ij}}\right).\end{split}$
(33g)

Each gradient in the above expression is

$\displaystyle\begin{split}\mathrm{vec}\left(\partialderivative{y_{t}}{u_{ij}}\right)&=\mathrm{vec}\left(\partialderivative{u_{ij}}A^{t+1}x_{0}+\sum_{l=0}^{t}A^{t-l}(B\gamma_{l}+w_{l})\right),\\\
&=\sum_{l=0}^{t}(I\otimes
A^{t-l}B)\mathrm{vec}\left(\partialderivative{\gamma_{l}}{u_{ij}}\right),\\\
&=(I\otimes A^{t-i}B)\mathrm{vec}\left(e_{ij}\right),\quad(t\geq
i)\end{split}$ (33h)
$\displaystyle\begin{split}&\mathrm{vec}\left(\partialderivative{s_{t}}{u_{ij}}\right)=\mathrm{vec}\left(\partialderivative{u_{ij}}\atan((\norm{x_{t}}-\alpha_{1})\alpha_{2})/\pi+0.5\right)\\\
&=\mathrm{vec}\left(\frac{1}{\pi\alpha_{2}^{2}}\frac{1}{1/\alpha_{2}^{2}+(\norm{x_{t}}-\alpha_{1})^{2}}\partialderivative{u_{ij}}\big{(}(\norm{x_{t}}-\alpha_{1})^{2}\big{)}\right)\\\
&\partialderivative{u_{ij}}\big{(}(\norm{x_{t}}-\alpha_{1})^{2}\big{)}=2(\norm{x_{t}}-\alpha_{1})\frac{x_{t}}{\norm{x_{t}}}\partialderivative{x_{t}}{u_{ij}}\end{split}$
(33i)
$\displaystyle\begin{split}&\mathrm{vec}\left(\partialderivative{z_{t}^{T}}{u_{ij}}\right)=\\\
&\mathrm{vec}\left(\partialderivative{u_{ij}}\begin{bmatrix}A^{t}x_{0}+\sum_{l=0}^{t-1}A^{t-1-l}(Bu_{l}+w_{l})\\\
u_{t}\end{bmatrix}^{T}\right)\\\ &=\begin{bmatrix}(I\otimes
A^{t-1-i}B)\mathrm{vec}\left(\partialderivative{u_{i}}{u_{ij}}\right),\quad(t\geq
i)\text{ else 0}\\\ \mathrm{vec}\left(e_{ij}\right),\quad(i=t),\text{ else
0}\end{bmatrix}^{T}\end{split}$ (33j)

Going to the second term, $\Delta$, in the estimator, we obtain

$\displaystyle\mathrm{vec}\left(\partialderivative{\Delta}{u_{ij}}\right)=-(\Delta\otimes\Delta)\mathrm{vec}\left(\partialderivative{u_{ij}}\sum_{t=0}^{T-1}z_{t}s_{t}z_{t}^{T}\right)$
(33k)
$\displaystyle=-(\Delta\otimes\Delta)\sum_{t=0}^{T-1}\Bigg{(}(z_{t}s_{t}^{T}\otimes
I)$ (33l) $\displaystyle\times\begin{bmatrix}(I\otimes
A^{t-1-i}B)\mathrm{vec}\left(e_{ij}\right),\leavevmode\nobreak\ (t\geq
j),\text{else 0}\\\ \mathrm{vec}\left(e_{ij}\right),\leavevmode\nobreak\
(j=t),\text{else 0}\end{bmatrix}$ (33m) $\displaystyle+(z_{t}\otimes
z_{t})\mathrm{vec}\left(\partialderivative{s_{t}}{\gamma_{i}}\right)+(I\otimes
Z_{T}s_{t})$ (33n) $\displaystyle\times\begin{bmatrix}(I\otimes
A^{t-1-i}B)\mathrm{vec}\left(e_{ij}\right),\leavevmode\nobreak\ (t\geq
j),\text{else 0}\\\ \mathrm{vec}\left(e_{ij}\right),\leavevmode\nobreak\
(i=t),\text{else 0}\end{bmatrix}^{T}\Bigg{)}$ (33o)

The gradient is then well-defined except if $\Delta$ were to be ill-defined
due to lack of invertibility; however, the prior $\Lambda_{0}$ is chosen to be
non-singular and obviates this possibility. The resulting expression contains
a rational and polynomial term, such that there exists a polynomial bounding
function. Overall then, for each term in the gradient, there exists an
integrable bounding function such that the dominated convergence can be
applied and the exchange of gradient (limit) and integral is valid.