# Performance Bounds for Neural Network Estimators: Applications in Fault
Detection

Navid Hashemi, Mahyar Fazlyab, Justin Ruths N. Hashemi and J. Ruths are with
the Department of Mechanical Engineering, The University of Texas at Dallas.
Email: (nxh150030<EMAIL_ADDRESS>M. Fazlyab is with the John Hopkins
Mathematical Institute for Data Sciences. Email<EMAIL_ADDRESS>
###### Abstract

We exploit recent results in quantifying the robustness of neural networks to
input variations to construct and tune a model-based anomaly detector, where
the data-driven estimator model is provided by an autoregressive neural
network. In tuning, we specifically provide upper bounds on the rate of false
alarms expected under normal operation. To accomplish this, we provide a
theory extension to allow for the propagation of multiple confidence
ellipsoids through a neural network. The ellipsoid that bounds the output of
the neural network under the input variation informs the sensitivity - and
thus the threshold tuning - of the detector. We demonstrate this approach on a
linear and nonlinear dynamical system.

##  I Introduction

The rise in interest in data-driven techniques is a response to the need for
better models of complex systems. Model-based fault and anomaly detection in
dynamical systems typically employs first-principle models (Newton’s and
Kirchhoff’s Laws, reaction kinetics, etc.) of systems to identify
discrepancies between the observed and predicted sensor measurements. However,
to use this general approach on large-scale operations, such as petro-chemical
refineries and power distribution grids, operators need scalable and efficient
techniques for creating models. Data-driven tools offer a compelling option
because models can be devised simply from past data, do not require intimate
knowledge of system parameters, and can be recomputed regularly to update the
dynamics that may change over time.

Artificial neural networks, in particular, provide relatively easy training
and generalization, simple architecture, good ability to approximate nonlinear
functions, and robust to inexact input data [1]. Neural networks can be used
to identify and control nonlinear dynamic systems because they can approximate
a wide range of nonlinear functions. Over the past two decades, fault
detection has employed neural networks in a variety of ways, including using
them as classifiers to categorize normal or various faulty behaviors [2] and
using them as a model of the system to provide an estimate for subsequent
fault detection [3, 4]. In the latter case, neural network estimators can be
implemented using an iterative approach (in which the past estimate(s) is part
of the input to produce the next estimate) or an autoregressive approach (in
which a finite history of past measurements is used to produce the next
estimate) [5]-[7]. The detection problem becomes particularly challenging when
the fault is unknown or cannot be modeled (e.g., we want to build a detector
that is sensitive to potentially unknown faults or anomalies, like attacks).
In this context, tuning of the detector is important to balance the
sensitivity of the detector with the prevalence of false alarms and while
data-driven techniques are not new to fault detection, tuning to-date has been
done empirically. In this paper we leverage recent results to quantify the
robustness of neural networks [8, 9] to develop a novel data-driven anomaly
detector with guarantees on the upper limit of the false alarm rate.

Quantifying the robustness of neural networks is originally motivated by their
vulnerability to adversarial attacks, i.e., carefully chosen small input
perturbations that can drastically change their output. To this end, a
plethora of tools have been developed to bound the output of neural networks
for a given range of inputs (e.g., the set of plausible attacks) [9]-[16]. Our
particular interest in this paper is the uncertainty propagation technique put
forth in [8], in which the authors develop a semi-definite program that
propagates an input ellipsoid (e.g., confidence region of a density function)
through the neural network to obtain _guaranteed_ ellipsoidal over-
approximation of the output set. In another context, Monte Carlo sampling, the
Unscented Transform (UT) and Extended Kalman Filtering (EKF) have been used to
take a set of samples from the input distribution, propagate them through the
neural network, and approximate the first and second moments of the output
distribution from them [17, 18]. These sampling methods scale to larger neural
networks but they lack formal guarantees.

##  II Background

Consider a discrete-time dynamical system whose state update can be described
by an unknown continuous function,
$\mathcal{F}:\mathbb{R}^{n}\to\mathbb{R}^{n}$,

$x_{k+1}=\mathcal{F}(x_{k}),$ (1)

which may arise as a freely evolving dynamical system or as a feedback control
system. Our observations (sensor measurements) of the system are linear
combinations of the states corrupted by additive zero-mean noise
$v_{k}\in\mathbb{R}^{p}$, with known covariance $\Sigma_{v}$,

$y_{k}=Hx_{k}+v_{k},$ (2)

where the sensor matrix $H\in\mathbb{R}^{p\times n}$ is known. Because the
system model is unknown, or possibly too complicated with too many parameters
to use in an effective way, we use a data-driven estimator with the assumption
that it is possible to construct an accurate nonlinear auto-regressive (NARX)
model of the output observations based on past measurement values (i.e., the
system is fully observable),

$y_{k+1}=f(\mathcal{Y}_{k,N}),$ (3)

where
$\mathcal{Y}_{k,N}=\textbf{vec}\left[\left\\{y_{k-i}\right\\}_{i=0}^{N}\right]\in\mathbb{R}^{p(N+1)}$.

Again, the mapping $f:\mathbb{R}^{n\times p}\to\mathbb{R}^{p}$ is unknown or
otherwise too complicated to form from first-principle models, so we
approximate this mapping through supervised training of a feedforward neural
network from time-series recordings of sensor measurements. At each time
$k\geq N$ the network is provided the vector of current and $N$ past
measurements $\mathcal{Y}_{k,N}$ as input training data and the next
measurement $y_{k+1}$ as the labeled output data of the neural network. The
network is trained under a large number of varying initial conditions.

The $\ell$-layer neural network, trained with input data $\mathcal{Y}_{k,N}$
and labeled data $y_{k+1}$, is

$\displaystyle\zeta_{k}$ $\displaystyle=\left[\ \mathcal{Y}_{k,N}\
\right]^{\top},$ (4) $\displaystyle z^{0}$ $\displaystyle=\zeta_{k},$
$\displaystyle z^{t+1}$ $\displaystyle=\phi(W^{t}z^{t}+b^{t}),\quad
t=0,1,\dots,\ell-1,$ $\displaystyle\hat{y}_{k+1}$
$\displaystyle=W^{\ell}z^{\ell}+b^{\ell},$

with ReLU activation function $\psi(s_{i})=\max{(0,s_{i})}$, $\forall
s_{i}\in\mathbb{R}$, and
$\phi([s_{1},s_{2},..,s_{d}])=[\psi(s_{1}),\psi(s_{2}),...,\psi(s_{d})]$. We
train the model to compute optimal weight matrices and bias vectors and we
call the output of the neural network, $\hat{y}_{k+1}$, as the
estimated/predicted output.

The primary role of the training process is to refine the neural network to
produce accurate estimates, i.e., predictions $\hat{y}_{k}$ that are close to
the actual measurements $y_{k}$. The purpose of this paper is not to refine
the training process to improve accuracy, but to quantify the robustness of
the prediction to input uncertainties for a given trained neural network
(e.g., uncertainty caused by sensor noise). Further, we aim to use the
robustness bounds as a way to quantify normal behavior from abnormal behavior.

###### Remark 1

The NARX architecture is one of several ways to construct a neural network
estimator of a dynamical system. An alternative choice is to feed the past
estimation $\hat{y}_{k}$ as an input to the prediction for $\hat{y}_{k+1}$.
The autoregressive approach we use here simplifies the training process and
will also simplify the robustness quantification since the bounds on the
neural network inputs are constant across time.

Figure 1: The proposed detector relies on the development of a prediction
ellipsoid (orange) $\mathcal{E}(\hat{\mu}_{k+1},\hat{\Sigma}_{k+1})$ that
bounds the variation of the estimate due to the noise inherent in the input
measurements $y_{k}$, $y_{k-1}$, …, $y_{k-N}$ and captured by the input
ellipsoids $\mathcal{E}(y_{k-i},\bar{\Sigma}_{v})$, $i=0,\dots,N$. The
framework to compute the prediction ellipsoid is presented in Section IV and
the detector that uses this ellipsoid to decide between normal and anomalous
behavior is presented in Section III.

##  III Neural Network-Based Anomaly Detector

In the same spirit as propagating the estimation covariance of a Luenberger
observer (e.g., Kalman Filter), evaluating the robustness of a neural network
estimator identifies the quality of the estimation - quantifying how much the
estimator can reduce the influence of the noise on the prediction. At the same
time, quantifying the inherent variation in the predictions due to uncertainty
(noise) also provides a bound on the variation observed during normal
behavior.

Consider the noisy measurement $y_{k}$ from (2) at a particular time step $k$,
which is composed of the non-noisy (deterministic) “ideal” measurement,

$y^{*}_{k}=Hx_{k},$ (5)

and a sample of the sensor noise, $v_{k}$. If a different realization of the
sensor noise $v_{k}$ was instead used, our goal is to bound how different the
prediction $\hat{y}_{k+1}$ would be under these two sensor realization
scenarios. Said differently, we would like to assess the robustness of the
prediction to the inherent perturbations caused by sensor noise.

Because the support of the noise could be large, or possibly infinite (e.g.,
Gaussian noise), it is practically useful to truncate the noise distribution
at a desired confidence set, such that the probability $\bar{p}$ that a noise
sample is drawn from within the confidence set is a desired value (typically
close to one). An ellipsoidal confidence set can be constructed using a scaled
version of the covariance of the noise distribution as the shape matrix and
characterized by

$\Pr[v_{k}^{\top}\Sigma_{v}^{-1}v_{k}\leq\alpha]=\Pr[v_{k}^{\top}(\underbrace{\alpha\Sigma_{v}}_{\bar{\Sigma}_{v}})^{-1}v_{k}\leq
1]=\bar{p}.$ (6)

This defines the $\bar{p}$-confidence ellipsoid on the sensor noise,
$\mathcal{E}(0,\bar{\Sigma}_{v})$, with an ellipsoid with center $\mu$ and
shape matrix $\Sigma$ defined as

$\mathcal{E}(\mu,\Sigma):=\left\\{\xi\ \big{|}\
(\xi-\mu)^{\top}\Sigma^{-1}(\xi-\mu)\leq 1\right\\},$ (7)

and the size of the ellipsoid to match the desired confidence level $\bar{p}$
is chosen by selecting $\alpha=2\Gamma^{-1}(\frac{p}{2},\bar{p})$, using the
inverse regularized lower incomplete gamma function and $p$ is the number of
measurements at each time step [19].

One way to interpret this bound is that since $y_{k}=y^{*}_{k}+v_{k}$, the
actual measurement $y_{k}$ is contained within an ellipsoid
$\mathcal{E}(y_{k}^{*},\bar{\Sigma}_{v})$, i.e., centered at the ideal
measurement $y^{*}_{k}$ (see Fig. 1a). However, we can also see that
$y^{*}_{k}=y_{k}-v_{k}$ which says the opposite - that the ideal measurement
is contained within the ellipsoid centered at the actual measurement,
$y^{*}_{k}\in\mathcal{E}(y_{k},\bar{\Sigma}_{v})$ (see Fig. 1b). The latter
interpretation is practically useful because we have access to the actual
measurements, but not the ideal (noise-free) measurements.

The primary contribution of this paper is to construct and tune a model-based
anomaly detector that employs a neural network to produce an estimate. Using
Fig. 1 as a guide, this objective can be decomposed into two pieces

  1. 1.

the framework to quantify the robustness of the neural network under input
perturbations ($N$ blue input confidence ellipsoids),
$\mathcal{E}(y_{k-N},\bar{\Sigma}_{v}),\cdots,\mathcal{E}(y_{k},\bar{\Sigma}_{v})$
to produce an ellipsoidal bound on the predictions (orange ellipsoid) which we
will denote by $\mathcal{E}(\hat{\mu}_{k+1},\hat{\Sigma}_{k+1})$; and

  2. 2.

the definition of a detector that uses this (orange) ellipsoid to bound normal
behavior through the use of the geometric sum.

We address the detector definition first in the following subsection before
presenting the mathematical framework to create the prediction bound in
Section IV.

###  III-A Detector Definition

Given the ellipsoidal confidence sets that quantify the perturbation of the
inputs to the neural network estimator, in Section IV we will compute an
ellipsoidal bound on the output of the estimate, which we call the prediction
bound $\mathcal{E}(\hat{\mu}_{k+1},\hat{\Sigma}_{k+1})$. This bound guarantees
that any $N+1$ past measurement vectors which stay within their respective
confidence sets are mapped by the neural network to a point within the
prediction bound ellipsoid.

From Fig. 1, since the ideal past measurements $y_{k-N}^{*}$, …, $y_{k}^{*}$
are within the input confidence ellipsoids, we know the neural network output
estimate $\hat{y}_{k+1}^{*}$ will be located within the prediction bound. The
accuracy of the estimator can be interpreted as the distance between this
prediction based on ideal measurements $\hat{y}_{k+1}^{*}$ and the next ideal
measurement $y_{k+1}^{*}$. In our framework presented here, we will assume
that the trained model is accurate such that this distance is relatively small
and can be neglected; however, an ellipsoidal bound on the estimation accuracy
could be easily integrated into our detector definition.

As discussed above, the actual noisy measurement at $k+1$ is, with probability
$\bar{p}$ within the confidence ellipsoid centered at the ideal measurement,
$y_{k+1}\in\mathcal{E}(y^{*}_{k+1},\bar{\Sigma}_{v})$. Although we do not know
the ideal measurement, we know it is within the prediction bound. This leads
us to suggest a detector that evaluates normal behavior as measurements based
on the geometric sum between the prediction bound and the confidence
ellipsoid.

###### Definition 1

A model-based anomaly detector that employs an autoregressive neural network
to produce an estimate from past measurements raises alarms using the
following logic:

$\left\\{\ \begin{aligned}
y_{k+1}\in\mathcal{E}(\hat{\mu}_{k+1},\hat{\Sigma}_{k+1})\oplus\mathcal{E}(0,\bar{\Sigma}_{v})\
&\to\ \ \text{no alarm},\\\ \text{otherwise}\hskip 105.2751pt&\to\ \
\text{alarm},\end{aligned}\right.$ (8)

where the prediction bound is characterized by a center $\hat{\mu}_{k+1}$ and
shape matrix $\hat{\Sigma}_{k+1}$; the noise $\bar{p}$-confidence ellipsoid is
characterized by the shape matrix $\bar{\Sigma}_{v}$; and $\oplus$ denotes the
geometric (Minkowski) sum of two sets
$\mathcal{S}_{1}\oplus\mathcal{S}_{2}=\\{s_{1}+s_{2}\ |\
s_{1}\in\mathcal{S}_{1},\ s_{2}\in\mathcal{S}_{2}\\}$.

###### Remark 2

Since the geometric sum of two ellipsoids is, in general, not an ellipsoid, it
is best to verify this inclusion directly rather than to first approximate the
sum as an ellipsoid. There are a number of computational techniques to verify
the inclusion of a point in the geometric sum of ellipsoidal sets such as
linear matrix inequalities [20] or geometric methods to compute the exact
geometric sum [21].

The purpose of a detector is to alert the operator of a system to potential
anomalies, such as faults or attacks. The number of false alarms, i.e., alarms
raised during normal operation, is a key way to assess the performance of a
detector. Tuning the detector to balance sensitivity and false alarms is a key
step towards making a detector usable in practice.

###### Proposition 1

The detector defined in Definition 1 has a false alarm rate that is upper
bounded by $1-\bar{p}^{N+2}$.

###### Proof:

The fact that the sensor noise is independent allows us to quantify
probabilities easily. The confidence ellipsoids are used $N+1$ times for times
$k-N$, …, $k-1$, $k$ (the input ellipsoids) and once for time $k+1$ in the
geometric sum. Thus since $v_{k}$ is independent, the probability that all
these $v_{k-N}$, …, $v_{k+1}$ fall within their confidence sets is
$\bar{p}^{N+2}$. Because the prediction ellipsoid is an outer bound and
because there could be noise realizations that lie outside their confidence
ellipsoids but still remain inside the prediction bound, this probability is a
lower bound. Thus under normal operation, we expect the probability of
generating a false alarm to be at most $1-\bar{p}^{N+2}$. ∎

##  IV Robustness of Prediction

$\displaystyle
E_{i}=\left[\begin{matrix}\begin{array}[]{c|c|c}0_{n_{i}\times\left(\sum_{j=1}^{i-1}n_{j}\right)}\
\ I_{n_{i}}\ \
0_{n_{i}\times\left(\sum_{j=i+1}^{q}n_{j}\right)}&0_{n_{i}\times\left(\sum_{t=2}^{\ell+1}N_{j}\right)}&0_{n_{i}\times
1}\\\ \hline\cr 0_{1\times
n_{\Gamma}}&0_{1\times\left(\sum_{t=2}^{\ell+1}N_{j}\right)}&1\end{array}\end{matrix}\right]$
($\star$)

Our detector design leverages the ability to compute a bound on the
predictions made by the neural network estimator under perturbations to the
input. To accomplish this bound, we use the formulation described in [8],
which provides an ellipsoidal bound on the output of a neural network given a
single ellipsoidally bounded input vector. Here, as seen in Fig. 1, our
approach uses multiple ellipsoidally bounded input vectors and, therefore,
requires an extension to the method in [8]. The fundamental difference in
these two approaches in the context of our problem is at what stage the
confidence level truncation is applied. In our proposed multi-ellipsoid input
approach we construct a confidence ellipsoid for each individual measurement.
To use the single ellipsoid input approach of [8] directly, we would define an
overall confidence ellipsoidal set on the entire stacked input vector. The
multiple ellipsoid approach has a few key benefits: (1) while the noise is
independent, the center of the input ellipsoids $y_{k-N}$, …, $y_{k}$ are
correlated due to the deterministic dynamics of the system, which complicates
concatenating the measurements; (2) the input ellipsoids all have the same
shape matrix whereas the single input ellipsoid shape matrix would be time
dependent; and (3) in Appendix A we show that the multiple ellipsoid input
approach results in a less conservative ellipsoidal prediction bound.

Here we follow the general approach in [8, 9], introducing the updates needed
for multiple ellipsoidally bounded inputs. Consider the following general
neural network with $\ell$ hidden layers,

$\displaystyle\Gamma_{k}$ $\displaystyle=[\gamma_{1,k}^{\top},\ \
\gamma_{2,k}^{\top},\ ...,{\gamma_{q,k}}^{\top}]^{\top}$ (9) $\displaystyle
z^{0}$ $\displaystyle=\Gamma_{k}$ $\displaystyle z^{t+1}$
$\displaystyle=\phi(W^{t}z^{t}+b^{t}),\quad t=0,1,...,\ell-1$
$\displaystyle\pi(z^{0})$ $\displaystyle=W^{\ell}z^{\ell}+b^{\ell},$

with $\gamma_{i,k}\in\mathbb{R}^{n_{i}}$,
$\Gamma_{k}\in\mathbb{R}^{n_{\Gamma}}$ ($\sum_{j=1}^{q}n_{i}=n_{\Gamma}$),
$\pi(z^{0})\in\mathbb{R}^{n_{\pi}}$ and $z^{t}\in\mathbb{R}^{N_{t}}$. By
concatenating all the post-activation values as $\mathbf{z}=\left[\
{z^{0}}^{\top},\
{z^{1}}^{\top},...,{z^{\ell}}^{\top}\right]^{\top}\in\mathbb{R}^{N_{\mathbf{z}}}$
($N_{\mathbf{z}}=\sum_{t=0}^{\ell}N_{t}$), and defining the “selector
matrices” $S^{t}$ such that $z^{t}=S^{t}\mathbf{z}$, we can rewrite the
network as

$\displaystyle\Gamma_{k}$ $\displaystyle=[\gamma_{1,k}^{\top},\ \
\gamma_{2,k}^{\top},\ ...,{\gamma_{q,k}}^{\top}]^{\top}$ (10) $\displaystyle
z^{0}$ $\displaystyle=S^{0}\mathbf{z}=\Gamma_{k}$ $\displaystyle B\mathbf{z}$
$\displaystyle=\phi(A\mathbf{z}+b)$ $\displaystyle\pi(z^{0})$
$\displaystyle=W^{\ell}S^{\ell}\mathbf{z}+b^{\ell},$

where

$\displaystyle A=\begin{bmatrix}W^{0}&0&\dots&0&0\\\ 0&W^{1}&\dots&0&0\\\
\vdots&\vdots&\ddots&\vdots&\vdots\\\
0&0&\dots&W^{\ell-1}&0\end{bmatrix},\quad b=\begin{bmatrix}b^{0}\\\ b^{1}\\\
\vdots\\\ b^{\ell-1}\end{bmatrix},$ $\displaystyle
B=\begin{bmatrix}0&I_{n_{1}}&\dots&0\\\ \vdots&\vdots&\ddots&\vdots\\\
0&0&\dots&I_{n_{\ell}}\end{bmatrix}.$ (11)

Suppose each input $\gamma_{i,k}$ takes values inside the ellipsoid
$\mathcal{E}(\mu_{\gamma_{i,k}},\Sigma_{i,k})$. Our goal is to bound the
resulting output $\pi(z^{0})$ of the neural network by an ellipsoid. To this
end, we first abstract the neural network via Quadratic Constraints (QC) [9]
and then use the S-procedure to propagate the $q$ input ellipsoids through the
network. We begin with the following definition.

######  Definition 2 (QC for functions [22])

Let $\phi:\mathbb{R}^{d}\to\mathbb{R}^{d}$ and suppose
$\mathcal{Q}_{\phi}\subset\mathbb{S}^{2d+1}$ is the set of all symmetric
indefinite matrices $Q$ such that

$\begin{bmatrix}x\\\ \phi(x)\\\ 1\end{bmatrix}^{\top}Q\begin{bmatrix}x\\\
\phi(x)\\\ 1\end{bmatrix}\geq 0,\quad\forall\ x\in\mathcal{X},$ (12)

where $\mathcal{X}\subset\mathbb{R}^{d}$ is a nonempty set. Then we say $\phi$
satisfies the QC defined by $\mathcal{Q}_{\phi}$ on $\mathcal{X}$.

###### Lemma 1

Consider the neural network described in (10).

  1. 1.

Suppose $\gamma_{i,k}\in\mathcal{E}(\mu_{\gamma_{i,k}},\Sigma_{i,k})$. Then
for any $\tau_{i}\geq 0$ we have

$\begin{bmatrix}\mathbf{z}\\\
1\end{bmatrix}^{\top}\left(\sum_{i=1}^{q}\tau_{i}M_{i}\right)\begin{bmatrix}\mathbf{z}\\\
1\end{bmatrix}\geq 0,$ (13)

where

$M_{i}\\!=\\!E_{i}^{\top}\begin{bmatrix}-\Sigma_{\gamma_{i,k}}^{-1}&\Sigma_{\gamma_{i,k}}^{-1}\mu_{\gamma_{i,k}}\\\
\mu_{\gamma_{i,k}}^{\top}\Sigma_{\gamma_{i,k}}^{-1}&-\mu_{\gamma_{i,k}}^{\top}\Sigma_{\gamma_{i,k}}^{-1}\mu_{\gamma_{i,k}}+1\end{bmatrix}E_{i},$

and $E_{i}$ is defined as in ($\star$ ‣ IV).

  2. 2.

Let $U\in\mathbb{S}^{n_{\pi}}$, $V\in\mathbb{R}^{n_{\pi}}$. Suppose
$\mathbf{z}$ satisfies the the quadratic inequality

$\begin{bmatrix}\mathbf{z}\\\
1\end{bmatrix}^{\top}M_{out}\begin{bmatrix}\mathbf{z}\\\ 1\end{bmatrix}\leq
0,$ (14)

where

$M_{out}\\!=\\!\begin{bmatrix}W^{\ell}S^{\ell}&b^{\ell}\\\
0&1\end{bmatrix}^{\top}\begin{bmatrix}U^{2}&UV\\\
V^{\top}U&V^{\top}V\\!-\\!1\end{bmatrix}\begin{bmatrix}W^{\ell}S^{\ell}&b^{\ell}\\\
0&1\end{bmatrix}.$

Then we have $\pi(z^{0})\in\mathcal{E}(-U^{-1}V,U^{-2})$.

###### Proof:

(1) The proof is a slight modification of [8]. We know that each input
$\gamma_{i,k}$ is bounded by the ellipsoid
$\mathcal{E}(\mu_{\gamma_{i,k}},\Sigma_{\gamma_{i,k}})$ , which means

$(\gamma_{i,k}-\mu_{\gamma_{i,k}})^{\top}\Sigma_{\gamma_{i,k}}^{-1}(\gamma_{i,k}-\mu_{\gamma_{i,k}})<1.$
(15)

This can be rewritten as

$\begin{bmatrix}\gamma_{i,k}\\\
1\end{bmatrix}^{\top}\begin{bmatrix}-\Sigma_{\gamma_{i,k}}^{-1}&\Sigma_{\gamma_{i,k}}^{-1}\mu_{\gamma_{i,k}}\\\
\mu_{\gamma_{i,k}}^{\top}\Sigma_{\gamma_{i,k}}^{-1}&-\mu_{\gamma_{i,k}}^{\top}\Sigma_{\gamma_{i,k}}^{-1}\mu_{\gamma_{i,k}}+1\end{bmatrix}\begin{bmatrix}\gamma_{i,k}\\\
1\end{bmatrix}\geq 0.$

Using the selector matrix $E_{i}$ defined in ($\star$ ‣ IV), we can rewrite
the preceding inequality as

$\begin{bmatrix}\mathbf{z}\\\
1\end{bmatrix}^{\top}\underbrace{E_{i}^{\top}\begin{bmatrix}-\Sigma_{\gamma_{i,k}}^{-1}&\Sigma_{\gamma_{i,k}}^{-1}\mu_{\gamma_{i,k}}\\\
\mu_{\gamma_{i,k}}^{\top}\Sigma_{\gamma_{i,k}}^{-1}&-\mu_{\gamma_{i,k}}^{\top}\Sigma_{\gamma_{i,k}}^{-1}\mu_{\gamma_{i,k}}+1\end{bmatrix}E_{i}}_{M_{i}}\begin{bmatrix}\mathbf{z}\\\
1\end{bmatrix}\geq 0,$

This implies for any $\tau_{i}\geq 0$, the inequality (13) holds. (2) For the
proof of the second part, we rewrite $\pi(z^{0})\in\mathcal{E}(-U^{-1}V,\
U^{-2})$ as

$\begin{bmatrix}\pi(z^{0})\\\ 1\end{bmatrix}^{\top}\begin{bmatrix}U^{2}&UV\\\
V^{\top}U&V^{\top}V-1\end{bmatrix}\begin{bmatrix}\pi(z^{0})\\\
1\end{bmatrix}\leq 0.$ (16)

This can be mapped to the space of vector $\mathbf{z}$ as,

$\begin{bmatrix}\mathbf{z}\\\
1\end{bmatrix}^{\top}\begin{bmatrix}W^{\ell}S^{\ell}&b^{\ell}\\\
0&1\end{bmatrix}^{\top}\begin{bmatrix}U^{2}&UV\\\
V^{\top}U&V^{\top}V\\!-\\!1\end{bmatrix}\begin{bmatrix}W^{\ell}S^{\ell}&b^{\ell}\\\
0&1\end{bmatrix}\begin{bmatrix}\mathbf{z}\\\ 1\end{bmatrix}\leq 0,$

which characterizes $M_{out}$ provided in the lemma. ∎

The quadratic constraint for the activation layers, denoted by $\mathcal{Q}$
is presented in [9] (Lemma 4) and remains unchanged in our formulation. Now
consider the implicit equation $B\mathbf{z}=\phi(A\mathbf{z}+b)$ in (10)
describing the neural network, where $\phi$ satisfies the quadratic constraint
defined by $\mathcal{Q}$. By Definition 2 this implies that for any
$Q\in\mathcal{Q}$,

$\begin{bmatrix}\mathbf{z}\\\
1\end{bmatrix}^{\top}\underbrace{\begin{bmatrix}A&b\\\ B&0\\\
0&1\end{bmatrix}^{\top}Q\begin{bmatrix}A&b\\\ B&0\\\
0&1\end{bmatrix}}_{M_{mid}(Q)}\begin{bmatrix}\mathbf{z}\\\ 1\end{bmatrix}\geq
0.$ (17)

Returning to our main goal, bounding the output $\pi(z^{0})$ by an ellipsoid
given ellipsoidal bounds on the inputs $\gamma_{i,k}$, we know that the
quadratic inequalities (13) and (17) hold and we would like to conclude the
quadratic inequality (14). To this end we use the S-procedure. Therefore, the
multiple ellipsoidally bounded input version of Theorem 1 (Output covering
ellipsoid) of [8] can be presented as follows.

###### Theorem 1

Consider the multi-layer neural network described by (9). Suppose
$\gamma_{i}\in\mathcal{E}(\mu_{\gamma_{i,k}},\Sigma_{\gamma_{i,k}})$ and ReLU
activation function $\phi$ satisfies the quadratic constraints defined by
(17). Let $U\in\mathbb{S}^{n_{\pi}}$, $V\in\mathbb{R}^{n_{\pi}}$ be two
matrices that satisfy $M(\tau,Q,U,V)\preceq 0$ for some $Q\in\mathcal{Q}$ and
$\tau\in\mathbb{R}_{+}^{q}$, where

$M\\!=\\!\left[\begin{matrix}\begin{array}[]{c|c}\sum_{i=1}^{q}\tau_{i}M_{i}\\!+\\!M_{mid}(Q)\\!-\\!ee^{\top}&\begin{matrix}0_{(N_{\mathbf{z}}-N_{\ell})\times
n_{\pi}}\\\ \\\ W^{{\ell}^{\top}U}\\\ \\\
{b^{\ell}}^{\top}U+V^{\top}\end{matrix}\\\
\hline\cr\begin{matrix}0_{n_{\pi}\times(N_{\mathbf{z}}-N_{\ell})}\ \ \
UW^{\ell}\ \ \
Ub^{\ell}+V\end{matrix}&-I_{n_{\pi}}\end{array}\end{matrix}\right].$ (18)

Then we have
$\pi(z^{0})\in\mathcal{E}_{\pi(z^{0})}=\mathcal{E}(-U^{-1}V,U^{-2})$.

###### Proof:

According to the first part of Lemma 1, in order to guarantee
$\pi(z^{0})\in\mathcal{E}_{\pi(z^{0})}=\mathcal{E}(-U^{-1}V,U^{-2})$, we need
to provide a sufficient condition to satisfy the inequality,

$\begin{bmatrix}\mathbf{z}\\\
1\end{bmatrix}^{\top}M_{out}\begin{bmatrix}\mathbf{z}\\\ 1\end{bmatrix}\leq
0.$ (19)

Suppose the following matrix inequality holds,

$\sum_{i=1}^{q}\tau_{i}M_{i}+M_{mid}(Q)+M_{out}\preceq 0,$ (20)

for some $\tau_{i}\geq 0$, $Q\in\mathcal{Q}$. Based on (20) we can conclude,

$\sum_{i=1}^{q}\underbrace{\begin{bmatrix}\mathbf{z}\\\
1\end{bmatrix}^{\top}\\!\\!\\!\tau_{i}M_{i}\begin{bmatrix}\mathbf{z}\\\
1\end{bmatrix}}_{\text{$\geq 0$}}+\underbrace{\begin{bmatrix}\mathbf{z}\\\
1\end{bmatrix}^{\top}\\!\\!\\!M_{mid}\begin{bmatrix}\mathbf{z}\\\
1\end{bmatrix}}_{\text{$\geq 0$}}+\begin{bmatrix}\mathbf{z}\\\
1\end{bmatrix}^{\top}\\!\\!\\!M_{out}\begin{bmatrix}\mathbf{z}\\\
1\end{bmatrix}\ \leq 0,$ (21)

which provides what we need from $M_{out}$. The matrix inequality in (20) is
not linear in $U$ and $V$. However, it can be made linear through the Schur
complement. First, rewrite $M_{out}$ as

$\displaystyle M_{out}$
$\displaystyle=\begin{bmatrix}W^{\ell}S^{\ell}&b^{\ell}\\\
0&1\end{bmatrix}^{\top}\begin{bmatrix}U&V\end{bmatrix}^{\top}\begin{bmatrix}U&V\end{bmatrix}\begin{bmatrix}W^{\ell}S^{\ell}&b^{\ell}\\\
0&1\end{bmatrix}-$ $\displaystyle\begin{bmatrix}W^{\ell}S^{\ell}&b^{\ell}\\\
0&1\end{bmatrix}^{\top}\begin{bmatrix}0_{1\times
n_{\pi}}&1\end{bmatrix}^{\top}\begin{bmatrix}0_{1\times
n_{\pi}}&1\end{bmatrix}\begin{bmatrix}W^{\ell}S^{\ell}&b^{\ell}\\\
0&1\end{bmatrix}$ $\displaystyle:=F^{\top}F-ee^{\top},$ (22)

where

$\displaystyle
F=\begin{bmatrix}0_{n_{\pi\times(N_{\mathbf{z}}-N_{\ell})}}&UW^{\ell}&Ub^{\ell}+V\end{bmatrix},$
(23) $\displaystyle e=\underbrace{[\ 0,\ 0,\ 0,\ ..,\
1]^{\top}}_{\text{$N_{\mathbf{z}}+1$}}.$

Through the application of the Schur complement on $M_{out}$, the matrix
inequality in (20) can be written as (18), which is linear in
$\left(V,U,Q\right)$. ∎

Using the result of Theorem 1, the tightest bound on $\pi(z^{0})$ is obtained
from the following semidefinite program,

$\left\\{\begin{aligned} \mathrm{min}_{U,V,Q}&\ -\mathrm{logdet}(U)\
\text{or}\ -\mathrm{tr}(U)\\\ \text{s.t.}\quad&M(\tau,Q,U,V)\preceq 0,\\\
&\tau\geq 0,\ Q\in\mathcal{Q}\text{ from
\cite[cite]{[\@@bibref{}{fazlyab2019safety}{}{}]} (Lemma
4)},\end{aligned}\right.$ (24)

where the objective function can be any metric related to the ellipsoid
volume.

###  IV-A Application to Prediction Bound Propagation

We now apply the general multiple ellipsoid input robustness bound developed
in Theorem 1 to find the prediction bound
$\mathcal{E}(\hat{\mu}_{k+1},\hat{\Sigma}_{k+1})$ with input ellipsoids
$\mathcal{E}(y_{k-N},\bar{\Sigma}_{v})$, $\cdots$,
$\mathcal{E}(y_{k},\bar{\Sigma}_{v})$.

###### Proposition 2

Let $\gamma_{i,k}\in\mathcal{E}(y_{k-i+1},\bar{\Sigma}_{v})$ for
$i=1,2,\cdots,q=N+1$ and define $\hat{\Sigma}_{k+1}=U^{-2}$,
$\hat{\mu}_{k+1}=-U^{-1}V$, and $\mathcal{C}_{k}$ as the actual (non-
ellipsoidal) prediction set of the neural network. Then
$\pi(z^{0})\in\mathcal{C}_{k}$ is bounded by
$\mathcal{E}(\hat{\mu}_{k+1},\hat{\Sigma}_{k+1})$ where the tightest bound
comes from the convex optimization (24).

##  V Numerical Experiments

We consider a linear beam and slider and a nonlinear water tank cascade to
demonstrate our proposed detector.

(a) Beam and slider

(b) $k=20$ Prediction bound

(c) Robustness due to training

Figure 2: (a) The linear beam and slider system. (b) The effectiveness of the
prediction bound at $k=20$. (c) A neural network trained with noisy data has
better robustness characteristics than one trained with noiseless data and is
quantified well by our prediction bounds.

###  V-A Linear Beam and Slider

In the beam and slider (see Fig. 2(a)), a sensor measures the position of the
slider which slides along the beam with a rate equal to its distance from the
origin. The beam rotates with constant angular velocity of
$3\pi\frac{rad}{sec}$. The sensor noise is Gaussian distributed,
$v_{k}\sim\mathcal{N}(0,\Sigma_{v})$, and we truncate it with a $\bar{p}=95\%$
confidence ellipse $\mathcal{E}(y_{k},\bar{\Sigma}_{v})$ with

$\Sigma_{v}=\begin{bmatrix}0.0214&0.0112\\\ 0.0112&0.0217\end{bmatrix}\
\text{and}\ \bar{\Sigma}_{v}=\begin{bmatrix}0.1282&0.0671\\\
0.0671&0.1300\end{bmatrix}.$

We train a NARX feedforward neural network with two hidden layers (10 neurons
in the first layer, 2 neurons in the second) with ReLU activation functions
from noisy sensor data. The data was generated through simulation of the beam
and slider system based on the linear discrete time model,

$\left\\{\begin{aligned}
x_{k+1}&=0.8\begin{bmatrix}[r]\text{cos}(\beta)&-\text{sin}(\beta)\\\
\text{sin}(\beta)&\text{cos}(\beta)\end{bmatrix}x_{k},\quad\beta=\frac{3\pi}{5},\\\
y_{k}&=x_{k}+v_{k}.\end{aligned}\right.$ (25)

We use the current measurement and one past measurement (i.e., $N=1$) to build
the NARX model (we still call it NARX because the neural network returns a
nonlinear approximation of this linear model).

Figure 2(b) demonstrates a snapshot of the efficacy of our methods under
normal operation at $k=20$. The prediction bound
$\mathcal{E}(\hat{x}_{21},\hat{\Sigma}_{21})$ (orange ellipse) can be quite
tight on the actual prediction set $\mathcal{C}_{21}$ (gray region), which is
the set of all possible neural network output estimates under all possible
combinations from the input ellipsoids. As designed, the estimate
$\hat{y}_{21}$ (is not necessarily the center $\hat{\mu}_{21}$) and ideal
measurement $y^{*}_{21}$ (assuming good accuracy of the estimation) are within
the prediction bound. The actual measurement $y_{21}$ is within the geometric
sum of the prediction bound (orange) and confidence ellipsoid (blue); in this
case it is also within the prediction bound, as the noise realization is
relatively small. Under normal operation we received 0.7% false alarms
compared with the upper bound on the false alarm rate
$\bar{\mathcal{A}}_{s}=1-\bar{p}^{3}=14.26\%$. This conservatism below the
upper bound illustrates that the distribution within the prediction boundary,
and the subsequent geometric sum, is an important factor for improving our
ability to predict the false alarm rate and an important direction for future
work.

Despite the conservatism, we demonstrate that it is an effective tool for
detecting anomalies in behavior. The first fault we consider is a vibration
generated on the shaft because of the rotor. This imposes a new additive
periodic displacement on the shaft and consequently the slider in a fixed
direction

$\delta x_{k}=\begin{bmatrix}\delta_{1}\\\
\delta_{1}\end{bmatrix}\text{sin}(k),$ (26)

where $\delta_{1}=0.3$ and $k$ is the time index. Under this fault the alarm
rate raises to $27.17\%$.

The second scenario we consider is a sensor anomaly in which the sensor
measurement is displaced by $\delta_{2}=0.3$,

$y_{k}=x_{k}+v_{k}+\begin{bmatrix}\delta_{2}\\\ \delta_{2}\end{bmatrix}.$ (27)

Under this fault the alarm rate raises to $25.35\%$.

The prediction bound provides a quantification of the robustness of a neural
network to perturbed input. It is intuitive that a neural network trained by
noisy data should have enhanced robustness compared with a neural network
trained only by ideal (no noise) sensor measurements since the model
parameters are learned to filter the noise more effectively. The prediction
bound allows us to characterize this intuitive relationship quantitatively. In
Fig. 2(c), we provide the prediction bounds for this system at exactly
identical conditions, with the only change that one model is trained with the
actual noisy data and the other uses the underlying ideal (no noise) outputs.
The larger prediction ellipsoid bound for the neural network trained without
noise (dashed ellipse) quantifies the value of training with noisy data.

(a) Tank System

(b) $k=302$ Prediction Bound

Figure 3: (a) The nonlinear water tank system. (b) The effectiveness of the
prediction bound at $k=302$.

###  V-B Nonlinear Water Tank Cascade

We now consider a two tank cascade (see Fig. 3(a)) in which flow out of the
tanks is driven by gravity and yields a nonlinear dynamic,

$\left\\{\begin{aligned}
\dot{h_{1}}(t)&=Q_{in}-c_{d}A_{d}\sqrt{2gh_{1}(t)},\\\
\dot{h_{2}}(t)&=c_{d}A_{d}\sqrt{2gh_{1}(t)}-c_{d}A_{d}\sqrt{2gh_{2}(t)},\\\
y_{k}&=h(t_{k})+v_{k}.\end{aligned}\right.$ (28)

where $t_{k}=k\Delta t$ with inter-sampling period $\Delta t=0.02$ sec,
$Q_{in}=15$, $c_{d}=0.9$, $A_{d}=1$ cm2, and $g$ is the gravitational
constant. For simplicity we use the same noise distribution and confidence
level $\bar{p}$ as the prior example.

We train the neural network (two hidden layers: 20 neurons in the first layer;
5 neurons in the second layer) with the current noisy measurement and three
(i.e., $N=3$) past noisy measurements generated by simulating the above
equations.

Figure 3(b) provides a snapshot of the performance of our approach under
normal operation, at $k=302$. Under normal operation the detector generates
zero false alarms compared with the upper bound of
$\bar{\mathcal{A}}_{s}=1-\bar{p}^{N+2}=22.62\%$.

Here we simulate a fault in the drainage of the lower tank (e.g., debris
caught in the drain) such that 20% of the lower drain area is blocked. In this
case, our proposed detector generates alarms at a rate of $58.1\%$ in steady
state.

##  VI Conclusion

When first principle models of dynamical systems are too difficult to
construct there is a growing tendency to turn to data-driven techniques. In
this paper, we leverage recent results on the robustness of neural network
predictions under input perturbations to compute bounds on the estimates
produced from noisy measurements when a neural network is used to approximate
an autoregressive model of the output measurements. We use this bound to
define normal behavior under typical measurement noise and use this as a
boundary to detect abnormal behavior.

This detector performs well and provides tight detection for the extreme
points of normal behavior. This study has revealed that despite the tightness
on these extreme points, a more accurate picture of the distribution of normal
behavior is needed to produce tighter alarm rates. Nonetheless, we demonstrate
the proposed detector is able to identify anomalies of a variety of types in a
systematic way.

## References

  * [1] M. H. Hassoun _et al._ , _Fundamentals of artificial neural networks_. MIT press, 1995.
  * [2] W. Liu, “An extended kalman filter and neural network cascade fault diagnosis strategy for the glutamic acid fermentation process,” _Artificial intelligence in engineering_ , vol. 13, no. 2, pp. 131–140, 1999.
  * [3] A. T. Vemuri and M. M. Polycarpou, “Neural-network-based robust fault diagnosis in robotic systems,” _IEEE Transactions on neural networks_ , vol. 8, no. 6, pp. 1410–1420, 1997.
  * [4] M. Wlas, Z. Krzeminski, J. Guzinski, H. Abu-Rub, and H. A. Toliyat, “Artificial-neural-network-based sensorless nonlinear control of induction motors,” _IEEE Transactions on Energy Conversion_ , vol. 20, no. 3, pp. 520–528, 2005.
  * [5] P. M. Frank and B. Köppen-Seliger, “New developments using ai in fault diagnosis,” _Engineering Applications of Artificial Intelligence_ , vol. 10, no. 1, pp. 3–14, 1997.
  * [6] Z. Li, L. Ma, and K. Khorasani, “Fault detection in reaction wheel of a satellite using observer-based dynamic neural networks,” in _International Symposium on Neural Networks_. Springer, 2005, pp. 584–590.
  * [7] A. Abbaspour, P. Aboutalebi, K. K. Yen, and A. Sargolzaei, “Neural adaptive observer-based sensor and actuator fault detection in nonlinear systems: Application in uav,” _ISA transactions_ , vol. 67, pp. 317–329, 2017.
  * [8] M. Fazlyab, M. Morari, and G. J. Pappas, “Probabilistic verification and reachability analysis of neural networks via semidefinite programming,” _arXiv preprint arXiv:1910.04249_ , 2019.
  * [9] ——, “Safety verification and robustness analysis of neural networks via quadratic constraints and semidefinite programming,” _arXiv preprint arXiv:1903.01287_ , 2019.
  * [10] S. Dutta, S. Jha, S. Sanakaranarayanan, and A. Tiwari, “Output range analysis for deep neural networks,” _arXiv preprint arXiv:1709.09130_ , 2017.
  * [11] A. Lomuscio and L. Maganti, “An approach to reachability analysis for feed-forward relu neural networks,” _arXiv preprint arXiv:1706.07351_ , 2017\.
  * [12] E. Wong and Z. Kolter, “Provable defenses against adversarial examples via the convex outer adversarial polytope,” in _International Conference on Machine Learning_ , 2018, pp. 5286–5295.
  * [13] X. Huang, M. Kwiatkowska, S. Wang, and M. Wu, “Safety verification of deep neural networks,” in _International Conference on Computer Aided Verification_. Springer, 2017, pp. 3–29.
  * [14] A. Raghunathan, J. Steinhardt, and P. S. Liang, “Semidefinite relaxations for certifying robustness to adversarial examples,” in _Advances in Neural Information Processing Systems_ , 2018, pp. 10 877–10 887.
  * [15] G. Katz, C. Barrett, D. L. Dill, K. Julian, and M. J. Kochenderfer, “Reluplex: An efficient smt solver for verifying deep neural networks,” in _International Conference on Computer Aided Verification_. Springer, 2017, pp. 97–117.
  * [16] N. Hashemi, J. Ruths, and M. Fazlyab, “Certifying incremental quadratic constraints for neural networks,” _arXiv preprint arXiv:2012.05981_ , 2020\.
  * [17] A. H. Abdelaziz, S. Watanabe, J. R. Hershey, E. Vincent, and D. Kolossa, “Uncertainty propagation through deep neural networks,” 2015.
  * [18] J. S. Titensky, H. Jananthan, and J. Kepner, “Uncertainty propagation in deep neural networks using extended kalman filtering,” _arXiv preprint arXiv:1809.06009_ , 2018.
  * [19] C. Murguia and J. Ruths, “On model-based detectors for linear time-invariant stochastic systems under sensor attacks,” _IET Control Theory & Applications_, vol. 13, no. 8, pp. 1051–1061, 2019.
  * [20] S. Boyd, S. P. Boyd, and L. Vandenberghe, _Convex optimization_. Cambridge university press, 2004.
  * [21] N. Hashemi and J. Ruths, “Co-design for performance and security: Geometric tools,” _arXiv preprint arXiv:2006.08739_ , 2020.
  * [22] M. Fazlyab, M. Morari, and G. J. Pappas, “Safety verification and robustness analysis of neural networks via quadratic constraints and semidefinite programming,” _IEEE Transactions on Automatic Control_ , pp. 1–1, 2020.

##  Appendix A comparison of multiple input approach with single input
approach

Consider the multiple ellipsoidal bounded inputs $\gamma_{i,k}$ introduced in
Lemma 1. The quadratic constraint for this ellipsoid results to matrix $M_{i}$
and the sum $\sum_{i=1}^{q}\tau_{i}M_{i}$ can be written as (suppressing the
$k$ index for readability)

$\begin{bmatrix}-\tau_{1}\Sigma_{\gamma_{1}}^{-1}\\!\\!&0&\cdots\\!\\!&\tau_{1}\Sigma_{\gamma_{i}}^{-1}\mu_{\gamma_{1}}\\\
0&-\tau_{2}\Sigma_{\gamma_{2}}^{-1}\\!\\!&\cdots\\!\\!&\tau_{2}\Sigma_{\gamma_{2}}^{-1}\mu_{\gamma_{2}}\\\
\vdots&\vdots&\ddots\\!\\!&\vdots\\\
\tau_{1}\mu_{\gamma_{1}}^{\top}\Sigma_{\gamma_{1}}^{-1}\\!\\!&\tau_{2}\mu_{\gamma_{2}}^{\top}\Sigma_{\gamma_{2}}^{-1}\\!\\!&\dots\\!\\!&-\sum_{i=1}^{q}\tau_{i}\\!\left(\mu_{\gamma_{i}}^{\top}\Sigma_{\gamma_{i}}^{-1}\mu_{\gamma_{i}}\\!\\!-\\!1\right)\end{bmatrix}.$
(29)

On the other hand if we follow [8] and define an overall ellipsoidal bound on
the input vector, this bound should satisfy,

$\sum_{i=1}^{q}(\gamma_{i,k}-\mu_{\gamma_{i,k}})^{\top}\Sigma_{\gamma_{i,k}}^{-1}(\gamma_{i,k}-\mu_{\gamma_{i,k}})\leq
q$ (30)

(compare to (15)) where the equality happens where all the points
$\gamma_{i,k}$ are selected from the boundary of the input ellipsoidal bounds.
Constructing $\sum_{i=1}^{q}\tau_{i}M_{i}$ in this case has only one term so
$q=1$ and $\tau_{1}=\tau$ leading to

$\tau\begin{bmatrix}-\Sigma_{\gamma_{1}}^{-1}&0&\cdots&\Sigma_{\gamma_{i}}^{-1}\mu_{\gamma_{1}}\\\
0&-\Sigma_{\gamma_{2}}^{-1}&\cdots&\Sigma_{\gamma_{2}}^{-1}\mu_{\gamma_{2}}\\\
\vdots&\vdots&\ddots&\vdots\\\
\mu_{\gamma_{1}}^{\top}\Sigma_{\gamma_{1}}^{-1}&\mu_{\gamma_{2}}^{\top}\Sigma_{\gamma_{2}}^{-1}&\dots&-\sum_{i=1}^{q}\mu_{\gamma_{i}}^{\top}\Sigma_{\gamma_{i}}^{-1}\mu_{\gamma_{i}}+q\end{bmatrix}.$
(31)

Clearly (31) is a specific form of (29), in which the single decision variable
$\tau$ in (31) is replaced with $q$ decision variables
$\\{\tau_{i}\\}_{i=1}^{q}$ in (29). Hence any solution for (31) can be
expressed by $\eqref{eq:new}$ justifying that the multiple input ellipsoid
approach can do no worse than the single input ellipsoid approach (i.e., it is
as good or better than the single input approach).


A Scalable Gradient-Free Method for Bayesian
Experimental Design with Implicit Models

Jiaxin Zhang Sirui Bi Guannan Zhang Computer Science and Mathematics Division
Oak Ridge National Laboratory<EMAIL_ADDRESS>Computational Sciences and
Engineering Division Oak Ridge National Laboratory<EMAIL_ADDRESS>Computer
Science and Mathematics Division Oak Ridge National Laboratory<EMAIL_ADDRESS>
###### Abstract

Bayesian experimental design (BED) is to answer the question that how to
choose designs that maximize the information gathering. For implicit models,
where the likelihood is intractable but sampling is possible, conventional BED
methods have difficulties in efficiently estimating the posterior distribution
and maximizing the mutual information (MI) between data and parameters. Recent
work proposed the use of gradient ascent to maximize a lower bound on MI to
deal with these issues. However, the approach requires a sampling path to
compute the pathwise gradient of the MI lower bound with respect to the design
variables, and such a pathwise gradient is usually inaccessible for implicit
models. In this paper, we propose a novel approach that leverages recent
advances in stochastic approximate gradient ascent incorporated with a
smoothed variational MI estimator for efficient and robust BED. Without the
necessity of pathwise gradients, our approach allows the design process to be
achieved through a unified procedure with an approximate gradient for implicit
models. Several experiments show that our approach outperforms baseline
methods, and significantly improves the scalability of BED in high-dimensional
problems.

##  1 Introduction

Experimental design plays an essential role in all scientific disciplines. Our
ultimate goal is to determine designs that maximize the information gathered
through the experiments so that improve our understanding on model comparison
or parameter estimations. A broadly used approach is Bayesian experimental
design (BED) (Chaloner and Verdinelli, 1995) that aims to find optimal design
$\bm{\xi}^{*}$ to maximize a utility function $I(\bm{\xi})$, which is
typically defined by the mutual information (MI) between data and model
parameters

$I(\bm{\xi})=\mathbb{E}_{p(\bm{\theta},\bm{y}|\bm{\xi})}\left[\log
p(\bm{\theta}|\bm{y},\bm{\xi})-\log p(\bm{\theta})\right]$ (1)

where $\bm{\xi}$ represent the experimental design, $p(\bm{\theta})$ is the
prior distribution of model parameters and $p(\bm{\theta}|\bm{y},\bm{\xi})$ is
the posterior distribution of $\bm{\theta}$ given data $\bm{y}$ with design
$\bm{\xi}$. However, finding $\bm{\xi}^{*}$ by maximizing $I(\bm{\xi})$ is a
challenging task in practice because evaluating $I(\bm{\xi})$ needs a nested
Monte Carlo simulator, which is computationally intensive, particularly in
high-dimensional settings (Drovandi et al., 2018).

Most of the existing BED studies focus on the explicit models (Chaloner and
Verdinelli, 1995; Sebastiani and Wynn, 2000; Foster et al., 2019, 2020) where
the likelihood is analytically known but in nature and science, the more
common scenario is the implicit models (Kleinegesse and Gutmann, 2019;
Overstall et al., 2018), where the likelihood is intractable and not evaluated
directly but sampling is possible. In other words, the model is specified
based on a stochastic data generating simulator and typically has no access to
the analytical form and the gradients of the joint density
$p(\bm{\theta},\bm{y}|\bm{\xi})$ and marginal density $p(\bm{y}|\bm{\xi})$.
The resulting BED scheme typically shares a two-stage feature: build a
pointwise estimator of $I(\bm{\xi})$ and then feed this “black-box” estimator
to a separate outer-level optimizer to find the optimal design $\bm{\xi}^{*}$.
This framework substantially increases the overall computational cost and is
difficult to scale the BED process to a high dimensional design space. Recent
studies (Kleinegesse and Gutmann, 2020; Foster et al., 2020; Harbisher et al.,
2019) alleviate the challenges by introducing stochastic gradient-based
approaches but they rely on the models with tractable likelihood functions or
assume the gradient can be reasonably approximated by pathwise gradient
estimators with sampling path (Kleinegesse and Gutmann, 2020), which requires
that we can sample from the data distribution $p(\bm{y}|\bm{\theta},\bm{\xi})$
by sampling from a base distribution $p(\bm{\epsilon})$ and then transforming
the samples through a specific function
$g(\bm{\epsilon};\bm{\theta},\bm{\xi})$, which is called the sampling path,
i.e., $\bm{y}\sim
p(\bm{y}|\bm{\theta},\bm{\xi})\Longleftrightarrow\bm{y}=g(\bm{\epsilon};\bm{\theta},\bm{\xi}),\bm{\epsilon}\sim
p(\bm{\epsilon})$. This is unlike the scope of our paper that focuses on the
BED for implicit models without gradients.

###  1.1 Related Work

Foster et al. (2019) recently proposed to use a lower bound of MI for BED.
This study relies on variational approximations to the likelihood and
posterior but it is a two-stage approach where the optimal designs were
determined by a separate Bayesian optimization (BO). Unfortunately, this
approach has a limitation in scaling to high-dimensional design problems. A
follow-up study developed by Foster et al. (2020) aims to address the
scalability issue by introducing a unified stochastic gradient-based approach.
However, they assumed the models with the tractable explicit likelihood or
gradient approximations are available. In the scope of BED for implicit
models, Ao and Li (2020) proposed an approximate KLD based BED method for
models with intractable likelihoods; Kleinegesse and Gutmann (2019);
Kleinegesse et al. (2020) have recently considered the use of MI combined with
likelihood-free inference by ratio estimation to approximate the posterior
distribution but this method is often computationally intensive. The authors
rectify this in a follow-up study (Kleinegesse and Gutmann, 2020) that
leverages mutual information neural estimation (MINE) (Belghazi et al., 2018)
to jointly determine the optimal design and the posterior distribution.
However, for the experimental design without gradients, this method falls back
to the gradient-free methods, i.e., BO at the expense of reduced scalability,
to solve the optimization problems within the two-stage framework. Moreover,
the variational BED methods (Kleinegesse and Gutmann, 2020; Foster et al.,
2019, 2020) may yield an unstable design with a high variance on the posterior
distribution due to the use of variational MI estimators that exhibit a high
variance argued by Song and Ermon (2019).

Contributions In this paper, we propose a novel Bayesian experimental design
approach for implicit models without available gradients. This scalable method
can address the aforementioned technical challenges, particularly in high
dimensional problems.

  * •

We propose a general unified framework that leverages stochastic approximate
gradient without the requirement or assumption of pathwise gradients and
sampling paths for implicit models.

  * •

We propose to use a smoothed MI lower bound to conduct robust MI estimation
and optimization, which allows the variance of the design and posterior
distribution to be much lower than existing approaches.

  * •

We show the superior performance of the proposed approach through several
experiments and demonstrate that the approach enables the optimization to be
performed by a stochastic gradient ascent algorithm and thus well scaled to
considerable high dimensional design problems.

##  2 Bayesian Experimental Design

The Bayesian experimental design (BED) framework aims at choosing an
experimental design $\bm{\xi}$ to maximize the information gained about some
parameters of interest $\bm{\theta}$ from the outcome $\bm{y}$ of the
experiment. Typically, the BED framework begins with a Bayesian model of the
experimental process, including a prior distribution $p(\bm{\theta})$ and a
likelihood $p(\bm{y}|\bm{\theta},\bm{\xi})$. The information gained about
$\bm{\theta}$ from running the experiment with design $\bm{\xi}$ and observed
outcome $\bm{y}$ can be interpreted by the reduction in entropy from the prior
to posterior

$\textup{IG}(\bm{y},\bm{\xi})=\mathcal{Q}[p(\bm{\theta})]-\mathcal{Q}[p(\bm{\theta}|\bm{y},\bm{\xi})].$
(2)

To define a metric to quantify the utility of the design $\bm{\xi}$ before
running experiments, an expected information gain (EIG), $I(\bm{\xi})$ is
often used by

$I(\bm{\xi})=\mathbb{E}_{p(\bm{y}|\bm{\xi})}[\mathcal{Q}[p(\bm{\theta})]-\mathcal{Q}[p(\bm{\theta}|\bm{y},\bm{\xi})]].$
(3)

Eq. (3) can also be interpreted as a mutual information (MI) between
$\bm{\theta}$ and $\bm{y}$ with a specified $\bm{\xi}$,

$I_{\rm
MI}(\bm{\xi})=\mathbb{E}_{p(\bm{\theta})p(\bm{y}|\bm{\theta},\bm{\xi})}\left[\log\frac{p(\bm{y}|\bm{\theta},\bm{\xi})}{p(\bm{y}|\bm{\xi})}\right].$
(4)

The BED problem is therefore defined as

$\bm{\xi}^{*}=\operatorname*{arg\,max}_{\bm{\xi}\in\bm{\Xi}}I_{\rm
MI}(\bm{\xi}),$ (5)

where $\bm{\Xi}$ is the feasible design domain. The most challenging task in
the BED framework is how to efficiently and accurately estimate $I_{\rm
MI}(\bm{\xi})$ in Eq. (4) and optimize $I_{\rm MI}(\bm{\xi})$ via Eq. (5) to
obtain the optimal design $\bm{\xi}^{*}$. The following section will discuss
the two core tasks in the BED framework, that is MI estimation and
optimization.

###  2.1 Mutual Information Estimation

Mutual information (MI) estimation plays a critical role in many important
problems, not only the BED framework but also other machine learning tasks
such as reinforcement learning (Pathak et al., 2017) and representation
learning (Chen et al., 2016; Oord et al., 2018). However, estimating mutual
information from samples is always challenging (McAllester and Stratos, 2020).
This is because that, in general, neither $p(\bm{\theta}|\bm{y},\bm{\xi})$ and
$p(\bm{y}|\bm{\xi})$ have analytical closed-form so that classical Monte Carlo
methods are intractable to compute the integral, specifically in high
dimensions. One potential approach is to use a nested Monte Carlo (NMC)
estimator (Myung et al., 2013; Rainforth et al., 2018), which is given by

$I_{\rm
NMC}(\bm{\xi})=\frac{1}{N}\sum_{i=1}^{N}\log\left[\frac{p(\bm{y}_{i}|\bm{\theta}_{i,0},\bm{\xi})}{\frac{1}{m}\sum_{j=1}^{M}p(\bm{y}_{i}|\bm{\theta}_{i,j},\bm{\xi})}\right],$
(6)

where $\bm{\theta}_{i,j}\sim p(\bm{\theta})$ are independent and identically
distributed (i.i.d) samples, and $\bm{y}_{i}\sim
p(\bm{y}|\bm{\theta}=\bm{\theta}_{i,0},\bm{\xi})$. However, NMC estimator
requires both the inner and outer integrals and therefore converges slowly at
an overall rate of $\mathcal{O}(T^{-1/3})$ (Rainforth et al., 2018) in the
total computational cost $T=\mathcal{O}(NM)$.

Recently, there has been an increasingly interest in MI estimation with
variational methods (Barber and Agakov, 2003; Nguyen et al., 2010; Foster et
al., 2019), which can benefit from a natural incorporation with deep learning
algorithms (Alemi et al., 2016; Oord et al., 2018; Belghazi et al., 2018;
Poole et al., 2019).

Although the variational approaches to MI estimation have been widely used,
there are still several limitations pointed out by Song and Ermon (2019), for
example, a high variance of the estimators in MINE (Belghazi et al., 2018).
This will bring additional challenges for the subsequent MI optimization. To
mitigate the issues, Song and Ermon (2019) proposed a unified framework over
variational estimators that treat variational MI estimation as an optimization
problem over density ratios. This is achieved by utilizing the role of
partition function estimation with variance reduction techniques. This
improved MI estimation, named by a smoothed MI lower bound estimator, $I_{\rm
SMILE}$, has been demonstrated by compared to existing variational methods in
Barber and Agakov (2003); Poole et al. (2019).

###  2.2 Mutual Information Optimization

The BED problem is to find the optimal design that maximizes the MI
estimation. Given a point-by-point base MI estimator, a variety of different
approaches can be used for the subsequent optimization over designs (Amzal et
al., 2006; Müller, 2005; Rainforth, 2017). At a high-level, most existing
methods belong to a two-stage procedure where an MI estimator of $I_{\rm
MI}(\bm{\xi})$ is first made and then followed by a separate optimization
algorithm which is used to select the candidate design $\bm{\xi}$ to evaluate
next. However, this two-stage framework can be computationally intensive
because it needs an extra level of nesting to the optimization process in Eq.
(5); in other words, $I_{\rm MI}(\bm{\xi})$ must be separately estimated for
each $\bm{\xi}$, which significantly increase the overall computational cost.

Recent studies (Foster et al., 2019; Kleinegesse and Gutmann, 2019) tend to
use Bayesian optimization (BO) (Snoek et al., 2012) for this purpose due to
its sample efficiency, robustness and capability to naturally deal with noisy
MI estimator (Foster et al., 2019). However, it is difficult for BO to scale
the overall BED process to high-dimensional design settings (Foster et al.,
2020; Kleinegesse and Gutmann, 2020). Broadly speaking, BO is prohibitively
slow and approaches the performance of “random search” if the dimensionality
of design space is $\Omega$>100 (Snoek et al., 2015). The scalability of BO is
still the major challenge, even though some recent studies (Wang et al., 2017;
Li et al., 2018; Mutny and Krause, 2018; Rana et al., 2017; Li et al., 2016;
Eriksson et al., 2019) are involved to make improvements.

To open the door of high-dimensional BED, Foster et al. (2019) proposed a
unified gradient-based BED method but they focus on models with tractable
likelihood functions or assume that the gradients are available. Kleinegesse
and Gutmann (2020) proposed a MINEBED method that leverages the neural MI
estimation to jointly determine the optimal design and the posterior using a
gradient-based optimization algorithm. However, the MINEBED only works for the
experimental design with available pathwise gradient estimators. As for the
implicit models without gradient, Kleinegesse and Gutmann (2020) still
consider a two-stage framework that uses BO as the outer-level optimizer at
the expense of reduced scalability.

##  3 The SAGABED Approach

We here show how to perform a stochastic approximate gradient ascent (SAGA)
method to design optimal experiments for implicit models without gradients and
address the computational challenges in scaling to high-dimensional problems.

###  3.1 Smoothed MI Estimator

Belghazi et al. (2018) have recently proposed to estimate MI by gradient
ascent over neural networks and argued that the lower bound can be tightened
by optimizing the neural network parameters. This MI estimator is typically
named by MINE-$f$ or $f$-GAN KL Nowozin et al. (2016)

$\displaystyle I_{\rm MINE}(\bm{\xi},\bm{\psi})=$
$\displaystyle\;\mathbb{E}_{p(\bm{\theta},\bm{y}|\bm{\xi})}[\mathcal{T}_{\bm{\psi}}(\bm{\theta},\bm{y})]$
(7)
$\displaystyle-\log\left(\mathbb{E}_{p(\bm{\theta})p(\bm{y}|\bm{\xi})}\left[e^{\mathcal{T}_{\bm{\psi}}(\bm{\theta},\bm{y})}\right]\right),$

where $\mathcal{T}_{\bm{\psi}}(\bm{\theta},\bm{y})$ is a neural network that
is parametrized by $\bm{\psi}$ with model parameter $\bm{\theta}$ and data
$\bm{y}$ as input. Incorporating neural network parameters $\bm{\psi}$ with
design parameters $\bm{\xi}$, the BED problem can be formulated by maximizing
the overall objective

$\bm{\xi}^{*}=\operatorname*{arg\,max}_{\bm{\xi}}\max_{\bm{\psi}}\left\\{I_{\rm
MINE}(\bm{\xi},\bm{\psi})\right\\}.$ (8)

The optimal design $\bm{\xi}^{*}$ can be obtained by maximizing the MI
estimator in Eq. (7) through a joint gradient-based algorithm (Kleinegesse and
Gutmann, 2020; Foster et al., 2020) or a separate gradient-free updating
scheme (Foster et al., 2019; Kleinegesse and Gutmann, 2020) of $\bm{\xi}$ and
$\bm{\psi}$. The accuracy, efficiency and robustness for estimating and
optimizing MI therefore become very important to the BED tasks. Unfortunately,
$I_{\rm MINE}$ exhibits variance that could grow exponentially with the ground
truth MI and leads to poor bias-variance trade-offs in practice (Song and
Ermon, 2019). The high variance of $I_{\rm MINE}$ may weaken the robustness
the final optimal design $\bm{\xi}^{*}$ and thus cause a higher variance of
the posterior distributions. We propose to use a smoothed mutual information
lower-bound estimator $I_{\rm SMILE}$ with hyperparameter $\tau$ (Song and
Ermon, 2019)

$\displaystyle I_{\rm SMILE}(\bm{\xi},\bm{\psi})$ (9) $\displaystyle=$
$\displaystyle\;\mathbb{E}_{p(\bm{\theta},\bm{y}|\bm{\xi})}[\mathcal{T}_{\bm{\psi}}(\bm{\theta},\bm{y})]$
$\displaystyle-\log\left(\mathbb{E}_{p(\bm{\theta})p(\bm{y}|\bm{\xi})}[{\rm
clip}(e^{\mathcal{T}_{\bm{\psi}}(\bm{\theta},\bm{y})},e^{-\tau},e^{\tau})]\right),$

where the clip function is defined as

$\textup{clip}(u,v,w)=\max(\min(u,w),v).$ (10)

The choice of $\tau$ affects the bias-variance trade-off: when
$\tau\rightarrow\infty$, $I_{\rm SMILE}$ converges to $I_{\rm MINE}$; with a
smaller $\tau$, the variance is reduced at the cost of increasing bias (Song
and Ermon, 2019). The improved MI estimation via variance reduction techniques
is benefit to the optimizing process in Eq. (8) and thus yields a robust
optimal design $\bm{\xi}^{*}$.

###  3.2 Stochastic Gradient Approximate

In the context of BED for implicit models, the pathwise gradients of the MI
lower bound are unavailable. A potential choice for approximating the gradient
is to use Gaussian Smoothing (GS) (Nesterov and Spokoiny, 2017) method to make
the function smooth. The smoothed loss is defined by

$f_{\sigma}(\bm{\xi})=\mathbb{E}_{\bm{\epsilon}\sim\mathcal{N}(0,\mathbf{I}_{d})}\left[f(\bm{\xi}+\sigma\bm{\epsilon})\right],$
(11)

where $\mathcal{N}(0,\mathbf{I}_{d})$ is the $d$-dimensional standard Gaussian
distribution, and $\sigma>0$ is the smoothing radius. The standard GS
represents the $\nabla f_{\sigma}(\bm{\xi})$ as an $d$-dimensional integral
and estimate it by drawing $M$ random samples
$\\{\bm{\epsilon}_{i}\\}_{i=1}^{M}$ from $\mathcal{N}(0,\mathbf{I}_{d})$,
i.e.,

$\displaystyle\nabla f_{\sigma}(\bm{\xi})$
$\displaystyle=\frac{1}{\sigma}\mathbb{E}_{\bm{\epsilon}\sim\mathcal{N}(0,\mathbf{I}_{d})}\left[f(\bm{\xi}+\sigma\bm{\epsilon})\,\bm{\epsilon}\right]$
(12)
$\displaystyle\approx\frac{1}{M\sigma}\sum_{i=1}^{M}f(\bm{\xi}+\sigma\bm{\epsilon}_{i})\bm{\epsilon}_{i}.$

This method using random MC sampling is also called Evolution Strategies (ES)
which have been successfully applied in a variety of high-dimensional
optimization problems, including RL tasks (Salimans et al., 2017; Choromanski
et al., 2018) due to the advantages in scalability and parallelization
capability. However, gradient estimator of ES tends to have a higher variance
for high-dimensional space such that it requires a large number of samples to
be robust (Nesterov and Spokoiny, 2017). To reduce the variance of the
gradient estimator and incorporate with the smoothed MI estimator $I_{\rm
SMILE}$, we propose to use the Guided ES (GES) algorithm that leverages the
recent advances in variance-reduced ES (Maheswaranathan et al., 2018).
Specifically, GES generates a subspace by keeping track of the previous $k$
surrogate gradients during optimization, and leverage this prior information
by changing the distribution of $\bm{\epsilon}_{i}$ in Eq. (12) to
$\mathcal{N}(0,\bm{\Sigma})$ with

$\bm{\Sigma}=(\alpha/n)\cdot\mathbf{I}_{n}+(1-\alpha)/k\cdot UU^{T},$ (13)

where $k$ and $n$ are the subspace and parameter dimensions respectively, $U$
denotes an $n\times k$ orthonormal basis for the subspace, and $\alpha$ is a
hyperparameter that trades off variance between the subspace and full
parameter space. This improved search distribution allows a low-variance
estimate of the descent direction $\nabla f_{\sigma}^{G}(\bm{\xi})$, which can
then be passed to a stochastic gradient ascent (SGA) optimizer. The GES
algorithm is further improved by utilizing historical estimated gradients to
build a low-dimensional subspace for sampling search directions and adaptively
update the importance of this subspace through a exploitation and exploration
trade-off (Liu et al., 2020). To this end, the variance of the GES gradient
estimator $\hat{\nabla}f_{\sigma}^{G}(\bm{\xi})$ is low and it can be
naturally incorporated with $I_{\rm SMILE}$ to find a robust optimal design
$\bm{\xi}^{*}$.

###  3.3 SAGABED Algorithm for Implicit Models

To address the grand challenges in high-dimensional BED for implicit models
without gradients, we propose a novel and scalable approach, named by SAGABED,
which integrates a smoothed MI lower bound and an approximate gradient
estimator. The details of SAGABED algorithm is provided in Algorithm 1.

Algorithm 1 : The SAGABED algorithm

1: Require: neural network architectures, learning rates $\ell_{\psi}$ and
$\ell_{\xi}$, $\tau$ in $I_{\rm SMILE}$, total prior samples $n$, total
iterations $T$, implicit model $\mathcal{M}$

2: Process:

3: Initialize a design $\bm{\xi}_{0}$ by random sampling

4: Initialize neural network parameter $\bm{\psi}_{0}$

5: for $t=0:T-1$ do

6: Draw $n$ samples from the prior distribution of the model parameters
$\bm{\theta}$: $\bm{\theta}^{(1)},...,\bm{\theta}^{(n)}\sim p(\bm{\theta})$

7: Compute the corresponding data samples $\bm{y}^{(i)}$, $i=1,...,n$ using
the current design $\bm{\xi}_{t}$ and an implicit model $\mathcal{M}$

8: Evaluate the smoothed MI lower bound $I_{\rm SMILE}$ by Eq. (9) at the
current design $\bm{\xi}_{t}$ and network parameters $\bm{\psi}_{t}$

9: Compute the approximate gradient estimator $\nabla_{\bm{\xi}}^{*}I_{\rm
SMILE}(\bm{\xi}_{t},\bm{\psi}_{t})$ using the GES algorithm

10: Evaluate the gradient of the $I_{\rm SMILE}$ with respect to the network
parameters $\nabla_{\bm{\psi}}I_{\rm SMILE}(\bm{\xi},\bm{\psi})$

11: Update design $\bm{\xi}_{t}$ via gradient ascent:
$\bm{\xi}_{t+1}=\bm{\xi}_{t}+\ell_{\xi}\nabla_{\bm{\xi}}^{*}I_{\rm
SMILE}(\bm{\xi}_{t},\bm{\psi}_{t})$

12: Update neural network parameters $\bm{\psi}_{t}$ via gradient ascent:
$\bm{\psi}_{t+1}=\bm{\psi}_{t}+\ell_{\psi}\nabla_{\bm{\psi}}I_{\rm
SMILE}(\bm{\xi}_{t},\bm{\psi}_{t})$

13: end for

In the following, we discuss some important features of the SAGABED approach,
specifically for high-dimensional design problems.

  * •

Unified framework vs two-stage framework
Without the requirement of pathwise gradients for implicit models, we utilize
the stochastic approximate gradients and construct a unified framework that
allows the design process to be performed by a simultaneous optimization with
respect to both the variational and design parameters. The existing two-stage
framework that builds a pointwise MI estimator before feeding this estimator
to an outer-level optimizer is often computationally intensive.

  * •

Scalability, portability, and parallelization
we proposes the stochastic approximate gradient ascent procedure that
naturally avoids to the scalability issue in gradient-free methods including
BO. The proposed framework can be easily incorporated with other MI estimators
and implicit models because we only need the forward simulation value to
approximate the gradient using guided ES algorithm. Compared to gradient-based
methods, our computational cost is slightly higher but we benefit from the
parallelization capability of Guided ES methods and thus improve the
computational efficiency.

  * •

Robust estimation with a lower variance
The smoothed MI lower bound used in this method allows us to perform a robust
MI estimation and optimization for Bayesian experimental design. The resulting
low variance of the optimal design and posterior samples enable a more
preciously estimate of the model parameters.

After determining the optimal design $\bm{\xi}^{*}$ by maximizing the MI lower
bound, we can obtain an estimate of the posterior
$p(\bm{\theta}|\bm{y},\bm{\xi}^{*})$ given the learned neural network
$\mathcal{T}_{\bm{\psi}^{*}}(\bm{\theta},\bm{y})$ and prior distribution

$p(\bm{\theta}|\bm{y},\bm{\xi})={\rm
clip}(e^{\mathcal{T}_{\bm{\psi}}(\bm{\theta},\bm{y})-1},e^{-\tau},e^{\tau})p(\bm{\theta}).$
(14)

The relationship in Eq. (14) allows to easily generate posterior samples
$\bm{\theta}_{i}\sim p(\bm{\theta}|\bm{y},\bm{\xi}^{*})$ using Markov chain
Monte Carlo (MCMC) sampling or categorical sampling (Kleinegesse and Gutmann,
2020) since the posterior density can be quickly evaluated via Eq. (14).

##  4 Experiments

We here demonstrate our method in three examples: a linear model with several
noises, a real-world pharmacokinetic model (Ryan et al., 2014), and a
scientific quantum control problem (McMichael, 2020). We compare SAGABED
against two baselines: the two-stage framework using BO and the gradient-based
framework (if gradient is available). The first two have sampling path
gradients, allowing us to compare with gradient-based methods, while the third
one does not.

###  4.1  Noisy Linear Regression

We first show our approach through a classical linear model with noisy sources
that has been used by Kleinegesse and Gutmann (2020) and Foster et al. (2020).
We assume the model to be

$\bm{y}=\theta_{1}\bm{1}+\theta_{2}\bm{\xi}+\bm{\epsilon}+\bm{\nu},$ (15)

where $\bm{y}$ is response variables,
$\bm{\theta}=[\theta_{1},\theta_{2}]^{T}$ are model parameters,
$\bm{\epsilon}\sim\mathcal{N}(0,1)$ and $\bm{\nu}\sim\Gamma(2,2)$ are noise
terms. The design problem is to make $D$ measurements to better estimate the
model parameters $\bm{\theta}$ by constructing a design vector
$\bm{\xi}=[\xi_{1},...,\xi_{D}]^{T}$ which consists of individual experimental
design. Using the linear model in Eq (15), we can obtain the corresponding
independent measurement $\bm{y}_{i}$, which gives a data vector
$\bm{y}=[y_{1},...,y_{D}]^{T}$.

This toy example enables us to use numerical integration to approximate the
posterior and MI but here we assume it is an implicit model for testing our
method. Since the sampling path is given by Eq (15), we can easily compute the
pathwise gradients with respect to the designs
$\nabla_{\bm{\xi}}I(\bm{\xi},\bm{\psi})$ in Kleinegesse and Gutmann (2020) and
thus find the optimal designs $\bm{\xi}^{*}=[\xi_{1}^{*},...,\xi_{D}^{*}]$
using gradient-based approach, which is chosen as a baseline for comparison.

We start from a simple case with only one measurement, i.e., $D=1$. The
initial design is randomly drawn from the design domain $\bm{\xi}\in[-10,10]$
and 10,000 samples of model parameters $\bm{\theta}$ are generated from the
the prior distribution $p(\bm{\theta})=\mathcal{N}(0,3^{2})$. To estimate the
MI lower bound, we use a neural network
$\mathcal{T}_{\bm{\psi}}(\bm{\theta},\bm{y})$ with one layer of 100 neurons,
with a ReLU activation function, as well as Adam optimizer with learning rates
$\ell_{\psi}=10^{-4}$ and $\ell_{\xi}=10^{-2}$. The choice of $\tau$ in
$I_{\rm SMILE}$ affects the bias-variance trade-off, so we compared
$\tau=1,5,10$ and selected the $\tau=5$ with a smaller variance but an
acceptable bias for all examples.

Figure 1 shows the MI lower bound as a function of neural network training
epochs. Except for the proposed SAGABED method that integrates $I_{\rm SMILE}$
and Guided ES algorithm, there are two baselines for comparison: $I_{\rm
MINE}$ with SGD (blue curve) and $I_{\rm MINE}$ with BO (green curve). When
$D=1$, three methods show a fast convergence of the MI lower bound that is
around 2.6, close to a reference MI value computed by the nested MC method and
the computing details can be found in the supplement material. The final
optimal design found by three methods are all at the boundary, i.e.,
$\xi^{*}=$ 10 or -10, which is intuitive due to a larger signal-to-noise ratio
at designs for the linear model.

For more complex cases, e.g., high-dimensional design problems, we adjust the
neural network architectures by hyperparameter optimization and determine to
use one layer with 150 neurons for $D=10$ and five layers with 50 neurons in
each layer for $D=$50 and 100. We use the same activate function, learning
rate, and optimizer as before. As shown in Figure 1, the final MI lower bound
estimated by SAGABED has a bias to the reference MI values but negligibly
small compared to the baselines. The gradient-based method shows a comparative
performance as SAGABED but it is not stable and tends to collapse during
training, e.g., $D=100$ in Figure 1, in the meantime, it might be difficult to
escape from a local minimum, which may yield suboptimal results, particularly
in high-dimensional design space. The method using BO is limited by its
computational and scalable issues, and shows a large bias compared to the
reference MI value at $D=50$ and 100. We also demonstrate superior performance
of the SAGABED on estimation variance that is much smaller than the baselines
across all tasks shown in Figure 1.

Figure 1: The MI lower bound as a function of neural network training epochs
for $D$=1, 10, 50 and 100 measurements in noisy linear model. The dotted lines
are reference MI values at optimized design $t^{*}$ computed by the nested MC
approach.

After finding the optimal design $\bm{\xi}^{*}$, we can compute the
corresponding data $\bm{y}^{*}$ by performing a real-world experiment; then we
can estimate the model parameters $\bm{\theta}$. Here we assume to have a true
model parameter $\bm{\theta}_{\rm true}=[1,4]$ and then use it to generate
$\bm{y}^{*}$. We can compute the posterior density using the trained neural
network $\mathcal{T}_{\bm{\psi}^{*}}(\bm{\theta},\bm{y})$ and Eq. (14). We
therefore obtain the posterior samples from MCMC sampling and use them to
estimate the model parameters. Table 1 shows the estimating mean and standard
deviation of the posterior samples of the model parameters
$\bm{\theta}=(\hat{\theta}_{1},\hat{\theta}_{2})$. The estimating error in
$D=1$ is very large because only one measurement is naturally infeasible for
accurate estimation. As more measurements are taken, the posterior
distribution is narrower and more accurate for $D=50$ than for $D=10$.
However, $I_{\rm MINE}$ with BO due to the scalability issue and computational
challenge still shows a relatively high variance at $D=100$, while SAGABED and
gradient-based method estimate the parameters more precisely and accurately.
The resulting difference in parameter estimation exactly maps the difference
in the MI lower bound estimator.

Table 1: Estimating mean and standard deviation of the posterior samples of the model parameters $\bm{\theta}$ using optimal designs $\bm{d}^{*}$ and real data observation $\bm{y}^{*}$ (use $\bm{\theta}_{\textup{true}}=[1,4]$ to generate $\bm{y}^{*}$) Method | D=1 | D=10 | D=50 | D=100
---|---|---|---|---
| $\hat{\theta}_{1}$ | $\hat{\theta}_{2}$ | $\hat{\theta}_{1}$ | $\hat{\theta}_{2}$ | $\hat{\theta}_{1}$ | $\hat{\theta}_{2}$ | $\hat{\theta}_{1}$ | $\hat{\theta}_{2}$
MINE+SGD | -1.39$\pm$2.54 | 6.03$\pm$0.93 | 0.51$\pm$0.44 | 2.99$\pm$0.67 | 1.20$\pm$0.18 | 3.79$\pm$0.23 | 0.97$\pm$0.05 | 4.04$\pm$0.04
MINE+BO | -1.42$\pm$0.81 | 2.98$\pm$1.19 | 1.22$\pm$0.58 | 4.93$\pm$0.91 | 0.71$\pm$0.25 | 3.66$\pm$0.40 | 1.35$\pm$0.21 | 4.79$\pm$0.26
SMILE+GES | 2.76$\pm$1.36 | 5.74$\pm$3.08 | 0.83$\pm$0.56 | 4.69$\pm$0.58 | 1.11$\pm$0.13 | 4.25$\pm$0.19 | 1.02$\pm$0.04 | 3.98$\pm$0.03

###  4.2  Pharmacokinetic Model

Next we illustrate our approach for the design problem of determining the
optimal measurements of blood sampling times for a pharmacokinetic (PK) study,
which basically involves the administration of a drug to an individual and
then investigate the drug distribution, absorption and elimination at certain
times to analyze the underlying kinetics. The computational model introduced
by Ryan et al. (2014) is used to simulate the drug concentration at a specific
time. There are three governed model parameters: the volume of distribution
$V$, the absorption rate $k_{a}$ and the elimination rate $k_{e}$. The PK
model is defined by a latent variable for the drug concentration at design
time $t$,

$z(t)=\frac{D_{v}}{V}\frac{k_{a}}{k_{a}-k_{e}}\left(e^{-k_{e}t}-e^{-k_{a}t}\right)(1+\epsilon_{1t})+\epsilon_{2t},$
(16)

where $k_{a}>k_{e}$ is a constraint, $D_{v}=400$ is a single fixed dose at the
beginning of the experiment, $t\in[0,24]$ is the design domain,
$\epsilon_{1t}\sim\mathcal{N}(0,0.01)$ and
$\epsilon_{2t}\sim\mathcal{N}(0,0.1)$ are multiplicative noisy sources that
are often observed in the drug data. Note that Eq. (16) is a simplified model
and we assumes that the design parameters $\bm{\theta}=(V,k_{a},k_{e})$ are
the same for the group of patients. As suggested by Ryan et al. (2014), the
prior distribution for $\bm{\theta}$ is given by

$\log(\bm{\theta})\sim\mathcal{N}\left(\left[\begin{matrix}\log 20\\\ \log
1\\\ \log 0.1\end{matrix}\right],\left[\begin{matrix}0.05&0&0\\\ 0&0.05&0\\\
0&0&0.05\end{matrix}\right]\right)$ (17)

In this example, we assume to take $D$ measurements that mean we have a group
of $D$ patients but take only one blood sample from each patient at time $t$.
Thus the design parameters are $\bm{\xi}=[t_{1},...,t_{D}]^{T}$ and the
corresponding observations are $\bm{y}=[y_{1},...,y_{D}]^{T}$. It is noted
that the PK model in Eq. (16) allows us to analytically derive the pathwise
gradients in terms of the sampling path (Kleinegesse and Gutmann, 2020)

$\frac{\partial z}{\partial
t}=\frac{D}{V}\frac{k_{a}}{k_{a}-k_{e}}\left(k_{a}e^{-k_{a}t}-k_{e}e^{-k_{e}t}\right)(1+\epsilon_{1t}),$
(18)

The available gradients in Eq (18) enable to find the optimal design of blood
sampling times using gradient-based methods, which can be considered as a
baseline method.

Instead of the simpler case by setting $D=1$ in the linear example, we
initially set $D=10$ here and randomly generate 10 design time in $t\in[0,24]$
hours as the initial state. Then 10,000 model parameters samples
$\bm{\theta}^{(i)}$ are drawn from the prior distribution in Eq. (17). Using
these samples, we can collect the corresponding data samples $\bm{y}^{(i)}$
using the current design $t$. For neural network model
$\mathcal{T}_{\bm{\psi}}(\bm{\theta},\bm{y})$, we use one hidden layer of 100
neurons with ReLU activation function as well as Adam optimizer with learning
rates $\ell_{\psi}=10^{-4}$ and $\ell_{\xi}=10^{-2}$. Hyperparameters tuning
details can be found in the supplementary material.

Figure 2: MI lower bound as a function of neural network training epochs for
different number of measurements: $D$=10, 50, 100 and 500 in the study of PK
model.

Figure 2 shows the MI lower bound as a function of the neural network training
epochs in different sampling times $D=10,50,100$ and $500$ for the PK study.
We compare the SAGABED (refers SMILE+GES in Figure 2) with the baselines that
use MINE with SGD and BO optimizer respectively. When $D=10$, the MI lower
bound using SAGABED converges to around 3.6 that is close to a reference MI
value at the optimal design $t^{*}$. The baselines show a comparative
performance compared to SAGABED. As more measurements $D=50$ are performed,
all three methods converge to a higher MI lower bound; this is as expected
since more information are obtained from more data of model parameters.
SAGABED shows a negligibly small bias compared to the reference MI value and
has a low variance, which outperforms the other two baselines.

We now focus on the high dimensional cases, e.g., $D=100$ and 500, that are
more challenging settings. Here we use a 5 layered neural network with 50
neurons for each layer instead of the neural network architectures in the low
dimensional cases. While SAGABED shows a larger bias of the final MI lower
bound compared to the reference value due to increasing dimensions, it is
still reasonably acceptable and higher than the baselines. In the high
dimensional setting, the gradient-based method increases fast initially but it
is easy to be trapped into local minima; the BO method has no this issue but
does not scale to such high dimensionality. Also, the MI lower bound estimated
by $I_{\rm MINE}$ shows a high variance whereas $I_{\rm SMILE}$ has a lower
variance across all different tasks. That is an important feature that may
potentially reduce the variance of final optimized design and also avoid
additional efforts on optimizing hyperparameters, e.g., learning rates.

Figure 3: Marginal posterior distributions of the model parameters for $D=10$
(solid curves) and $D=50$ (dotted curves). Kernel density estimate is used to
approximate the posterior samples. SAGABED shows a much narrower posterior
distribution of the model parameters.

Once the optimal design $t^{*}$ is obtained, we can conduct the experiment to
generate the real-world data $\bm{y}^{*}$ by assuming a true model parameter
$\theta_{\rm true}=[20,2.0,0.2]$. Then we compute the posterior density
$p(\bm{\theta}|t,^{*},y^{*})$ using the learned neural network
$\mathcal{T}_{\bm{\psi}^{*}}(\bm{\theta},\bm{y})$, and obtain the posterior
samples using MCMC sampling. Figure 3 shows the marginal posterior
distribution of each model parameters. Using $I_{\rm SMILE}$, we achieve a
more accurate and much narrower estimate of all the model parameters with a
lower variance than the baseline methods.

###  4.3 Tuning for Quantum Control

Typically, the reliable capability to manipulate qbit states is critical to
quantum technology. For instance, radio-frequency pulses can be used to change
of state of spin-up and spin-down, and patients benefit from MRI scans that
use these approaches. In this example, we aim to conduct SAGABED to simulate a
tuning process so that we can actively control the desired duration and
frequency of pulses to flip electron spins in a reliable scheme. Specifically,
the model parameters $\bm{\theta}=(\theta_{1},\theta_{2})$ are Rabi frequency
and true center of resonance relative to the reference frequency, and the
design parameters $\bm{\xi}=(t,\Delta f)$ are the duration time of microwave
pulse and the detuning relative to a reference frequency. We here use the
physical simulator developed by McMichael (2020) to simulate the tuning
process for quantum control. The physical simulator used here can be
interpreted as an implicit model where the gradient can not be computed
exactly. We assume that our measurements are subject to Gaussian noise,
$\mathcal{N}(0,1)$. The simulation domain is discretized by 101$\times$101
grids and a finer grid would lead to a more accurate prediction but at a
higher computational cost. The design domain is set to $t\in[0,1]$ and $\Delta
f\ \in[-10,10]$ and we therefore assume uniform prior distributions
$p(t)=U(0,1)$ and $p(\Delta_{f})=U(-10,10)$.

We use a one-layer neural network with 200 neurons and a ReLU activate
function as well as Adam optimizer with learning rates $\ell_{\psi}=10^{-3}$
and $\ell_{\xi}=10^{-2}$. We still start from a simple case, i.e., only one
measurement $N=1$, and then gradually increase the number of measurements to
$N=5,10,50,100$, and 500 finally. We randomly draw 10,000 uniform prior
samples $\bm{\theta}^{(i)}$ and then simulate 10,000 corresponding samples of
the outcome $\bm{y}$ that is the Rabi counts. The optimal design results are
shown in the left column in Figure 4. Once we have done the neural network
training, we assume the true model parameters $\bm{\theta}_{\rm
true}=[3.85,1.67]$ as the true Rabi and detuning frequency. As similar to the
previous examples, we can conduct a real-world measurement of the tuning
control process given the optimal design $\bm{\xi}^{*}$ and learned neural
network $\mathcal{T}_{\bm{\psi}^{*}}(\bm{\theta},\bm{y})$. As a result, the
posterior samples are obtained and illustrated by the right column in Figure
4.

Our main purpose in this example is to demonstrate the superior performance of
SAGABED on the robustness of the optimal design and their effect on the
variance of the posterior samples. As a comparison, BED with the BO method
shows a similar performance on the low-dimensional cases but a high variance
when the dimensionality increases to 100. Even more design are collected, the
posterior distribution in the BO method is difficult to be narrowed, as shown
in Figure 4. On the contrary, SAGABED demonstrates a faster narrowing rate to
the true model parameters and it displays a much lower variance of the
posterior samples in the high-dimensional design problems. This also
illustrates the significant advantages of SAGABED in addressing the
scalability challenges compared to the BO method.

Figure 4: Performance comparison between SAGABED and BO method for tuning
quantum pulse. The contour image (left column) shows the model photon counts
for optically detected spin manipulation for pulse duration and amounts of
detuning from the spin’s natural resonance frequency. The right column
displays the evolution of the posterior distribution with the number of
designed measurements. The red points are the true model parameters.

##  5 Conclusion

In this paper, we develop a general unified framework that utilizes the
stochastic approximate gradient for BED with implicit models. Without the
requirement or assumption of pathwise gradients, our approach allows the
optimization to be carried out by stochastic gradient ascent algorithms and
therefore scaled to substantial high dimensional design problems. Several
experiments demonstrate that our approach outperforms the baseline methods,
and significantly improves the scalability of BED in high dimensional
settings.

The future work will focus on the extension of our proposed framework to
sequential Bayesian experimental design (SBED) that involves an iterative
update of the prior and posterior distribution. We plan to improve the
computational efficiency of the SAGABED, which may play a more critical role
in the SBED scenario. We are also interested in improving the stochastic
approximate gradient methods and enhancing the global exploration capability
to escape the sub-optimal design caused by the issue of local minima in
current SGA algorithms, e.g., standard and guided ES (Zhang et al., 2020,
2021).

#### Acknowledgements

This work was supported by the U.S. Department of Energy, Office of Science,
Office of Advanced Scientific Computing Research (ASCR), Applied Mathematics
program; and by the Artificial Intelligence Initiative at the Oak Ridge
National Laboratory (ORNL). This work used resources of the Oak Ridge
Leadership Computing Facility, which is supported by the Office of Science of
the U.S. Department of Energy under Contract No. DE-AC05-00OR22725.

## References

  * Alemi et al. (2016) Alexander A Alemi, Ian Fischer, Joshua V Dillon, and Kevin Murphy.  Deep variational information bottleneck.  _arXiv preprint arXiv:1612.00410_ , 2016.
  * Amzal et al. (2006) Billy Amzal, Frédéric Y Bois, Eric Parent, and Christian P Robert.  Bayesian-optimal design via interacting particle systems.  _Journal of the American Statistical association_ , 101(474):773–785, 2006.
  * Ao and Li (2020) Ziqiao Ao and Jinglai Li.  An approximate $\\{$KLD$\\}$ based experimental design for models with intractable likelihoods.  _arXiv preprint arXiv:2004.00715_ , 2020.
  * Barber and Agakov (2003) David Barber and Felix V Agakov.  The im algorithm: a variational approach to information maximization.  In _Advances in neural information processing systems_ , page None, 2003.
  * Belghazi et al. (2018) Mohamed Ishmael Belghazi, Aristide Baratin, Sai Rajeswar, Sherjil Ozair, Yoshua Bengio, Aaron Courville, and R Devon Hjelm.  Mine: mutual information neural estimation.  _arXiv preprint arXiv:1801.04062_ , 2018.
  * Chaloner and Verdinelli (1995) Kathryn Chaloner and Isabella Verdinelli.  Bayesian experimental design: A review.  _Statistical Science_ , pages 273–304, 1995.
  * Chen et al. (2016) Xi Chen, Yan Duan, Rein Houthooft, John Schulman, Ilya Sutskever, and Pieter Abbeel.  Infogan: Interpretable representation learning by information maximizing generative adversarial nets.  In _Advances in neural information processing systems_ , pages 2172–2180, 2016.
  * Choromanski et al. (2018) Krzysztof Choromanski, Mark Rowland, Vikas Sindhwani, Richard E Turner, and Adrian Weller.  Structured evolution with compact architectures for scalable policy optimization.  _arXiv preprint arXiv:1804.02395_ , 2018.
  * Drovandi et al. (2018) Christopher C Drovandi, Minh-Ngoc Tran, et al.  Improving the efficiency of fully Bayesian optimal design of experiments using randomised quasi-Monte Carlo.  _Bayesian Analysis_ , 13(1):139–162, 2018.
  * Eriksson et al. (2019) David Eriksson, Michael Pearce, Jacob Gardner, Ryan D Turner, and Matthias Poloczek.  Scalable global optimization via local Bayesian optimization.  In _Advances in Neural Information Processing Systems_ , pages 5496–5507, 2019.
  * Foster et al. (2019) Adam Foster, Martin Jankowiak, Elias Bingham, Paul Horsfall, Yee Whye Teh, Thomas Rainforth, and Noah Goodman.  Variational Bayesian optimal experimental design.  In _Advances in Neural Information Processing Systems_ , pages 14036–14047, 2019.
  * Foster et al. (2020) Adam Foster, Martin Jankowiak, Matthew O’Meara, Yee Whye Teh, and Tom Rainforth.  A unified stochastic gradient approach to designing Bayesian-optimal experiments.  In _International Conference on Artificial Intelligence and Statistics_ , pages 2959–2969. PMLR, 2020.
  * Harbisher et al. (2019) Sophie Harbisher, Colin S Gillespie, and Dennis Prangle.  Bayesian optimal design using stochastic gradient optimisation and Fisher information gain.  _arXiv preprint arXiv:1904.05703_ , 2019.
  * Kleinegesse and Gutmann (2019) Steven Kleinegesse and Michael U Gutmann.  Efficient Bayesian experimental design for implicit models.  In _The 22nd International Conference on Artificial Intelligence and Statistics_ , pages 476–485, 2019.
  * Kleinegesse and Gutmann (2020) Steven Kleinegesse and Michael U Gutmann.  Bayesian experimental design for implicit models by mutual information neural estimation.  _arXiv preprint arXiv:2002.08129_ , 2020.
  * Kleinegesse et al. (2020) Steven Kleinegesse, Christopher Drovandi, and Michael U Gutmann.  Sequential Bayesian experimental design for implicit models via mutual information.  _arXiv preprint arXiv:2003.09379_ , 2020.
  * Li et al. (2018) Cheng Li, Sunil Gupta, Santu Rana, Vu Nguyen, Svetha Venkatesh, and Alistair Shilton.  High dimensional Bayesian optimization using dropout.  _arXiv preprint arXiv:1802.05400_ , 2018.
  * Li et al. (2016) Chun-Liang Li, Kirthevasan Kandasamy, Barnabás Póczos, and Jeff Schneider.  High dimensional Bayesian optimization via restricted projection pursuit models.  In _Artificial Intelligence and Statistics_ , pages 884–892, 2016\.
  * Liu et al. (2020) Fei-Yu Liu, Zi-Niu Li, and Chao Qian.  Self-guided evolution strategies with historical estimated gradients.  _IJCAI_ , 2020.
  * Maheswaranathan et al. (2018) Niru Maheswaranathan, Luke Metz, George Tucker, Dami Choi, and Jascha Sohl-Dickstein.  Guided evolutionary strategies: escaping the curse of dimensionality in random search.  2018\.
  * McAllester and Stratos (2020) David McAllester and Karl Stratos.  Formal limitations on the measurement of mutual information.  In _International Conference on Artificial Intelligence and Statistics_ , pages 875–884, 2020.
  * McMichael (2020) R. D. McMichael.  Optimal Bayesian experimental design, 2020.  URL https://github.com/usnistgov/optbayesexpt.
  * Müller (2005) Peter Müller.  Simulation based optimal design.  _Handbook of Statistics_ , 25:509–518, 2005.
  * Mutny and Krause (2018) Mojmir Mutny and Andreas Krause.  Efficient high dimensional Bayesian optimization with additivity and quadrature fourier features.  In _Advances in Neural Information Processing Systems_ , pages 9005–9016, 2018.
  * Myung et al. (2013) Jay I Myung, Daniel R Cavagnaro, and Mark A Pitt.  A tutorial on adaptive design optimization.  _Journal of mathematical psychology_ , 57(3-4):53–67, 2013.
  * Nesterov and Spokoiny (2017) Yurii Nesterov and Vladimir Spokoiny.  Random gradient-free minimization of convex functions.  _Foundations of Computational Mathematics_ , 17(2):527–566, 2017.
  * Nguyen et al. (2010) XuanLong Nguyen, Martin J Wainwright, and Michael I Jordan.  Estimating divergence functionals and the likelihood ratio by convex risk minimization.  _IEEE Transactions on Information Theory_ , 56(11):5847–5861, 2010.
  * Nowozin et al. (2016) Sebastian Nowozin, Botond Cseke, and Ryota Tomioka.  f-gan: Training generative neural samplers using variational divergence minimization.  In _Advances in neural information processing systems_ , pages 271–279, 2016.
  * Oord et al. (2018) Aaron van den Oord, Yazhe Li, and Oriol Vinyals.  Representation learning with contrastive predictive coding.  _arXiv preprint arXiv:1807.03748_ , 2018.
  * Overstall et al. (2018) Antony Overstall, James McGree, et al.  Bayesian design of experiments for intractable likelihood models using coupled auxiliary models and multivariate emulation.  _Bayesian Analysis_ , 2018.
  * Pathak et al. (2017) Deepak Pathak, Pulkit Agrawal, Alexei A Efros, and Trevor Darrell.  Curiosity-driven exploration by self-supervised prediction.  In _Proceedings of the IEEE Conference on Computer Vision and Pattern Recognition Workshops_ , pages 16–17, 2017.
  * Poole et al. (2019) Ben Poole, Sherjil Ozair, Aaron van den Oord, Alexander A Alemi, and George Tucker.  On variational bounds of mutual information.  _arXiv preprint arXiv:1905.06922_ , 2019.
  * Rainforth (2017) Thomas William Gamlen Rainforth.  _Automating inference, learning, and design using probabilistic programming_.  PhD thesis, University of Oxford, 2017.
  * Rainforth et al. (2018) Tom Rainforth, Rob Cornish, Hongseok Yang, Andrew Warrington, and Frank Wood.  On nesting Monte Carlo estimators.  In _International Conference on Machine Learning_ , pages 4267–4276. PMLR, 2018.
  * Rana et al. (2017) Santu Rana, Cheng Li, Sunil Gupta, Vu Nguyen, and Svetha Venkatesh.  High dimensional Bayesian optimization with elastic Gaussian process.  In _International Conference on Machine Learning_ , pages 2883–2891, 2017.
  * Ryan et al. (2014) Elizabeth G Ryan, Christopher C Drovandi, M Helen Thompson, and Anthony N Pettitt.  Towards Bayesian experimental design for nonlinear models that require a large number of sampling times.  _Computational Statistics & Data Analysis_, 70:45–60, 2014.
  * Salimans et al. (2017) Tim Salimans, Jonathan Ho, Xi Chen, Szymon Sidor, and Ilya Sutskever.  Evolution strategies as a scalable alternative to reinforcement learning.  _arXiv preprint arXiv:1703.03864_ , 2017.
  * Sebastiani and Wynn (2000) Paola Sebastiani and Henry P Wynn.  Maximum entropy sampling and optimal Bayesian experimental design.  _Journal of the Royal Statistical Society: Series B (Statistical Methodology)_ , 62(1):145–157, 2000.
  * Snoek et al. (2012) Jasper Snoek, Hugo Larochelle, and Ryan P Adams.  Practical Bayesian optimization of machine learning algorithms.  In _Advances in neural information processing systems_ , pages 2951–2959, 2012.
  * Snoek et al. (2015) Jasper Snoek, Oren Rippel, Kevin Swersky, Ryan Kiros, Nadathur Satish, Narayanan Sundaram, Mostofa Patwary, Mr Prabhat, and Ryan Adams.  Scalable Bayesian optimization using deep neural networks.  In _International conference on machine learning_ , pages 2171–2180, 2015.
  * Song and Ermon (2019) Jiaming Song and Stefano Ermon.  Understanding the limitations of variational mutual information estimators.  _arXiv preprint arXiv:1910.06222_ , 2019.
  * Wang et al. (2017) Zi Wang, Chengtao Li, Stefanie Jegelka, and Pushmeet Kohli.  Batched high-dimensional Bayesian optimization via structural kernel learning.  _arXiv preprint arXiv:1703.01973_ , 2017.
  * Zhang et al. (2020) Jiaxin Zhang, Hoang Tran, Dan Lu, and Guannan Zhang.  A novel evolution strategy with directional Gaussian smoothing for blackbox optimization.  _arXiv preprint arXiv:2002.03001_ , 2020.
  * Zhang et al. (2021) Jiaxin Zhang, Sirui Bi, and Guannan Zhang.  A directional Gaussian smoothing optimization method for computational inverse design in nanophotonics.  _Materials & Design_, 197:109213, 2021.