# Active Discrimination Learning for Gaussian Process Models

Elham Yousefi1, Luc Pronzato2, Markus Hainy1, Werner G. Müller1, Henry P.
Wynn3

(1JKU Linz, 2CNRS, 3London School of Economics)

###### Abstract

The paper covers the design and analysis of experiments to discriminate
between two Gaussian process models, such as those widely used in computer
experiments, kriging, sensor location and machine learning. Two frameworks are
considered. First, we study sequential constructions, where successive design
(observation) points are selected, either as additional points to an existing
design or from the beginning of observation. The selection relies on the
maximisation of the difference between the symmetric Kullback Leibler
divergences for the two models, which depends on the observations, or on the
mean squared error of both models, which does not. Then, we consider static
criteria, such as the familiar log-likelihood ratios and the Fréchet distance
between the covariance functions of the two models. Other distance-based
criteria, simpler to compute than previous ones, are also introduced, for
which, considering the framework of approximate design, a necessary condition
for the optimality of a design measure is provided. The paper includes a study
of the mathematical links between different criteria and numerical
illustrations are provided.

Keywords: model discrimination; Gaussian random field; kriging

##  1 Introduction

The term ‘active learning’ (cf. Hino, (2020) for a recent review) has replaced
the traditional (sequential or adaptive) ‘design of experiments’ in the
computer science literature, typically when the response is approximated by
Gaussian process regression (GPR, cf. Sauer et al., (2022)). It refers to
selecting the most suitable inputs to achieve the maximum of information from
the outputs, usually with the aim of improving prediction accuracy. A good
overview is given in Chapter 6 of Gramacy, (2020).

Frequently the aim of an experiment – in the broad sense of any data
acquisition exercise – may rather be the discrimination between two or more
potential explanatory models. When data can be sequentially collected during
the experimental process, the literature goes back to the classic procedure of
Hunter and Reiner, (1965) and has generated ongoing research (see e.g. Schwaab
et al., (2008), Olofsson et al., (2018) and Heirung et al., (2019)). When the
design needs to be fixed before the experiment and thus no intermediate data
will be available, the literature is less developed. While in the classical
(non)linear regression case the criterion of T-optimality (cf. Atkinson and
Fedorov, (1975)) and the numerous papers extending it was a major step, a
similar breakthrough for Gaussian process regression is lacking.

With this paper we would like to investigate various sequential/adaptive and
non-sequential design schemes for GPRs and their relative properties. When the
observations associated with the already collected points are available, one
may base the criterion on the predictions and prediction errors (Section 3.1).
On the one hand, one natural choice will be to put the next design point where
the symmetric Kullback-Leibler divergence between those two predictive
(normal) distributions differs most. On the other hand, when the associated
observations are not available, the incremental construction of the designs
could be based on the mean squared error (MSE) for both models, assuming in
turn that either of the two models is the true one (Section 3.2). The static
construction of a set of optimal designs of given size for nominal model
parameters is the last mode we have considered (Section 4). Our first choice
is to use the difference between the expected values of the log likelihood
ratios, assuming in turn that either of the two models is the true one. This
is actually a function of the symmetric Kullback-Leibler divergence, which
also arises from Bayesian considerations. In a similar spirit, the Fréchet
distance between two covariance matrices provides another natural criterion.
Some further novel but simple approaches are considered in this paper as well.
In particular we are interested whether complex likelihood-based criteria like
the Kullback-Leibler-divergence can be effectively replaced by simpler ones
based directly on the respective covariance kernels. The construction of
optimal design measures for model discrimination (approximate design theory)
is considered in Section 5.

Eventually, to compare the discriminatory power of the resulting designs from
different criteria, one can compute the correct classification (hit) rates
after selecting the model with the higher likelihood value. A numerical
illustration is provided in Section 6 for two Matérn kernels with different
smoothness.

##  2 Notation

One of the most popular design criteria for discriminating between rival
models is T-optimality (Atkinson and Fedorov,, 1975). This criterion is only
applicable when the observations are independent and normally distributed with
a constant variance. López-Fidalgo et al., (2007) generalised the normality
assumption and developed an optimal discriminating design criterion to choose
among non-normal models. The criterion is based on the log-likelihood ratio
test under the assumption of independent observations. We denote by
$\varphi_{0}(y,x,\theta_{0})$ and $\varphi_{1}(y,x,\theta_{1})$ the two rival
probability density functions for one observation $y$ at point $x$. The
following system of hypotheses might be considered:

$\displaystyle H_{0}:\varphi(y,x)=\varphi_{0}(y,x,\theta_{0})$ $\displaystyle
H_{1}:\varphi(y,x)=\varphi_{1}(y,x,\theta_{1})$

where $\varphi_{1}(y,x,\theta_{1})$ is assumed to be the true model. A common
test statistic is the log-likelihood ratio given as

$\displaystyle
L=-\log\dfrac{\varphi_{0}(y,x,\theta_{0})}{\varphi_{1}(y,x,\theta_{1})}=\log\dfrac{\varphi_{1}(y,x,\theta_{1})}{\varphi_{0}(y,x,\theta_{0})},$

where the null hypothesis is rejected when
$\varphi_{1}(y,x,\theta_{1})>\varphi_{0}(y,x,\theta_{0})$ or equivalently when
$L>0$. The power of the test refers to the expected value of the log-
likelihood ratio criterion under the alternative hypothesis $H_{1}$. We have

$\displaystyle\mathsf{E}_{H_{1}}(L)=\mathsf{E}_{1}(L)$ $\displaystyle=$
$\displaystyle\int\varphi_{1}(y,x,\theta_{1})\log\left\\{\dfrac{\varphi_{1}(y,x,\theta_{1})}{\varphi_{0}(y,x,\theta_{0})}\right\\}\mathrm{d}y$
(1) $\displaystyle=$ $\displaystyle D_{KL}(\varphi_{1}\rVert\varphi_{0}),$

where $D_{KL}(\varphi_{1}\rVert\varphi_{0})$ is the Kullback–Leibler distance
between the true and the alternative model (Kullback and Leibler,, 1951).

Interchanging the two models in the null and the alternative hypothesis, the
power of the test would be

$\displaystyle\mathsf{E}_{0}(-L)$ $\displaystyle=$ $\displaystyle
D_{KL}(\varphi_{0}\rVert\varphi_{1}).$ (2)

If it is not clear in advance which of the two models is the true model, one
might consider to search for a design optimising a convex combination of (1)
and (2), most commonly using weights $1/2$ for each model. This would be
equivalent to maximising the symmetric Kullback-Leibler distance

$\displaystyle D_{KL}(\varphi_{0},\varphi_{1})$ $\displaystyle=$
$\displaystyle\frac{1}{2}\left[D_{KL}(\varphi_{0}\rVert\varphi_{1})+D_{KL}(\varphi_{1}\rVert\varphi_{0})\right].$

In this paper we will consider random fields, i.e. we will allow for
correlated observations. When the random field is Gaussian, we might still
base the design strategy on the log-likelihood ratio criterion to choose among
two rival models.

For a positive definite kernel $K(x,x^{\prime})$ and an $n$-point design
$\mathbf{X}_{n}=(x_{1},\ldots,x_{n})$, $\mathbf{k}_{n}(x)$ is the
$n$-dimensional vector $(K(x,x_{1}),\ldots,K(x,x_{n}))^{\top}$ and
$\mathbf{K}_{n}$ is the $n\times n$ (kernel) matrix with elements
$\\{\mathbf{K}_{n}\\}_{i,j}=K(x_{i},x_{j})$. Although $x$ is not bold, it may
correspond to a point in a (compact) set ${\mathscr{X}}\subset\mathds{R}^{d}$.
Assume that $Y(x)$ corresponds to the realisation of a random field $Z_{x}$,
indexed by $x$ in ${\mathscr{X}}$, with zero mean $\mathsf{E}\\{Z_{x}\\}=0$
for all $x$ and covariance
$\mathsf{E}\\{Z_{x}Z_{x^{\prime}}\\}=K(x,x^{\prime})$ for all
$(x,x^{\prime})\in{\mathscr{X}}^{2}$. Our prediction of a future observation
$Y(x)$ based on observations
$\mathbf{Y}_{n}=(Y(x_{1}),\ldots,Y(x_{n}))^{\top}$ corresponds to the best
linear unbiased predictor (BLUP)
$\widehat{\eta}_{n}(x)=\mathbf{k}_{n}^{\top}(x)\mathbf{K}_{n}^{-1}\mathbf{Y}_{n}$.
The associated prediction error is $e_{n}(x)=Y(x)-\widehat{\eta}_{n}(x)$ and
we have

$\displaystyle\mathsf{E}\\{e_{n}^{2}(x)\\}=\rho_{n}^{2}(x)=K(x,x)-\mathbf{k}_{n}^{\top}(x)\mathbf{K}_{n}^{-1}\mathbf{k}_{n}(x)\,.$

The index $n$ will often be omitted when there is no ambiguity, and in that
case $\mathbf{k}_{i}(x)=\mathbf{k}_{n,i}(x)$,
$\mathbf{K}_{i}=\mathbf{K}_{n,i}$, $e_{i}(x)=e_{n,i}(x)$,
$\rho_{i}^{2}(x)=\rho_{n,i}^{2}(x)$ will refer instead to model $i$, with
$i\in\\{0,1\\}$. We shall need to distinguish between the cases where the
truth is model $0$ or model $1$, and following Stein, (1999, p. 58) we denote
by $\mathsf{E}_{i}$ the expectation computed with model $i$ assumed to be
true. We reserve the notation $\rho_{i}^{2}(x)$ to the case where the
expectation is computed with the true model; i.e.,

$\displaystyle\rho_{i}^{2}(x)=\mathsf{E}_{i}\\{e_{i}^{2}(x)\\}\,.$

Hence we have
$\rho_{0}^{2}(x)=\mathsf{E}_{0}\\{e_{0}^{2}(x)\\}=K_{0}(x,x)-\mathbf{k}_{0}^{\top}(x)\mathbf{K}_{0}^{-1}\mathbf{k}_{0}(x)$
and calculation gives

$\displaystyle\mathsf{E}_{0}\\{e_{1}^{2}(x)\\}$ $\displaystyle=$
$\displaystyle
K_{0}(x,x)+\mathbf{k}_{1}^{\top}(x)\mathbf{K}_{1}^{-1}\mathbf{K}_{0}\mathbf{K}_{1}^{-1}\mathbf{k}_{1}(x)-2\,\mathbf{k}_{1}^{\top}(x)\mathbf{K}_{1}^{-1}\mathbf{k}_{0}(x)\,,$
(3) $\displaystyle\mathsf{E}_{0}\\{[e_{1}(x)-e_{0}(x)]^{2}\\}$
$\displaystyle=$
$\displaystyle\mathsf{E}_{0}\\{e_{1}^{2}(x)\\}-\mathsf{E}_{0}\\{e_{0}^{2}(x)\\}\,,$

with an obvious permutation of indices 0 and 1 when assuming the model 1 is
true to compute $\mathsf{E}_{1}\\{\cdot\\}$.

If model 0 is correct, the prediction error is larger when we use model 1 for
prediction than if we use the BLUP (i.e., model 0). Stein, (1999, p. 58) shows
that the relation

$\displaystyle\frac{\mathsf{E}_{0}\\{e_{1}^{2}(x)\\}}{\mathsf{E}_{0}\\{e_{0}^{2}(x)\\}}=1+\frac{\mathsf{E}_{0}\\{[e_{1}(x)-e_{0}(x)]^{2}\\}}{\mathsf{E}_{0}\\{e_{0}^{2}(x)\\}}$

shown above is valid more generally for models with linear trends. Also of
interest is the assumed mean squared error (MSE)
$\mathsf{E}_{1}\\{e_{1}^{2}(x)\\}$ when we use model 1 for assessing the
prediction error (because we think it is correct) while the truth is model 0,
and in particular the ratio

$\displaystyle\frac{\mathsf{E}_{1}\\{e_{1}^{2}(x)\\}}{\mathsf{E}_{0}\\{e_{1}^{2}(x)\\}}=\frac{K_{1}(x,x)-\mathbf{k}_{1}^{\top}(x)\mathbf{K}_{1}^{-1}\mathbf{k}_{1}(x)}{\mathsf{E}_{0}\\{e_{1}^{2}(x)\\}}\,,$

which may be larger or smaller than one.

Another important issue concerns the choice of covariance parameters in
$K_{0}$ and $K_{1}$. Denote
$K_{i}(x,x^{\prime})=\sigma_{i}^{2}\,C_{i,\theta_{i}}(x,x^{\prime})$, $i=0,1$,
$(x,x^{\prime})\in{\mathscr{X}}^{2}$, where the $\sigma_{i}^{2}$ define the
variance, the $\theta_{i}$ may correspond to correlation lengths in a
translation invariant model and are thus scalar in the isotropic case, and
$C(x,x^{\prime})$ defines a correlation.

##  3 Prediction-based discrimination

For the incremental construction of a design for model discrimination, points
are added conditionally on previous design points. We can distinguish the case
where the observations associated with those previous points are available and
can thus be used to construct a sequence of predictions (sequential, i.e.,
conditional, construction) from the unconditional case where observations are
not used.

###  3.1 Sequential (conditional) design

Consider stage $n$, where $n$ design points $\mathbf{X}_{n}$ and $n$
observations $\mathbf{Y}_{n}$ are available. Assuming that the random field is
Gaussian, when model $i$ is true we have
$Y(x)\sim{\mathscr{N}}(\widehat{\eta}_{n,i}(x),\rho_{n,i}^{2}(x))$. A rather
natural choice is to choose the next design point $x_{n+1}$ where the
symmetric Kullback-Leibler divergence between those two normal distributions
differs most; that is,

$\displaystyle x_{n+1}$ $\displaystyle\in$
$\displaystyle\mathrm{Arg}\max_{x\in{\mathscr{X}}}\frac{\rho_{n,0}^{2}(x)}{\rho_{n,1}^{2}(x)}+\frac{\rho_{n,1}^{2}(x)}{\rho_{n,0}^{2}(x)}+[\widehat{\eta}_{n,1}(x)-\widehat{\eta}_{n,0}(x)]^{2}\,\left[\frac{1}{\rho_{n,0}^{2}(x)}+\frac{1}{\rho_{n,1}^{2}(x)}\right]\,.$
(4)

Other variants could be considered as well, such as

$\displaystyle x_{n+1}$ $\displaystyle\in$
$\displaystyle\mathrm{Arg}\max_{x\in{\mathscr{X}}}\,[\widehat{\eta}_{n,1}(x)-\widehat{\eta}_{n,0}(x)]^{2}\,,$
$\displaystyle x_{n+1}$ $\displaystyle\in$
$\displaystyle\mathrm{Arg}\max_{x\in{\mathscr{X}}}\,\frac{[\widehat{\eta}_{n,1}(x)-\widehat{\eta}_{n,0}(x)]^{2}}{\rho_{n,0}^{2}(x)+\rho_{n,1}^{2}(x)}\,,$
$\displaystyle x_{n+1}$ $\displaystyle\in$
$\displaystyle\mathrm{Arg}\max_{x\in{\mathscr{X}}}\,[\widehat{\eta}_{n,1}(x)-\widehat{\eta}_{n,0}(x)]^{2}\,\left[\frac{1}{\rho_{n,0}^{2}(x)}+\frac{1}{\rho_{n,1}^{2}(x)}\right]\,.$

They will not be considered in the rest of the paper.

If necessary one can use plug-in estimates $\widehat{\sigma}_{n,i}^{2}$ and
$\widehat{\theta}_{n,i}$ of $\sigma_{i}^{2}$ and $\theta_{i}$, for instance
maximum likelihood (ML) or leave-one-out estimates based on $\mathbf{X}_{n}$
and $\mathbf{Y}_{n}$, when we choose $x_{n+1}$. Note that the value of
$\sigma^{2}$ does not affect the BLUP
$\widehat{\eta}_{n}(x)=\mathbf{k}_{n}^{\top}\mathbf{K}_{n}^{-1}\mathbf{Y}_{n}$.
In the paper we do not address the issues related to the estimation of
$\sigma^{2}$ or of the correlation length or smoothness parameters of the
kernel; one may refer to Karvonen et al., (2020) and the recent papers
Karvonen, (2022); Karvonen and Oates, (2022) for a detailed investigation. The
connection between the notion of microergodicity, related to the consistency
of the maximum-likelihood estimator, and discrimination through a KL
divergence criterion is nevertheless considered in Example 1 below.

###  3.2 Incremental (unconditional) design

Consider stage $n$, where $n$ design points $\mathbf{X}_{n}$ are available. We
base the choice of the next point on the difference between the MSEs for both
models, assuming that one or the other is true. For instance, assuming that
model 0 is true, the difference between the MSEs is
$\mathsf{E}_{0}\\{e_{1}^{2}(x)\\}-\mathsf{E}_{0}\\{e_{0}^{2}(x)\\}=\mathsf{E}_{0}\\{[e_{1}(x)-e_{0}(x)]^{2}\\}=\mathsf{E}_{0}\\{[\widehat{\eta}_{n,1}(x)-\widehat{\eta}_{n,0}(x)]^{2}\\}$.

A first, un-normalised, version is thus

$\displaystyle\phi_{A}(x)$ $\displaystyle=$
$\displaystyle\mathsf{E}_{0}\\{[e_{1}(x)-e_{0}(x)]^{2}\\}+\mathsf{E}_{1}\\{[e_{1}(x)-e_{0}(x)]^{2}\\}\,$
(5) $\displaystyle=$
$\displaystyle\mathsf{E}_{0}\\{e_{1}^{2}(x)\\}-\mathsf{E}_{0}\\{e_{0}^{2}(x)\\}+\mathsf{E}_{1}\\{e_{0}^{2}(x)\\}-\mathsf{E}_{1}\\{e_{1}^{2}(x)\\}\,.$

A normalisation seems in order here too, such as

$\displaystyle\phi_{B}(x)=\frac{\mathsf{E}_{0}\\{[e_{1}(x)-e_{0}(x)]^{2}\\}}{\rho_{n,0}^{2}(x)}+\frac{\mathsf{E}_{1}\\{[e_{1}(x)-e_{0}(x)]^{2}\\}}{\rho_{n,1}^{2}(x)}=\frac{\mathsf{E}_{0}\\{e_{1}^{2}(x)\\}}{\mathsf{E}_{0}\\{e_{0}^{2}(x)\\}}+\frac{\mathsf{E}_{1}\\{e_{0}^{2}(x)\\}}{\mathsf{E}_{1}\\{e_{1}^{2}(x)\\}}-2\,.$
(6)

A third criterion is based on the variation of the symmetric Kullback-Leibler
divergence (11) of Section 4 when adding an $(n+1)$-th point $x$ to
$\mathbf{X}_{n}$. Direct calculation, using

$\displaystyle\mathbf{K}_{n+1,i}=\left(\begin{array}[]{cc}\mathbf{K}_{n,i}&\mathbf{k}_{n,i}(x)\\\
\mathbf{k}_{n,i}^{\top}(x)&K_{i}(x,x)\\\ \end{array}\right)\,,\ i=0,1\,,$

and the expression of the inverse of a block matrix, gives

$\displaystyle\Phi_{KL\,[K0,K1]}(\mathbf{X}_{n}\cup\\{x\\})=\Phi_{KL\,[K0,K1]}(\mathbf{X}_{n})+\frac{1}{2}\,\left[\frac{\mathsf{E}_{1}\\{e_{0}^{2}(x)\\}}{\mathsf{E}_{0}\\{e_{0}^{2}(x)\\}}+\frac{\mathsf{E}_{0}\\{e_{1}^{2}(x)\\}}{\mathsf{E}_{1}\\{e_{1}^{2}(x)\\}}\right]-1\,.$

We thus define

$\displaystyle\phi_{KL}(x)=\frac{1}{2}\,\left[\frac{\mathsf{E}_{1}\\{e_{0}^{2}(x)\\}}{\mathsf{E}_{0}\\{e_{0}^{2}(x)\\}}+\frac{\mathsf{E}_{0}\\{e_{1}^{2}(x)\\}}{\mathsf{E}_{1}\\{e_{1}^{2}(x)\\}}\right]-1\,,$
(8)

to be maximised with respect to $x\in{\mathscr{X}}$.

Although the $\sigma_{i}^{2}$ do not affect predictions,
$\mathsf{E}_{i}\\{e_{j}^{2}(x)\\}$ is proportional to $\sigma_{i}^{2}$. Unless
specific information is available, it seems reasonable to assume that
$\sigma_{0}^{2}=\sigma_{1}^{2}=1$. Other parameters $\theta_{i}$ should be
chosen to make the two kernels the most similar, which seems easier to
consider in the approach presented in Section 4, see (12). In the rest of this
section we suppose that the parameters of both kernels are fixed.

The un-normalised version $\phi_{A}(x)$ given by (5) could be used to derive a
one-step (non-incremental) criterion, in the same spirit as those of Section
4, through integration with respect to $x$ for a given measure $\mu$ on
${\mathscr{X}}$. Indeed, we have

$\displaystyle\mathsf{E}_{0}\\{[e_{1}(x)-e_{0}(x)]^{2}\\}=\mathbf{k}_{0}^{\top}(x)\mathbf{K}_{0}^{-1}\mathbf{k}_{0}(\mathbf{x})+\mathbf{k}_{1}^{\top}(x)\mathbf{K}_{1}^{-1}\mathbf{K}_{0}\mathbf{K}_{1}^{-1}\mathbf{k}_{1}(x)-2\,\mathbf{k}_{1}^{\top}(x)\mathbf{K}_{1}^{-1}\mathbf{k}_{0}(x)\,,$

so that

$\displaystyle\int_{\mathscr{X}}\mathsf{E}_{0}\\{[e_{1}(x)-e_{0}(x)]^{2}\\}\,\mathrm{d}\mu(x)=\mathrm{trace}\left[\mathbf{K}_{0}^{-1}\mathbf{A}_{0}(\mu)+\mathbf{K}_{1}^{-1}\mathbf{K}_{0}\mathbf{K}_{1}^{-1}\mathbf{A}_{1}(\mu)-2\,\mathbf{K}_{1}^{-1}\mathbf{A}_{0,1}(\mu)\right]\,,$

where
$\mathbf{A}_{i}(\mu)=\int_{\mathscr{X}}\mathbf{k}_{i}(x)\mathbf{k}_{i}^{\top}(x)\,\mathrm{d}\mu(x)$,
$i=0,1$, and
$\mathbf{A}_{0,1}(\mu)=\int_{\mathscr{X}}\mathbf{k}_{0}(x)\mathbf{k}_{1}^{\top}(x)\,\mathrm{d}\mu(x)$.
Similarly,

$\displaystyle\int_{\mathscr{X}}\mathsf{E}_{1}\\{[e_{1}(x)-e_{0}(x)]^{2}\\}\,\mathrm{d}\mu(x)=\mathrm{trace}\left[\mathbf{K}_{1}^{-1}\mathbf{A}_{1}(\mu)+\mathbf{K}_{0}^{-1}\mathbf{K}_{1}\mathbf{K}_{0}^{-1}\mathbf{A}_{0}(\mu)-2\,\mathbf{K}_{0}^{-1}\mathbf{A}_{0,1}(\mu)\right]\,.$

The matrices $\mathbf{A}_{i}(\mu)$ and $\mathbf{A}_{0,1}(\mu)$ can be
calculated explicitly for some kernels and measures $\mu$. This happens in
particular when ${\mathscr{X}}=[0,1]^{d}$, the two kernels $K_{i}$ are
separable, i.e., products of one-dimensional kernels on $[0,1]$, and $\mu$ is
uniform on ${\mathscr{X}}$.

#### Example 1: exponential covariance, no microergodic parameters.

We consider Example 6 in Stein, (1999, p. 74) and take
$K_{i}(x,x^{\prime})=\mathsf{e}^{-\alpha_{i}|x-x^{\prime}|}/\alpha_{i}$,
$i=0,1$. The example focuses on two difficulties: first, the two kernels only
differ by their parameter values; second, the particular relation between the
variance and correlation length makes the parameters $\alpha_{i}$ not
microergodic and they cannot be estimated consistently from observations on a
bounded interval; see Stein, (1999, Chap. 6). It is interesting to investigate
the behaviour of the criteria (5), (6) and (8) in this particular situation.

We suppose that $n$ observations are made at $x_{i}=(i-1)/(n-1)$,
$i=1,\ldots,n\geq 2$. We denote $\delta=\delta_{n}=1/[2(n-1)]$ the half-
distance between two design points. The particular Markovian property of
random processes with kernels $K_{i}$ simplifies the analysis. The prediction
and MSE at a given $x\in(0,1)$ only depend on the position of $x$ relative to
its two closest neighbouring design points; moreover, all other points have no
influence. Therefore, due to the regular repartition of the $x_{i}$, we only
need to consider the behaviour in one (any) interval
$\mathds{I}_{i}=[a_{i},b_{i}]=[x_{i},x_{i+1}]$.

We always have $\phi_{A}(x)\to 0$ as $x\to x_{i}\in{\mathscr{X}}_{n}$.
Numerical calculation shows that for $\delta_{n}$ small enough,
$\phi_{A}(\cdot)$ has a unique maximum in $\mathds{I}_{i}$ at the centre
$C_{i}=(x_{i}+x_{i+1})/2$. The next design point $x_{n+1}$ that maximises
$\phi_{A}(\cdot)$ is then taken at $C_{i}$ for one of the $n-1$ intervals, and
we get

$\displaystyle\phi_{A}(C_{i})=\frac{1}{4}\,\frac{(\alpha_{1}-\alpha_{0})^{2}(\alpha_{1}+\alpha_{0})^{3}}{\alpha_{0}\alpha_{1}}\,\delta_{n}^{4}+\mathcal{O}(\delta_{n}^{5})\,,\
n\to\infty\,.$

Similar results apply to the case where the design $\mathbf{X}_{n}$ contains
the endpoints 0 and 1 and its covering radius
$\operatorname{\mathsf{CR}}(\mathbf{X}_{n})=\max_{x\in[0,1]}\min_{i=1,\ldots,n}|x-x_{i}|$
tends to zero, the points $x_{i}$ being not necessarily equally spaced:
$C_{i}$ is then the centre of the largest interval $[x_{i},x_{i+1}]$ and
$\delta_{n}=\operatorname{\mathsf{CR}}(\mathbf{X}_{n})$.

When $\delta_{n}$ is large compared to the correlation lengths $1/\alpha_{0}$
and $1/\alpha_{1}$, there exist two maxima, symmetric with respect to $C_{i}$,
that get closer to the extremities of $\mathds{I}_{i}$ as $\alpha_{1}$
increases, and $C_{i}$ corresponds to a local minimum of $\phi_{A}(\cdot)$.
This happens for instance when $\alpha_{0}\,\delta_{n}=1$ and
$\alpha_{1}\,\delta_{n}\gtrsim 2.600455$.

A similar behaviour is observed for $\phi_{B}(x)$ and $\phi_{KL}(x)$: for
small enough $\delta_{n}$ they both have a unique maximum in $\mathds{I}_{i}$
at $C_{i}$, with now

$\displaystyle\phi_{B}(C_{i})$ $\displaystyle=$
$\displaystyle\frac{1}{4}\,\frac{(\alpha_{1}-\alpha_{0})^{2}(\alpha_{1}+\alpha_{0})^{3}}{\alpha_{0}\alpha_{1}}\,\delta_{n}^{3}+\mathcal{O}(\delta_{n}^{4})\,,\
n\to\infty\,,$ $\displaystyle\phi_{KL}(C_{i})$ $\displaystyle=$
$\displaystyle\frac{1}{8}\,\frac{(\alpha_{1}-\alpha_{0})^{2}(\alpha_{1}+\alpha_{0})^{3}}{\alpha_{0}\alpha_{1}}\,\delta_{n}^{3}+\mathcal{O}(\delta_{n}^{4})\,,\
n\to\infty\,.$

Also, $\phi_{B}(x)\to 0$ and $\phi_{KL}(x)\to 0$ as $x\to
x_{i}\in\mathbf{X}_{n}$. For large values of $\delta_{n}$ compared to the
correlation lengths $1/\alpha_{0}$ and $1/\alpha_{1}$, there exist two maxima
in $\mathds{I}_{i}$, symmetric with respect to $C_{i}$. When
$\alpha_{0}\,\delta_{n}=1$, this happens for instance when
$\alpha_{1}\,\delta_{n}\gtrsim 2.020178$ for $\phi_{B}(\cdot)$ and when
$\alpha_{1}\,\delta_{n}\gtrsim 7.251623$ for $\phi_{KL}(\cdot)$. However, in
the second case the function is practically flat between the two maxima.

The left panel of Figure 1 presents $\phi_{A}(x)$, $\phi_{B}(x)$ and
$\phi_{KL}(x)$ for $x\in[x_{1},x_{2}]=[0,0.1]$ when $n=11$ ($\delta_{n}=0.05$)
and $\alpha_{0}=1$, $\alpha_{1}=10$. The right panel is for
$\alpha_{0}\,\delta_{n}=1$, $\alpha_{1}\,\delta_{n}=10$.

Figure 1: $\phi_{A}(x)$, $\phi_{B}(x)$ and $\phi_{KL}(x)$,
$x\in[x_{1},x_{2}]$, for $n=11$ ($\delta_{n}=0.05$) in Example 1. Left:
$\alpha_{0}=1$, $\alpha_{1}=10$; Right: $\alpha_{0}=20$, $\alpha_{1}=200$.

This behaviour of $\phi_{KL}(C_{i})$ for small $\delta_{n}$ sheds light on the
fact that $\alpha$ is not estimable in this model. Indeed, consider a sequence
of embedded $n_{k}$-point designs $\mathbf{X}_{n_{k}}$, initialised with the
design $\mathbf{X}_{n}=\mathbf{X}_{n_{0}}$ considered above and with
$n_{k}=2^{k}\,(n_{0}-1)+1$, all these designs having the form
$x_{i}=(i-1)/(n_{k}-1)$, $i=1,\ldots,n_{k}$. Then,
$\operatorname{\mathsf{CR}}(\mathbf{X}_{n_{k}})=\operatorname{\mathsf{CR}}(\mathbf{X}_{j})=\delta_{j}=1/[2(n_{k}-1)]$
for $j=n_{k},\ldots,n_{k+1}-1=2\,n_{k}-2$. For $k$ large enough, the increase
in Kullback-Leibler divergence (11) from $\mathbf{X}_{n_{k}}$ to
$\mathbf{X}_{n_{k+1}}$ is thus bounded by $c/(n_{k}-1)^{2}$ for some $c>0$, so
that the expected log-likelihood ratio
$\mathsf{E}_{0}\\{L_{n_{k}}\\}-\mathsf{E}_{1}\\{L_{n_{k}}\\}$ remains bounded
as $k\to\infty$.

More generally, denote by $0\leq x_{1}\leq x_{2}\leq\cdots\leq x_{n}\leq 1$
the ordered points of an $n$-point design $\mathbf{X}_{n}$ in $[0,1]$, $n\geq
3$. Let $i^{*}\geq 3$ be such that
$|x_{i^{*}-2}-x_{i^{*}}|=\min_{i=3,\ldots,n}|x_{i-2}-x_{i}|$. Then necessarily
$|x_{i^{*}-2}-x_{i^{*}}|\leq 1/(\lceil n/2\rceil-1)$. Indeed, consider the
following iterative modification of $\mathbf{X}_{n}$ that cannot decrease
$\min_{i=3,\ldots,n}|x_{i-2}-x_{i}|$: first, move $x_{1}$ to zero, then move
$x_{2}$ to $x_{1}$; leave $x_{3}$ unchanged, but move $x_{4}$ to $x_{3}$, etc.
For $n$ even, the design $\mathbf{X}_{n}^{\prime}$ obtained is the duplication
of an $(n/2)$-points design; for $n$ odd, only the right-most point $x_{n}$
remains single. In the fist case, the minimum distance between points of
$\mathbf{X}_{n}^{\prime}$ is at most $1/(n/2-1)$, in the second case it is at
most $1/(\lceil n/2\rceil-1)$. We then define
$\mathbf{X}_{n-1}=\mathbf{X}_{n}\setminus\\{x_{i^{*}-1}\\}$. For $n$ large
enough, the increase in Kullback-Leibler divergence (11) from
$\mathbf{X}_{n-1}$ to $\mathbf{X}_{n}$ is thus bounded by $c/(\lceil
n/2\rceil-1)^{3}$ for some $c>0$ depending on $\alpha_{0}$ and $\alpha_{1}$.
Starting from some design $\mathbf{X}_{n_{0}}$, we thus have, for $n_{0}$
large enough,

$\displaystyle\Phi_{KL\,[K0,K1]}(\mathbf{X}_{n})-\Phi_{KL\,[K0,K1]}(\mathbf{X}_{n_{0}})\leq
c\sum_{k=n_{0}+1}^{n}\frac{1}{(\lceil k/2\rceil-1)^{3}}\,,$

which implies $\lim_{n\to\infty}\Phi_{KL\,[K0,K1]}(\mathbf{X}_{n})\leq B$ for
some $B<\infty$. Assuming, without any loss of generality, that model 0 is
correct, we have $0\leq\mathsf{E}_{0}\\{L_{n}\\}\leq B$ (we get
$0\leq\mathsf{E}_{1}\\{-L_{n}\\}\leq B$ when we assume that model 1 is
correct), implying in particular that $L_{n}$ does not tend to infinity a.s.
and the ML estimator of $\alpha$ is not strongly consistent.

#### Example 2: exponential covariance, microergodic parameters.

Consider now two exponential covariance models with identical variances (which
we take equal to one without any loss of generality):
$K_{i}(x,x^{\prime})=\mathsf{e}^{-\alpha_{i}|x-x^{\prime}|}$, $i=0,1$.

Again, $\phi_{A}(x)\to 0$ as $x\to x_{i}\in\mathbf{X}_{n}$ and
$\phi_{A}(\cdot)$ has a unique maximum at $C_{i}$ for small enough
$\delta_{n}$, with now

$\displaystyle\phi_{A}(C_{i})=\frac{1}{2}\,(\alpha_{1}^{2}-\alpha_{0}^{2})^{2}\,\delta_{n}^{4}+\mathcal{O}(\delta_{n}^{5})\,,\
n\to\infty\,.$

There are two maxima for $\phi_{A}(\cdot)$ in $\mathds{I}_{i}$, symmetric with
respect to $C_{i}$ for large $\delta_{n}$: when $\alpha_{0}\,\delta_{n}=1$,
this happens for instance when $\alpha_{1}\,\delta_{n}\gtrsim 2.558545$.
Nothing is changed for $\phi_{B}(\cdot)$ compared to Example 1 as the
variances cancel in the ratios that define $\phi_{B}(\cdot)$, see (3) and (6).
The situation is quite different for $\phi_{KL}(\cdot)$, with

$\displaystyle\phi_{KL}(C_{i})=\frac{1}{2}\,\frac{(\alpha_{1}-\alpha_{0})^{2}}{\alpha_{0}\alpha_{1}}+\mathcal{O}(\delta_{n})\,,\
n\to\infty\,,$

indicating that it is indeed possible to distinguish between the two models
much more efficiently with this criterion than with the two others.
Interestingly enough, the best choice for next design point is not at $C_{i}$
but always as close as possible to one of the endpoints $a_{i}$ or $b_{i}$,
with however a criterion value similar to that in the centre $C_{i}$ when
$\delta_{n}$ is small enough, as $\lim_{x\to
x_{i}}\phi_{KL}(x)=(\alpha_{1}-\alpha_{0})^{2}/(2\,\alpha_{0}\alpha_{1})$.
Here, the same sequence of embedded designs as in Example 1 ensures that
$\mathsf{E}_{0}\\{L_{n_{k}}\\}-\mathsf{E}_{1}\\{L_{n_{k}}\\}\to\infty$ as
$k\to\infty$. Figure 2 presents $\phi_{A}(x)$, $\phi_{B}(x)$ and
$\phi_{KL}(x)$ in the same configuration as in Figure 1 but for the kernels
$K_{i}(x,x^{\prime})=\mathsf{e}^{-\alpha_{i}|x-x^{\prime}|}$, $i=0,1$.

Figure 2: $\phi_{A}(x)$, $\phi_{B}(x)$ and $\phi_{KL}(x)$,
$x\in[x_{1},x_{2}]$, for $n=11$ ($\delta_{n}=0.05$) in Example 2. Left:
$\alpha_{0}=1$, $\alpha_{1}=10$; Right: $\alpha_{0}=20$, $\alpha_{1}=200$.

#### Example 3: Matérn kernels.

Take $K_{0}$ and $K_{1}$ as the 3/2 and 5/2 Matérn kernels, respectively:

$\displaystyle K_{0,\theta}(x,x^{\prime})$ $\displaystyle=$
$\displaystyle(1+\sqrt{3}\theta\,|x-x^{\prime}|)\,\exp(-\sqrt{3}\theta\,|x-x^{\prime}|)\
\mbox{ (Mat\'{e}rn 3/2)}\,,$ (9) $\displaystyle K_{1,\theta}(x,x^{\prime})$
$\displaystyle=$
$\displaystyle[1+\sqrt{5}\theta\,|x-x^{\prime}|+5\theta^{2}\,|x-x^{\prime}|^{2}/3]\,\exp(-\sqrt{5}\theta\,|x-x^{\prime}|)\
\mbox{ (Mat\'{e}rn 5/2)}\,.$ (10)

We take $\theta=\theta_{0}=1$ in $K_{0,\theta}$ and adjust $\theta=\theta_{1}$
in $K_{1,\theta}$ to minimise
$\phi_{2\,[K_{0,\theta_{0}},K_{1,\theta_{1}}]}(\mu)$ defined by Eq. (14) in
Section 4 with $\mu$ the uniform measure on $[0,1]$, which gives
$\theta_{1}\simeq 1.1275$. The left panel of Figure 3 shows
$K_{0,\theta_{0}=1}(x,0)$ and $K_{1,\theta}(x,0)$ for $\theta=1$ and
$\theta=\theta_{1}$ when $x\in[0,1]$. The right panel presents $\phi_{B}(x)$
and $\phi_{KL}(x)$ for the same $n=11$-point equally spaced design
$\mathbf{X}_{n}$ as in Example 1 and $x\in[0,1]$ for $K_{0,1}$ and
$K_{1,1.1275}$ (the value of $\phi_{A}(x)$ does not exceed $0.65\,10^{-4}$ and
is not shown). The behaviours of $\phi_{B}(x)$ and $\phi_{KL}(x)$ are now
different in different intervals $[x_{i},x_{i+1}]$ (they remain symmetric with
respect to $1/2$, however), the maximum of $\phi_{KL}(x)$ is obtained at the
central point $x_{5}$. The behaviour of $\phi_{KL}(\cdot)$ could be related to
the fact that discriminating between $K_{0}$ and $K_{1}$ amounts to estimating
the smoothness of the realisation, which requires that some design points are
close to each other.

Figure 3: Left: $K_{0,1}(x,0)$, $K_{1,1}(x,0)$ and $K_{1,1.1275}(x,0)$,
$x\in[0,1]$. Right: $\phi_{B}(x)$ and $\phi_{KL}(x)$ for $x\in[0,1]$ and the
same $11$-point equally spaced design
$\mathbf{X}_{n}=\\{0,1/10,2/10,\ldots,1\\}$ as in Example 1, with $K_{0,1}$
and $K_{1,1.1275}$.

##  4 Distance-based discrimination

We will now consider criteria which are directly based on the discrepancies of
the covariance kernels. Ideally those should be simpler to compute and still
exhibit reasonable efficiencies and some similar properties. The starting
point is again the use of the log-likelihood ratio criterion to choose among
the two models. Assuming that the random field is Gaussian, the probability
densities of observations $\mathbf{Y}_{n}$ for the two models are

$\displaystyle\varphi_{n,i}(\mathbf{Y}_{n})=\frac{1}{(2\pi)^{n/2}\,\det^{1/2}\mathbf{K}_{n,i}}\,\exp\left[-\frac{1}{2}\,\mathbf{Y}_{n}^{\top}\mathbf{K}_{n,i}^{-1}\mathbf{Y}_{n}\right]\,,\
i=0,1\,.$

The expected value of the log-likelihood ratio
$L_{n}=\log\varphi(\mathbf{Y}_{n}|0)-\log\varphi(\mathbf{Y}_{n}|1)$ under
model 0 is

$\displaystyle\mathsf{E}_{0}\\{L_{n}\\}=\frac{1}{2}\log\det(\mathbf{K}_{n,1}\mathbf{K}_{n,0}^{-1})-\frac{n}{2}+\frac{1}{2}\,\mathrm{trace}(\mathbf{K}_{n,0}\mathbf{K}_{n,1}^{-1})$

and similarly

$\displaystyle\mathsf{E}_{1}\\{L_{n}\\}=\frac{1}{2}\log\det(\mathbf{K}_{n,1}\mathbf{K}_{n,0}^{-1})+\frac{n}{2}-\frac{1}{2}\,\mathrm{trace}(\mathbf{K}_{n,1}\mathbf{K}_{n,0}^{-1})\,.$

A good discriminating design should make the difference
$\mathsf{E}_{0}\\{L_{n}\\}-\mathsf{E}_{1}\\{L_{n}\\}$ as large as possible;
that is, we should choose $\mathbf{X}_{n}$ that maximises

$\displaystyle\Phi_{KL\,[K_{0},K_{1}]}(\mathbf{X}_{n})=\mathsf{E}_{0}\\{L_{n}\\}-\mathsf{E}_{1}\\{L_{n}\\}$
$\displaystyle=$
$\displaystyle\frac{1}{2}\left[\mathrm{trace}(\mathbf{K}_{n,0}\mathbf{K}_{n,1}^{-1})+\mathrm{trace}(\mathbf{K}_{n,1}\mathbf{K}_{n,0}^{-1})\right]-n\,$
(11) $\displaystyle=$ $\displaystyle
2\,D_{KL}(\varphi_{n,0},\varphi_{n,1})\,,$

i.e. twice the symmetric Kullback-Leibler divergence between the normal
distributions with densities $\varphi_{n,0}$ and $\varphi_{n,1}$.

We may enforce the normalisation $\sigma_{0}^{2}=\sigma_{1}^{2}=1$ and choose
the $\theta_{i}$ to make the two kernels most similar in the sense of the
criterion $\Phi(\cdot,\cdot)$ considered; that is, maximise

$\displaystyle\min_{\theta_{0}\in\Theta_{0},\,\theta_{1}\in\Theta_{1}}\Phi_{KL\,[K_{0},K_{1}]}(\mathbf{X}_{n})\,.$
(12)

The choice of $\Theta_{0}$ and $\Theta_{1}$ is important; in particular,
unconstrained minimisation over the $\theta_{i}$ could make both kernels
completely flat or on the opposite close to Dirac distributions. It may thus
be preferable to fix $\theta_{0}$ and minimise over $\theta_{1}$ without
constraints. Also, the Kullback-Leibler distance is sensitive to kernel
matrices being near singularity, which might happen if design points are very
close to each other. Pronzato et al., (2019) suggest a family of criteria
based on matrix distances derived from Bregman divergences between functions
of covariance matrices from Kiefer’s $\varphi_{p}$-class of functions
(Kiefer,, 1974). If $p\in(0,1)$, these criteria are rather insensitive to
eigenvalues close or equal to zero. Alternatively, they suggest criteria
computed as Bregman divergences between squared volumes of random
$k$-dimensional simplices for $k\in\\{2,\ldots,d-1\\}$, which have similar
properties.

The index $n$ is omitted in the following and we consider fixed parameters for
both kernels. The Fréchet-distance criterion

$\displaystyle\Phi_{F\,[K_{0},K_{1}]}(\mathbf{X}_{n})=\mathrm{trace}\left[\mathbf{K}_{0}+\mathbf{K}_{1}-2\,(\mathbf{K}_{0}\mathbf{K}_{1})^{1/2}\right],$
(13)

related to the Kantorovich (Wasserstein) distance, seems of particular
interest due to the absence of matrix inversion. The expression is puzzling
since the two matrices do not necessarily commute, but the paper Dowson and
Landau, (1982) is illuminating.

Other matrix “entry-wise” distances will be considered, in particular the one
based on the (squared) Frobenius norm,

$\displaystyle\Phi_{2\,[K_{0},K_{1}]}(\mathbf{X}_{n})=\mathrm{trace}\left(\mathbf{K}_{0}^{2}+\mathbf{K}_{1}^{2}-2\,\mathbf{K}_{0}\mathbf{K}_{1}\right)=\mathrm{trace}\left[(\mathbf{K}_{0}-\mathbf{K}_{1})^{2}\right],$

which corresponds to the substitution of $\mathbf{K}_{i}^{2}$ for
$\mathbf{K}_{i}$ in (13) for $i=0,1$. Denote more generally

$\displaystyle\Phi_{p\,[K_{0},K_{1}]}(\mathbf{X}_{n})=\|\mathbf{K}_{1}-\mathbf{K}_{0}\|_{p}^{p}=\sum_{i,j=1}^{n}|\\{\mathbf{K}_{1}-\mathbf{K}_{0}\\}_{i,j}|^{p}=\mathbf{1}_{n}^{\top}|\mathbf{K}_{1}-\mathbf{K}_{0}|^{\odot
p}\mathbf{1}_{n}\,,\ p>0\,,$

where $\mathbf{1}_{n}$ is the $n$-dimensional vector with all components equal
to 1, the absolute value is applied entry-wise and ⊙p denotes power $p$
applied entry-wise.

Figure 4 shows the values of the criteria $\Phi_{i\,[K_{0,1},K_{1,\theta}]}$,
$i=1,2$, $\Phi_{F\,[K_{0,1},K_{1,\theta}]}$ and
$\Phi_{KL\,[K_{0,1},K_{1,\theta}]}$ as functions of $\theta$ for the two
kernels $K_{0,\theta}$ and $K_{1,\theta}$ given by (9) and (10) and the same
regular design as in Example 1: $x_{i}=(i-1)/(n-1)$, $i=1,\ldots,11$. The
criteria are re-scaled so that their maximum equals one on the interval
considered for $\theta$. Note the similarity between
$\Phi_{2\,[K_{0,1},K_{1,\theta}]}(\mathbf{X}_{n})$ and
$\Phi_{F\,[K_{0,1},K_{1,\theta}]}(\mathbf{X}_{n})$ and the closeness between
the distance-minimising $\theta$ for $\Phi_{1}$, $\Phi_{2}$ and $\Phi_{F}$.
Also note the good agreement with the value $\theta_{1}\simeq 1.1275$ that
minimises $\phi_{2\,[K_{0,1},K_{1,\theta_{1}}]}(\mu)$ from Eq. (14), see
Example 3. The optimal $\theta$ for
$\Phi_{KL\,[K_{0,1},K_{1,\theta}]}(\mathbf{X}_{n})$ is much different,
however, showing that the criteria do not necessarily agree between them.

Figure 4: $\Phi_{i\,[K_{0,1},K_{1,\theta}]}(\mathbf{X}_{n})$, $i=1,2$,
$\Phi_{F\,[K_{0,1},K_{1,\theta}]}(\mathbf{X}_{n})$ and
$\Phi_{KL\,[K_{0,1},K_{1,\theta}]}(\mathbf{X}_{n})$ as functions of
$\theta\in[0.75,3]$ for the same $11$-point equally spaced design
$\mathbf{X}_{n}$ as in Example 1 and $K_{0,\theta}$, $K_{1,\theta}$ given by
(9) and (10), respectively.

An interesting feature of the family of criteria
$\Phi_{p\,[K_{0},K_{1}]}(\cdot)$, $p>0$, is that they extend straightforwardly
to a design measure version. Indeed, defining $\xi_{n}$ as the empirical
measure on the points in $\mathbf{X}_{n}$,
$\xi_{n}=(1/n)\,\sum_{i=1}^{n}\delta_{x_{i}}$, we can write

$\displaystyle\Phi_{p\,[K_{0},K_{1}]}(\mathbf{X}_{n})=n^{2}\,\phi_{p\,[K_{0},K_{1}]}(\xi_{n})\,,$

where we define, for any design (probability) measure on ${\mathscr{X}}$,

$\displaystyle\phi_{p}(\xi)=\phi_{p\,[K_{0},K_{1}]}(\xi)$ $\displaystyle=$
$\displaystyle\int_{{\mathscr{X}}^{2}}|K_{1}(x,x^{\prime})-K_{0}(x,x^{\prime})|^{p}\,\mathrm{d}\xi(x)\mathrm{d}\xi(x^{\prime})\,.$
(14)

Denote by $F_{p\,[K_{0},K_{1}]}(\xi;\nu)$ the directional derivative of
$\phi_{p\,[K_{0},K_{1}]}(\cdot)$ at $\xi$ in the direction $\nu$,

$\displaystyle F_{p\,[K_{0},K_{1}]}(\xi;\nu)=\lim_{\alpha\rightarrow
0^{+}}\frac{\phi_{p\,[K_{0},K_{1}]}[(1-\alpha)\xi+\alpha\nu]-\phi_{p\,[K_{0},K_{1}]}(\xi)}{\alpha}\,.$

Direct calculation gives

$\displaystyle
F_{p\,[K_{0},K_{1}]}(\xi;\nu)=2\left[\int_{{\mathscr{X}}^{2}}|K_{1}(x,x^{\prime})-K_{0}(x,x^{\prime})|^{p}\,\mathrm{d}\nu(x)\mathrm{d}\xi(x^{\prime})-\phi_{p\,[K_{0},K_{1}]}(\xi)\right]\,,$

and thus in particular

$\displaystyle
F_{p\,[K_{0},K_{1}]}(\xi;\delta_{x})=2\left[\int_{\mathscr{X}}|K_{1}(x,x^{\prime})-K_{0}(x,x^{\prime})|^{p}\,\mathrm{d}\xi(x^{\prime})-\phi_{p\,[K_{0},K_{1}]}(\xi)\right]\,.$

One can easily check that the criterion is neither concave nor convex in
general (as the matrix $|\mathbf{K}_{1}-\mathbf{K}_{0}|^{\odot p}$ can have
both positive and negative eigenvalues), but we nevertheless have a necessary
condition for optimality.

######  Theorem 1.

If the probability measure $\xi^{*}$ on ${\mathscr{X}}$ maximises
$\phi_{p\,[K_{0},K_{1}]}(\xi)$, then

$\displaystyle\forall x\in{\mathscr{X}},\
\int_{\mathscr{X}}|K_{1}(x,x^{\prime})-K_{0}(x,x^{\prime})|^{p}\,\mathrm{d}\xi^{*}(x^{\prime})\leq\phi_{p\,[K_{0},K_{1}]}(\xi^{*})\,.$

Moreover,
$\int_{\mathscr{X}}|K_{1}(x,x^{\prime})-K_{0}(x,x^{\prime})|^{p}\,\mathrm{d}\xi^{*}(x^{\prime})=\phi_{p\,[K_{0},K_{1}]}(\xi^{*})$
for $\xi^{*}$-almost every $x\in{\mathscr{X}}$.

This suggests the following simple incremental construction: at iteration $n$,
with $\mathbf{X}_{n}$ the current design and $\xi_{n}$ the associated
empirical measure, choose
$x_{n+1}\in\mathrm{Arg}\max_{x\in{\mathscr{X}}}F_{p\,[K_{0},K_{1}]}(\xi_{n};\delta_{x})=\mathbf{1}_{n}^{\top}|\mathbf{k}_{n,0}(x)-\mathbf{k}_{n,1}(x)|^{\odot
p}$. It will be used in the numerical example of Section 6.2.

##  5 Optimal design measures

In this section we explain why the determination of optimal design measures
maximising $\phi_{p}(\xi)$ is generally difficult, even when limiting
ourselves to the satisfaction of the necessary condition in Theorem 1. At the
same time, we can characterise measures that are approximately optimal for
large $p$.

We assume that the two kernels are isotropic, i.e., such that
$K_{i}(x,x^{\prime})=\Psi_{i}(\|x-x^{\prime}\|)$, $i=0,1$, and that the
functions $\Psi_{i}$ are differentiable except possibly at 0 where they only
admit a right derivative. We define $\psi(t)=|\Psi_{1}(t)-\Psi_{0}(t)|$,
$t\in\mathds{R}^{+}$, and assume that the kernels have been normalised so that
$K_{0}(x,x)=K_{1}(x,x)$; that is, $\psi(0)=0$. Also, we only consider the case
where the function $\psi(\cdot)$ has a unique global maximum on
$\mathds{R}^{+}$. This assumption is not very restrictive. Consider again the
two Matérn kernels (9) and (10). Figure 5 shows the evolution of $\psi^{2}(t)$
for $K_{0}=K_{0,1}$ and $K_{1}=K_{1,\theta_{1}}$ with two different values of
$\theta_{1}$: $\theta_{1}=1$ and $\theta_{1}\simeq 1.1275$; the latter
minimises $\phi_{2\,[K_{0,1},K_{1,\theta}]}(\mu)$ for $\mu$ being the uniform
measure on $[0,1]$.

Figure 5: $\psi^{2}(t)$ for $K_{0}=K_{0,1}$ and $K_{1}=K_{1,\theta_{1}}$ with
two different values of $\theta_{1}$.

In the following, we shall consider normalised functions $\psi(\cdot)$, such
that $\max_{t\in\mathds{R}^{+}}\psi(t)=1$. We denote by $\Delta$ the (unique)
value such that $\psi(\Delta)=1$. On Figure 5, $\Delta\simeq 0.7$ when
$K_{1}=K_{1,1}$.

###  5.1 A simplified problem with an explicit optimal solution

Consider the extreme case where $\psi=\psi_{*}$ defined by

$\displaystyle\psi_{*}(t)=\left\\{\begin{array}[]{ll}1&\mbox{if }t=\Delta,\\\
0&\mbox{otherwise.}\end{array}\right.$ (17)

Note that $\psi_{*}^{p}(t)=\psi_{*}(t)$ for any $p>0$; we can thus restrict
our attention to $p=1$ for the maximisation of $\phi_{p}(\xi)$ defined by
(14); that is, we consider

$\displaystyle\phi_{1}(\xi)=\int_{{\mathscr{X}}^{2}}\psi_{*}(\|x-x^{\prime}\|)\,\mathrm{d}\xi(x)\mathrm{d}\xi(x^{\prime})\,.$

######  Theorem 2.

When $\psi=\psi_{*}$ and ${\mathscr{X}}\subset\mathds{R}^{d}$ is large enough
to contain a regular $d$ simplex with edge length $\Delta$, any measure
$\xi^{*}$ allocating weight $1/(d+1)$ at each vertex of such a simplex
maximises $\phi_{1}(\xi)$, and $\phi_{1}(\xi^{*})=d/(d+1)$.

###### Proof.

Since $\phi_{1}(\xi)=0$ when $\xi$ is continuous with respect to the Lebesgue
measure on ${\mathscr{X}}$, we can restrict our attention to measures without
any continuous component. Assume that $\xi=\sum_{i=1}^{n}w_{i}\delta_{x_{i}}$,
with $w_{i}\geq 0$ for all $i$ and $\sum_{i=1}^{n}w_{i}=1$, $n\in\mathds{N}$.
Consider the graph $\mathcal{G}(\xi)$ having the $x_{i}$ as vertices, with an
edge $(i,j)$ connecting $x_{i}$ and $x_{j}$ if and only if
$\|x_{i}-x_{j}\|=\Delta$. We have

$\displaystyle\phi_{1}(\xi)=\sum_{(i,j)\in\mathcal{G}(\xi)}w_{i}w_{j},$

and Theorem 1 of Motzkin and Straus, (1965) implies that $\phi_{1}(\xi)$ is
maximum when $\xi$ is uniform on the maximal complete subgraph of
$\mathcal{G}(\xi)$. The maximal achievable order is $d+1$, obtained when the
$x_{i}$ are the vertices of a regular simplex in ${\mathscr{X}}$ with edge
length $\Delta$. Motzkin and Straus, (1965) also indicate in their Theorem 1
that $\phi_{1}(\xi^{*})=1-1/(d+1)$. This is easily recovered knowing that
$\mathcal{G}(\xi^{*})$ is fully connected with order $d+1$. Indeed, we then
have

$\displaystyle\phi_{1}(\xi)=\sum_{i=1}^{d+1}w_{i}\sum_{\stackrel{{\scriptstyle
j=1}}{{j\neq i}}}^{d+1}w_{j}=\sum_{\stackrel{{\scriptstyle i,j=1}}{{j\neq
i}}}^{d+1}w_{i}w_{j}=1-\sum_{i=1}^{d+1}w_{i}^{2},$

which is maximum when all $w_{i}$ equal $1/(d+1)$. ∎

###  5.2 Optimal designs for $\psi(t)=|\Psi_{1}(t)-\Psi_{0}(t)|$

The optimal designs of Theorem 2 are natural candidates for being optimal when
we return to the case of interest $\psi(t)=|\Psi_{1}(t)-\Psi_{0}(t)|$. In the
light of Theorem 1, for a given probability measure $\xi$ on ${\mathscr{X}}$,
we consider the function

$\displaystyle\delta_{\xi}(x)=\int_{\mathscr{X}}\psi^{p}(\|x-x^{\prime}\|)\,\mathrm{d}\xi(x^{\prime})-\phi_{p}(\xi),$

which must satisfy $\delta_{\xi}(x)\leq 0$ for all $x\in{\mathscr{X}}$ when
$\xi$ is optimal. For an optimal measure $\xi^{*}$ as in Theorem 2, with
support $x_{1},\ldots,x_{d+1}$ forming a regular $d$-simplex, we have

$\displaystyle\delta_{\xi^{*}}(x)=\frac{1}{d+1}\,\left[\sum_{i=1}^{d+1}\psi^{p}(\|x-x_{i}\|)-d\right].$

One can readily check that $\delta_{\xi^{*}}(x_{i})=0$ for all $i$ (as
$\psi(\|x_{i}-x_{j}\|)=\psi(\Delta)=1$ for $i\neq j$ and $\psi(0)=0$).
Moreover, since $\psi(\cdot)$ is differentiable everywhere except possibly at
zero, when $p>1$ the gradient of $\delta_{\xi^{*}}(x)$ equals zero at each
$x_{i}$. However, these $d+1$ stationary points may sometimes correspond to
local minima — a situation when of course $\xi^{*}$ is not optimal. The left
panel of Figure 6 shows an illustration ($d=2$) for $p=1.5$,
$K_{0}(x,x^{\prime})=\exp(-\|x-x^{\prime}\|)$ and $K_{1}$ being the Matérn 5/2
kernel $K_{1,1}$. The measure $\xi^{*}$ is supported at the vertices of the
equilateral triangle $(0,0),(\Delta,0),(\Delta/2,\sqrt{(}3)\Delta/2)$
(indicated in blue on the figure), with $\Delta\simeq 0.53$ (the value where
$\psi(\cdot)$ is maximum). Here the $x_{i}$ correspond to local minima of
$\delta_{\xi^{*}}(x)$, $\psi(\cdot)$ is not differentiable at zero but $p>1$
so that $\delta_{\xi^{*}}(\cdot)$ is differentiable.

When $p\to\infty$, $\psi^{p}(\cdot)$ approaches the (discontinuous) function
$\psi_{*}(\cdot)$, suggesting that $\xi^{*}$ may become close to being optimal
for $\phi_{p}$ when $p$ is large enough. However, when ${\mathscr{X}}$ is
large, $\xi^{*}$ is never truly optimal, no matter how large $p$ is. Indeed,
suppose that ${\mathscr{X}}$ contains a point $x_{*}$ corresponding to the
symmetric of a vertex $x_{k}$ of the simplex defining the support of $\xi^{*}$
with respect to the opposite face of that simplex. Direct calculation gives

$\displaystyle
L=\|x_{k}-x_{*}\|=2\,\Delta\,\left(\frac{d+1}{2\,d}\right)^{1/2}.$

The right panel of Figure 6 shows an illustration for $K_{0}$ and $K_{1}$
being the Matérn 3/2 and Matérn 5/2 kernels $K_{0,1}$ and $K_{1,1}$,
respectively. The measure $\xi^{*}$ is supported at the vertices of the
equilateral triangle with vertices
$(0,0),(\Delta,0),(\Delta/2,\sqrt{(}3)\Delta/2)$ with now $\Delta\simeq 0.7$.
At the point $x_{*}$, symmetric to $x_{k}$, indicated in red on the figure, we
have

$\displaystyle\delta_{\xi^{*}}(x_{*})$ $\displaystyle=$
$\displaystyle\frac{1}{d+1}\,\left[\sum_{\stackrel{{\scriptstyle i=1}}{{i\neq
k}}}^{d+1}\psi^{p}(\|x_{*}-x_{i}\|)+\psi^{p}(\|x_{*}-x_{k}\|)-d\right]$ (18)
$\displaystyle=$ $\displaystyle\frac{1}{d+1}\,\psi^{p}(L)>0\,,$

where the second equality follows from $\|x_{*}-x_{i}\|=\Delta$ for all $i\neq
k$, implying that $\xi^{*}$ is not optimal. Another, more direct, proof of the
non-optimality of $\xi^{*}$ is to consider the measure $\widehat{\xi}$ that
sets weights $1/(d+1)$ at all $x_{i}\neq x_{k}$ and weights $1/[2(d+1)]$ at
$x_{k}$ and its symmetric $x_{*}$. Direct calculation gives

$\displaystyle\phi_{p}(\widehat{\xi})=\frac{d}{d+1}\,\left(1-\frac{1}{d+1}\right)+\frac{2}{2\,(d+1)}\left[\frac{d}{d+1}+\frac{1}{2\,(d+1)}\,\psi^{p}(L)\right].$

The first term on the right-hand side comes from the $d$ vertices $x_{i}$,
$i\neq k$, each one having weight $1/(d+1)$ and being at distance $\Delta$ of
all other vertices, those having total weight $1-1/(d+1)$. The second term
comes from the two symmetric points $x_{k}$ and $x_{*}$, each one with weight
$1/[2(d+1)]$. Each of these two points is at distance $\Delta$ from $d$
vertices with weights $1/(d+1)$ and at distance $L$ of the other opposite
point with weight $1/[2(d+1)]$. We get after simplification

$\displaystyle\phi_{p}(\widehat{\xi})=\frac{d}{d+1}+\frac{\psi^{p}(L)}{2\,(d+1)^{2}}>\phi_{p}(\xi^{*})=\frac{d}{d+1},$

showing that $\xi^{*}$ is not optimal. Note that, for symmetry reasons, the
design $\widehat{\xi}$ is not optimal for large enough ${\mathscr{X}}$. The
determination of a truly optimal design seems very difficult. In the
simplified problem of Section 5.1, where the criterion is based on the
function $\psi_{*}$ defined by (17), the measures $\xi^{*}$ and
$\widehat{\xi}$ supported on $d+1$ and $d+2$ points, respectively, have the
same criterion value $\phi_{p}(\xi^{*})=\phi_{p}(\widehat{\xi})=d/(d+1)$ for
all $p>0$.

Although $\xi^{*}$ is not optimal, since $\psi(\|x_{*}-x_{k}\|)<1$ (as
$\psi(t)$ takes its maximum value 1 for $t=\Delta$), (18) suggests that
$\xi^{*}$ may be only marginally suboptimal when $p$ is large enough.
Moreover, as the right panel of Figure 6 illustrates, a design $\xi^{*}$
supported on a regular simplex is optimal provided that ${\mathscr{X}}$ is
small enough and $p$ is large enough to make $\delta_{\xi^{*}}(x)$ concave at
each $x_{i}$ (for symmetry reasons, we only need to check concavity at one
vertex). In fact, $p>2$ is sufficient. Indeed, assuming that $p>2$ and that
$\psi(\cdot)$ is twice differentiable everywhere, with second-order derivative
$\psi^{\prime\prime}(\cdot)$, except possibly at zero, direct calculation
gives

$\displaystyle\frac{d^{2}\delta_{\xi^{*}}(x)}{dxdx^{\top}}\bigg{|}_{x=x_{1}}=\frac{1}{d+1}\,\frac{p\,\psi^{p-1}(\Delta)\psi^{\prime\prime}(\Delta)}{\Delta^{2}}\sum_{i=2}^{d+1}(x_{1}-x_{i})(x_{1}-x_{i})^{\top}\,,$

which is negative-definite (since $\psi^{\prime\prime}(\Delta)<0$,
$\psi(\cdot)$ being maximal at $\Delta$). The right panel of Figure 6 gives an
illustration. Note that $p<2$ on the left panel, and the $x_{i}$ correspond to
local minimas of $\delta_{\xi^{*}}(\cdot)$. Figure 7 shows a plot of
$\delta_{\xi^{*}}(x)$ for $p=2$ and $K_{0}$ and $K_{1}$ being the Matérn 3/2
and Matérn 5/2 kernels $K_{0,1}$ and $K_{1,1.07}$, respectively, suggesting
that the form of optimal designs may be in general quite complicated.

Figure 6: Surface plot of $\delta_{\xi^{*}}(x)$ ($x\in\mathds{R}^{2}$), the
support of $\xi^{*}$ corresponds to the vertices of the equilateral triangle
in blue. Left: $K_{0}(x,x^{\prime})=\exp(-\|x-x^{\prime}\|)$ and
$K_{1}=K_{1,1}$ ($\Delta\simeq 0.53$), $p=1.5$; Right: $K_{0}=K_{0,1}$,
$K_{1}=K_{1,1}$ ($\Delta\simeq 0.7$), $p=10$; the red point $x_{*}$ is the
symmetric of the origin $(0,0)$ with respect to the opposite side of the
triangle. Figure 7: Surface plot of $\delta_{\xi^{*}}(x)$
($x\in\mathds{R}^{2}$), the support of $\xi^{*}$ corresponds to the vertices
of the equilateral triangle in blue: $K_{0}=K_{0,1}$, $K_{1}=K_{1,1.07}$
($\Delta\simeq 1.92$), $p=2$.

##  6 A numerical example

###  6.1 Exact designs

In this section, we consider numerical evaluations of designs resulting from
the prediction-based and distance-based criteria. Here, the rival models are
the isotropic versions of the covariance kernels used in Example 3 (Section
3.2) for the design space ${\mathscr{X}}=[0,10]^{2}$, discretised at $n=25$
equally spaced points in each dimension. For an agreement on the setting of
correlation lengths in both kernels, we have applied a minimisation procedure.
Specifically, we have taken $\theta=\theta_{0}=1$ in
$K_{0,\theta}(x,x^{\prime})$ and adjusted the parameter in the second kernel
minimising each of the distance-based criteria for the design
$\mathbf{X}_{625}$ corresponding to the full grid. This resulted in
$\theta_{1}=1.0047$, $1.0285$, $1.0955$ and $1.3403$, respectively, for
$\Phi_{F},\Phi_{1},\Phi_{2}$ and $\Phi_{KL}$. We have finally chosen
$\theta_{1}=1.07$, which seems to be compatible with the above values.

The left panel in Figure 8 shows the plot of the two Matérn covariance
functions at the assumed parameter values. This plot illustrates the
similarity of the kernels which we aim to discriminate. The right panel in the
figure refers to the plot of the absolute difference between the covariance
kernels. The red line corresponds to the distance where the absolute
difference between them is maximal. This is denoted by $\Delta$, which is
equal to $\Delta=1.92$ in this case.

Figure 8: Left: Plot of the Matérn covariance functions at the assumed
parameter setting. Right:
$\psi(t)=|K_{0,\theta_{0}}(t,0)-K_{0,\theta_{1}}(t,0)|,(\theta_{0}=1,\theta_{1}=1.07)$.

The sequential approach is the only case where the observations
$\mathbf{Y}_{n}$ corresponding to the previous design points $\mathbf{X}_{n}$
are used in the design construction. We use this information to estimate the
parameter setting at each step. The (box)plots of the maximum likelihood (ML)
estimates $\hat{\theta}_{0}$ and $\hat{\theta}_{1}$ of the inverse correlation
lengths $\theta_{0}$ and $\theta_{1}$ of $K_{0,\theta}(x,x^{\prime})$ and
$K_{1,\theta}(x,x^{\prime})$, respectively, are presented in Figure 9. This
refers to the case where the first kernel, Matérn 3/2, is the data generator.
The $\hat{\theta}_{0}$ estimates converge to their null value, $\theta_{0}=1$,
drawn as a red dashed line in the left panel of Figure 9, as expected due to
the consistency of the ML estimator in this case. For the second kernel to be
similar to the first one (i.e., less smooth), the $\hat{\theta}_{1}$ estimates
have increased (see the right panel). The decrease of the correlation length
causes the covariance kernel to drop faster as a function of distance. We
defer from presenting the opposite case (where the Matérn 5/2 is the data
generator), which is similar.

Figure 9: Maximum likelihood estimates of the correlation lengths in Matérn
kernels.

Apart from the methods applied in Section 4, we have considered some other
static approaches for discrimination. $D_{s}$-optimal design is a natural
candidate that can be applied in the distance-based fashion. For
$D_{s}$-optimality, we require the general form of the Matérn covariance
kernel, which is based on the modified Bessel function of the second kind
(denoted by $C_{\nu}$). It is given by

$\displaystyle\mathbf{K}_{\nu}(r)=\dfrac{2^{1-\nu}}{\Gamma(\nu)}\left(\sqrt{2\nu}~{}r\theta\right)^{\nu}C_{\nu}\left(\sqrt{2\nu}~{}r\theta\right).$
(19)

Smoothness, $\nu$, is considered as the parameter of interest, while the
correlation length $\theta$ is assumed as nuisance. The first off-diagonal
element in the $2\times 2$ information matrix, associated with the estimation
of parameters $\boldsymbol{\theta}=(\theta,\nu)$, is

$\displaystyle
M(\mathbf{X}_{n},\boldsymbol{\theta})_{12}=\dfrac{1}{2}\mathrm{trace}\left\\{\mathbf{K}_{\nu}^{-1}\dfrac{\partial\mathbf{K}_{\nu}}{\partial\theta}\mathbf{K}_{\nu}^{-1}\dfrac{\partial\mathbf{K}_{\nu}}{\partial\nu}\right\\},$
(20)

see, e.g., Eq. (6.19) in Müller, (2007). The other elements in the information
matrix are calculated similarly. We have used the supplementary material of
Lee et al., (2018) to compute the partial derivatives of the Matérn covariance
kernel. Finally, the $D_{s}$-criterion is

$\displaystyle\Phi_{D_{s}}=|M(\mathbf{X}_{n},\boldsymbol{\theta})|/|M(\mathbf{X}_{n},\boldsymbol{\theta})_{11}|,$
(21)

where $M(\mathbf{X}_{n},\boldsymbol{\theta})_{11}$ is the element of the
information matrix corresponding to the nuisance parameter (i.e., in
$M(\mathbf{X}_{n},\boldsymbol{\theta})_{11}$ both partial derivatives are
calculated with respect to $\theta$). In the examples to follow we consider
local $D_{s}$-optimal design; that is, the parameters $\theta$ and $\nu$ are
set at given values.

From a Bayesian perspective, models can be discriminated optimally when the
difference between the expected entropies of the prior and the posterior model
probabilities is maximised. This criterion underlies a famous sequential
procedure put forward by Box and Hill, (1967) and Hill and Hunter, (1969).
Since such criteria typically cannot be computed analytically, several bounds
were derived. The upper bound proposed by Box and Hill, (1967) is equivalent
to the symmetric Kullback-Leibler divergence $\Phi_{KL}$. Hoffmann, (2017)
derives a lower bound based on a lower bound for the Kullback-Leibler
divergence between a mixture of two normals, which is given by Eq. (24) and is
denoted by $\Phi_{\Gamma}$. Here, we assume equal prior probabilities. A more
detailed account of Bayesian design criteria and their bounds is given in
Appendix A.

Table 1: Comparison of average hit rates in different methods for the first
numerical example.

| Average hit rate
---|---
Design size | $5$ | $6$ | $7$ | $8$ | $9$ | $10$ | $20$ | $30$ | $40$ | $50$
Sequential (4) | $0.500$ | $0.535$ | $0.540$ | $0.595$ | $0.570$ | $0.640$ | $0.695$ | $0.715$ | $0.740$ | $0.770$
$\phi_{A}$ | $0.505$ | $0.500$ | $0.530$ | $0.525$ | $0.505$ | $0.510$ | $0.520$ | $0.535$ | $0.585$ | $0.635$
$\phi_{B}$ | $0.520$ | $0.545$ | $0.575$ | $0.585$ | $0.615$ | $0.650$ | $0.785$ | $0.875$ | $0.900$ | $0.910$
$\phi_{KL}$ | $0.520$ | $0.545$ | $0.575$ | $0.585$ | $0.615$ | $0.650$ | $0.785$ | $0.870$ | $0.915$ | $0.925$
$\Phi_{F}$ | $0.580$ | 0.625 | 0.620 | $0.625$ | $0.670$ | 0.715 | $0.795$ | 0.900 | $0.925$ | $0.950$
$\Phi_{1}$ | $0.525$ | $0.520$ | $0.555$ | $0.540$ | $0.550$ | $0.610$ | $0.725$ | $0.890$ | $0.910$ | $0.920$
$\Phi_{2}$ | $0.525$ | $0.520$ | $0.555$ | $0.540$ | $0.550$ | $0.610$ | $0.715$ | $0.860$ | $0.890$ | $0.910$
$\Phi_{KL}$ | $0.580$ | 0.625 | 0.620 | $0.625$ | $0.670$ | 0.715 | $0.795$ | $0.895$ | $0.925$ | 0.955
$\Phi_{\Gamma}$ | 0.595 | 0.625 | $0.610$ | 0.645 | 0.675 | $0.700$ | $0.795$ | $0.895$ | 0.935 | $0.940$
$\Phi_{D_{s}}$ | $0.540$ | $0.575$ | $0.590$ | $0.620$ | $0.650$ | $0.675$ | 0.805 | $0.850$ | $0.855$ | $0.925$

Table 1 collects simulation results for the given example. We have included
the sequential procedure (4) as a benchmark for orientation. For all other
approaches the true parameter values are used in the covariance kernels.
Concerning static (distance-based) designs based on maximisation of
$\Phi_{F},\Phi_{1},\Phi_{2},\Phi_{KL},\Phi_{\Gamma},\Phi_{D_{s}}$, for each
design size considered we first built a an incremental design and then used a
classical exchange-type algorithm to improve it. These designs are thus not
necessarily nested, i.e., $\mathbf{X}_{n}\not\subset\mathbf{X}_{n^{\prime}}$
for $n<n^{\prime}$.

Each design of size $n$ was then evaluated by generating $N=100$ independent
sets of $n$ observations generated with the assumed true model, evaluating the
likelihood functions for these sets of observations for both models, and then
deciding for each set of observations which model has the higher likelihood
value. The hit rate is the fraction of sets of observations where the assumed
true model has the higher likelihood value. The procedure was repeated by
assuming the other model to be the true one. The two hit rates are then
averaged and stated in Table 1, which contains the results for all the
criteria and design sizes we considered. For the special case of the
sequential construction (4), the design path depends on the observations
generated at the previously selected design points; that is, unlike for the
other criteria, for a given design size $n$ each random run produces a
different design. To compute the hit rates for a particular $n$ we used
$N=100$ independent runs of the experiment.

The hit rates reported in Table 1 reflect the discriminatory power of the
corresponding designs. One can observe that $\Phi_{F}$ and as expected
$\Phi_{KL}$ are outperforming in terms of hit rates. The Bayesian lower bound
criterion $\Phi_{\Gamma}$ is similar to the symmetric $\Phi_{KL}$. The
sequential design strategy (4) does not behave as well as the outperforming
ones. It is, however, the realistic scenario that one might consider in
applications as it does not assume knowledge of the kernel parameters. The
effect of this knowledge can thus be partially calibrated for by comparing the
first line against the other criteria.

###  6.2 Optimal design measure for $\phi_{p}$

Theorem 1 also allows the use of approximate designs as it presents a
necessary condition for optimality of the family of criteria
$\phi_{p},~{}p>0$. This is more extensively discussed in the previous section.
Here we present the numerical results for two specific cases of $p=2$ and
$p=10$. To reach a design which might be numerically optimal (or at least
nearly optimal), we have applied the Fedorov-Wynn algorithm (Fedorov,, 1971;
Wynn,, 1970) on a dense regular grid of candidate points.

Numerical results show that for very small $p$ (e.g., $p=1$) explicit optimal
measures are hard to derive. The left panel in Figure 10 presents the measure
$\xi_{2}^{*}$ obtained for $\phi_{2}$. To construct $\xi_{2}^{*}$, we have
first calculated an optimal design on a dense grid by applying $1000$
iterations of the Fedorov-Wynn algorithm (see the comment following Theorem
1); the design measure obtained is supported on 9 grid points. We then applied
a continuous optimisation algorithm (library NLopt (Johnson,, 2021) through
its R-interface nloptr) initialised at this $9$-point design. The 9 support
points of the resulting design measure $\xi_{2}^{*}$ are independent of the
grid size; they receive unequal weights, proportional to the disk areas on
Figure 10-left. Any translation or rotation of $\xi_{2}^{*}$ yields the same
value of $\phi_{2}$.

As the order $p$ increases, we eventually reach an optimal measure with only
three support points and equal weights. The right panel in Figure 10
corresponds to the optimal design measure computed for $\phi_{10}$. This has,
similarly as before, resulted from application of a continuous optimisation
initialised at an optimal 3-point design calculated with the Fedorov-Wynn
algorithm on a grid. This optimal design measure $\xi_{10}^{*}$ has three
support points, drawn as blue dots, with equal weights 1/3 represented by the
areas of the red disks. The blue line segments between every two locations
have length $\Delta\simeq 1.92$, reflecting the ideal interpoint distance (see
the right panel of Figure 8), in agreement with corresponding discussions in
Section 5. Also here the optimal designs are rotationally and translationally
invariant, and thus any design of such type is optimal as long as the design
region is large enough to fit it.

Figure 10: Left: The optimal measure for $\phi_{2}$. Right: The optimal
measure for $\phi_{10}$.

##  7 Conclusions

In this paper we have considered the design problem for the discrimination of
Gaussian process regression models. This problem differs considerably from the
well-treated one in standard regression models and thus offers a multitude of
challenges. While the KL-divergence is a straightforward criterion, it comes
with the price of being computationally demanding and lacking convenient
simplifications such as design measures. We have therefore introduced a family
of criteria that allow such a simplification at least in special cases and
have investigated its properties. We have also compared the performance of
these and other potential criteria on several examples and see that KL-
divergence can be effectively replaced by simpler criteria without much loss
in efficiency. In particular designs based on the Fréchet-distance between
covariance kernels seem to be competitive. Results from the approximate design
computations indicate that for classical isotropic kernels, designs with $d+1$
support points placed at the vertices of a simplex of suitable size are
optimal for distance-based criteria $\phi_{p}$ when $p$ is large enough.

## Acknowledgments

This work was partly supported by project INDEX (INcremental Design of
EXperiments) ANR-18-CE91-0007 of the French National Research Agency (ANR) and
I3903-N32 of the Austrian Science Fund (FWF).

## References

  * Atkinson and Fedorov, (1975) Atkinson, A. C. and Fedorov, V. V. (1975).  The design of experiments for discriminating between two rival models.  Biometrika, 62(1):57–70.
  * Box and Hill, (1967) Box, G. E. P. and Hill, W. J. (1967).  Discrimination among mechanistic models.  Technometrics, 9(1):57–71.
  * Dowson and Landau, (1982) Dowson, D. C. and Landau, B. V. (1982).  The Fréchet distance between multivariate normal distributions.  Journal of Multivariate Analysis, 12(3):450–455.
  * Fedorov, (1971) Fedorov, V. V. (1971).  The design of experiments in the multiresponse case.  Theory of Probability & Its Applications, 16(2):323–332.
  * Gramacy, (2020) Gramacy, R. B. (2020).  Surrogates : Gaussian Process Modeling, Design, and Optimization for the Applied Sciences.  Chapman and Hall/CRC.
  * Heirung et al., (2019) Heirung, T. A. N., Santos, T. L. M., and Mesbah, A. (2019).  Model predictive control with active learning for stochastic systems with structural model uncertainty: Online model discrimination.  Computers & Chemical Engineering, 128:128–140.
  * Hershey and Olsen, (2007) Hershey, J. R. and Olsen, P. A. (2007).  Approximating the Kullback Leibler divergence between Gaussian mixture models.  In 2007 IEEE International Conference on Acoustics, Speech and Signal Processing – ICASSP ’07, volume 4, pages IV–317–IV–320.
  * Hill and Hunter, (1969) Hill, W. J. and Hunter, W. G. (1969).  A Note on Designs for Model Discrimination: Variance Unknown Case.  Technometrics, 11(2):396–400.
  * Hino, (2020) Hino, H. (2020).  Active learning: Problem settings and recent developments.  Technical report, arXiv:2012.04225.
  * Hoffmann, (2017) Hoffmann, C. (2017).  Numerical Aspects of Uncertainty in the Design of Optimal Experiments for Model Discrimination.  PhD thesis, Ruprecht-Karls-Universität Heidelberg.
  * Hunter and Reiner, (1965) Hunter, W. and Reiner, A. (1965).  Designs for discriminating between two rival models.  Technometrics, 7(3):307–323.
  * Johnson, (2021) Johnson, S. G. (2021).  The NLopt nonlinear-optimization package, http://github.com/stevengj/nlopt.
  * Karvonen, (2022) Karvonen, T. (2022).  Asymptotic bounds for smoothness parameter estimates in Gaussian process interpolation.  arXiv preprint arXiv:2203.05400.
  * Karvonen and Oates, (2022) Karvonen, T. and Oates, C. (2022).  Maximum likelihood estimation in Gaussian process regression is ill-posed.  arXiv preprint arXiv:2203.09179.
  * Karvonen et al., (2020) Karvonen, T., Wynne, G., Tronarp, F., Oates, C., and Särkkä, S. (2020).  Maximum likelihood estimation and uncertainty quantification for Gaussian process approximation of deterministic functions.  SIAM/ASA Journal on Uncertainty Quantification, 8(3):926–958.
  * Kiefer, (1974) Kiefer, J. (1974).  General equivalence theory for optimum designs (approximate theory).  The Annals of Statistics, 2(5):849–879.
  * Kullback and Leibler, (1951) Kullback, S. and Leibler, R. A. (1951).  On information and sufficiency.  The Annals of Mathematical Statistics, 22(1):79–86.
  * Lee et al., (2018) Lee, X. J., Hainy, M., McKeone, J. P., Drovandi, C. C., and Pettitt, A. N. (2018).  ABC model selection for spatial extremes models applied to South Australian maximum temperature data.  Computational Statistics & Data Analysis, 128:128–144.
  * López-Fidalgo et al., (2007) López-Fidalgo, J., Tommasi, C., and Trandafir, P. C. (2007).  An optimal experimental design criterion for discriminating between non-normal models.  Journal of the Royal Statistical Society: Series B (Statistical Methodology), 69(2):231–242.
  * Motzkin and Straus, (1965) Motzkin, T. S. and Straus, E. G. (1965).  Maxima for graphs and a new proof of a theorem of Turán.  Canadian Journal of Mathematics, 17:533–540.
  * Müller, (2007) Müller, W. G. (2007).  Collecting Spatial Data: Optimum Design of Experiments for Random Fields.  Springer Berlin, Heidelberg, 3 edition.
  * Olofsson et al., (2018) Olofsson, S., Deisenroth, M. P., and Misener, R. (2018).  Design of experiments for model discrimination using Gaussian process surrogate models.  In Eden, M. R., Ierapetritou, M. G., and Towler, G. P., editors, 13th International Symposium on Process Systems Engineering (PSE 2018), volume 44 of Computer Aided Chemical Engineering, pages 847–852. Elsevier.
  * Pronzato et al., (2019) Pronzato, L., Wynn, H. P., and Zhigljavsky, A. (2019).  Bregman divergences based on optimal design criteria and simplicial measures of dispersion.  Statistical Papers, 60(2):545–564.
  * Sauer et al., (2022) Sauer, A., Gramacy, R. B., and Higdon, D. (2022).  Active learning for deep Gaussian process surrogates.  Technometrics, 0(0):1–15.
  * Schwaab et al., (2008) Schwaab, M., Luiz Monteiro, J., and Carlos Pinto, J. (2008).  Sequential experimental design for model discrimination: Taking into account the posterior covariance matrix of differences between model predictions.  Chemical Engineering Science, 63(9):2408–2419.
  * Stein, (1999) Stein, M. (1999).  Interpolation of Spatial Data: Some Theory for Kriging.  Springer Series in Statistics. Springer, Heidelberg.
  * Wynn, (1970) Wynn, H. P. (1970).  The sequential generation of $D$-optimum experimental designs.  The Annals of Mathematical Statistics, 41(5):1655–1664.

## Appendix

##  Appendix A Notes on Box-Hill-Hunter Bayesian criteria for model
discrimination between Gaussian random fields

Chapter 5 of Hoffmann, (2017) contains an overview of Bayesian design criteria
for model discrimination and some useful bounds on them. We assume there are
$M$ models $m_{0},\ldots,m_{M-1}$. The most common Bayesian design criterion
for model discrimination has the following form:

$\Phi_{\Lambda}(\mathbf{X}_{k})=-\sum_{i=0}^{M-1}p(m_{i})\log(p(m_{i}))+\int_{\mathbf{Y}_{k}\in\mathcal{Y}}p(\mathbf{Y}_{k})\sum_{i=0}^{M-1}p(m_{i}|\mathbf{Y}_{k})\log(p(m_{i}|\mathbf{Y}_{k}))\>\mathrm{d}\mathbf{Y}_{k},$
(22)

where the data $\mathbf{Y}_{k}=(Y_{1}(x_{1}),\ldots,Y_{k}(x_{k}))^{\top}$ are
observed at the design $\mathbf{X}_{k}=(x_{1},\ldots,x_{k})$, $p(m_{i})$
denotes the prior and $p(m_{i}|\mathbf{Y}_{k})$ the posterior model
probability of model $m_{i}$ and $p(\mathbf{Y}_{k})$ is the marginal
distribution of $\mathbf{Y}_{k}$ with respect to the models. Hence, this
criterion is the (expected) difference of the model entropy and the
conditional model entropy (conditional on the observations). The posterior
model probability $p(m_{i}|\mathbf{Y}_{k})$ is defined by

$p(m_{i}|\mathbf{Y}_{k})\propto p(\mathbf{Y}_{k}|m_{i})p(m_{i}),$

where $p(\mathbf{Y}_{k}|m_{i})$ is the likelihood of model $m_{i}$
(marginalised over the parameters), and $p(\mathbf{Y}_{k})$ is given by

$p(\mathbf{Y}_{k})=\sum_{i=0}^{M-1}p(\mathbf{Y}_{k}|m_{i})p(m_{i}).$

The first term in (22) does not depend on the design and can therefore be
ignored.

A common alternative formulation of criterion (22) is the one adopted by Box
and Hill, (1967) and Hill and Hunter, (1969), which will henceforth be called
Box-Hill-Hunter (BHH) criterion:

$\Phi_{\Lambda}(\mathbf{X}_{k})=\sum_{i=0}^{M-1}p(m_{i})\int_{\mathbf{Y}_{k}\in\mathcal{Y}}p(\mathbf{Y}_{k}|m_{i})\log\left(\frac{p(\mathbf{Y}_{k}|m_{i})}{p(\mathbf{Y}_{k})}\right)\>\mathrm{d}\mathbf{Y}_{k}.$
(23)

In our case, if we assume point priors for the kernel parameters, we have

$p(\mathbf{Y}_{k}|m_{i})=\varphi(\mathbf{Y}_{k}|\boldsymbol{\eta}_{k,i},\mathbf{K}_{k,i}),$

where
$\boldsymbol{\eta}_{k,i}=(\eta_{1,i}(x_{1}),\ldots,\eta_{k,i}(x_{k}))^{\top}$
is the mean vector of model $i$ at design $\mathbf{X}_{k}$, $\mathbf{K}_{k,i}$
is the $k\times k$ kernel matrix of model $i$ with elements given by
$\\{\mathbf{K}_{k,i}\\}_{j,l}=K_{i}(x_{j},x_{l})$, and
$\varphi(\cdot|\boldsymbol{\eta},\mathbf{K})$ is the normal pdf with mean
vector $\boldsymbol{\eta}$ and variance-covariance matrix $\mathbf{K}$.

For example, for a static design involving $n$ design points, we set $k=n$ and
assume that $\boldsymbol{\eta}_{n,i}=\mathbf{0}$ for each design
$\mathbf{X}_{n}$. The model probabilities $p(m_{i})$ would just be the prior
model probabilities before having collected any observations.

In a sequential design setting, where $n$ observations $\mathbf{Y}_{n}$ have
already been observed at locations $\mathbf{X}_{n}$ and we want to find the
optimal design point $x$ where to collect our next observation, we have $k=1$
and set $\boldsymbol{\eta}_{k,i}$ to the conditional mean
$\widehat{\eta}_{n,i}(x)=\mathbf{k}_{n,i}(x)^{\top}\mathbf{K}_{n,i}^{-1}\mathbf{Y}_{n}$
and $\mathbf{K}_{k,i}$ to the conditional variance
$\rho_{n,i}^{2}(x)=K_{i}(x,x)-\mathbf{k}_{n,i}(x)^{\top}\mathbf{K}_{n,i}^{-1}\,\mathbf{k}_{n,i}(x)$,
where $\mathbf{k}_{n,i}(x)^{\top}=(K_{i}(x,x_{1}),\ldots,K_{i}(x,x_{n}))$, see
Section 3.1. The prior model probabilities would have to be set to the
posterior model probabilities given the already observed data:

$p(m_{i})=p(m_{i}|\mathbf{Y}_{n})\propto\varphi(\mathbf{Y}_{n}|\,\mathbf{0},\mathbf{K}_{n,i})\>p(m_{i}).$

It follows that $p(\mathbf{Y}_{k})$ is a mixture of normal distributions. The
criterion representations (22) and (23) cannot be computed directly. However,
several bounds have been developed for the criterion, the most famous being
the classic upper bound derived by Box and Hill, (1967).

###  A.1 Upper bound

The upper bound has the following form (see also Hoffmann, (2017, Thm. 5.2, p.
168)):

$\Phi_{U}(\mathbf{X}_{k})=\frac{1}{2}\sum_{i=0}^{M-1}\sum_{j=0}^{M-1}p(m_{i})p(m_{j})\left\\{\left\|\boldsymbol{\eta}_{k,i}-\boldsymbol{\eta}_{k,j}\right\|^{2}_{\mathbf{K}_{k,j}^{-1}}+\mathrm{trace}\left(\mathbf{K}_{k,i}\mathbf{K}_{k,j}^{-1}\right)-n\right\\}.$

For $M=2$, the formula simplifies to

$\displaystyle\Phi_{U}(\mathbf{X}_{k})$ $\displaystyle=$
$\displaystyle\frac{1}{2}p(m_{0})p(m_{1})\biggl{\\{}\left\|\boldsymbol{\eta}_{k,0}-\boldsymbol{\eta}_{k,1}\right\|^{2}_{\mathbf{K}_{k,0}^{-1}}+\left\|\boldsymbol{\eta}_{k,0}-\boldsymbol{\eta}_{k,1}\right\|^{2}_{\mathbf{K}_{k,1}^{-1}}$
$\displaystyle\qquad\qquad\qquad+\mathrm{trace}\left(\mathbf{K}_{k,0}\mathbf{K}_{k,1}^{-1}\right)+\mathrm{trace}\left(\mathbf{K}_{k,1}\mathbf{K}_{k,0}^{-1}\right)-2n\biggr{\\}}.$

This is equivalent to the symmetric Kullback-Leibler divergence that we use as
the criterion $\Phi_{KL}$ (with $p(m_{0})=p(m_{1})=1/2$ and
$\boldsymbol{\eta}_{k,0}=\boldsymbol{\eta}_{k,1}=\mathbf{0}$) .

###  A.2 Lower bound

Hershey and Olsen, (2007, Sec. 7) derive a lower bound for the Kullback-
Leibler divergence between a mixture of two normals, see also Hoffmann, (2017,
Thm. 5.4 and Cor. 5.5, pp. 173–174). This result is then used by Hoffmann,
(2017) to find a lower bound for the BHH criterion
$\Phi_{\Lambda}(\mathbf{X}_{k})$ (Hoffmann,, 2017, Thm. 5.9, p. 178). This
lower bound is given by

$\Phi_{\Gamma}(\mathbf{X}_{k})=-\sum_{i=0}^{M-1}p(m_{i})\log\left\\{\sum_{j=0}^{M-1}p(m_{j})\exp\left(-\frac{1}{2}\boldsymbol{\Gamma}(\mathbf{X}_{k})_{ij}\right)\right\\},$

where

$\boldsymbol{\Gamma}(\mathbf{X}_{k})_{ij}=\left\|\boldsymbol{\eta}_{k,i}-\boldsymbol{\eta}_{k,j}\right\|^{2}_{\mathbf{K}_{k,j}^{-1}}+\mathrm{trace}\left(\mathbf{K}_{k,i}\mathbf{K}_{k,j}^{-1}\right)-\log\det\left(\mathbf{K}_{k,i}\mathbf{K}_{k,j}^{-1}\right)-n.$

For $M=2$, as is the relevant case for our setup we get

$\displaystyle\Phi_{\Gamma}(\mathbf{X}_{k})$ $\displaystyle=$
$\displaystyle-p(m_{0})\log\biggl{\\{}p(m_{0})$ (24)
$\displaystyle\qquad\qquad\quad+p(m_{1})\exp\biggl{(}-\frac{1}{2}\bigl{[}\left\|\boldsymbol{\eta}_{k,0}-\boldsymbol{\eta}_{k,1}\right\|^{2}_{\mathbf{K}_{k,1}^{-1}}+\mathrm{trace}\left(\mathbf{K}_{k,0}\mathbf{K}_{k,1}^{-1}\right)$
$\displaystyle\qquad\qquad\qquad\qquad\qquad\qquad-\log\det\left(\mathbf{K}_{k,0}\mathbf{K}_{k,1}^{-1}\right)-n\bigr{]}\biggr{)}\biggr{\\}}$
$\displaystyle-p(m_{1})\log\biggl{\\{}p(m_{1})$
$\displaystyle\qquad\qquad\quad+p(m_{0})\exp\biggl{(}-\frac{1}{2}\bigl{[}\left\|\boldsymbol{\eta}_{k,0}-\boldsymbol{\eta}_{k,1}\right\|^{2}_{\mathbf{K}_{k,0}^{-1}}+\mathrm{trace}\left(\mathbf{K}_{k,1}\mathbf{K}_{k,0}^{-1}\right)$
$\displaystyle\qquad\qquad\qquad\qquad\qquad\qquad-\log\det\left(\mathbf{K}_{k,1}\mathbf{K}_{k,0}^{-1}\right)-n\bigr{]}\biggr{)}\biggr{\\}}$

where
$\varphi_{i}(\cdot)=\varphi(\cdot|\boldsymbol{\eta}_{k,i},\mathbf{K}_{k,i})$,
which we are also using to compute designs in Section 6.1 (again with
$p(m_{0})=p(m_{1})=1/2$ and
$\boldsymbol{\eta}_{k,0}=\boldsymbol{\eta}_{k,1}=\mathbf{0}$).