text:
\section{Introduction}

Selecting a good design for an experiment can be crucial to extracting useful information and controlling costs.
Applications include medical interventions \citep{Amzal:2006}, epidemic modelling \citep{Cook:2008}, pharmacokinetics \citep{Ryan:2014, Overstall:2017} and ecology \citep{Gillespie:2019}.
In modern applications it is increasingly feasible to take a large number of measurements --
e.g.~placing sensors \citep{Krause:2009} or making observations in a numerical integration problem \citep{Oates:2020} --
or make other complex designs -- e.g.~selecting a time series of chemical input levels to a synthetic biology experimental \citep{Bandiera:2018}.
Therefore finding high dimensional designs is an increasingly relevant task.

We focus on the Bayesian approach to optimal experimental design, which takes into account existing knowledge and uncertainty about the process being studied before the experiment is undertaken.
In this framework an experimenter must select a design.
They then receive some utility based on the outcome of the experiment.
The aim is to select the design optimising expected utility given the experimenter's prior beliefs.

Most utility functions in Bayesian experimental design require posterior calculations, such evaluating the evidence.
Optimisation of the design requires these calculations to be repeated for a large number of potential designs and simulated datasets.
This can be expensive and prohibit scaling up these methods to models with many parameters, or designs with many unknowns to select.
Some sophisticated approaches have been developed, including converting the optimisation problem into Markov chain Monte Carlo on the joint space of designs, parameters, and observations \citep{Muller:1999}, and using variational inference to learn approximate surrogate posteriors \citep{Foster:2019}, but these still effectively solve inference problems and hence remain costly.
(We comment on other related approaches in Section \ref{sec:related}.)

This paper presents approaches which avoid the need for posterior inference.
First we consider a utility function based on the trace of the Fisher information, which is often available in a simple closed form.
\cite{Walker:2016} presents an information theoretic justification for this utility.
One contribution of this paper is to show that it also emerges naturally from the decision theoretic framework of \cite{Bernardo:1979}.
Our derivation is based on judging the quality of a parameter estimate through the Hyv\"arinen score \citep{Hyvarinen:2005}, rather than through the logarithmic score as in \cite{Bernardo:1979}.
We demonstrate that it is straightforward to optimise the resulting expected utility using recent developments in stochastic gradient optimisation \citep{Kingma:2015} and automatic differentiation \citep{Baydin:2017}.
Compared to existing methods, this approach is fast and scales easily to higher dimensional designs.
However, a drawback is that the method can converge to poor local maxima.
We show how a second stage of optimisation similar to that of \cite{Overstall:2017} can often be used to find the overall optimal design.

A more fundamental limitation of the above approach is that it sometimes produces poor designs in practice e.g.~requiring all observations to occur at a single time point.
We provide an explanation:
optimising this expected utility encourages designs giving accurate inference of one linear combination of the parameters but not necessarily others.
Furthermore this utility is not invariant to reparameterisation:
this can alter which parameter combinations are most rewarding to infer accurately.

The main contribution of this paper is to address both these issues by introducing an adversarial approach.
We propose a game theoretic framework in which, as before, the experimenter chooses a design to optimise their expected utility.
Now there is also an adversary who selects a linear transformation of the parameters.
We investigate the optimal designs in this framework under the game theoretic solution concept of subgame perfect equilibrium, and prove they are invariant to reparameterisation.
The presence of the adversary also encourages designs not to neglect informativeness of any parameter combination:
if it does then the adversary can choose a parameterisation concentrating on this weakness.
We show it is possible to find optimal designs in this game theoretic framework using generic \emph{gradient descent ascent} methods, which have been much studied recently in the machine learning literature \citep[e.g][]{Heusel:2017, Jin:2019}.

In the remainder of the paper,
Section \ref{sec:background} presents background on Bayesian experimental design.
Section \ref{sec:DTtheory} presents results on the decision theoretic framework of \cite{Bernardo:1979}, and Section \ref{sec:GTtheory}, on our proposed game theoretic extension.
Proofs and other technical material are presented in the appendices.
Section \ref{sec:opt} discusses details of gradient based optimisation for both approaches.

We illustrate our methods on a simple Poisson model where optimal designs can be derived analytically (Section \ref{sec:illustration}).
We provide a detailed simulation study on a pharmacokinetic model (Section \ref{sec:pk}) showing our adversarial approach produces a sensible design and is at least 10 times faster than competing methods.
We also present a geostatistical regression example, where hundreds of design choices can be optimised in under a minute (Section \ref{sec:geostat}).
Code for these examples is available at
\url{https://github.com/dennisprangle/AdversarialDesignCode}.
All examples were run on a desktop PC with 12 CPU cores.
Our conclusion, Section \ref{sec:discussion}, summarises our findings and recommendations to implement our methods.
It also discusses limitations of our work, including the requirement of tractable Fisher information, and future research directions to address these.

\subsection{Related work} \label{sec:related}

As discussed above, our theoretical contribution builds on
\cite{Bernardo:1979}'s decision theoretic framework for Bayesian experimental design,
and also on \cite{Walker:2016}'s justification for the use of a utility based on the trace of the Fisher information.
We also discuss \cite{Overstall:2020}'s comments on this utility later.

Gradient-based optimisation methods for experimental design have been explored previously.
\cite{Pronzato:1985} optimise the expected determinant of the Fisher information using analytically derived gradients.
\cite{Huan:2013, Huan:2014} optimise expected Shannon information using gradients
(either derived analytically or based on finite differences)
for a biased numerical approximation to the utility.
\cite{Foster:2019} use a first stage of variational inference to learn an approximate posterior, and then use this to produce a surrogate expected utility function which is optimised by gradient-based methods in a second stage.

\cite{Foster:2020} and \cite{Kleinegesse:2020} propose jointly optimising a design and some tuning parameters defining a lower bound on the expected Shannon information gain utility (this utility is discussed in Section \ref{sec:utilities}).
This can be framed as a single optimisation problem, allowing for easy implementation using stochastic gradient or Bayesian optimisation methods.
Producing a tight lower bound is equivalent to, or closely related to, being able to perform exact Bayesian inference for the optimal design.

\cite{Overstall:2017} propose an alternative coordinate ascent approach which loops over the components of the design, updating each in turn.
To perform an update, designs
are selected from the one dimensional search space (in which the current component only is updated)
and Monte Carlo estimates of expected utility calculated.
A Gaussian process is fitted to the expected utility estimates and used to propose an improved value for the design component under consideration.
This is accepted or rejected based on a Bayesian test of whether it improves expected utility, using a large number of simulations under the current and proposed designs.

Many other algorithms have been proposed for Bayesian experimental design.
An influential method of \cite{Muller:1999} performs optimal design using Markov Chain Monte Carlo.
\citealp{Amzal:2006} and \citealp{Kuck:2006} extend this approach to use sequential Monte Carlo.
In other approaches \cite{Ryan:2014} look at high dimensional designs with a low dimensional parametric form,
\cite{Price:2018} use evolutionary algorithms, and
\cite{Gillespie:2019} perform a sophisticated search on a discrete grid of designs.

\section{Background} \label{sec:background}

Optimal experimental design concerns the following problem.
An experimenter must select a design $\tau$.
The experiment produces data $y$ with likelihood $f(y|\theta; \tau)$,
where $\theta$ is a vector of unknown parameters.
The goal is to select the design which optimises some notion of the quality of the experiment,
typically based on its informativeness and its cost.
We mostly, but not exclusively, concentrate on the case where
$\tau \in \mathcal{T} \subseteq \mathbb{R}^d$
and $\theta \in \Theta \subseteq \mathbb{R}^p$
for some closed sets (under the Euclidean metric) $\mathcal{T}$ and $\Theta$,
and where $y$ is a vector of $n$ observations.
Designs $\tau$ of this form often represent times or locations for measurements to be taken.
In this case $\tau$ can be seen as a set of \emph{design points}, $\tau_1, \tau_2, \ldots, \tau_d$.

This section reviews relevant details of the statistical background to our paper.
First Sections \ref{sec:Bayes} and \ref{sec:FIM} give some necessary definitions on Bayesian statistics and Fisher information.
Section \ref{sec:BED} describes Bayesian experimental design, which is based on maximising expected utility.
Finally, Section \ref{sec:utilities} summarises some common utility functions.

\subsection{Bayesian framework} \label{sec:Bayes}

We work in the Bayesian framework and introduce a prior density $\pi(\theta)$ for $\theta$.
We will often make use of the posterior density and the prior predictive density (or \emph{evidence}) for $y$.
In our setting both depend on the experimental design $\tau$,
\begin{align}
\pi(\theta | y; \tau) &= \pi(\theta) f(y | \theta; \tau) / \pi(y; \tau), \label{eq:posterior} \\
\pi(y; \tau) &= \E_{\theta \sim \pi(\theta)}[f(y | \theta; \tau)] \label{eq:prior predictive}.
\end{align}
The prior and model define a joint density\footnote{
For discrete $y$ a density with respect to a product of Lebesgue and counting measures can be used.},
\begin{equation} \label{eq:joint}
\pi(\theta, y; \tau) = \pi(\theta) f(y | \theta; \tau)
= \pi(\theta | y; \tau) \pi(y; \tau).
\end{equation}

\subsection{Fisher information} \label{sec:FIM}

We will make frequent use of the Fisher information matrix for $\theta$,
\begin{equation} \label{eq:FisherInfo}
\mathcal{I}_\theta(\theta; \tau) = \E_{y \sim f(y | \theta; \tau)} [ u(y, \theta; \tau) u(y, \theta; \tau)^T ],
\end{equation}
which is based on the score function,
\begin{equation} \label{eq:score}
u(y,\theta;\tau) = \nabla_\theta \log f(y|\theta; \tau).
\end{equation}
The subscript in $\mathcal{I}_\theta$ denotes which variable is used for differentiation in \eqref{eq:score}.
Often this is obvious from the context so the subscript will be dropped.
We will also sometimes omit the dependence on $\tau$ where this is not relevant.
We will focus on models where $u$ and $\mathcal{I}$ are well defined.

\paragraph{Examples}

The Fisher information for exponential family models is based on the variance of the sufficient statistics (see e.g.~\citealp{Lehmann:2006}) which is often available in closed form.

Two examples we use in this paper are:
\begin{itemize}
\item
The Poisson distribution, $y \sim Poisson(\phi)$.
Here $\mathcal{I}(\phi) = 1/\phi$.
\item
The multivariate normal distribution with known variance,
$y \sim N(\mu, \Sigma)$.
Here $\mathcal{I}(\mu) = \Sigma^{-1}$.
\end{itemize}

\paragraph{Reparameterisation}

Consider a probability model with parameters $\theta$,
and a function $\phi(\theta)$
producing an alternative vector of parameters $\phi$
(which may be shorter or longer than $\theta$).
Let $J(\phi)$ be the Jacobian
i.e.~the matrix whose row $i$ column $j$ entry
is $\frac{\partial \phi_i}{\partial \theta_j}$.
Then it follows from \eqref{eq:FisherInfo} that
\begin{equation} \label{eq:reparam}
\mathcal{I}_\theta(\theta) = J(\phi)^T \mathcal{I}_\phi(\phi) J(\phi).
\end{equation}
An application of this result is to
the model $y \sim N(x(\theta,\tau), \Sigma)$.
Using \eqref{eq:reparam},
\begin{equation} \label{eq:mvn_info}
\mathcal{I}_\theta(\theta) = J^T \mathcal{I}_x(x) J = J(x)^T \Sigma^{-1} J(x).
\end{equation}

\subsection{Bayesian experimental design} \label{sec:BED}

The Bayesian approach to experimental design involves selecting a function $\mathcal{U} = \mathcal{U}(\tau, \theta, y)$,
giving the utility of selecting design $\tau$ given observations $y$ and true parameters $\theta$.
(As we shall see, many choices of $\mathcal{U}$ do not depend on all these possible inputs.)
We try to maximise the \emph{expected utility} of $\tau$ i.e. the prior predictive utility
\begin{equation} \label{eq:expected utility}
\mathcal{J}(\tau) = \E_{(\theta,y) \sim \pi(\theta, y; \tau)}[\mathcal{U}(\tau, \theta, y)].
\end{equation}
See \cite{Chaloner:1995}, \cite{Atkinson:2007} and \cite{Ryan:2016} for comprehensive surveys of Bayesian experimental design.

\subsection{Utility functions} \label{sec:utilities}

Ideally a utility function could be specified for each application,
perhaps by eliciting preferences over different $(\tau,\theta,y)$ combinations from the experimenter \citep[e.g.][]{Wolfson:1996}.
However this is rarely feasible in practice.
Instead several generic choices of utility have been proposed.

\paragraph{Shannon information gain (SIG)}

This is a popular and well-motivated utility choice,
\begin{equation} \label{eq:SIG}
\mathcal{U}_{\text{SIG}}(\tau, \theta, y) = \log \pi(\theta | y; \tau) - \log \pi(\theta),
\end{equation}
introduced by \cite{Lindley:1956},
which is particularly relevant later in this paper.
Designs maximising the expectation of $\mathcal{U}_{\text{SIG}}$
have an appealing information theoretic interpretation:
they maximise expected reduction in Shannon entropy from prior to posterior
Furthermore, \cite{Bernardo:1979} gave a decision theoretic derivation of $\mathcal{U}_{\text{SIG}}$ which we recap in Section \ref{sec:DTtheory}.
A helpful property of $\mathcal{U}_{\text{SIG}}$ is that it is \emph{reparameterisation invariant} i.e.~unchanged under a bijective transformation of $\theta$.
Thus the resulting designs are not affected by the choice of parameterisation.

In practice, optimising the expectation of \eqref{eq:SIG} is complicated by the need to evaluate the posterior density.
However, using \eqref{eq:posterior} the utility can be rewritten as
\begin{equation}
\mathcal{U}_{\text{SIG}}(\tau, \theta, y) = \log f(y | \theta; \tau) - \log \pi(y; \tau). \label{eq:SIG2}
\end{equation}
Optimisation now requires evaluation of a posterior summary: the \emph{log evidence} $\log \pi(y; \tau)$.
A common SIG estimate replaces $\pi(y; \tau)$ in \eqref{eq:SIG2} with a Monte Carlo approximation
\begin{equation} \label{eq:SIGMC}
\hat{\pi}(y; \tau) = \frac{1}{L} \sum_{\ell=1}^L f(y | \theta^{(\ell)}; \tau),
\end{equation}
where $\theta^{(\ell)}$ are independent prior samples.
A typical choice of $L$ is 1000 \citep{Overstall:2017},
which makes each utility evaluation somewhat computationally expensive.
Furthermore, a biased estimate of $\mathcal{U}_{\text{SIG}}$ is produced.
More efficient approaches are possible,
as outlined in Section \ref{sec:related}.
However the need for log evidence estimation remains a source of computational expense and approximation error.

\paragraph{Other utilities and Bayesian justification}

Many alternative utility functions to $\mathcal{U}_{\text{SIG}}$ are also used in practice e.g.~posterior precision and and mean squared error between the posterior mean and true parameters.
See \cite{Chaloner:1995} and \cite{Ryan:2016} for detailed reviews.
Like SIG, both the example utilities just given are also based on posterior calculations.
Indeed \cite{Ryan:2016} argue that for a utility to be ``fully Bayesian'', it must be a functional of the posterior,
and other utilities, such as scalar summaries of the Fisher information, are ``pseudo-Bayesian''.
One contribution of this paper is to instead use a decision theoretic justification for which utilities to use in a Bayesian framework, which, surprisingly, provides support to some apparently pseudo-Bayesian utilities
(see discussion in Section \ref{eq:FIGprop}),
and then to further develop this into a game theoretic approach.

\section{Theory: decision theory approach} \label{sec:DTtheory}

This section explores a decision theoretic approach to underlie Bayesian experimental design.
Section \ref{sec:DTframework} presents the framework, which is taken from \cite{Bernardo:1979}.
Section \ref{sec:scoring rules} describes some background on \emph{scoring rules}.
Section \ref{sec:DTres} presents a novel theoretical result showing how the framework supports maximisation of expected utility for various classes of utility function derived from scoring rules.
This includes Shannon information gain, as shown by \cite{Bernardo:1979},
but also the trace of the Fisher information, as proposed by \cite{Walker:2016}.
We conclude in Section \ref{eq:FIGprop} by discussing advantages and disadvantages of the latter utility, and how these reveal limitations of the decision theory framework, motivating our modification in next section.

\subsection{Decision theoretic framework} \label{sec:DTframework}

\cite{Bernardo:1979} proposed the following decision theoretic framework for Bayesian experimental design.
For a summary see the left half of Table \ref{tab:dec}.
The experimenter selects a design, and then nature generates parameters $\theta$ from the prior $\pi(\theta)$ and observations $y$ from the likelihood $f(y|\theta;\tau)$.
The experimenter only observes $y$ and must now choose $a$, a density for $\theta$.
They then receive a reward $\mathcal{R}(a,\theta)$.
In this section and Section \ref{sec:GTtheory} we will assume that:
\begin{enumerate}
\item[A1] The reward is the negative of a \emph{strictly proper scoring rule}, as defined in Section \ref{sec:scoring rules}.
\end{enumerate}
Let $\mathcal{S}(a,\theta) = -\mathcal{R}(a,\theta)$ denote the scoring rule
in the current framework.

\subsection{Scoring rules} \label{sec:scoring rules}

A \emph{scoring rule} $\mathcal{S}(q, \theta)$ measures the quality of a distribution -- in this paper represented by its density $q(\theta)$ -- to model an uncertain quantity, given a realised value $\theta$.
Low scores represent a good match.
A scoring rule is \emph{strictly proper} if, given any $p(\theta)$,
$\E_{\theta \sim p(\theta)}[\mathcal{S}(q, \theta)]$ is uniquely minimised by $q=p$.
For more background on scoring rules see for example \cite {Gneiting:2007} and \cite{Parry:2012}.

Given a scoring rule, two related quantities are
\begin{align*}
\mathcal{H}[p(\theta)] &= \E_{\theta \sim p(\theta)}[\mathcal{S}(p(\theta), \theta)] &&
\text{entropy of $p$} \\
\mathcal{D}[p(\theta), q(\theta)] &= \E_{\theta \sim p(\theta)}[\mathcal{S}(q(\theta), \theta) - \mathcal{S}(p(\theta), \theta)] \!&&
\text{divergence from $p$ to $q$}
\end{align*}
Appendix \ref{app:scoring} gives details of the entropy and divergence for two strictly proper scoring rules which will be used below, logarithmic score and Hyv\"arinen score \citep{Hyvarinen:2005}.
Note that Hyv\"arinen score uses only the derivatives of $\log q(\theta)$ so can be calculated from unnormalised densities.

\subsection{Results} \label{sec:DTres}

Result \ref{res:DTmain} is a general result characterising solutions of the decision theoretic framework.
It uses the following extra assumption:
\begin{enumerate}
\item[A2.]
Both
$\E_{(\theta) \sim \pi(\theta)}[ \, |S(\pi(\theta), \theta)| \, ]$
and, for any $\tau$,
$\E_{(\theta,y) \sim \pi(\theta, y; \tau)}[ \, |S(\pi(\theta|y;\tau), \theta)| \, ]$
are finite.
\end{enumerate}
\begin{result} \label{res:DTmain}
Assume A1 and A2.
Then the following are equivalent, in the sense of sharing the same set of optimal designs:
\begin{enumerate}
\item
The experimenter acts to maximise their expected reward.
\item
The experimenter selects $\tau$ to maximise the expectation, with respect to $\pi(\theta, y; \tau)$, of any of the following
\begin{align*}
\mathcal{U}_{\text{entropy}} &= -\mathcal{H}[\pi(\theta | y; \tau)], \\
\mathcal{U}_{\text{entropy diff}} &= \mathcal{H}[\pi(\theta)] - \mathcal{H}[\pi(\theta | y; \tau)], \\
\mathcal{U}_{\text{divergence}} &= \mathcal{D}[\pi(\theta | y; \tau), \pi(\theta)].
\end{align*}
\end{enumerate}
\end{result}
(Note that arguments of utilities are omitted in this section to simplify notation.)

\begin{proof}
See Appendix \ref{sec:proofDTmain}.
\end{proof}

The next result provides further equivalent utility choices for particular scoring rules.
\begin{result} \label{res:SIG/FIG}
For the logarithmic scoring rule, and assuming A2, maximising the expectation with respect to $\pi(\theta, y; \tau)$
of either $\mathcal{U}_{\text{divergence}}$ or $\mathcal{U}_{\text{SIG}}$, defined in \eqref{eq:SIG},
gives the same set of optimal designs.

\noindent For the Hyv\"arinen scoring rule the same is true for $\mathcal{U}_{\text{divergence}}$ and both of
\[
\mathcal{U}_{\text{FIG}} = ||\nabla_\tau \log f(y|\theta;\tau)||^2, \quad
\mathcal{U}_{\text{trace}} = \trace \mathcal{I}(\theta; \tau),
\]
assuming regularity condition A5\footnote{A5 implies A2, and also contains some other conditions.}, defined in Appendix \ref{app:scoring}.
\end{result}
Here $||\theta||$ represents the $L_2$ norm i.e.~$||x|| = \sqrt{x^T x}$.
\begin{proof}
See Appendix \ref{sec:proofSIGFIG}.
\end{proof}

The first part of Result \ref{res:SIG/FIG} provides a decision theoretic derivation of the Shannon information gain utility.
The argument used is essentially the same as that of \cite{Bernardo:1979}.
In an analogy to this, we refer to $\mathcal{U}_{\text{FIG}}$ in the second part of the result as \emph{Fisher information gain}\footnote{
Some earlier preprints of this paper used this name for $\mathcal{U}_{\text{trace}}$ instead.
Our usage here is a closer analogy to Shannon information gain, and matches that of \cite{Overstall:2020}.}.
We note that \cite{Walker:2016} proved directly that $\mathcal{U}_{\text{trace}}$ and the corresponding $\mathcal{U}_{\text{entropy diff}}$ and $\mathcal{U}_{\text{divergence}}$ all have the same expectations up to an additive constant,
and were therefore equivalent when used in experimental design.
The second part of Result \ref{res:SIG/FIG} shows that the same conclusion arises from a decision theoretic approach based on the Hyv\"arinen score.

\subsection{Fisher information gain properties} \label{eq:FIGprop}

Result \ref{res:SIG/FIG} supports maximising the expectation of $\mathcal{U}_{\text{FIG}}$, or equivalently of $\mathcal{U}_\text{trace}$,
\begin{equation} \label{eq:Jfig}
\mathcal{J}_{\text{FIG}}(\tau) = \E_{\theta \sim \pi(\theta)}[ \trace \mathcal{I}(\theta; \tau) ] = \trace \bar{\mathcal{I}}(\tau).
\end{equation}
where $\bar{\mathcal{I}}(\tau) = \E_{\theta \sim \pi(\theta)}[ \mathcal{I}(\theta; \tau) ]$.
In this section we discuss the properties of this objective.
We refer to this as the FIG approach to experimental design.
We will see that despite computational advantages it has several undesirable properties,
which illustrated later in our examples.
Section \ref{sec:GTtheory} addresses these issues by generalising the decision theoretic framework.

\paragraph{Bayesian justification}

Since the definition of $\mathcal{U}_{\text{trace}}$ does not involve $y$,
it is not a functional of the posterior,
and therefore is pseudo-Bayesian under the terminology of \cite{Ryan:2016}
(discussed at the end of Section \ref{sec:utilities}.)
However, as pointed out by \cite{Walker:2016}, $\mathcal{J}_{\text{FIG}}$ also results from using utilities which are functionals of the posterior -- e.g.~$\mathcal{U}_{\text{divergence}}$ in our notation --
so can be regarded as fully Bayesian.
This shows the divide between pseudo-Bayesian and fully Bayesian approaches can be hard to clearly define, and motivates using a decision or game theoretic approach to do so.

\paragraph{Computational advantages}

The Fisher information is often available in closed form which can easily be evaluated, and allows easy evaluation of gradients.
Then optimisation of $\mathcal{J}_{\text{FIG}}(\tau)$ is straightforward using standard stochastic optimisation methods, as described in Section \ref{sec:opt}.
In particular, the objective does not involve calculating the log evidence, which creates optimisation difficulties for $\mathcal{J}_{\text{SIG}}$, or other use of explicit posterior inference.
The case where evaluation of the Fisher information is more complicated is discussed in Appendix \ref{sec:estimationMore}.

\paragraph{Lack of reparameterisation invariance}

The optimum of $\mathcal{J}_{\text{FIG}}(\tau)$ can change under a reparameterisation of $\theta$.
We illustrate this with a simple example in Section \ref{sec:illustration}.
This property is undesirable as the optimal design is affected by the seemingly irrelevant choice of what parameterisation is used.
In contrast, the design optimising expected Shannon information gain is reparameterisation invariant.

\paragraph{Uninformative designs}

Optimal FIG designs are often highly informative for some linear combinations of parameters, but uninformative for others.
This is illustrated in Figure \ref{fig:pk_posterior} where the design results in a long narrow posterior which has uninformative margins for two parameters,
while a competing design is informative for all parameters.
Also FIG designs can have an excessive amount of \emph{replication}: repeated observations at the same time or location.
For example, in Figure \ref{fig:pk_SGD_traces} shows a design where observations are made at the same time.
These repeated observations provide increasingly accurate inference on one parameter (or parameter combination) at the expense of the others.

An informal theoretical explanation for this behaviour is as follows.
The optimal FIG design maximises the trace of $\bar{\mathcal{I}}(\tau)$
which equals the sum of its eigenvalues.
Often this sum is maximised when one eigenvalue is large and the others are small:
we provide an example in Section \ref{sec:illustration},
and also \citet{Overstall:2020} describes a linear model where this occurs.
In a Bayesian setting, $\mathcal{I}(\theta; \tau)$ is an approximation to posterior precision \cite[see e.g.][]{Vaart:2000}, so this corresponds to the typical posterior having one parameter combination which is accurately inferred at the expense of the others.
\cite{Overstall:2020} also points out that the optimal $\tau$ can produce a singular $\bar{\mathcal{I}}(\tau)$ matrix,
which can cause some statistical methods to break down,
and that a cause is that $\mathcal{J}_{\text{FIG}}$ only uses terms on the diagonal of $\bar{\mathcal{I}}(\tau)$, allowing off-diagonal terms to make it singular.

\section{Theory: game theory approach} \label{sec:GTtheory}

The previous section provided a theoretical framework supporting the use of $\mathcal{U}_\text{trace}$,
but found this produced designs with several undesirable features,
including dependence on the choice of parameterisation and uninformativeness for some linear combinations of parameters.
Here we modify the framework to produce a design which is robust to linear reparameterisation, which, as we shall see, encourages informativeness for all parameter combinations.

Below, Section \ref{sec:GTframework} outlines our framework.
Section \ref{sec:GTdefs} reviews some game theory definitions,
which are used in Section \ref{sec:GTresults}
to characterise optimal designs supported by the framework,
including an adversarial variation on FIG which we refer to as ADV.
Section \ref{sec:properties} describes the properties of ADV and its advantages over FIG.

\subsection{Game theoretic framework} \label{sec:GTframework}

We propose the following game theoretic modification to the framework of Section \ref{sec:DTtheory},
which is summarised in the right half of Table \ref{tab:dec}.
As before, initially the experimenter selects a design $\tau$.
We introduce an \emph{adversary} who now selects the parameterisation.
To so they selects an invertible matrix $A$ defining parameters $\phi = A^{-1} \theta$.
The experimenter must then select $a$, which is now a density for $\phi$,
and their reward is $\mathcal{R}(a,\phi)$.
The adversary's reward is $-\mathcal{R}(a,\phi)$:
they aim to find the parameterisation for which the design does worst.

We continue to use assumption A1,
which now gives a scoring rule of $\mathcal{S}(a,\phi) = -\mathcal{R}(a,\phi)$.
We also introduce another assumption, which is discussed in Appendix \ref{app:det}:
\begin{enumerate}
\item[A3]
The adversary is restricted to selecting $A$ with $\det A = 1$.
\end{enumerate}

\begin{table}[bthp]
\begin{tabular}{p{0mm}p{60mm}|p{0mm}p{70mm}}
& Decision theoretic && Game theoretic \\
\toprule[1pt]
\textbullet & \textbf{E} selects design $\tau$ &
\textbullet & \textbf{E} selects design $\tau$ \\
&& \textbullet & \textcolor{blue}{\textbf{A} selects matrix $A$} \\
\textbullet & \textbf{N} generates parameters $\theta$ from prior &
\textbullet & \textbf{N} generates parameters $\theta$ from prior \\
& (unobserved until reward allocated) && (unobserved until reward allocated) \\
\textbullet & \textbf{N} generates data $y$ from model &
\textbullet & \textbf{N} generates data $y$ from model \\
\textbullet & \textbf{E} selects $a$, estimated density for $\theta$ &
\textbullet & \textcolor{blue}{\textbf{E} selects $a$, estimated density for $\phi = A^{-1} \theta$} \\
\textbullet & \textbf{E} receives reward $\mathcal{R}(a,\theta)$ &
\textbullet & \textcolor{blue}{\textbf{E} receives reward $\mathcal{R}(a,\phi)$} \\
&&
\textbullet & \textcolor{blue}{\textbf{A} receives reward $-\mathcal{R}(a,\phi)$} \\
\bottomrule[1pt]
\end{tabular}
\caption{Two frameworks for experimental design:
the approach of \cite{Bernardo:1979} (left),
and our modification (right), with changes shown in blue.
The abbreviations are for \textbf{E}xperimenter, \textbf{A}dversary and \textbf{N}ature.
Players observe all actions except $\theta$.
}
\label{tab:dec}
\end{table}

\subsection{Game theory definitions} \label{sec:GTdefs}

We refer to $\tau,A,\theta,y,a$ as \emph{actions} of our game.
The game specifies a mapping from $\tau,A,\theta,y,a$ to real-valued rewards for the \emph{players}: the experimenter and adversary.
The actions $\theta, y$ are random samples from specified distributions.
The remaining actions are chosen by the players according to \emph{decision rules},
$\tau, A(\tau), a(\tau, A, y)$.

In the remainder of this section we explore the set of \emph{subgame perfect equilibria} (SPEs): players act to optimise their expected reward under the assumption that later decisions will also do so.
In our setting, SPE is a condition on decision rules.
The first condition is that $a(\tau,A,y)$ must output some $a$ maximising the experimenter's expected reward given $(\tau,A,y)$.
The second condition is that $A(\tau)$ must output $A$ maximising the adversary's expected reward given $\tau$ and $a(\tau,A,y)$.
Finally $\tau$ must maximise the experimenter's expected reward given $A(\tau)$ and $a(\tau,A,y)$.
Note that there may be multiple or no decision rule combinations meeting these conditions.

SPEs are often argued \citep[see e.g.][]{Jin:2019}
to be an appropriate solution concept for games in which players act sequentially.
The alternative solution concept of \emph{Nash equilibria} is more appropriate for games with simultaneous actions.
For more background on game theory and solution concepts see e.g.~\cite{Osborne:1994}.

\subsection{Results} \label{sec:GTresults}

Our first result is that under logarithmic score, the game theoretic framework is essentially equivalent to the decision theoretic framework.

\begin{result} \label{res:GTlogscore}
The following sets are equal:
\begin{itemize}
\item Designs which are selected in SPEs of the game theoretic framework with logarithmic score and assumptions A1,A3.
\item Optimal designs in the decision theoretic framework with assumption A1 and logarithmic score.
\end{itemize}
\end{result}

\begin{proof}
See Appendix \ref{sec:proofGTlogscore}.
\end{proof}

The next result shows that for Hyv\"arinen score, the game theoretic framework addresses the issue of sensitivity to reparameterisation which occurred for the decision theoretic framework.
The result focuses on the case of a \emph{linear reparameterisation}
which we define as reparameterising $\theta$ to $B \theta$ for some invertible $B$.

\begin{result} \label{res:GTreparam}
The set of designs which are selected in SPEs of the game theoretic framework with Hyv\"arinen score and assumptions A1,A3 is unchanged by a linear reparameterisation.
\end{result}

\begin{proof}
See Appendix \ref{sec:proofGTreparam}.
\end{proof}

Our main result characterises which designs are selected in SPEs under Hyv\"arinen score.
First we introduce the following technical assumption:
\begin{enumerate}
\item[A4]
There exists some $\tau$ such that $\E_{\theta \sim \pi(\theta)}[\det \mathcal{I}(\theta; \tau)] > 0$.
\end{enumerate}
Informally this means there is a design under which the model is guaranteed to provide some information on every parameter, or linear combination of parameters.
We also use a regularity condition A5 (see Appendix \ref{app:scoring}).

\begin{result} \label{res:GTmain}
Consider the game theoretic framework with Hyv\"arinen score under assumptions A1, A3-A5.
Then the set of action pairs which are selected in SPEs
are those which solve $\min_\tau \max_A \mathcal{K}(\tau,A)$
where
\begin{equation} \label{eq:K}
\mathcal{K}(\tau,A)
= -\E_{\theta \sim \pi(\theta)} \trace [ A^T \mathcal{I}(\theta; \tau) A ].
\end{equation}
Also, the set of designs which are selected in SPEs are those maximising
\begin{equation} \label{eq:Jadv}
\mathcal{J}_{\text{ADV}}(\tau) = \det \bar{\mathcal{I}}(\tau),
\end{equation}
where $\bar{\mathcal{I}}(\tau) = \E_{\theta \sim \pi(\theta)}[ \mathcal{I}(\theta; \tau) ]$.
\end{result}

\begin{proof}
See Appendix \ref{sec:proofGTmain}.
\end{proof}

\subsection{Adversarial objective properties} \label{sec:properties}

Using the decision theoretic framework with Hyv\"arinen score produced the objective
$\mathcal{J}_{\text{FIG}}(\tau) = \trace \bar{\mathcal{I}}(\tau)$.
We have shown that the game theoretic framework instead produces the objective
$\mathcal{J}_{\text{ADV}}(\tau) = \det \bar{\mathcal{I}}(\tau)$.
This has improved properties, as we discuss now.

\paragraph{Linear reparameterisation invariance}

By Result \ref{res:GTreparam}, the set of optimal designs from $\mathcal{J}_{\text{ADV}}(\tau)$ is invariant to linear reparameterisation, unlike $\mathcal{J}_{\text{FIG}}(\tau)$.
This is also easy to show directly.
Consider a linear reparameterisation $\phi = B \theta$ with $\det B \neq 0$.
From \eqref{eq:reparam}, the Fisher information  is $\mathcal{I}_\phi(\phi; \tau) = B^{-T} \mathcal{I}_\theta(\theta; \tau) B^{-1}$.
Thus
$\bar{\mathcal{I}}_\phi(\tau) = B^{-T} \bar{\mathcal{I}}_\theta(\tau) B^{-1}$
and the reparameterised objective $\mathcal{J}_{\phi, \text{ADV}}(\tau)$ equals $\mathcal{J}_{\theta, \text{ADV}}(\tau) (\det B)^{-2}$:
the original objective multiplied by a positive constant.
Hence the set of optimal designs is unchanged.

\paragraph{Equivalent $A$ matrices}

From the cyclic property of trace,
\[
\mathcal{K}(\tau,A)
= -\E_{\theta \sim \pi(\theta)} \trace [ A A^T \mathcal{I}(\theta; \tau) ].
\]
Hence any two $A$ matrices producing the same $A A^T$ are equivalent
in that they give the same $\mathcal{K}(\cdot,A)$ function.
In particular, given any $A$, a Cholesky factor of $A A^T$ is equivalent.
So when we perform minimax optimisation, we can restrict $A$ to be a Cholesky matrix
(i.e.~lower triangular with positive diagonal entries)
with determinant 1.
This will help set up our optimisation algorithm later -- see Section \ref{sec:Arep}.

\paragraph{Informative designs}

As discussed in Section \ref{eq:FIGprop}, $\mathcal{J}_{\text{FIG}}(\tau)$ sometimes produces designs in which one parameter combination was inferred well but others were not.
The objective $\mathcal{J}_{\text{ADV}}(\tau)$ penalises such designs.
This is because it is the product of the eigenvalues of $\bar{\mathcal{I}}(\tau)$.
Thus, compared to $\mathcal{J}_{\text{FIG}}(\tau)$, it is much less advantageous to make one eigenvalue large and the others small.
See Figure \ref{fig:pk_posterior} for an illustration of informativeness
-- the posterior under the ADV design is concentrated with respect to both parameters shown in this plot, unlike that under the FIG design.

A related property of the objective $\mathcal{J}_{\text{FIG}}(\tau)$ is that it sometimes produces singular or near-singular $\bar{\mathcal{I}}(\tau)$ matrices.
The objective $\mathcal{J}_{\text{ADV}}(\tau)$ avoids this by directly penalising matrices with low determinants.
One reason this is possible is that the determinant operation is affected by off-diagonal elements of $\bar{\mathcal{I}}(\tau)$, unlike trace.

\paragraph{Bayesian justification}

Unlike $\mathcal{J}_{\text{FIG}}(\tau)$, $\mathcal{J}_{\text{ADV}}(\tau)$ is not defined as the expected value of a utility function.
Therefore the ``pseudo-Bayesian'' or ``fully Bayesian'' definitions of \cite{Ryan:2016} discussed in Section \ref{sec:utilities} cannot be directly applied.
However, $\mathcal{J}_{\text{ADV}}(\tau)$ meets the spirit of the pseudo-Bayesian definition as it depends on Fisher information matrices rather than posteriors.
Nonetheless, we have shown that it emerges from a game theoretic framework based on enabling an experimenter to estimate the full posterior well.

\paragraph{Link to $D$-optimality}

Classical optimal design using the $D$-optimality criterion requires a design to maximise $\det \mathcal{I}(\theta; \tau)$ at a reference $\theta$ value.
This can be generalised to a Bayesian setting incorporating parameter uncertainty in several ways: see \cite{Atkinson:2007} for a list of objectives, which includes $\mathcal{J}_{\text{ADV}}(\tau)$.
One contribution of our work is to provide theoretical support and an optimisation method for this particular choice.

\paragraph{Computational tractability}

It is hard to optimise $\mathcal{J}_{\text{ADV}}(\tau)$ directly for continuous $\tau$, as it is not easy to obtain an unbiased gradient estimate due to the non-linearity of the determinant operator.
However optimisation based on \eqref{eq:K} is tractable, as described in the next section.

\section{Optimisation} \label{sec:opt}

Theorem \ref{res:GTmain} of the previous section
motivates performing optimal design by finding minimax solutions of
$\mathcal{K}(\tau,A) = -\E_{\theta \sim \pi(\theta)} \trace [ A^T \mathcal{I}(\theta; \tau) A ]$.
Since the function $\mathcal{K}$ is an expectation,
this allows straightforward computation of unbiased gradient estimates
when $\mathcal{I}(\theta; \tau)$ can be evaluated.
(See Appendix \ref{sec:estimationMore} for discussion of the case where this is not possible.)
This section describes how as a consequence it is possible to find candidate minimax solutions using generic gradient based optimisation methods.
Throughout we assume $\nabla_\tau \mathcal{I}(\theta; \tau)$ and $\nabla_\tau \E_{\theta \sim \pi(\theta)} \mathcal{I}(\theta; \tau)$ exist.

Section \ref{sec:GDA} gives background on gradient based optimisation methods for minimax problems.
Section \ref{sec:Arep} discusses how to deal with the constraint $\det A = 1$.
Section \ref{sec:Kgrad} describes calculation of unbiased gradient estimates.
Section \ref{sec:main alg} presents our optimisation algorithm,
and describes various implementation details.

\subsection{Gradient descent ascent} \label{sec:GDA}

Algorithm \ref{alg:GDA}, \emph{gradient descent ascent}, attempts to solve $\min_x \max_y f(x,y)$.
It iteratively updates $(x_t, y_t)$ based on unbiased gradient estimates
$g_{x,t} \approx -\nabla_x f(x_t,y_t), g_{y,t} \approx \nabla_y f(x_t,y_t)$.
For an overview and history of GDA see \cite{Lin:2019}.
GDA generalises \emph{stochastic gradient descent} (SGD),
which is the special case where $f$ is a function of $x$ alone,
and only $x_t$ updates are needed.

\begin{algorithm}[htbp] \caption{Gradient descent ascent (GDA)} \label{alg:GDA}
\begin{algorithmic}[1]
\STATE Input: Initial values $x_1, y_1$, update subroutines $h_x, h_y$
\FOR{$t = 1,2,\ldots$}
\STATE Compute $g_{x,t}, g_{y,t}$, unbiased estimates of $-\nabla_x f(x_t,y_t)$ and $\nabla_y f(x_t,y_t)$.
\STATE Update estimates using $x_{t+1} = x_t + h_x(g_{x,t})$, $y_{t+1} = y_t + h_y(g_{y,t})$.
\ENDFOR
\end{algorithmic}
\end{algorithm}

A simple update rule for GDA,
which we will refer to as the \emph{default} update rule,
is $z_{t+1} = z_t + \alpha_{z,t} g_{z,t}$ (for $z \in \{ x,y \}$),
given predefined learning rate sequences $\alpha_{x,t}, \alpha_{y,t}$.
Stochastic approximation theory \citep[see e.g.][]{Kushner:2003} suggests using learning rates such that $\sum_{t=1}^\infty \alpha_{z,t} = \infty$,
$\sum_{t=1}^\infty \alpha_{z,t}^2 < \infty$
(for $z \in \{x,y\}$).
This ensures that SGD using this update rule and unbiased gradient estimates converges to a local minimum, under appropriate regularity conditions.

More sophisticated update rules have been developed,
effectively tuning the learning rates adaptively and using different learning rates for each component of the $x$ and $y$ vectors.
We use the popular Adam update rule \citep{Kingma:2015}.
(We use the default Adam tuning parameters, which typically produce good empirical performance but not asymptotic convergence.
Convergence guarantees are possible by setting some parameters to decay: see \citealp{Kingma:2015})

Convergence of GDA under any learning rule is more complicated than SGD,
and is an area of active research.
One issue is that the dynamics can produce limit cycles as well as limit points.
Algorithms to avoid this have been suggested, including a \emph{two time-scale update rule} \citep{Borkar:1997, Heusel:2017} in which $\lim_{t \to \infty} \alpha_{y,t} / \alpha_{x,t} = \infty$
(or a similar condition under the Adam learning rule.)
Another issue is to characterise the limit points of GDA as an appropriate local generalisation of minimax solutions \citep{Heusel:2017, Jin:2019, Lin:2019}.

However, for our experimental design application, we find empirically that standard GDA methods with the Adam update rule suffice to produce sensible designs.
Therefore we recommend using GDA with this update rule,
and checking for convergence using diagnostics, multiple runs from different initial values, and post-processing (all detailed later in this section).

\subsection{Representation of $A$} \label{sec:Arep}

We wish to solve $\min_\tau \max_A \mathcal{K}(\tau,A)$
under the constraint that $\det A = 1$.
As discussed in Section \ref{sec:properties}, it is sufficient to search for $A$ over Cholesky matrices
(i.e.~lower triangular with positive diagonal entries)
with determinant 1.
Such matrices can be represented as
\begin{equation} \label{eq:A(eta)}
A(\eta) =
\begin{pmatrix}
\exp(\eta_{11}) & 0           & \ldots & 0 \\
\eta_{21}       & \exp(\eta_{22}) & \ldots & 0 \\
\vdots      & \vdots      & \ddots & \vdots \\
\eta_{p1}       & \eta_{p2}       & \ldots & \exp(-\sum_{i=1}^{p-1} \eta_{ii})
\end{pmatrix}.
\end{equation}
This maps an unconstrained real vector $\eta$ of $\eta_{ij}$ variables to the space of matrices of interest.
We can now solve $\min_\tau \max_\eta \mathcal{K}(\tau,A(\eta))$ using GDA.
We initialise $\eta = 0$ so that initially $A=I$, corresponding to the adversary making no reparameterisation.

\subsection{Gradient estimation} \label{sec:Kgrad}

In this section we will consider estimating $\nabla_z \mathcal{K}(\tau, A(\eta))$ for $z \in \{\tau, \eta\}$
in the common case where it is easy to evaluate the Fisher information and related gradients.
See Appendix \ref{sec:estimationMore} for discussion of cases where these are intractable.

In a few cases it is possible to directly evaluate
$\bar{\mathcal{I}}(\tau) = \E_{\theta \sim \pi(\theta)}[ \mathcal{I}(\theta; \tau) ]$.
See Sections \ref{sec:illustration} and \ref{sec:geostat} for example.
Then, using \eqref{eq:K},
we have $\mathcal{K}(\tau, A) = - \trace[ A^T \bar{\mathcal{I}}(\tau) A ]$,
and gradients can be evaluated using automatic differentiation \citep{Baydin:2017} or derivation of direct expressions for them.
Our code performs automatic differentiation using the PyTorch library \citep{Paszke:2019}.

More commonly it is necessary to derive unbiased gradient estimates.
From the definition of $\mathcal{K}(\tau, A)$, \eqref{eq:K},
an unbiased estimate is
\begin{equation} \label{eq:Kest}
\hat{\mathcal{K}}(\tau, A) = -\trace \left[A^T \left\{ \frac{1}{K} \sum_{k=1}^K \mathcal{I}(\theta^{(k)}; \tau) \right\} A \right].
\end{equation}
We can also prove the following result.
\begin{result} \label{res:unbiased}
Under appropriate regularity conditions on $\trace [A^T \mathcal{I}(\theta; \tau) A]$
(see Appendix \ref{sec:regularity}),
then $\nabla_z \hat{\mathcal{K}}(\tau, A(\eta))$ is an unbiased estimate of
$\nabla_z \mathcal{K}(\tau, A(\eta))$ for $z \in \{\tau, \eta\}$.
\end{result}

\begin{proof}
See Appendix \ref{sec:proofUnbiased}.
\end{proof}
We calculate gradients of $\hat{\mathcal{K}}(\tau, A)$ using automatic differentiation in PyTorch.
Note that using a larger $K$ in \eqref{eq:Kest}
reduces the variance but increases computational cost.
We explore this trade-off in Section \ref{sec:pk}.

\subsection{Main algorithm} \label{sec:main alg}

Algorithm \ref{alg:GDA_BED} applies GDA to our experimental design setting.
Note that we perform $R$ replications of GDA in parallel from different initial conditions.
This is typically more efficient than repeating the algorithm $R$ times in serial.
One reason is that the calculations in step 3 and 4 are amenable to parallelisation.
(We ran our experiments on a CPU, but our PyTorch code can easily exploit GPU parallelisation, allowing for further speed improvements.)
Another reason is that the same simulations in step 3 can be reused for all replications.
Algorithm \ref{alg:GDA_BED} can also be used for SGD to optimise the FIG objective by keeping $\eta=0$ and only updating $\tau$.

\begin{algorithm}[htbp] \caption{Gradient descent ascent for Bayesian experimental design} \label{alg:GDA_BED}
\begin{algorithmic}[1]
\STATE Input:
Number of samples to use in \eqref{eq:Kest} $K$,
number of parallel replications to perform $R$,
initial values $\tau_1^i, \eta^i_1$ (for $i = 1,2,\ldots,R$),
update subroutines $h_\tau, h_\eta$
\FOR{$t = 1,2,\ldots$}
\STATE Sample $\theta^{(k)}$ from the prior for $k = 1, 2, \ldots, K$.
\STATE Compute $g^i_{\tau,t}, g^i_{\eta,t}$, unbiased estimates of $-\nabla_\tau \mathcal{K}(\tau_t^i, A(\eta_t^i))$ and $\nabla_\eta \mathcal{K}(\tau_t^i, A(\eta_t^i))$ using automatic differentiation\footnotemark of \eqref{eq:Kest} for all $i$.
\STATE Update estimates using $\tau^i_{t+1} = \tau^i_t + h_\tau(g^i_{\tau,t})$, $\eta^i_{t+1} = \eta^i_t + h_\eta(g^i_{\eta,t})$ for all $i$.
\ENDFOR
\end{algorithmic}
\end{algorithm}
\footnotetext{
In practice we implement this step by differentiating
$\sum_{i=1}^R \hat{\mathcal{K}}(\tau_t^i, A(\eta_t^i))$
with respect to all $\tau, \eta$ variables.
This is more efficient in backwards mode automatic differentiation, as implemented in PyTorch.}

In the remainder of this section various implementation details are discussed.

\subsubsection{Diagnostics}

Step 4 of Algorithm \ref{alg:GDA_BED} involves calculating $\hat{\mathcal{K}}(\tau,A)$ from \eqref{eq:Kest},
whose gradients are then found using automatic differentiation.
This Monte Carlo estimate of $\mathcal{K}$ can be used as a diagnostic of the algorithm's performance.
As we are performing minimax optimisation, it will typically rise and fall before reaching an equilibrium.
Also the values of $\hat{\mathcal{K}}$ from different parallel runs can be highly correlated since they are based on the same $\theta^{(k)}$ samples.
Both phenomena can be seen in the bottom left graph of Figure \ref{fig:pk_GDA_traces}.
(The presence of correlation is illustrated by the right hand side of the graph seems to show a single thick line.
In fact these are multiple lines with different values which are highly correlated with each other.)

An alternative diagnostic is to estimate $\mathcal{J}_{\text{ADV}}$ as defined in \eqref{eq:Jadv}.
A Monte Carlo estimate is
\[
\hat{\mathcal{J}}(\tau) = \det \left( \frac{1}{J} \sum_{j=1}^J \mathcal{I}(\tilde{\theta}^{(j)}; \tau) \right).
\]
Our algorithm can optionally calculate $\hat{\mathcal{J}}$ for designs produced during its execution.
To do so we initially sample $\tilde{\theta}^{(j)}$ values from the prior for $j = 1,2,\ldots,J$ (we take $J=1000$).
Then $\hat{\mathcal{J}}$ is calculated for each design $\tau_t^i$ produced during the algorithm.
Unlike $\hat{\mathcal{K}}$, this adds an extra computational cost.
However this diagnostic has the advantage that it directly estimates the objective $\mathcal{J}_{\text{ADV}}$ so larger values roughly correspond to better performing designs.
For example, the bottom right graph of Figure \ref{fig:pk_GDA_traces} directly traces the performance of designs during optimisation.
we
In the case of SGD optimisation we can return $\hat{\mathcal{K}}(\tau,I)$ as a diagnostic.
This directly estimates $-\mathcal{J}_{\text{FIG}}(\tau)$ but does so using a different $\theta$ sample each time, which adds some variability to the diagnostic.
To remedy this we can also calculate
\[
\hat{\mathcal{J}}_{\text{FIG}}(\tau) = \trace \left( \frac{1}{J} \sum_{j=1}^J \mathcal{I}(\tilde{\theta}^{(j)}; \tau) \right)
\]
for a fixed sample of $\tilde{\theta}^{(j)}$ values, as above.

\subsubsection{Termination}

We run Algorithm \ref{alg:GDA_BED} for a fixed number of iterations or fixed runtime.
Alternatively it could be run until a convergence condition is met for one of the diagnostics above,
or for the size of the increments to $\tau$ or $A$.

\subsubsection{Optimisation under constraints}

We often wish to find the optimal design under a constraint: $\tau \in \mathcal{T} \subset \mathbb{R}^d$.
In the applications in this paper this was simple to achieve by simple pragmatic approaches.
In Section \ref{sec:illustration} we represented $\tau$ as the transformation of an unbounded vector.
In Section \ref{sec:pk} the designs remained in $\mathcal{T}$ under unconstrained optimisation, so no modification was needed.
In Section \ref{sec:geostat} we added a large penalty to rewards for $\tau \not \in \mathcal{T}$, whose gradient moves designs back towards $\mathcal{T}$.

In more complex settings these methods may not suffice.
A more sophisticated alternative would be to compose each GDA update with a projection operation into $\mathcal{T}$ \citep{Kushner:2003}.

\subsubsection{Local optima}

GDA can often converge to multiple possible locally optimal designs.
To attempt to find the global optimum we run GDA multiple times from different initial values of $\tau$ (keeping $\eta$ initialised as a zero vector),
and compare their $\hat{\mathcal{J}}$ diagnostics.

In some settings we can also use a \emph{point exchange algorithm}.
Suppose we must select multiple \emph{design points} from some region.
Optimal designs often have a high degree of \emph{replication}: they consist of a small number of clusters of repeated observations \citep{Gotwalt:2009, Overstall:2017, Binois:2019}.
While gradient based optimisation can find the optimal cluster locations well,
there may be a large number of local optima, differing by the number of points in each cluster.
See Figure \ref{fig:pk_GDA_designs} for an example.
A point exchange algorithm takes cluster locations as input and uses discrete optimisation to find optimal cluster sizes.
We use a simple approach detailed in Appendix \ref{sec:pointExchange}.
There is scope for developing more sophisticated approaches in future work.

\section{Poisson model illustration} \label{sec:illustration}

This section provides a simple illustrative example of the properties of the SIG, FIG and ADV approaches to experimental design.
The setting is that an experimenter must divide a unit of time between two experiments making Poisson observations with different rates.
More precisely, the design is $\tau \in [0,1]$.
There are 2 independent observations $y_1 \sim Poisson(\tau \theta_1 \omega_1),
y_2 \sim Poisson([1-\tau] \theta_2 \omega_2)$.
We assume $\omega_1 > \omega_2 > 0$ and that $\theta_1, \theta_2$ have independent $Gamma(2,1)$ priors.

\paragraph{Analytic results}

Optimal designs for this example can be derived analytically.
We summarise the results here: see Appendix \ref{app:Poisson} for derivations.
The optimal design is $\tau=1$ under FIG and $\tau=1/2$ under SIG or ADV.
Under the FIG design, $y_2$ is always zero, so this design is uninformative for $\theta_2$
(i.e.~the posterior always equals the prior).
However the SIG/ADV design is informative for both parameters.
Also, the SIG and ADV designs are invariant to linear reparameterisations.
However for FIG, a linear reparameterisation can change the optimal design to $\tau_1=0$ (or make all values of $\tau$ optimal).

\paragraph{Numerical optimisation}

We perform a numerical analysis of this example with $\omega_1=2, \omega_2=1$.
Since the design $\tau \in [0,1]$ is bounded
we optimise a transformation, $\lambda = \text{logit}(\tau)$,
to allow the use of unconstrained optimisation.
Following \eqref{eq:A(eta)} we take adversary's action to be
\[
A(\eta) =
\begin{pmatrix}
\exp(\eta_{11}) & 0           \\
\eta_{21}       & \exp(-\eta_{11})
\end{pmatrix}.
\]

Figure \ref{fig:illustration} shows a GDA vector field using default update rules with a particular choice of learning rates.
The figure shows $\lambda$ and $\eta_{11}$ for $\eta_{21}$ fixed at zero.
(In practice other $\eta_{21}$ values quickly converge to zero.)
The vector field illustrates spiral trajectories converging to a fixed point.
Additionally, it shows that for any fixed value of $\eta_{11}$, $\frac{\partial \mathcal{K}}{\partial \lambda}$ has a fixed sign.
This illustrates that for fixed $A$ (i.e.~the FIG setting),
SGD optimisation produces a design converging to either $\tau=0$ or $1$.

The figure also shows GDA trajectories as the adversary's learning rate is varied.
In all cases, the trajectories converge on the fixed point.
However, convergence is much faster as the adversary's learning rate is increased.

\begin{figure}[htp] \label{fig:illustration}
\begin{center}
\includegraphics[width=0.7\textwidth]{poisson_vector_field.pdf}
\end{center}
\caption{GDA vector field and trajectories for the Poisson example.
Paths are shown for experimenter learning rate $10^{-2}$
and adversary learning rate $10^{-3}$ (blue), $10^{-4}$ (green), $10^{-5}$ (red).
Points indicate every 500 steps of optimisation.
The vector field corresponds to the $10^{-3}$ adversary learning rate.}
\end{figure}

\section{Pharmacokinetic example} \label{sec:pk}

This section contains simulations studies on a pharmacokinetic model.
The goal is to investigate the performance of our FIG and ADV approaches and compare them with existing methods for SIG.

\subsection{Model}

Pharmacokinetics studies the time course of drug absorption,
distribution, metabolism, and excretion.
Concentration level has to be observed from via samples of fluid
(such as blood, plasma or urine) from the subject taken at preplanned time points.
The number of observations are constrained by budget and resources,
as well as patient comfort and wellbeing.
We assume that observed concentration, $y_j$, at time $\tau_i$ (in hours) is distributed as
\[
y_j \sim N(x(\theta, \tau_i), \sigma^2),
\qquad \text{where}\
x(\theta, \tau_i) = \frac{D\theta_2 (\exp[-\theta_1 \tau_i] - \exp[-\theta_2 \tau_i]),}{\theta_3(\theta_2 - \theta_1)},
\]
and $D=400$.
Concentrations at different times are assumed to be independent.
We follow \cite{Ryan:2014} and \cite{Overstall:2017} who investigate Bayesian experimental design for a similar model.
However we remove a multiplicative noise term for simplicity.
Following this earlier work we assume independent log normal prior distributions
\[
\theta_1 \sim LN(\log 0.1, 0.05), \theta_2 \sim LN(\log 1, 0.05), \theta_3 \sim LN(\log 20, 0.05),
\]
and aim to find $15$ observation times in $[0,24]$. Also we treat $\sigma = 0.1$ as known.
Most previous work requires consecutive observations to be at least $15$ minutes apart.
We do not enforce this condition as its implementation would vary between methods making it more difficult to draw fair conclusions.

\subsection{Methods}

The Fisher information matrix for this model is available in closed form.
See Appendix \ref{app:pkFIG} for details.
This allows the optimal design under ADV to be found using Algorithm \ref{alg:GDA_BED},
and the optimal design under FIG using a SGD variant (only updating $\tau$).
In both cases we used the Adam update rule.
When implementing these algorithms, we found no need to use constrained optimisation as the designs remained in the interval $[0,24]$ in any case.
Both algorithms performed 100 replications in parallel, based on 100 initial designs sampled from a uniform distribution over the search space.
We considered several choices of $K$, the number samples to estimate $\mathcal{K}$ in \eqref{eq:Kest}, and selected $K=1$.
See Appendix \ref{app:Kcomparison} for details.

We compare our results to two methods of finding the optimal SIG design:
the approximate coordinate exchange (ACE) algorithm of \cite{Overstall:2017},
and the prior contrastive estimation (PCE) algorithm of \cite{Foster:2020}.
(PCE is one of several methods in \cite{Foster:2020}, and is not their overall recommendation. However we found PCE converged more quickly than the alternatives for this example, so we use it as lower bound on the speed of their methods.)
Implementation details of these methods are given in Appendices \ref{app:pkACE}--\ref{app:pkPCE}.
PCE allowed 100 replications to be run in parallel.
We ran replications of ACE in serial, noting that the ACE code utilises multiple cores during its execution in any case.
As ACE took much longer to run we used only 30 replications.

\subsection{Results}

Table \ref{tab:timings} shows mean times for each method
i.e.~run time divided by number of designs returned,
demonstrating the speed advantage of GDA and SGD.
Below we comment in more detail on each algorithm.

\begin{table}[hbtp]
\begin{tabular}{c|ccccccc}
                      & GDA & GDA+PE & SGD & SGD+PE & ACE  & PCE \\
\hline
Mean time (seconds)   & 1.4 & 2.2    & 1.4 & 1.5    & 8012 & 36  \\
Number of repetitions & 100 & 100    & 100 & 100    & 30   & 100
\end{tabular}
\caption{Mean times to run optimisation methods on the pharmacokinetic example.
(n.b.~PE is ``point exchange'').
Both ACE and PCE have not fully converged after the time quoted, so the figures given are lower bounds on the run time.
}
\label{tab:timings}
\end{table}

\paragraph{ADV}

The top half of Figure \ref{fig:pk_GDA_traces} shows trajectories of $\tau$ and $A$ during GDA optimisation for a single replication.
The design points eventually settle into clusters of repeated observations at three times.
The bottom half shows estimated $\mathcal{K}$ and $\mathcal{J}$ objectives over 100 replications.
Although all runs have converged after $100,000$ iterations, the objective values are not identical.

Figure \ref{fig:pk_GDA_designs} shows final designs for all replications.
The design points typically converge to 3 clusters,
around times $1.1$, $3.4$ and $14.0$.
However the cluster sizes vary between runs, as different local optima are found.
This explains why runs converge to slightly different objective values.
Employing point exchange reduces the variation in cluster sizes.
(There is some variability in exact cluster locations.
Further investigation showed that this is mainly due to cluster locations changing slightly depending on cluster sizes.)
PE finds two candidates for an optimal design, with 6 or 7 points near time 1.
(To find an overall optimum one could now estimate $\hat{J}$ for both using a larger sample size.)
However GDA runs produced a maximum of 6 points for this cluster.
This highlights the importance of PE in escaping local optima.

Finally, Figure \ref{fig:pk_posterior} illustrates the posterior from a ADV design, and shows that it is concentrated around the true parameters,
unlike that for a FIG design.

\begin{figure}[htp]
\begin{center}
\includegraphics[width=0.45\textwidth]{design_GDA_trace.pdf}
\includegraphics[width=0.45\textwidth]{A_GDA_trace.pdf} \\
\includegraphics[width=0.45\textwidth]{K_GDA_trace.pdf}
\includegraphics[width=0.45\textwidth]{J_GDA_trace.pdf}
\end{center}
\caption{Trace plots for GDA on the pharmacokinetic example.
The top row shows output for a single run of GDA.
The bottom row shows output for 100 runs.}
\label{fig:pk_GDA_traces}
\end{figure}

\begin{figure}[htp]
\begin{center}
\includegraphics[width=0.45\textwidth]{GDA_designs.pdf}
\includegraphics[width=0.45\textwidth]{GDA_and_gda_pe_designs.pdf} \\
\end{center}
\caption{Designs output from 100 optimisation runs for the pharmacokinetic example.
The horizontal axis shows the index of each point in the sorted design
i.e.~observation times are shown in increasing order from left to right.
Points are plotted as semi-transparent, so rare results are light, and common results are dark.
The left plot shows GDA output after 100,000 iterations.
In the right plot point exchange is also applied.
}
\label{fig:pk_GDA_designs}
\end{figure}

\begin{figure}[htp]
\begin{center}
\includegraphics[width=0.7\textwidth]{PK_posterior.pdf}
\end{center}
\caption{Bivariate posteriors for $\theta_1, \theta_3$ in the pharmacokinetic example.
Samples are plotted from: the prior (green circles),
the posterior with FIG design (red triangles),
the posterior with ADV design (blue pluses).
The ADV and FIG designs used 100,000 iterations of GDA or SGD,
followed by point exchange optimisation.
The true parameter values, sampled from the prior, are $\theta_1=0.14,  \theta_3=12.29$, consistent with the blue pluses.
Datasets for use in inference were then sampled from the model,
using the same observation noise for both ADV and FIG.
Posterior samples are importance sampling output.}
\label{fig:pk_posterior}
\end{figure}

\paragraph{FIG}

Full details of SGD optimisation results are given in Appendix \ref{app:pkSGD}.
Briefly, SGD designs always converge to a single cluster around time 12.
Intuitively this is a poor design:
the three $\theta$ parameters of the $x(\theta,\tau)$ function cannot all be identified from repeated observations at a single time point $\tau$.
Figure \ref{fig:pk_posterior} illustrates that this design indeed performs poorly by presenting a typical posterior.
The narrowness of the FIG design posterior show that it is informative about some function of $\theta_1$ and $\theta_3$.
However the posterior is mostly uninformative about $\theta_1$ and $\theta_3$ marginally:
it stretches almost the full length of the prior distribution,
unlike the ADV design posterior which is highly concentrated.

\paragraph{SIG}

Running ACE under the default tuning choices took a few hours, while our approach took only a few minutes to run 100 repetitions.
Furthermore, the ACE results do not seem to have converged to the optimal design in the time.
See Appendix \ref{app:pkACE} for more details.
We ran PCE for roughly 10 times longer than our approach and found the results concentrated near 3 observation times similar to those for ADV.
However this runtime appears to only be a lower bound on the time needed for convergence, as running PCE for longer increased the concentration of the design points.
We also explored alternative methods from \cite{Foster:2020} but were not able to improve on the PCE results.
See Appendix \ref{app:pkPCE} for more details.

\subsection{Conclusion}

We have shown that our optimisation method to find ADV and FIG designs is faster than SIG optimisation methods by at least a factor of 10.
The true advantage may be greater as the SIG methods do not appear to have fully converged in the time stated.

The FIG design can give uninformative marginal posteriors for some parameters.
The ADV design avoids this drawback, and appears to be similar to the SIG design.
Multiple runs of ADV produced designs with similar cluster locations but varying cluster sizes.
Post-processing these using point exchange found two candidates for optimal cluster sizes.
This included one which was different from any cluster sizes found without post-processing, illustrating the importance of this step.

\section{Geostatistical regression example} \label{sec:geostat}

This section considers an example requiring hundreds of design choices,
to illustrate how our method can scale up to higher dimensional applications.

\subsection{Model}

Consider the following geostatistical regression model.
Here a design $\tau$ is a $d \times 2$ matrix whose rows specify measurement locations.
The model assumes normal observations with a linear trend and squared exponential covariance function with a nugget effect, giving
\begin{align*}
y &\sim N(x(\theta, \tau), \Sigma(\tau)), \quad
&x_i &= \theta_1 \tau_{i1} + \theta_2 \tau_{i2}, \\
\Sigma &= \sigma_1^2 I + \sigma_2^2 R(\tau), \quad
&R_{ij} &= \exp \left[-\sum_{k=1}^2 (\tau_{ik} - \tau_{jk})^2 / \ell^2 \right].
\end{align*}
For simplicity we assume that $\sigma^2_1, \sigma^2_2$ (observation variance components) and $\ell$ (covariance length scale) are known.
Hence the unknown parameters are $\theta_1$ and $\theta_2$ (trends).

\subsection{Methods}

Using \eqref{eq:mvn_info} the Fisher information matrix is
\[
\mathcal{I}(\tau) = \tau^T \Sigma(\tau)^{-1} \tau.
\]
Note that this does not depend on $\theta$.
Hence $\mathcal{I}(\tau)$ is also the expected Fisher information matrix
and we do not need to use Monte Carlo to estimate it.

We performed simulation studies
for $\ell = 0.01,0.02,0.04,0.08$ with $\sigma_1=1$ and $\sigma_2=3$
to search for $500$ design points restricted to a unit square centred at the origin.
The design was initialised as independent uniform draws.
Each run used Algorithm \ref{alg:GDA_BED} with the Adam update rule for 1000 iterations,
which was enough for convergence (see Appendix \ref{app:geostats}).
We implemented constrained optimisation by adding a $L_1$ penalty to designs outside the unit square.

As this example aims to illustrate the time required by ADV, we did not investigate repeated runs from different initial designs, or post-processing using point exchange.
In any case, it seems unlikely the latter would help as there is little evidence of replicated observations in the results.

\subsection{Results}

Optimisation took on average 19.4 seconds.
Figure \ref{fig:geodesigns} shows the resulting designs.
For small $\ell$ values, the design points cluster in the corners.
For larger values, the designs are spread across the region
with varying spatial structures.
For all runs $A$ remained very close to $I$ throughout optimisation,
reflecting the symmetry of $\theta_1$ and $\theta_2$ in the model.
Hence FIG would also work well for this example.

\begin{figure}[htbp]
  \centering
  \setlength\tabcolsep{0pt}
  \begin{tabular}{cccc}
    $\ell=0.01$ & $\ell = 0.02$ & $\ell = 0.04$ & $\ell = 0.08$ \\
    \begin{minipage}{0.25\textwidth} \includegraphics[width =\textwidth]{geo_design_0.pdf} \end{minipage} &
    \begin{minipage}{0.25\textwidth} \includegraphics[width =\textwidth]{geo_design_1.pdf} \end{minipage} &
    \begin{minipage}{0.25\textwidth} \includegraphics[width =\textwidth]{geo_design_2.pdf} \end{minipage} &
    \begin{minipage}{0.25\textwidth} \includegraphics[width =\textwidth]{geo_design_3.pdf} \end{minipage}
  \end{tabular}
  \caption{Geostatistical regression model designs returned for various choices of $\ell$.
The design space is a unit square centred at the origin.}
  \label{fig:geodesigns}
\end{figure}

\section{Discussion} \label{sec:discussion}

We have presented a gradient descent ascent algorithm Bayesian optimal design using an objective function based on the Fisher information.
This provides improvements in speed and scalability to higher dimensional designs by avoiding the need for posterior inference.
We also provide a novel game theoretic justification for our objective,
and provide theoretical insights into the choice of utilities for Bayesian optimal design from decision/game theoretic principles,
extending the work of \cite{Bernardo:1979} and \cite{Walker:2016}.

In our simulation studies our approach finds locally optimal design faster than other state of the art methods, by a factor of at least 10.
To assess whether multiple locally optimal designs exist we recommend performing multiple runs of our algorithm from random initial designs, which can easily be done in a parallel version of our algorithm.
In situations where locally optimal designs involve clusters of repeated observations,
we recommend performing post-processing using a point exchange algorithm.
Although we did not observe any, GDA dynamics can converge to limit cycles rather than the desired solutions.
To check this we recommend checking trace plots of $\tau$ values produced during optimisation for cyclic behaviour.

\subsection{Limitations and future research}

\paragraph{Intractable Fisher information}

Our methods have assumed that the Fisher information, and associated gradients, can easily be evaluated.
When the Fisher information cannot be evaluated in this way but the score function \eqref{eq:score} can, it is possible to produce an unbiased estimate of the Fisher information.
Appendix \ref{sec:estimationMore} describes this and how it can be used for experimental design.

However for many models the score function cannot be evaluated due to the presence of latent variables, such as nuisance parameters.
Unbiased estimation of the Fisher information is possible if one can sample from the posterior for the latent variables.
See Appendix \ref{sec:estimationMore} for details.
An interesting direction for future research is to approximate this posterior during the optimisation process using machine learning approaches, similar to \cite{Foster:2020} and \cite{Kleinegesse:2020}.
This would help unlock applications to experimental design for model choice and to optimise predictive performance \citep{Ryan:2016}, both cases where model parameters are treated as nuisance parameters.

\paragraph{Discrete designs}

Gradient based optimisation is only available for a continuous space of possible designs.
It would be of interest to develop analogous methods to this paper for discrete designs.
These could be based on discrete optimisation algorithms,
or involve relaxation of the discrete problem to a continuous approximation.

\paragraph{Discrete parameters}

A limitation of our work is that it does not apply to discrete parameters,
since Hyv\"arinen score and Fisher information are only defined for continuous parameters.
Discrete analogues have been proposed \citep{Dawid:2012, Shao:2019},
which would be interesting to investigate.

\paragraph{Alternative optimisers}

Variations on GDA with better convergence guarantees,
such as two time scaled update rules \citep{Heusel:2017},
are an active topic of research and could be used for the objective in this paper.
Another possibility for future work is modifying generic gradient based optimisation methods to avoid local optima in optimal design problems
e.g.~using tempering methods, or non-local updates such as line search, as used by \cite{Overstall:2017}.

\paragraph{Variations to game theoretic framework}

Several details of our game theoretic framework could be altered.
One possibility is to consider alternative scoring rules.
Many alternative proper scoring rules exist beyond logarithmic and Hyv\"arinen \citep{Parry:2012}.
Alternatively, non-proper scoring rules could be used which emphasise a particularly important aspect of inference to the task at hand e.g.~the tails of the distribution \citep{Loaiza:2019}.

Another variation is to allow more freedom to the adversary.
This could include making more general, non-linear, reparameterisations,
or being able to condition their reparameterisation on the observations $y$.

\section*{Acknowledgements}

The authors thank Kevin Wilson, Chris Oates, John Matthews, Phil Dawid, Kyle Cranmer, Iain Murray and anonymous referees for helpful comments, and gratefully acknowledge an EPSRC studentship supporting Sophie Harbisher.

\bibliographystyle{apalike}

redpajama_set_name: RedPajamaArXiv