Bayesian experimental design without posterior calculations: an adversarial approach

TL;DR

This paper introduces an adversarial method for Bayesian experimental design that avoids costly posterior calculations by optimizing Fisher information, addressing invariance issues and improving design efficiency.

Contribution

It proposes a novel adversarial approach that enhances Bayesian experimental design by eliminating the need for posterior computations and ensuring invariance to reparameterisation.

Findings

The adversarial approach avoids posterior calculations.

It maintains invariance to reparameterisation.

It enables efficient gradient-based optimization.

Abstract

Most computational approaches to Bayesian experimental design require making posterior calculations repeatedly for a large number of potential designs and/or simulated datasets. This can be expensive and prohibit scaling up these methods to models with many parameters, or designs with many unknowns to select. We introduce an efficient alternative approach without posterior calculations, based on optimising the expected trace of the Fisher information, as discussed by Walker (2016). We illustrate drawbacks of this approach, including lack of invariance to reparameterisation and encouraging designs in which one parameter combination is inferred accurately but not any others. We show these can be avoided by using an adversarial approach: the experimenter must select their design while a critic attempts to select the least favourable parameterisation. We present theoretical properties of this approach and show it can be used with gradient based optimisation methods to find designs efficiently in practice.

Tables

Table 1. Table 1: 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.

	GDA	GDA+PE	SGD	SGD+PE	ACE	PCE
Mean time (seconds)	1.4	2.2	1.4	1.5	8012	36
Number of repetitions	100	100	100	100	30	100

Table 2. Table 2: Summary of two scoring rules and related quantities. Here p ​ ( θ ) 𝑝 𝜃 p(\theta) is the true density of an unobserved quantity and q ​ ( θ ) 𝑞 𝜃 q(\theta) is a generic density. Note that Δ Δ \Delta is the Laplacian operator (sum of second partial derivatives) with respect to θ 𝜃 \theta . For a derivation of the Hyvärinen divergence see Appendix A of Hyvärinen, ( 2005 ) . The other derivations are straightforward.

	Logarithmic score	Hyvärinen score
Scoring rule	$-\log q\left(\theta \right)$	$2\Delta \log q\left(\theta \right)+{\Vert \nabla \log q\left(\theta \right)\Vert }^2$
Entropy	$-{𝔼}_{\theta \sim p\left(\theta \right)}\left[\log p\left(\theta \right)\right]$	$-{𝔼}_{\theta \sim p\left(\theta \right)}\left[{\Vert \nabla \log p\left(\theta \right)\Vert }^2\right]$
Divergence	${𝔼}_{\theta \sim p\left(\theta \right)}\left[\log p\left(\theta \right)-\log q\left(\theta \right)\right]$	${𝔼}_{\theta \sim p\left(\theta \right)}\left[{\Vert \nabla \log p\left(\theta \right)-\nabla \log q\left(\theta \right)\Vert }^2\right]$

Table 3. Table 3: Two frameworks for experimental design: the approach of Bernardo, ( 1979 ) (left), and our modification (right), with changes shown in blue. The abbreviations are for E xperimenter, C ritic and N ature. Players observe all actions except θ 𝜃 \theta .

	Decision theoretic		Game theoretic
•	E selects design $\tau$	•	E selects design $\tau$
		•	C selects matrix $A$
•	N generates parameters $\theta$ from prior	•	N generates parameters $\theta$ from prior
	(unobserved until reward allocated)		(unobserved until reward allocated)
•	N generates data $y$ from model	•	N generates data $y$ from model
•	E selects $a$, estimated density for $\theta$	•	E selects $a$, estimated density for $\varphi =A^{-1}\theta$
•	E receives reward $ℛ(a,\theta )$	•	E receives reward $ℛ(a,\varphi )$
		•	C receives reward $-ℛ(a,\varphi )$

Equations

$$\pi (\theta |y;\tau )$$

$$\pi (y;\tau )$$

$$\pi (\theta ,y;\tau )=\pi (\theta )f(y|\theta ;\tau )=\pi (\theta |y;\tau )\pi (y;\tau ).$$

$${\mathcal{I}}_{\theta }(\theta ;\tau )=E_{y\sim f(y|\theta ;\tau )}[u(y,\theta ;\tau )u(y,\theta ;\tau {)}^T],$$

$$u(y,\theta ;\tau )={\nabla }_{\theta }\log f(y|\theta ;\tau ).$$

$${\mathcal{I}}_{\theta }(\theta )=J(\varphi {)}^T{\mathcal{I}}_{\varphi }(\varphi )J(\varphi ).$$

$${\mathcal{I}}_{\theta }(\theta )=J(x{)}^T{\mathcal{I}}_x(x)J(x)=J(x{)}^T{\Sigma }^{-1}J(x).$$

$$\mathcal{J}(\tau )=E_{(\theta ,y)\sim \pi (\theta ,y;\tau )}[\mathcal{U}(\tau ,\theta ,y)].$$

$${\mathcal{U}}_{\text{SIG}}(\tau ,\theta ,y)=\log \pi (\theta |y;\tau )-\log \pi (\theta ),$$

$${\mathcal{U}}_{\text{SIG}}(\tau ,\theta ,y)=\log f(y|\theta ;\tau )-\log \pi (y;\tau ).$$

$$\hat{\pi }(y;\tau )=\frac{1}{L}\sum _{\ell =1}^{L}f(y|{\theta }^{(\ell )};\tau ),$$

$$\mathcal{H}[p(\theta )]$$

$$\mathcal{D}[p(\theta ),q(\theta )]$$

$${\mathcal{U}}_{\text{entropy}}$$

$${\mathcal{U}}_{\text{entropy diff}}$$

$${\mathcal{U}}_{\text{divergence}}$$

$${\mathcal{U}}_{\text{FIG}}=||{\nabla }_{\tau }\log f(y|\theta ;\tau )|{|}^2,{\mathcal{U}}_{\text{trace}}=\mathrm{tr}\mathcal{I}(\theta ;\tau ),$$

$${\mathcal{J}}_{\text{FIG}}(\tau )=E_{\theta \sim \pi (\theta )}[\mathrm{tr}\mathcal{I}(\theta ;\tau )]=\mathrm{tr}\bar{\mathcal{I}}(\tau ).$$

$$\mathcal{K}(\tau ,A)=-E_{\theta \sim \pi (\theta )}\mathrm{tr}[A^T\mathcal{I}(\theta ;\tau )A].$$

$${\mathcal{J}}_{\text{ADV}}(\tau )=\det \bar{\mathcal{I}}(\tau ),$$

$$\mathcal{K}(\tau ,A)=-E_{\theta \sim \pi (\theta )}\mathrm{tr}[A{A}^T\mathcal{I}(\theta ;\tau )].$$

$$A(\eta )=\left(\begin{array}{cccc}\exp ({\eta }_{11})& 0& \dots & 0\\ {\eta }_{21}& \exp ({\eta }_{22})& \dots & 0\\ ⋮& ⋮& \ddots & ⋮\\ {\eta }_{p1}& {\eta }_{p2}& \dots & \exp (-{\sum }_{i=1}^{p-1}{\eta }_{ii})\end{array}\right).$$

$$\hat{\mathcal{K}}(\tau ,A)=-\mathrm{tr}\left[A^T\left\{\frac{1}{K}\sum _{k=1}^K\mathcal{I}({\theta }^{(k)};\tau )\right\}A\right].$$

$${\hat{\mathcal{J}}}_{\text{ADV}}(\tau )=\det \left(\frac{1}{J}\sum _{j=1}^J\mathcal{I}({\tilde{\theta }}^{(j)};\tau )\right).$$

$${\hat{\mathcal{J}}}_{\text{FIG}}(\tau )=\mathrm{tr}\left(\frac{1}{J}\sum _{j=1}^J\mathcal{I}({\tilde{\theta }}^{(j)};\tau )\right)$$

$$A(\eta )=\left(\begin{array}{cc}\exp ({\eta }_{11})& 0\\ {\eta }_{21}& \exp (-{\eta }_{11})\end{array}\right).$$

$$y_i\sim N(x(\theta ,{\tau }_i),{\sigma }^2),\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)},$$

$${\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),$$

$$y_i\sim N(x(\theta ,{\tau }_i),{\sigma }_1^2+{\sigma }_2^{2}x(\theta ,{\tau }_i{)}^2),$$

$$y$$

$$\Sigma$$

$$-E_{\theta \sim \pi (\theta |y;\tau )}[\mathcal{S}(a,\theta )].$$

$${\mathcal{U}}_{\text{divergence}}={\mathcal{U}}_{\text{entropy}}+E_{\theta \sim \pi (\theta |y;\tau )}[S(\pi (\theta ),\theta )].$$

Full text

1 Introduction

Selecting a good design for an experiment can be crucial to extracting useful information and controlling costs. Applications include medical interventions (Amzal et al.,, 2006), epidemic modelling (Cook et al.,, 2008), pharmacokinetics (Ryan et al.,, 2014; Overstall and Woods,, 2017) and ecology (Gillespie and Boys,, 2019). In modern applications it is increasingly feasible to take a large number of measurements – e.g. placing sensors (Krause et al.,, 2009) or making observations in a numerical integration problem (Oates et al.,, 2020) – or make other complex designs – e.g. selecting a time series of chemical input levels to a synthetic biology experiment (Bandiera et al.,, 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 as 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 (Müller,, 1999), and using variational inference to learn approximate surrogate posteriors (Foster et al.,, 2019, 2020), but these still effectively solve inference problems and hence remain costly. (We comment on other related approaches in Section 1.2.)

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. 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 Bernardo, (1979). Our derivation is based on judging the quality of a parameter estimate through the Hyvärinen score (Hyvärinen,, 2005), rather than through the logarithmic score as in Bernardo, (1979). We demonstrate that it is straightforward to optimise the resulting expected utility using recent developments in stochastic gradient optimisation (Kingma and Ba,, 2015) and automatic differentiation (Baydin et al.,, 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 Overstall and Woods, (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.

We 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 a critic 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 critic also encourages designs not to neglect the posterior accuracy of any parameter combination: if one did then the critic could choose a parameterisation concentrating on this weakness. We show it is possible to find optimal designs in this game theoretic framework using generic gradient descent ascent methods, which have been much studied recently in the machine learning literature (e.g. Heusel et al.,, 2017; Jin et al.,, 2020).

Below, we summarise our contributions in Section 1.1 and related work in Section 1.2. In the remainder of the paper, Section 2 presents background on Bayesian experimental design. Section 3 presents results on the decision theoretic framework of Bernardo, (1979), and Section 4, on our proposed game theoretic extension. Proofs and other technical material are presented in the appendices. Section 5 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 6). We provide a detailed simulation study on a pharmacokinetic model (Section 7), 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 8). Code for these examples is available at https://github.com/dennisprangle/AdversarialDesignCode. All examples were run on a desktop PC with 12 CPU cores. Our conclusion, Section 9, summarises our findings and recommendations to implement our methods. It also discusses limitations of our work, and future research directions to address these, including a discussion of the intractable Fisher information case (detailed further in the appendices).

1.1 Contributions

Our main contribution is a faster approach to Bayesian experimental design using the Fisher information. Where the Fisher information is available in closed form, our approach outperforms existing methods (Section 7) and scales easily to designs with hundreds of design choices (Section 8).

In addition, we contribute by extending the decision theoretic framework for Bayesian experimental design of Bernardo, (1979) to allow other proper scoring rules in addition to logarithmic score (Section 3.1). We also show that using the Hyvärinen score in the decision theoretic framework results in a utility based on the trace of the Fisher information, referred to as Fisher information gain or FIG in the paper (Section 3.3). This provides a decision-theoretic justification to a measure that was previously suggested by Walker, (2016) based on information theoretic arguments. We explore the limitations of FIG (Section 3.4) and address them by extending the decision theoretic approach to a game theoretic approach, and show the relevant solution concept gives more intuitively informative designs (Section 4). We also show that the game theoretic solution can easily be obtained in practice using gradient-based minimax optimisation.

Finally, we provide new insight into two long-standing questions in Bayesian experimental design. Firstly: what counts as a Bayesian utility function? (See Sections 3.4 and 4.4 under the heading “Bayesian justification”.) Secondly: which of several possible Bayesian generalisations of $D$-optimality should be preferred? (See Section 4.4 under the heading “Link to $D$-optimality”.)

1.2 Related work

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

Gradient-based optimisation methods for experimental design have been explored previously. Pronzato and Walter, (1985) optimise the expected determinant of the Fisher information using analytically derived gradients. Huan and Marzouk, (2013, 2014) optimise expected Shannon information using gradients (either derived analytically or based on finite differences) for a biased numerical approximation to the utility. Foster et al., (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.

Foster et al., (2020) and Kleinegesse and Gutmann, (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 2.4). 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.

Overstall and Woods, (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 only the current component 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 Müller, (1999) performs optimal design using Markov chain Monte Carlo. Amzal et al.,, 2006 and Kück et al.,, 2006 extend this approach to use sequential Monte Carlo. In other approaches, Ryan et al., (2014) look at high dimensional designs with a low dimensional parametric form, Price et al., (2018) use evolutionary algorithms, and Gillespie and Boys, (2019) search a discrete grid of designs.

2 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 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 design points, $\tau_{1},\tau_{2},\ldots,\tau_{d}$.

This section reviews relevant details of the statistical background. First Sections 2.1 and 2.2 give some necessary definitions on Bayesian statistics and Fisher information. Section 2.3 describes Bayesian experimental design, which is based on maximising expected utility. Finally, Section 2.4 summarises some common utility functions.

2.1 Bayesian framework

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 evidence) for $y$. In our setting both depend on the experimental design $\tau$,

[TABLE]

The prior and model define a joint density111 For discrete $y$ a density with respect to a product of Lebesgue and counting measures can be used.,

[TABLE]

2.2 Fisher information

We will make frequent use of the Fisher information matrix (FIM) for $\theta$,

[TABLE]

which is based on the score function,

[TABLE]

The subscript in $\mathcal{I}_{\theta}$ denotes which variable is used for differentiation in (5). 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.

Examples

The FIM for exponential family models is based on the variance of the sufficient statistics (see e.g. Lehmann and Casella,, 2006) which is often available in closed form. Two examples we use in this paper are:

• •

The Poisson distribution, $y\sim Poisson(\phi)$. Here $\mathcal{I}(\phi)=1/\phi$.

• •

The multivariate normal distribution with known variance, $y\sim N(\mu,\Sigma)$. Here $\mathcal{I}(\mu)=\Sigma^{-1}$.

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 (4) that

[TABLE]

An application of this result is to the model $y\sim N(x(\theta,\tau),\Sigma)$. Using (6),

[TABLE]

2.3 Bayesian experimental design

The Bayesian approach to experimental design involves selecting a function $\mathcal{U}=\mathcal{U}(\tau,\theta,y)$, giving the utility of choosing 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 expected utility of $\tau$ i.e. the prior predictive utility

[TABLE]

See Chaloner and Verdinelli, (1995), Atkinson et al., (2007) and Ryan et al., (2016) for comprehensive surveys of Bayesian experimental design.

2.4 Utility functions

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

Shannon information gain (SIG)

This is a popular and well-motivated utility choice,

[TABLE]

introduced by 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, Bernardo, (1979) gave a decision theoretic derivation of $\mathcal{U}_{\text{SIG}}$ which we recap in Section 3. A helpful property of $\mathcal{U}_{\text{SIG}}$ is that it is 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 (9) is complicated by the need to evaluate the posterior density. However, using (1) the utility can be rewritten as

[TABLE]

Optimisation now requires evaluation of a posterior summary: the log evidence $\log\pi(y;\tau)$. A common SIG estimate replaces $\pi(y;\tau)$ in (10) with a Monte Carlo approximation

[TABLE]

where $\theta^{(\ell)}$ are independent prior samples. A typical choice of $L$ is 1000 (Overstall and Woods,, 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 1.2. However the need for log evidence estimation remains a source of computational expense and approximation error.

Other utilities and Bayesian justification

Many alternative utility functions to $\mathcal{U}_{\text{SIG}}$ are also used in practice e.g. posterior precision and mean squared error between the posterior mean and true parameters. See Chaloner and Verdinelli, (1995) and Ryan et al., (2016) for detailed reviews. Like SIG, both utilities just mentioned are also based on posterior calculations. Indeed Ryan et al., (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 FIM, are “pseudo-Bayesian”. (Throughout we use this definition of “pseudo-Bayesian”, but note some authors use the term differently e.g. Overstall,, 2020.) 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 3.4), and then to further develop this into a game theoretic approach.

3 Theory: decision theory approach

This section explores a decision theoretic approach to underlie Bayesian experimental design. Section 3.1 presents the framework, which is taken from Bernardo, (1979). Section 3.2 describes some background on scoring rules. Section 3.3 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 Bernardo, (1979), but also the trace of the FIM, as proposed by Walker, (2016). We conclude in Section 3.4 by discussing advantages and disadvantages of the latter utility, and how these reveal limitations of the decision theory framework, motivating our modification in Section 4.

3.1 Decision theoretic framework

Bernardo, (1979) proposed the following decision theoretic framework for Bayesian experimental design. 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$, receiving reward $\mathcal{R}(a,\theta)$. In this section and Section 4 we will assume that:

1. A1

The reward is the negative of a strictly proper scoring rule, as defined in Section 3.2.

Throughout this section $\mathcal{S}(a,\theta)=-\mathcal{R}(a,\theta)$ denotes the scoring rule.

3.2 Scoring rules

A 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 strictly proper if, given any $p(\theta)$, $\operatorname{\mathbb{E}}_{\theta\sim p(\theta)}[\mathcal{S}(q,\theta)]$ is uniquely minimised by $q=p$. For more background on scoring rules see for example Gneiting and Raftery, (2007) and Parry et al., (2012).

Given a scoring rule, two related quantities are

[TABLE]

Appendix A gives details of the entropy and divergence for two strictly proper scoring rules which will be used below: logarithmic score and Hyvärinen score (Hyvärinen,, 2005). Hyvärinen score uses only the derivatives of $\log q(\theta)$, so it can be calculated from unnormalised densities.

3.3 Results

Result 1 is a general result characterising solutions of the decision theoretic framework. It uses the following extra assumption:

1. A2

Both $\operatorname{\mathbb{E}}_{\theta\sim\pi(\theta)}[\,|S(\pi(\theta),\theta)|\,]$ and, for any $\tau$, $\operatorname{\mathbb{E}}_{(\theta,y)\sim\pi(\theta,y;\tau)}[\,|S(\pi(\theta|y;\tau),\theta)|\,]$ are finite.

Result 1**.**

Assume A1 and A2. Then the following are equivalent, in the sense of sharing the same set of optimal designs:

The experimenter acts to maximise their expected reward. 2. 2.

The experimenter selects $\tau$ to maximise the expectation, with respect to $\pi(\theta,y;\tau)$, of any of the following utilities:

[TABLE]

(Note that arguments of utilities are omitted in this section to simplify notation.)

The next result provides further equivalent utility choices for particular scoring rules.

Result 2**.**

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 (9), gives the same set of optimal designs.

For the Hyvärinen scoring rule the same is true for $\mathcal{U}_{\text{divergence}}$ and both of

[TABLE]

assuming regularity condition A5222A5 implies A2, and also contains some other conditions., defined in Appendix A.

Here $||\theta||$ represents the $L_{2}$ norm i.e. $||x||=\sqrt{x^{T}x}$.

The first part of Result 2 provides a decision theoretic derivation of the Shannon information gain utility. The argument used is essentially the same as that of Bernardo, (1979). In an analogy to this, we refer to $\mathcal{U}_{\text{FIG}}$ in the second part of the result as Fisher information gain333 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 Overstall, (2020).. We note that 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 2 shows that the same conclusion arises from a decision theoretic approach based on the Hyvärinen score.

3.4 Fisher information gain properties

Result 2 supports maximising the expectation of $\mathcal{U}_{\text{FIG}}$, or equivalently of $\mathcal{U}_{\text{trace}}$,

[TABLE]

where $\bar{\mathcal{I}}(\tau)=\operatorname{\mathbb{E}}_{\theta\sim\pi(\theta)}[\mathcal{I}(\theta;\tau)]$. We refer to this as the FIG approach to experimental design, and here we discuss its properties. We will see that despite computational advantages it has several undesirable properties, illustrated later in our examples. Section 4 addresses these issues by generalising the decision theoretic framework.

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 Ryan et al., (2016) (discussed at the end of Section 2.4.) However, as pointed out by 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 it 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 (/game) theoretic approach to do so.

Computational advantages

The FIM is often available in a 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 5. In particular, the objective does not involve calculating the log evidence, which creates optimisation difficulties for $\mathcal{J}_{\text{SIG}}$, or any other use of explicit posterior inference. The case where evaluation of the FIM is more complicated is discussed in Appendix G.

Lack of reparameterisation invariance

The optimum of $\mathcal{J}_{\text{FIG}}(\tau)$ can change under a reparameterisation of $\theta$: see Section 6 for a simple 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.

Intuitively uninformative designs

By construction, optimal FIG designs maximise a particular definition of informativeness. However these designs often have properties which seem “uninformative” in an intuitive sense.

One issue is that FIG designs often result in posterior distributions which are very diffuse for some linear combinations of parameters, although they may be concentrated for others. This is illustrated in Figure 5, where the design results in a long narrow posterior in which the marginals for two parameters are very similar to the prior marginals. In contrast, a competing design produces posterior marginals which are concentrated around the true values for both parameters.

Also, FIG designs can have an excessive amount of replication: repeated observations at the same time or location. For example, Appendix I.5 shows a FIG design with all observations at the same time. Such repeated observations provide increasingly accurate inference on one parameter (or parameter combination) at the expense of the others.

An informal 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 6, and also Overstall, (2020) describes a linear model where this occurs. In a Bayesian setting, $\mathcal{I}(\theta;\tau)$ is an approximation to posterior precision (see e.g. Van der Vaart,, 2000), so this corresponds to the typical posterior having one parameter combination which is accurately inferred at the expense of the others. Overstall, (2020) also points out that the optimal $\tau$ can produce a singular $\bar{\mathcal{I}}(\tau)$ matrix, causing some statistical methods to break down, and a reason for this 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.

4 Theory: game theory approach

The previous section provided a theoretical framework supporting the use of $\mathcal{U}_{\text{trace}}$, but found it produced designs with several undesirable features, including dependence on the choice of parameterisation and a diffuse marginal posterior 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 marginal posterior concentration for all parameter combinations.

Below, Section 4.1 outlines our framework. Section 4.2 reviews some game theory definitions, which are used in Section 4.3 to characterise optimal designs supported by the framework, including an adversarial variation on FIG which we refer to as ADV. Section 4.4 describes the properties of ADV and its advantages over FIG.

4.1 Game theoretic framework

We propose the following game theoretic framework. Initially, the experimenter selects a design $\tau$. We introduce a critic who now selects the parameterisation by choosing an invertible matrix $A$ defining parameters $\phi=A^{-1}\theta$. The experimenter must then select $a$, a density for $\phi$, receiving a reward $\mathcal{R}(a,\phi)$. The critic’s reward is $-\mathcal{R}(a,\phi)$: they aim to find the parameterisation for which the design does worst. This completes the specification of a game, as discussed further in Section 4.2.

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 D:

1. A3

The critic is restricted to selecting $A$ with $\det A=1$.

Table 3 in Appendix B summarises and contrasts both the decision-theoretic and game-theoretic frameworks.

4.2 Game theory definitions

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

In the remainder of this section we explore the set of 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 critic’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 existence and uniqueness of decision rules meeting these conditions is not automatic.

SPEs are often argued (see e.g. Jin et al.,, 2020) to be an appropriate solution concept for games in which players act sequentially. The alternative solution concept of Nash equilibria is more appropriate for games with simultaneous actions. For more background on game theory and solution concepts see e.g. Osborne and Rubinstein, (1994).

4.3 Results

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

Result 3**.**

The following sets are equal:

• •

Designs which are selected in SPEs of the game theoretic framework with logarithmic score and assumptions A1,A3.

• •

Optimal designs in the decision theoretic framework with assumption A1 and logarithmic score.

The next result shows that for Hyvärinen 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 linear reparameterisation which we define as reparameterising $\theta$ to $B\theta$ for some invertible $B$.

Result 4**.**

The set of designs which are selected in SPEs of the game theoretic framework with Hyvärinen score and assumptions A1,A3 is unchanged by a linear reparameterisation.

Finally, our main result characterises which designs are selected in SPEs under Hyvärinen score. It uses the following technical assumption:

1. A4

There exists some $\tau$ such that $\operatorname{\mathbb{E}}_{\theta\sim\pi(\theta)}[\det\mathcal{I}(\theta;\tau)]>0$.

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 A).

Result 5**.**

Consider the game theoretic framework with Hyvärinen score under assumptions A1, A3-A5. Then the set of $\tau,A$ pairs which are selected in SPEs are those which solve $\min_{\tau}\max_{A}\mathcal{K}(\tau,A)$ where

[TABLE]

Also, the set of designs which are selected in SPEs are those maximising

[TABLE]

where $\bar{\mathcal{I}}(\tau)=\operatorname{\mathbb{E}}_{\theta\sim\pi(\theta)}[\mathcal{I}(\theta;\tau)]$.

4.4 Adversarial objective properties

Using the decision theoretic framework with Hyvärinen score produced the objective $\mathcal{J}_{\text{FIG}}(\tau)=\operatorname{tr}\bar{\mathcal{I}}(\tau)$. Result 5 shows 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 in this subsection.

Result 5 also shows that maximisation of $\mathcal{J}_{\text{ADV}}(\tau)$ is equivalent to minimax optimisation of $\mathcal{K}(\tau,A)$, as defined in (13). This is helpful in discussing the properties of $\mathcal{J}_{\text{ADV}}(\tau)$, and particularly useful in defining a practical optimisation scheme, as detailed below and in Section 5.

Linear reparameterisation invariance

By Result 4, 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 (6), the FIM 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.

Equivalent $A$ matrices

From the cyclic property of trace,

[TABLE]

Hence any two $A$ matrices producing the same $AA^{T}$ are equivalent in that they give the same $\mathcal{K}(\cdot,A)$ function. In particular, given any $A$, a Cholesky factor of $AA^{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 5.2.

Intuitively informative designs

As discussed in Section 3.4, $\mathcal{J}_{\text{FIG}}(\tau)$ sometimes produces designs in which one parameter combination is inferred accurately but others are 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 5 for an illustration – 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 the trace.

Bayesian justification

Unlike $\mathcal{J}_{\text{FIG}}(\tau)$, $\mathcal{J}_{\text{ADV}}(\tau)$ is not defined as the expectation of a utility function. Therefore the “pseudo-Bayesian” or “fully Bayesian” definitions of Ryan et al., (2016) discussed in Section 2.4 cannot be directly applied. However, $\mathcal{J}_{\text{ADV}}(\tau)$ meets the spirit of the pseudo-Bayesian definition as it depends on FIMs rather than posteriors. Nonetheless, we have shown it emerges from a game theoretic framework based on enabling an experimenter to estimate the posterior well.

Link to $D$-optimality

Classical optimal design using $D$-optimality 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 Atkinson et al., (2007), Table 18.1, for a list of 5 possible objectives, including $\mathcal{J}_{\text{ADV}}(\tau)$ (and the objective of Pronzato and Walter, (1985), mentioned in Section 1.2, which differs from $\mathcal{J}_{\text{ADV}}(\tau)$ by swapping the order of determinant and expectation.) One contribution of our work is to provide theoretical support and an optimisation method for this particular choice.

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 (13) is tractable, as described in the next section.

5 Optimisation

Result 5 of the previous section motivates performing optimal design by finding minimax solutions of $\mathcal{K}(\tau,A)=-\operatorname{\mathbb{E}}_{\theta\sim\pi(\theta)}\operatorname{tr}[A^{T}\mathcal{I}(\theta;\tau)A]$. Since this is an expectation, computation of unbiased gradient estimates is straightforward when $\mathcal{I}(\theta;\tau)$ can be evaluated. (See Appendix G 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}\operatorname{\mathbb{E}}_{\theta\sim\pi(\theta)}\mathcal{I}(\theta;\tau)$ exist.

Section 5.1 gives background on gradient based optimisation methods for minimax problems. Section 5.2 discusses how to deal with the constraint $\det A=1$. Section 5.3 describes calculation of unbiased gradient estimates. Section 5.4 presents our optimisation algorithm, and describes various implementation details.

5.1 Gradient descent ascent

Algorithm 1, gradient descent ascent, attempts to solve $\min_{x}\max_{y}f(x,y)$. It iteratively updates $(x_{t},y_{t})$ based on $g_{x,t}$ and $g_{y,t}$, unbiased gradient estimates of $-\nabla_{x}f(x_{t},y_{t})$ and $\nabla_{y}f(x_{t},y_{t})$. For an overview and history of GDA see Lin et al., (2020). GDA generalises stochastic gradient descent (SGD), which is the special case where $f$ is a function of $x$ alone, and only $x_{t}$ updates are needed.

A simple update rule for GDA, which we will refer to as the 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 (see e.g. Kushner and Yin,, 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 (Kingma and Ba,, 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 Kingma and Ba,, 2015.)

Convergence of GDA under any update 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 two time-scale update rule (Borkar,, 1997; Heusel et al.,, 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 (Heusel et al.,, 2017; Jin et al.,, 2020; Lin et al.,, 2020).

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. Details of these are contained in the following subsections.

5.2 Representation of $A$

We wish to solve $\min_{\tau}\max_{A}\mathcal{K}(\tau,A)$ under the constraint that $\det A=1$. As discussed in Section 4.4, 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

[TABLE]

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 critic making no reparameterisation.

5.3 Gradient estimation

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 FIM and related gradients. See Appendix G for discussion of cases where these are intractable.

In a few cases it is possible to directly evaluate $\bar{\mathcal{I}}(\tau)=\operatorname{\mathbb{E}}_{\theta\sim\pi(\theta)}[\mathcal{I}(\theta;\tau)]$. See Sections 6 and 8 for example. Then, using (13), we have $\mathcal{K}(\tau,A)=-\operatorname{tr}[A^{T}\bar{\mathcal{I}}(\tau)A]$, and gradients can be evaluated using automatic differentiation (Baydin et al.,, 2017) or derivation of direct expressions for them. Our code performs automatic differentiation using the PyTorch library (Paszke et al.,, 2019).

More commonly it is necessary to derive unbiased gradient estimates. From the definition of $\mathcal{K}(\tau,A)$, (13), an unbiased estimate is

[TABLE]

We also have the following result.

Result 6**.**

Under appropriate regularity conditions on $\operatorname{tr}[A^{T}\mathcal{I}(\theta;\tau)A]$ (see Appendix F), 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\}$.

We calculate gradients of $\hat{\mathcal{K}}(\tau,A)$ using automatic differentiation in PyTorch. Note that using a larger $K$ in (16) reduces the variance but increases computational cost. We explore this trade-off in Section 7.

5.4 Main algorithm

Algorithm 2 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 2 can also be used for SGD to optimise the FIG objective by keeping $\eta=0$ and only updating $\tau$.

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

5.4.1 Diagnostics

Step 4 of Algorithm 2 involves calculating $\hat{\mathcal{K}}(\tau,A)$ from (16), 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 2. (The presence of correlation is illustrated by the fact that the right hand side of the graph seems to show a single thick line. In fact there 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 (14) by

[TABLE]

Our code can optionally calculate $\hat{\mathcal{J}}_{\text{ADV}}$ 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}}_{\text{ADV}}$ 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 correspond to better performing designs. For example, the bottom right graph of Figure 2 directly traces the performance of designs during optimisation.

For 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

[TABLE]

for a fixed sample of $\tilde{\theta}^{(j)}$ values, as above.

5.4.2 Termination

We run Algorithm 2 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$.

5.4.3 Optimisation under constraints

We often wish to find the optimal design under a constraint: $\tau\in\mathcal{T}\subset\mathbb{R}^{d}$. In our applications we achieved this using simple pragmatic approaches. In Section 6 we represent (scalar) $\tau$ as the transformation of an unconstrained variable. For most examples in Section 7 the designs remain in $\mathcal{T}$ under unconstrained optimisation, so no modification to this is needed. We address constraints for analyses that require a minimum time between observations in Section 7 or in Section 8 by adding a large penalty to $\mathcal{K}$ for $\tau\not\in\mathcal{T}$, whose gradient moves designs 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}$ (Kushner and Yin,, 2003).

5.4.4 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 point exchange algorithm. Suppose we must select multiple design points from some region. Optimal designs often have a high degree of replication: they consist of a small number of clusters of repeated observations (Gotwalt et al.,, 2009; Overstall and Woods,, 2017; Binois et al.,, 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 3 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 E. There is scope for developing more sophisticated approaches in future work.

6 Poisson model 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.

Analytic results

Optimal designs for this example can be derived analytically. We summarise the results here: see Appendix H 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 produces no information on $\theta_{2}$ i.e. the posterior always equals the prior. The SIG/ADV design avoids this undesirable property. Also, the SIG/ADV design is invariant to linear reparameterisations. However for FIG, linear reparameterisation can change the optimal design to $\tau_{1}=0$, or make all values of $\tau$ optimal.

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 (15) we take the critic’s action to be

[TABLE]

Figure 1 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 limit 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 critic learning rate is varied. In all cases, the trajectories converge on the limit point. However, convergence is much faster as the critic’s learning rate is increased.

7 Pharmacokinetic example

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

7.1 Model

Pharmacokinetics studies the time course of drug absorption, distribution, metabolism, and excretion. Concentration level is observed via samples of fluid – such as blood, plasma or urine – from the subject taken at preplanned time points. The number of observations is constrained by budget and resources, as well as patient comfort and wellbeing. We assume observed concentration, $y_{i}$, at time $\tau_{i}$ (in hours) is distributed as

[TABLE]

and $D=400$. Concentrations at different times are assumed to be independent. We assume independent log normal priors

[TABLE]

and aim to find $15$ observation times in $[0,24]$. Also we treat $\sigma^{2}=0.1$ as known.

A similar model to the above was used by Ryan et al., (2014) and Overstall and Woods, (2017). However we make two modifications to create a simple setting for comparing different methods. First, we omit a multiplicative noise term. Secondly, we do not enforce a 15 minute minimum time gap between consecutive observations, as its implementation would vary between methods making it more difficult to draw fair conclusions. In Section 7.5 we remove these modifications and show results from our approach for the realistic model used in previous work.

7.2 Methods

The FIM for this model is available in closed form – see Appendix I.1. This allows the optimal ADV design to be found using GDA (Algorithm 2), and the optimal FIG design 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. We performed 100 replications in parallel with each algorithm, based on 100 initial designs sampled from a uniform distribution over $[0,24]^{15}$. We considered several choices of $K$, the number samples to estimate $\mathcal{K}$ in (16), and selected $K=1$. See Appendix I.2 for details.

We compare our results to two methods of finding the optimal SIG design: the approximate coordinate exchange (ACE) algorithm of Overstall and Woods, (2017), and the prior contrastive estimation (PCE) algorithm of Foster et al., (2020). (PCE is one of several methods in Foster et al., 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 I.3 and I.4. PCE allowed 100 replications to be run in parallel. We ran replications of ACE serially, 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.

7.3 Results

Table 1 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’s results.

ADV

The top half of Figure 2 shows trajectories of $\tau$ and $A$ during a single replication of GDA optimisation. 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 3 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{\mathcal{J}}_{\text{ADV}}$ 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 to find global optima.

Figure 4 (left) shows $\hat{\mathcal{J}}_{\text{ADV}}$ values achieved by the final designs from GDA and the SIG methods, as well as for designs sampled from a uniform distribution over $[0,24]^{15}$. These estimated were calculated using equation (17). On average GDA outperforms the SIG methods and further improvement is achieved by also applying point exchange.

Finally, Figure 5 illustrates a posterior found from one particular ADV design, and shows that it is concentrated around the true parameters.

FIG

Full details of SGD optimisation results are given in Appendix I.5. 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 5 illustrates that this design indeed performs poorly by presenting a typical posterior. The narrowness of the FIG design posterior illustrates that the posterior is highly concentrated for some function of $\theta_{1}$ and $\theta_{3}$. However the posterior is diffuse for $\theta_{1}$ and $\theta_{3}$ marginally: it stretches almost the full length of the prior distribution, unlike the highly concentrated ADV design posterior.

SIG

Running ACE under the default tuning choices took a few hours, much slower than GDA which 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 I.3 for more details. (We also explored varying one tuning choice in ACE: the number of Monte Carlo samples for utility estimates used to fit the Gaussian process. However we found only marginal improvements over the recommended default tuning for ACE.)

We ran PCE for roughly 10 times longer than GDA and found the results concentrated near 3 observation times similar to those for GDA. 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 Foster et al., (2020) but were not able to improve on the PCE results. See Appendix I.4 for more details.

Figure 4 (right) shows estimates of expected Shannon information gain achieved by the final designs from all SIG methods, as well as GDA designs and designs sampled from a uniform distribution over $[0,24]^{15}$. The calculation method is described in Appendix I.6. On average GDA outperforms the other methods, with a slight further improvement from using point exchange. This suggests that GDA designs produce good performance under the SIG objective, and is also further evidence that the ACE and PCE methods have not fully converged to the overall optimum.

7.4 Conclusion of comparisons

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 overly diffuse marginal posteriors for some parameters. The ADV design avoids this drawback, and appears to be similar to the SIG design, and indeed gives competitive performance under the SIG objective. Multiple runs of ADV produced designs with similar cluster locations but varying cluster sizes. Post-processing these designs using point exchange further improved the ADV objective reached, illustrating the importance of this step. To explain the improvement, note that point exchange found two candidates for optimal cluster sizes, including one which was different from any cluster sizes found without post-processing.

7.5 Realistic example

Here we implement our ADV approach on a more realistic version of the pharmacokinetic example. Firstly, we now assume multiplicative noise:

[TABLE]

with $\sigma_{1}^{2}=0.1,\sigma_{2}^{2}=0.01$. Secondly we require gaps i.e. a constraint that observation times are at least $0.25$ hours apart. These changes result in the model used by Ryan et al., (2014) and Overstall and Woods, (2017).

Results exist for the FIM of a model with multiplicative noise (e.g. see Malagò and Pistone,, 2015, equation 23). However a lengthy analytic derivation is required which can easily result in errors. Therefore we implement a method which can be used when the FIM is intractable: Algorithm 3, described in Appendix G. Appendix J gives further implementation details for this example and more comments on the results.

Figure 6 shows the resulting designs. Without the gaps constraint, there are 3 clusters of repeated design points. The cluster locations differ depending on whether or not multiplicative noise is used. When gaps are enforced, the clusters remain at the same locations, but consist of spaced out, rather than repeated, design points.

8 Geostatistical regression example

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

8.1 Model

Consider the following geostatistical regression model. Here a design $\tau$ is a $d\times 2$ matrix whose rows specify measurement locations. We assume normal observations with a linear trend and squared exponential covariance function with a nugget effect, giving

[TABLE]

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).

8.2 Methods

Using (7) the FIM is $\mathcal{I}(\tau)=\tau^{T}\Sigma(\tau)^{-1}\tau$, which does not depend on $\theta$. Hence $\mathcal{I}(\tau)$ is also the expected FIM 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 2 with the Adam update rule for 1000 iterations, which was enough for convergence (see Appendix K). 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, the latter seems unlikely to help as there is little evidence of replicated observations in the results.

8.3 Results

Optimisation took on average 19.4 seconds. Figure 7 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 the identity matrix throughout optimisation, reflecting the symmetry of $\theta_{1}$ and $\theta_{2}$ in the model. Hence FIG would also work well for this example.

9 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 Bernardo, (1979) and Walker, (2016).

In 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 starting from multiple random initial designs, which can easily be done in a parallel version of our algorithm. If locally optimal designs involve clusters of repeated observations, we recommend post-processing using a point exchange algorithm. Although we did not observe any in our examples, GDA dynamics can converge to limit cycles rather than the desired solutions. We recommend checking for cyclic behaviour using trace plots of $\tau$ values produced during optimisation.

9.1 Limitations and future research

Intractable Fisher information

Our methods have assumed that the FIM, and associated gradients, can easily be evaluated. When the FIM cannot be evaluated but the likelihood or score function (5) can, it is possible to produce an unbiased estimate of the FIM. Appendix G describes this and how it can be used for experimental design. Section 7.5 contained an application.

However sometimes the likelihood and score function cannot be evaluated. One common reason is the presence of latent variables, such as nuisance parameters. Appendix G describes difficulties of implementing our approach in this setting, and sketches methods to do so, which we plan to investigate in future research.

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.

Discrete parameters

A limitation of our work is that it does not apply to discrete parameters, since Hyvärinen score and FIM are only defined for continuous parameters. Discrete analogues have been proposed (Dawid et al.,, 2012; Shao et al.,, 2019), which would be interesting to investigate.

Variance reduction

Optimisation efficiency could be increased by reducing the variance of our Monte Carlo gradient estimates. For instance, a reviewer suggests the use of randomised quasi-Monte Carlo, as in Drovandi and Tran, (2018).

Alternative optimisers

Variations on GDA with better convergence guarantees, such as two time scaled update rules (Heusel et al.,, 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 Overstall and Woods, (2017).

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ärinen (Parry et al.,, 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 (Loaiza-Maya et al.,, 2021). Another variation is to allow more freedom to the critic. This could include the ability to make more general, non-linear, reparameterisations, or to condition their reparameterisation on the observations, $y$.

Acknowledgements

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

Appendix A Scoring rule details

Section 3.2 of the main paper describes scoring rules and related concepts. Here, Table 2 summarises two strictly proper scoring rules of particular interest in this paper: logarithmic score and Hyvärinen score (Hyvärinen,, 2005). It also gives their associated entropies and divergences. The Hyvärinen results rely on the following regularity conditions:

1. H1.

$p(\theta)$ and $q(\theta)$ are twice differentiable with respect to all $\theta_{i}$. 2. H2.

$\operatorname{\mathbb{E}}_{\theta\sim p(\theta)}[\,||\nabla\log p(\theta)||^{2}],\operatorname{\mathbb{E}}_{\theta\sim q(\theta)}[\,||\nabla\log q(\theta)||^{2}]$ are finite. 3. H3.

$\nabla p(\theta)\to 0,\nabla q(\theta)\to 0$ for $||\theta||\to\infty$.

It this section we use $\nabla$ without a subscript to represents gradient with respect to $\theta$. Also recall that $||\theta||$ represents the $L_{2}$ norm. Finally, note that several of the results of this paper use the following assumption.

1. A5

Conditions H1-H3 hold for $p(\theta)=\pi(\theta|y;\tau)$ and $q(\theta)=\pi(\theta)$, given any $y,\tau$.

Appendix B Framework summaries

Table 3 summarises and contrasts the decision-theoretic and game-theoretic frameworks from Sections 3 and 4.

Appendix C Proofs

C.1 Proof of Result 1

We will show the following which suffices to give the required result. As in the statement of Result 1, all expectations are with respect to $\pi(\theta,y;\tau)$.

The expected reward from design $\tau$ equals $\operatorname{\mathbb{E}}[\mathcal{U}_{\text{entropy}}]$, assuming the experimenter selects $a$ to maximise their expected utility. 2. 2.

$\operatorname{\mathbb{E}}[\mathcal{U}_{\text{entropy diff}}]$ and $\operatorname{\mathbb{E}}[\mathcal{U}_{\text{divergence}}]$ also equal this value up to an additive constant.

Proofs of both parts are below. The first is a consequence of A1. The proof of the second part is based on the common observation that the divergence contains a constant term which can be ignored under maximisation. Similar results have appeared in the experimental design literature previously for a logarithmic score function (e.g. Shewry and Wynn,, 1987). This result simply generalises them to a general score function.

Note that manipulations of expectations used below are valid by A2 and Fubini’s theorem.

Part 1

First fix some values of $\tau$ and $y$ and consider the experimenter’s choice of $a$. This must maximise the expected reward, which is, using A1,

[TABLE]

Since $\mathcal{S}$ is a strictly proper scoring rule, the optimal choice of $a$ is the posterior $\pi(\theta|y;\tau)$. The resulting expected reward is $\mathcal{U}_{\text{entropy}}(\tau,y)$. (We write $\mathcal{U}_{\text{entropy}}$ having only the arguments $\tau,y$ as it does not depend on $\theta$.)

Now suppose $y$ is randomly generated by nature. Then the expected reward of design $\tau$ is the expectation of $\mathcal{U}_{\text{entropy}}(\tau,y)$ with respect to $f(y|\theta,\tau)$, or equivalently with respect to $\pi(\theta,y;\tau)=\pi(\theta)f(y|\theta;\tau)$, as required.

Part 2

From our definitions,

[TABLE]

We are interested in expectations of utility with respect to $\pi(\theta,y;\tau)$, which equals $\pi(\theta|y;\tau)\pi(y;\tau)$ by (3). Thus

[TABLE]

(Note A2 guarantees the terms above are finite.)

So $\mathcal{U}_{\text{divergence}}$ and $\mathcal{U}_{\text{entropy diff}}$ produce the same expected utility. Furthermore the expected utility of $\mathcal{U}_{\text{entropy}}$ only differs by a finite additive constant, $\mathcal{H}[\pi(\theta)]$.

C.2 Proof of Result 2

It’s sufficient to show that these utilities have the same expectation as $\mathcal{U}_{\text{divergence}}$ with respect to $\pi(\theta,y;\tau)$. From Table 2 (and using A5 for the Hyvärinen score case):

[TABLE]

The quantities within the expectations are $\mathcal{U}_{\text{SIG}}$ and $\mathcal{U}_{\text{FIG}}$ respectively. Thus these have same expectations as $\mathcal{U}_{\text{divergence}}$ under $\pi(\theta,y;\tau)$. Finally note that from (4), $\mathcal{U}_{\text{trace}}=\operatorname{\mathbb{E}}_{y\sim f(y|\theta;\tau)}[\mathcal{U}_{\text{FIG}}]$. So this has the same expectation under $\pi(\theta,y;\tau)$ as $\mathcal{U}_{\text{FIG}}$. (These manipulations of expectations are valid using Fubini’s theorem and assumption A2 for the log score case, or A5 for the Hyvärinen score case.)

C.3 Proof of Result 3

As the logarithmic score is a proper scoring rule, $a(\tau,A,y)$ must output the posterior. Using the change of variables formula, the posterior for $\phi=A^{-1}\theta$ is

[TABLE]

where $\pi_{\Theta}(\theta|y;\tau)$ is the posterior for $\theta$. Using A3, $\det A=1$ so $\pi_{\Phi}(\phi|y;\tau)=\pi_{\Theta}(\theta|y;\tau)$. Hence the expected reward to the experimenter from $\tau$ given $a(\tau,A,y)$ is, regardless of the choice of $A$, $\operatorname{\mathbb{E}}_{\theta,y\sim\pi(\theta,y;\tau)}[-\log\pi_{\Theta}(\theta|y;\tau)]$, the same as in the decision theoretic framework. Thus $\tau$ must optimise the same objective in the decision theoretic setting, or in a SPE of the game theoretic setting.

C.4 Proof of Result 4

Let “game 1” be based on the original parameters $\theta$, and “game 2” use the alternative parameters $B\theta$.

C.4.1 Part 1: Relation between action sequences

Here we show that there are invertible mappings $A_{2}=f(A_{1}),a_{2}=g(a_{1})$ such that actions $(\tau,A_{1},y,a_{1})$ in game 1, and $(\tau,A_{2},y,a_{2})$ in game 2 have the same expected rewards up to a multiplicative constant that depends only on $B$. In game 1, $a_{1}$ is a density $a_{1}(\phi_{1})$ for $\phi_{1}=A_{1}^{-1}\theta$. The expected reward to the experimenter is then $r_{1}=-\operatorname{\mathbb{E}}_{\theta\sim\pi(\theta|y;\tau)}[\mathcal{S}(a_{1},A_{1}^{-1}\theta)]$ where $\mathcal{S}$ is Hyvärinen score.

First suppose $\det B>0$. In this case let $A_{2}=k^{-1/p}BA_{1}$ where $k=\det B$ and $p=\dim(\theta)$. This is a valid action in game 2 since $\det A_{2}=1$. In game 2 the experimenter must select a density on $\phi_{2}=A_{2}^{-1}B\theta=k^{1/p}\phi_{1}$. Let $a_{2}(\phi_{2})=a_{1}(k^{-1/p}\phi_{2})/k$. Now consider the actions $(\tau,A_{2},y,a_{2})$ in game 2. The expected reward to the experimenter is then $r_{2}=-\operatorname{\mathbb{E}}_{\theta\sim\pi(\theta|y;\tau)}[\mathcal{S}(a_{2},k^{1/p}A_{1}^{-1}\theta)]$. From the definition of Hyvärinen score, $r_{2}=k^{-2/p}r_{1}$. (As reparameterising by a scalar preserves Hyvärinen score up to proportionality.)

Now suppose $\det B<0$. (Recall $B$ is invertible so $\det B\neq 0$.) Let $F$ be a $p\times p$ identity matrix with the top left entry changed to $-1$. Now let $A_{2}=k^{-1/p}BA_{1}F$ where $k=|\det B|$. Again $\det A_{2}=1$ so this is a valid action in game 2. In game 2 the experimenter must select a density on $\phi_{2}=A_{2}^{-1}B\theta=k^{1/p}F\phi_{1}$. Let $a_{2}(\phi_{2})=a_{1}(k^{-1/p}F\phi_{2})/k$. Now consider the actions $(\tau,A_{2},y,a_{2})$ in game 2. The expected reward to the experimenter is then $r_{2}^{\prime}=-\operatorname{\mathbb{E}}_{\theta\sim\pi(\theta|y;\tau)}[\mathcal{S}(a_{2},k^{1/p}FA_{1}^{-1}\theta)]$. From the definition of Hyvärinen score, $r_{2}^{\prime}=k^{-2/p}r_{1}$. (As Hyvärinen score is preserved under reparameterisation by axis-aligned reflections.)

C.4.2 Part 2: Relation between subgame perfect equilibria

Consider a SPE of game 1 with decision rules $\tau$, $A(\tau)$, $a(\tau,A,y)$. Define $A^{\prime}(\tau)=f(A(\tau))$ and $a^{\prime}(\tau,A,y)=g(a(\tau,A,y))$. Then $\tau$, $A^{\prime}(\tau)$, $a^{\prime}(\tau,A,y)$ is a SPE of game 2. This follows since any counterexample could be transformed to a counterexample for the original SPE in game 1. Similarly a SPE of game 2 can be transformed to a SPE of game 1 using the inverse mappings. This suffices to prove the main result.

C.5 Proof of Result 5

C.5.1 Part 1: Objective function for $\tau,A$

Initially the proof follows that of Results 1 and 2. First consider the choice of $a(\tau,A,y)$ in a SPE. This must maximise the expected reward to the experimenter, which is, using A1, $-\operatorname{\mathbb{E}}_{\phi\sim\pi_{\Phi}(\phi|y;\tau)}[\mathcal{S}(a,\phi)]$. Since $\mathcal{S}$ is a strictly proper scoring rule, the optimal choice of $a$ is the posterior $\pi_{\Phi}(\phi|y;\tau)$. Hence the optimal $a(\tau,A,y)$ outputs $\pi_{\Phi}(\phi|y;\tau)$, and the resulting expected reward to the experimenter is the negative posterior entropy

[TABLE]

The expected reward to the experimenter for this $a$ decision rule under random $y$ is

[TABLE]

(Noting that A5 implies that the terms above are finite and allows manipulation of expectations using Fubini’s theorem.) Define $\mathcal{K}(\tau,A)$ to be the negative of the first term on the right hand side of (20). Then

[TABLE]

Note that the second line requires assumption A5. Since $\phi=A^{-1}\theta$, applying (6) gives

[TABLE]

Here and in the remainder of the proof, $\pi$ and $\mathcal{I}$ without subscripts represent densities and FIM with respect to $\theta$.

Now consider the second term in (20), $\mathcal{H}[\pi_{\Phi}(\phi)]$. From Table 2 and assumption A5

[TABLE]

where $\phi=A^{-1}\theta$. Using the change of variables formula

[TABLE]

Thus $\mathcal{H}[\pi_{\Phi}(\phi)]$ equals $\mathcal{H}[\pi_{\Theta}(\theta)]$ regardless of $A$.

So under the optimal decision rule $a(\tau,A,y)$, the expected reward to the experimenter is $-\mathcal{K}(\tau,A)$, up to an additive constant. Thus the remainder of the nested optimisation problem to be solved to identify SPEs is, as claimed,

[TABLE]

C.5.2 Part 2: Lower bound on $\det\bar{\mathcal{I}}(\tau^{*})$

Here we prove that in a SPE, $\det\bar{\mathcal{I}}(\tau)>0$. This will be needed in part 4. To show this result, we consider $\tau^{*},A^{*}$ satisfying (22) such that $\det\bar{\mathcal{I}}(\tau^{*})=0$, and derive a contradiction.

Firstly, consider the case $\mathcal{K}(\tau^{*},A^{*})=0$. By A4 there exists some $\tilde{\tau}$ such that $\det\bar{\mathcal{I}}(\tilde{\tau})>0$. Then

[TABLE]

So all the eigenvalues of ${A^{*}}^{T}\bar{\mathcal{I}}(\tilde{\tau})A^{*}$ are positive. Hence

[TABLE]

which contradicts $\tau^{*},A^{*}$ satisfying (22).

Secondly, consider the case $\mathcal{K}(\tau^{*},A^{*})=r<0$. Express $\bar{\mathcal{I}}(\tau^{*})$ as an eigendecomposition $P\Lambda P^{T}$. Let $z$ be the number of zero eigenvalues (so there are $p-z$ non-zero eigenvalues). We must have $z>0$ since $\det\bar{\mathcal{I}}(\tau^{*})=0$. Let $\tilde{\Lambda}(\lambda)$ be matrix resulting from replacing the zero eigenvalues in $\Lambda$ with $\lambda^{1/z}$ and the non-zero eigenvalues with $\lambda^{-1/(p-z)}$. Define $\tilde{A}(\lambda)=P\tilde{\Lambda}(\lambda)P^{T}$. Then $\det\tilde{A}(\lambda)=1$ so it is a valid action. Also $\lim_{\lambda\to\infty}\operatorname{tr}[\tilde{A}(\lambda)^{T}\bar{\mathcal{I}}(\tau^{*})\tilde{A}(\lambda)]=0$, so for $\lambda$ sufficiently large

[TABLE]

which contradicts $\tau^{*},A^{*}$ satisfying (22).

C.5.3 Part 3: Optimal $\mathcal{K}$ given $\tau$

Fix some $\tau$ such that $\det\bar{\mathcal{I}}(\tau)>0$. This part shows that $\max_{A}\mathcal{K}(\tau,A)=-p[\det\bar{\mathcal{I}}(\tau)]^{1/p}$.

From (21), $\mathcal{K}(\tau,A)=-\text{tr}B$ where $B=A^{T}\bar{\mathcal{I}}(\tau)A$. Let $\lambda_{1},\ldots,\lambda_{p}$ be the eigenvalues of $B$. These are all positive: they are non-negative since $B$ is positive definite, and non-zero since $\prod_{i=1}^{p}\lambda_{i}=\det\bar{\mathcal{I}}(\tau)>0$. Consider minimising $\operatorname{tr}B$ subject to $\det A=1$, which implies $\det B=\det\bar{\mathcal{I}}(\tau)$. Equivalently we must minimise $\sum_{i=1}^{p}\lambda_{i}$ subject to $\prod_{i=1}^{p}\lambda_{i}=\det\bar{\mathcal{I}}(\tau)$. It is straightforward that the solution is $\lambda_{i}=[\det\bar{\mathcal{I}}(\tau)]^{1/p}$ for all $i$. This can be realised by taking $A$ as a Cholesky factor of $[\det\bar{\mathcal{I}}(\tau)]^{1/p}\bar{\mathcal{I}}(\tau)^{-1}$. The claimed result follows from this.

C.5.4 Part 4: Conclusion

We will partition the possible $\tau$ values into two types. In both cases we will show that the set of those which appear in SPEs equals the set of those which attain the maximum possible $\det\bar{\mathcal{I}}(\tau)$ value.

First consider $\tau$ such that $\det\bar{\mathcal{I}}(\tau)=0$. Such a $\tau$ does not appear in a SPE (by part 2) and does not maximise $\det\bar{\mathcal{I}}(\tau)$ (from A4). Hence the two sets are empty for this case.

Next consider the case of $\tau$ such that $\det\bar{\mathcal{I}}(\tau)>0$. Part 1 shows that SPE actions for $\tau,A$ are solutions to $\min_{\tau}\max_{A}\mathcal{K}(\tau,A)$. Part 3 shows that, for this case, $\max_{A}\mathcal{K}(\tau,A)=-p[\det\bar{\mathcal{I}}(\tau)]^{1/p}$. Thus for this case, the set of designs which appear in SPEs equals those maximising $\det\bar{\mathcal{I}}(\tau)$, as required.

C.6 Proof of Result 6

From (13) we have

[TABLE]

The regularity conditions allow interchange of differentiation and expectation (see Appendix F) so

[TABLE]

An unbiased Monte Carlo estimate of this is

[TABLE]

as required.

Appendix D Restricting $\det A$

This appendix justifies introducing assumption A3 in the game theoretic setting, which enforces $\det A=1$. First recall that we required $A$ to be invertible. For non-invertible $A$, $\phi=A^{-1}\theta$ would not be defined, and this is essential in the description of the game theoretic framework. Hence we must have $\det A\neq 0$.

Some further restriction on $A$ is necessary for a SPE to exist, at least for the scoring rules we consider. First consider the logarithmic score case. From (19) we have

[TABLE]

Without a further restriction on $A$, the critic can always improve their expected reward, $-\operatorname{\mathbb{E}}_{\phi\sim\pi_{\Phi}(\phi|y;\tau)}[\log\pi_{\Phi}(\phi|y;\tau)]$, by replacing $A$ with $\lambda A$ for $\lambda<1$. (This makes the posterior for $\phi$ flatter, benefiting the critic.)

Next consider the Hyvärinen score case, and the critic’s expected reward in a SPE after the experimenter plays $\tau$. From Section C.5 (part 1) this is

[TABLE]

Again, without a further restriction on $A$, the critic can always improve their expected reward by replacing $A$ with $\lambda A$ for $\lambda<1$. (This assumes that the expected reward is not already zero, which is a consequence of assumption A4.)

The above suggests that an additional restriction on $A$ is needed beyond $\det A\neq 0$. Assumption A3, $\det A=1$, is a natural restriction to use since $\det A$ appears in (23).

Appendix E Point exchange algorithms

As discussed in Section 5.4 we often wish to find a design made up of multiple design points from some region, and the optimal design is a small number of clusters of repeated observations. Gradient based optimisation can usually locate the optimal cluster locations well (i.e. determine where the clusters should be located), and here we present a simple point exchange algorithm (Atkinson et al.,, 2007; Overstall and Woods,, 2017) to find approximately optimal cluster sizes (i.e. determine how many design points should be at each location).

The method seeks to optimise one of the diagnostics introduced in Section 5.4,

[TABLE]

These are Monte Carlo estimates of $\mathcal{J}_{\text{FIG}}(\tau)$ from (12) and $\mathcal{J}_{\text{ADV}}(\tau)$ from (14). The $\tilde{\theta}^{(j)}$ terms are independent samples from the prior, with the same samples being reused in all evaluations of $\hat{\mathcal{J}}_{\text{FIG}}$ or $\hat{\mathcal{J}}_{\text{ADV}}$. Our implementation uses $J=1000$. We then seek to maximise the deterministic function $\hat{\mathcal{J}}_{\text{FIG}}(\tau)$ or $\hat{\mathcal{J}}_{\text{ADV}}(\tau)$.

We perform several iterations of greedy optimisation. The design $\tau^{(0)}$ is initialised as the output of gradient based optimisation. Iteration $i$ considers every possible candidate design formed from $\tau^{(i-1)}$ by replacing one design point with another. The candidate design optimising $\hat{\mathcal{J}}_{\text{FIG}}$ or $\hat{\mathcal{J}}_{\text{ADV}}$ is used as $\tau^{(i)}$. The algorithm terminates if $\tau^{(i)}=\tau^{(i-1)}$ or if a maximum number of iterations is exceeded.

Note that this is a simple approximate optimisation method, and there is scope for more sophisticated algorithms to be developed for use here.

Appendix F Regularity conditions

Several of our derivations, such as Result 6, require weak regularity conditions on $g(\theta,\tau)$ to allow interchange of differentiation and expectation i.e.

[TABLE]

Sufficient regularity conditions are that for all $\tau$,

$\operatorname{\mathbb{E}}_{\theta\sim\pi(\theta)}[\,|g(\theta,\tau)|]$ is finite 2. 2.

$g(\theta,\tau)$ is differentiable in each component of $\tau$ 3. 3.

$\operatorname{\mathbb{E}}_{\theta\sim\pi(\theta)}[\,|\nabla_{\theta}g(\theta,\tau)|]$ is finite

where absolute value $|\cdot|$ acts elementwise. The required equality then follows using Fubini’s theorem.

Appendix G Intractable Fisher information

This appendix considers how to apply our approach when the FIM is not available in a closed form that can easily be evaluated.

First we consider two cases where extending our approach is relatively straightforward, and give an algorithm which does so. Section G.1 looks at the case where the score function is available, and Section G.2 at the case where the likelihood is available. In the latter we also discuss scope for future work to improve the efficiency of the algorithm’s implementation in code.

Secondly we discuss two more complex cases which we plan to consider in future work. At the request of a reviewer, we sketch methods to do so, to show that this is feasible. Section G.3 looks at the case where there are latent variables, such as nuisance parameters, and Section G.4 at the case where the likelihood is intractable.

Recap

Recall that our ADV approach is to find minimax solutions of

[TABLE]

using gradient descent ascent. As in the main paper we take $A=A(\eta)$, following (15). GDA requires unbiased estimates of $\nabla_{z}\mathcal{K}(\tau,A)$ for $z\in\{\tau,\eta\}$.

G.1 Estimation using the score function

First we consider the case where the score function $u(y,\theta;\tau)=\nabla_{\theta}\log f(y|\theta;\tau)$ can be evaluated. We show that it is possible to produce unbiased estimates of the required gradients, provided assumption A6, stated below, is met. This allows us to perform optimisation using GDA, as described below in Algorithm 3.

Estimation of objective

Recall that (4) defines

[TABLE]

Substituting this into (24) gives

[TABLE]

which can easily be estimated unbiasedly by Monte Carlo.

Estimation of gradients

Getting unbiased estimators of $\nabla_{z}\mathcal{K}(\tau,A)$ is more complicated. This is because on taking the gradient of (26), it is generally not possible to exchange the order of gradient and expectation operators, as the distribution for $y$ depends on $\tau$. Therefore we use the pathwise derivative approach (see Mohamed et al.,, 2020 for a review) based on the following assumption:

1. A6

The observations $y$ can be generated as a transformation $y(\epsilon,\theta,\tau)$ of a random variable $\epsilon$ with fixed density $p(\epsilon)$.

This gives

[TABLE]

Now it is possible to exchange expectation and differentiation, assuming appropriate regularity conditions (see Appendix F),

[TABLE]

Hence an unbiased Monte Carlo estimate of $\mathcal{K}(\tau,A)$ is

[TABLE]

where $(\theta^{(k)},\epsilon^{(k)})$ are independent draws from $\pi(\theta)p(\epsilon)$. Also, it follows from the above that $\nabla_{z}\breve{\mathcal{K}}(\tau,A)$ is an unbiased estimate of $\nabla_{z}\mathcal{K}(\tau,A)$ for $z\in\{\tau,\eta\}$.

In some cases A6 does not hold, as $y$ cannot be represented as a suitable transformation of $\epsilon$. This includes the case of discrete $y$ and other common distributions such as Beta and Gamma. Various techniques are available to deal with these cases such as taking a continuous approximation to discrete variables (Maddison et al.,, 2017) or using implicit differentiation (Figurnov et al.,, 2018; Jankowiak and Obermeyer,, 2018).

Algorithm

Algorithm 3 presents our approach when the score function can be evaluated, based on the use of pathwise derivative gradient estimates.

G.2 Estimation using the likelihood

Here we consider the case where the likelihood can be evaluated easily but the score function cannot. In this case whenever the score function is required it can simply be calculated by automatic differentiation of the log likelihood, and we can use Algorithm 3. Section 7.5 has an example using this approach, with further details in Appendix J.

Our PyTorch code for this case is restricted to $R=K=1$ i.e. number of parallel replications and Monte Carlo sample size must both equal one. Otherwise an iteration of the algorithm requires calculating likelihood gradients for several $(\theta,y)$ combinations, then manipulating these, and finally calculating gradients again: it is difficult to perform such a calculation efficiently in PyTorch. In future work it would be interesting to implement this method in other frameworks which provide more flexibility in the use of automatic differentiation, such as Jax (Bradbury et al.,, 2018).

G.3 Latent variables case

Next we discuss a case where the likelihood cannot be evaluated as the observations depend on latent variables $\psi$. The latent variables could represent:

• •

Nuisance parameters that we do not wish to learn from our experiment, but which affect the distribution of $y$ nonetheless.

• •

Unobserved states in a time series model.

• •

Model parameters when we are interested only in model choice.

Below we sketch an approach to experimental design in the presence of latent variables. This illustrates the difficulty of this case but that progress is feasible. We plan to investigate this further in future work.

Estimation of score

The observations now depend on latent variables $\psi$, so their density is $f(y|\theta,\psi;\tau)$. Let the conditional prior density of $\psi$ be $\pi(\psi|\theta)$. We would like to work with the marginalised likelihood

[TABLE]

and calculate the corresponding marginalised score function $\nabla_{\tau}\log f(y|\theta;\tau)$ and the resulting Fisher information matrix. Unfortunately the marginalised likelihood involves an intractable integral, which typically cannot be evaluated.

However an unbiased estimate of the marginalised score function exists, using Fisher’s identity (Cappé et al.,, 2006; Poyiadjis et al.,, 2011)

[TABLE]

where $\pi(\psi|\theta,y;\tau)$ is a posterior for the latent variables conditional on $\theta$.

Difficulties

The two difficulties in implementing the ADV approach with latent variables are:

Inference for $\pi(\psi|\theta,y;\tau)$. 2. 2.

Getting unbiased estimates of $\nabla_{\tau}u(y,\theta;\tau)$.

The latter is required in estimating $\nabla_{\tau}\mathcal{K}(\tau,A)$, and it is difficult because on taking the gradient of (28), it is not possible to exchange the order of gradient and expectation operators, as the posterior distribution for $\psi$ depends on $\tau$.

Below we sketch an approach which deals with both difficulties.

Estimation of $\psi$ posterior

For now, suppose $\tau$ is fixed. We would like to have access to $\pi(\psi|\theta,y;\tau)$ for any choices of $\theta$ and $y$. An approximate approach is to use a version of neural density estimation. The goal is to find a function $g(\varepsilon,\theta,y)$ which given $\varepsilon\sim N(0,I)$ outputs samples from a close estimate of $\pi(\psi|\theta,y;\tau)$. The function $g$ can be a neural network of suitable form, such as a normalising flow (Papamakarios et al.,, 2021). The neural network’s tuning parameters can be trained by using SGD to maximise the likelihood based on samples of $(\theta,\psi,y)$ from the prior and model (Papamakarios et al.,, 2021). See Foster et al., (2019, 2020); Kleinegesse and Gutmann, (2020) for approaches to Bayesian experimental design under the SIG objective which make use of neural density estimation or closely related methods.

Estimation of objective gradients

Suppose one has access to a function $g(\varepsilon,\theta,y)$ which given $\varepsilon\sim N(0,I)$ produces samples from $\pi(\psi|\theta,y;\tau^{*})$ where $\tau^{*}$ is the optimal design. Fisher’s identity (28) then becomes

[TABLE]

where $\psi$ is a shorthand for555 Note that $\epsilon$ and $\varepsilon$ represent different random vectors. $g(\varepsilon,\theta,y)$. Similarly, we will assume A6 and write $y$ as shorthand for $y(\epsilon,\theta,\tau)$. Substituting (29) into (25) gives

[TABLE]

Here that $\varepsilon_{1}$ and $\varepsilon_{2}$ are independent $N(0,I)$ random vectors, and $\psi_{1},\psi_{2}$ are shorthand for $g(\varepsilon_{1},\theta,y)$ and $g(\varepsilon_{2},\theta,y)$.

The final expression above allows Monte Carlo estimation of the required gradients of the objective $\mathcal{K}(\tau,A)$.

Optimisation algorithm

We can now sketch an algorithm for optimal design, which is to iterate the following steps:

Sample $\theta,\psi,\epsilon$ given the current $\tau$. 2. 2.

Update $g$ by a step of SGD based on training data $(\theta,\psi,y(\epsilon,\theta,\tau))$. 3. 3.

Sample $\varepsilon_{1},\varepsilon_{2}$. 4. 4.

Produce unbiased estimates of $\nabla_{\tau}\mathcal{K}$ and $\nabla_{\eta}\mathcal{K}$. 5. 5.

Update $\tau$ and $\eta$ by SGD steps.

Upon convergence of $\tau$ to a limit point $\tau^{*}$, $g$ will approximate the corresponding posterior for $\psi$. Hence $\tau^{*}$ will approximate an optimal design.

G.4 Intractable likelihood

Finally, consider the case where it is only possible to sample $y$ given $\theta,\tau$. In this case we do not have access to an unbiased estimator of the score function. However it is possible to train an estimate of the likelihood $\hat{f}(y|\theta)$ using neural density estimation methods. This estimate can be used to define a score estimate $\hat{u}(y,\theta)=\nabla_{\theta}\log\hat{f}(y|\theta)$, allowing the approach of the previous subsection to be used. The difference is that now $u$ is estimated purely from simulations, while the previous section also made use of the full likelihood via (28).

Appendix H Poisson example derivations

Here we derive the results stated in Section 6 on the Poisson example. Recall that the design is $\tau\in[0,1]$, and the observations are $y_{1}\sim Poisson(\tau\theta_{1}\omega_{1}),y_{2}\sim Poisson((1-\tau)\theta_{2}\omega_{2})$, and these are independent. We assume $\omega_{1}>\omega_{2}>0$ and that $\theta_{1},\theta_{2}$ have independent $Gamma(2,1)$ priors.

Posterior distribution

Using standard conjugacy relationships we have independent posterior distributions $\theta_{1}|y_{1}\sim Gamma(2+y_{1},1+\tau\omega_{1})$, $\theta_{2}|y_{2}\sim Gamma(2+y_{2},1+(1-\tau)\omega_{2})$. Note that if $\tau=0$ then $y_{1}$ is a point mass at zero, and so the $\theta_{1}$ posterior equals its prior. Similarly if $\tau=1$ then the $\theta_{2}$ posterior equals its prior.

Expected Fisher information matrix

For independent observations depending on different parameters, the FIM is diagonal with entries given by the scalar Fisher informations of the individual parameters (this is a straightforward consequence of the definition (4)). Now consider $y\sim Poisson(\varphi)$. The Fisher information with respect to $\varphi$ is $1/\varphi$. Suppose $\varphi=k\vartheta$. Using (6), the Fisher information with respect to $\vartheta$ is $k^{2}/\varphi=k/\vartheta$.

Hence the FIM for this example is:

[TABLE]

Note that $1/\theta_{i}$ has an inverse gamma (2,1) distribution which has expectation 1. Thus

[TABLE]

ADV design

We have $\det[\bar{\mathcal{I}}_{\theta}(\theta)]=\tau(1-\tau)\omega_{1}\omega_{2}$ which is maximised by the design $\tau=1/2$. Furthermore we have already shown in Result 4 that this design is invariant to linear reparameterisation.

FIG design

We have $\operatorname{tr}[\bar{\mathcal{I}}_{\theta}(\theta)]=\tau\omega_{1}+(1-\tau)\omega_{2}$. Since $\omega_{1}>\omega_{2}$, this is maximised by $\tau=1$.

Next, consider a linear reparameterisation $\phi=B\theta$. Then

[TABLE]

where $\begin{pmatrix}\omega^{\prime}_{1}\\ \omega^{\prime}_{2}\end{pmatrix}$ is the elementwise product of $\operatorname{diag}(B^{-1}B^{-T})$ and $\begin{pmatrix}\omega_{1}\\ \omega_{2}\end{pmatrix}$. The resulting optimal design is

• •

$\tau=1$ if $\omega^{\prime}_{1}>\omega^{\prime}_{2}$

• •

$\tau=0$ if $\omega^{\prime}_{1}<\omega^{\prime}_{2}$

• •

$\tau\in[0,1]$ if $\omega^{\prime}_{1}=\omega^{\prime}_{2}$

and clearly $B$ can be chosen to allow all these possibilities.

SIG design

Without loss of generality suppose $\omega_{1}=\omega_{2}=1$. This is equivalent to a suitable reparameterisation of $\theta$, and the SIG design is invariant to reparameterisation.

For any $\tau<1/2$, let $h=1/2-\tau$ and define $X_{1}\sim Poisson(\tau\theta)$, $X_{2}\sim Poisson(h\theta)$, $X_{3}\sim Poisson(h\theta)$. Now let $S_{i}=\sum_{j=1}^{i}X_{i}$. so that $S_{1}\sim Poisson(\tau\theta)$, $S_{2}\sim Poisson(\theta/2)$, $S_{3}\sim Poisson([1-\tau]\theta)$. Note that $X_{1},X_{2},X_{3}$ are not independent, but are conditionally independent given $\theta$.

Define $v(\mathcal{A})$ as the expected Shannon information gain for parameter $\theta$ upon observing the set of variables $\mathcal{A}$. Using sufficiency, we get $v(S_{1})=v(X_{1})$, $v(S_{2})=v(\{X_{1},X_{2}\})=v(\{X_{1},X_{3}\})$ and $v(S_{3})=v(\{X_{1},X_{2},X_{3}\})$. We prove below that

[TABLE]

so that $v(S_{2})-v(S_{1})>v(S_{3})-v(S_{2})$. Rearranging gives $2v(S_{2})>v(S_{3})+v(S_{1})$.

In our experiment we have two observations $Y_{1},Y_{2}$, which are conditionally independent given $\theta$. It follows from the definition of Shannon information gain that $v(\{Y_{1},Y_{2}\})=v(Y_{1})+v(Y_{2})$. Hence the expected Shannon information gain from the experiment is $v(S_{3})+v(S_{1})$ under design $\tau$ (or $1-\tau$), and $2v(S_{2})$ under design $1/2$. Thus we have shown that the latter is optimal.

Proof of (30)

Here we adapt the proof of proposition 2 from Krause and Guestrin, (2005) (on submodularity of mutual information). Define $\mathcal{X}$ to be a subset of $\{X_{1},X_{2},X_{3}\}$. From Bayes theorem and conditional independence we have that

[TABLE]

Taking the expectation with respect to $\theta$ and all $X$ variables gives

[TABLE]

where $H$ denotes Shannon entropy and the dash denotes conditioning on $\theta$. Applying this result to both sides of (30) gives

[TABLE]

To prove (30) it remains to demonstrate the inequality

[TABLE]

This holds due to the “information can’t hurt” principle (see e.g. Cover and Thomas,, 2012, Theorem 2.6.5), with the non-independence of the $X_{i}$ variables giving strict inequality.

Appendix I Further details of pharmacokinetic example

I.1 Fisher information matrix

The pharmacokinetic model is of the form $y\sim N(x(\theta,\tau),\sigma^{2}I)$ with $\sigma$ fixed and $\dim x=n$. Using (7),

[TABLE]

where $J$ is a Jacobian matrix, with row $i$ column $j$ entry $\frac{\partial x_{i}}{\partial\theta_{j}}$. It remains to state these terms, which are

[TABLE]

I.2 Choice of $K$

Section 7.2 mentions that, when performing the GDA analysis of the pharmacokinetic example, we compared several choices of $K$, the number samples to estimate $\mathcal{K}$ in (16). Figure 8 shows the results of using $K=1,10,100$ by plotting $\mathcal{J}_{\text{ADV}}$ estimates, as defined in (17). It shows that convergence is slowest for $K=100$. We use $K=1$ in the main paper as it appears slightly quicker than $K=10$.

I.3 Results using Overstall and Woods, (2017) approach

Here we give some further details of our implementation of the pharmacokinetic example using the ACE method of Overstall and Woods, (2017) for SIG optimisation. The analysis presented in the main paper uses the default tuning choices in Overstall and Woods, (2017). As mentioned in the main text, we experimented with one tuning choice but found little improvement. The implementation uses the acebayes package (Overstall et al.,, 2019), and our code can be found at https://github.com/dennisprangle/AdversarialDesignCode/R.

The ACE software offers a second phase of point exchange optimisation. Exploratory analysis showed that using this makes little change to the designs found for this example. Therefore to avoid the extra computational cost, we did not use it for our main analysis.

Figure 9 shows the designs produced for 30 runs from different initial designs. The design points are spread out over the entire $[0,24]$ interval, although not uniformly. This suggests that ACE has not completely converged, as the SIG optimisation results using approaches from Foster et al., (2020) (see next subsection) are much less spread out.

I.4 Results using Foster et al., (2020) approaches

The main paper presents results for SIG optimisation for the pharmacokinetic example using methods from Foster et al., (2020). An adaptation of their code to reproduce our analysis is available at https://github.com/dennisprangle/pyro/tree/sgboed-reproduce. Foster et al., (2020) propose several methods. First we discuss prior contrastive estimation (PCE), as we found it converged most quickly for this application, and is therefore the method we report in the main paper to give a lower bound on the speed for methods from this paper.

To investigate convergence of PCE we considered designs produced after 1 and 2 hours. Figure 10 shows these designs, for 100 replications of the optimisation algorithm, which were computed in parallel. The designs after 2 hours are slightly different, suggesting that the method has not yet fully converged. In the main paper we report results after 1 hour, taking this as a lower bound on the time required for convergence.

Figure 10 has 3 clusters of design points, around666 We report these only to the nearest integer, which is less precise than cluster locations for other methods elsewhere in the paper. This is because the design points here are more variable across replications. $t=1,4,14$. The points are more spread out around these clusters than the ADV results from the main paper, but this may be due to incomplete convergence, as discussed shortly.

Foster et al., (2020) propose several alternatives to PCE, and indeed PCE is not their overall recommendation. For instance it is only guaranteed to converge to the optimal design in the limit of a large number of Monte Carlo samples being used to estimate gradients in each iteration. However for this example we found the other methods were slower or less numerically stable, as we describe now.

We also attempted to implement the Barber-Agakov (BA) and adaptive contrastive estimation (ACE2777 We refer to this as ACE2 to distinguish it from the ACE method of Overstall and Woods, (2017) discussed elsewhere in the paper.) methods of Foster et al., (2020). Both methods involve a neural network to map observations $y$ to an approximate posterior for $\theta$, $N(\mu(y),\Sigma(y))$. We used a neural network with 2 fully connected layers of 16 hidden units and ReLU activation functions. A final linear layer produces 9 outputs: 3 are used for $\mu$ and the other 6 for a Cholesky factor for $\Sigma(y)$.

We were unable to implement ACE2 as training produced NaN outputs. So we discuss only the results of BA. We ran 100 replications in parallel. The left plot in Figure 11 shows the results of one well-behaved replication. The design points converge to 3 clusters near $t=0.7,4.5,15.4$ after roughly 5,000 seconds (1.4 hours). These are similar but not identical to the cluster locations for ADV reported in the main paper. However many other replications did not converge so clearly to these locations, as shown by the right plot in Figure 11. Indeed, in a few replications we noticed the design remained very close to its initial values.

The BA runtime for the well-behaved replication is a little longer than the time used for PCE in the main paper – which we ran for 1 hour (3,600 seconds) for 100 replications in parallel – but of a similar order of magnitude. However other replications did not appear to have converged after running for several hours. As discussed by Foster et al., (2020), BA is likely to be somewhat slower than PCE as it needs to optimise a large number of neural network weights in addition to the design.

We expect that with further effort BA and ACE2 can be more successfully implemented for this example. However, we believe that the PCE results reported in the main paper capture a lower bound on the runtime for any of these methods on this example.

I.5 SGD results

Here we consider using the FIG approach rather than ADV for the pharmacokinetic example. To do so, Figure 12 shows results of optimising $\mathcal{J}_{\text{FIG}}$ using the SGD version of Algorithm 2. The top left plot shows the evolution of the design during a single replication. The design points converge to repeated observations at a single time. The top right and bottom left plots display estimated $\mathcal{K}$ and $\mathcal{J}_{\text{FIG}}$ objectives over 100 repeated runs, which show rapid convergence. The bottom right plot shows that in 100 runs, the design points always converge to a single cluster, around time $12$. Point exchange was not implemented as multiple clusters were not observed.

I.6 Shannon information gain

Figure 4 shows expected Shannon information gain values as a diagnostic for pharmacokinetic example design quality. This subsection derives the calculations used for this diagnostic.

From (8) and (9), the expectation of the Shannon information gain utility is, in general

[TABLE]

where $f(y|\theta;\tau)$ is the likelihood and $\pi(y;\tau)$ is the evidence.

The pharmacokinetic model has the form $y=x(\theta,\tau)+\epsilon$ where $\epsilon\sim N(0,\sigma^{2}\mathrm{I})$. Let $\varphi(\epsilon)$ be the corresponding density. Then:

[TABLE]

where $H$ represents Shannon entropy. Using the standard result for the Shannon entropy of the multivariate normal (see e.g. Cover and Thomas,, 2012, Theorem 8.4.1) gives

[TABLE]

since $\dim y=15$ and $\sigma=0.1$. Hence

[TABLE]

The second term can be estimated using a nested Monte Carlo approximation, to give

[TABLE]

Here the $y^{(i)}$s are independent samples from the prior predictive under design $\tau$, and the $\theta^{(j)}$s are independent samples from the prior. To reduce cost and variance, we reuse the same $\theta$ and $\epsilon$ samples in all evaluations of $\hat{\mathcal{J}}_{\text{SIG}}(\tau)$. Following Overstall and Woods, (2017) we take $B=B^{\prime}=1000$.

Appendix J Realistic pharmacokinetic example

Here we describe the analysis used for the realistic version of the pharmacokinetic example described in Section 7.5 of the main paper. Section J.1 describes the changes in methodology required, and Section J.2 gives further discussion of the results, which were presented in the main paper.

J.1 Methodology

Gaps

We wish to impose the constraint that observations must be separated by at least $0.25$ hours. To do so it was sufficient to follow the approach described in Section 5.4 and simply add the following penalty term to our objective:

[TABLE]

where $\tau_{(i)}$ is the $i$th observation time after sorting into increasing order. Thus a penalty is added whenever there are gaps between observation times less than $0.25$ hours.

Multiplicative noise

We also wish to use the multiplicative noise model (18) i.e.

[TABLE]

This is a situation where the likelihood is available, so we use Algorithm 3 as described in Section G.2.

Let $\epsilon$ be a $N(0,I)$ random vector of length 30. Then observations from the multiplicative noise model can be expressed as

[TABLE]

as required by assumption A6, and needed to implement Algorithm 3.

We found analyses using multiplicative noise take longer to converge. The likely explanation for this is that under the multiplicative noise the gradient estimates are more variable, since the Fisher information matrix must now be estimated. Therefore we use $500,000$ iterations here, rather than the $100,000$ iterations used for the other pharmacokinetic analyses. These analyses took approximately half an hour to run, roughly 10 times as long as the analyses without multiplicative noise. Also, as discussed in Section G.2, in this time the multiplicative noise analyses perform only a single replication of experimental design, while the others perform 100 replications in parallel.

J.2 Discussion

Figure 6 in the main paper shows designs found from four analyses, under all combinations of with/without gaps and with/without multiplicative noise. The design with multiplicative noise and gaps is similar to those found in previous work using SIG: see Figure 2c in Overstall and Woods,, 2017.

Close inspection of Figure 6 shows that the design points for multiplicative noise without gaps are slightly more spaced out than those without multiplicative noise. Inspection of trace plots suggest this is due to lack of full convergence of the optimisation algorithm. As mentioned above, this is likely due to more variable gradient estimates.

Each row of Figure 6 shows a design from a single run of a GDA optimisation algorithm. Therefore each gives a local optimum of our ADV objective (or a close approximation to one). As in the other pharmacokinetic examples, we expect that there are multiple local optima which differ in the number of design points in each cluster. Selecting between these would require further work, such as applying point exchange when no gaps are used, or developing a similar post-processing method for the case with gaps.

Appendix K Further details of geostatistical regression example

Figure 13 shows trace plots of $\mathcal{J}_{\text{ADV}}$ estimates, confirming convergence after 1000 iterations. It also shows that the improvement in this objective is largest for small $\ell$.