# Uncovering Regions of Maximum Dissimilarity
on Random Process Data

Miguel de Carvalho
School of Mathematics
University of Edinburgh
EH93FD, Edinburgh, UK
<EMAIL_ADDRESS>
&Gabriel Martos
Departamento de Matemática y Estadística
Universidad Torcuato Di Tella
Buenos Aires, Argentina
<EMAIL_ADDRESS>

###### Abstract

The comparison of local characteristics of two random processes can shed light
on periods of time or space at which the processes differ the most. This paper
proposes a method that learns about regions with a certain volume, where the
marginal attributes of two processes are less similar. The proposed methods
are devised in full generality for the setting where the data of interest are
themselves stochastic processes, and thus the proposed method can be used for
pointing out the regions of maximum dissimilarity with a certain volume, in
the contexts of functional data, time series, and point processes. The
parameter functions underlying both stochastic processes of interest are
modeled via a basis representation, and Bayesian inference is conducted via an
integrated nested Laplace approximation. The numerical studies validate the
proposed methods, and we showcase their application with case studies on
criminology, finance, and medicine.

_K_ eywords Functional Parameters, Multi-objective Optimization, Pairs of
Random Processes, Kolmogorov metric, Set Function Optimization, Youden J
Statistic

##  1 Introduction

Everyday millions of data patterns flow around the world at unprecedented
speed, thus leading to an explosion on the demand for modeling stochastic
process data—such as time series, point processes, and functional data; each
of these types of data plays a key role in machine learning, as can be seen,
for instance, from the recent papers of Berrendero et al. (2020), Faouzi and
Janati (2020), and Xu et al. (2020). Hand in hand with this shock on demand
arrived a pressing need for the development of data-intensive methods,
techniques, and algorithms for learning and comparing random processes.

###  1.1 The Learning Problem of Interest

In this paper we deal with the following problem on the comparison of
stochastic processes:

> learning problem. For a pair of random processes, what is the region—with a
> given volume—where they statistically differ the most?

Specifically, the task of interest will entail tracking down regions with a
certain volume, where the marginal attributes of two stochastic processes
differ the most. Throughout, we will refer to this problem as that of
uncovering regions of maximum discrimination (RMD) on random process data.

Let’s translate the learning problem into mathematical parlance. Since the
target of interest consists of a set that fulfills an optimization
criterion—i.e. a period of time or region over which a marginal feature of two
stochastic processes differ the most—some of the main concepts in this paper
can be framed as an optimization problem over a set function. The canonical
problem in set function optimization is of the type,

$\begin{array}[]{rl}\underset{A\subseteq\mathscr{A}}{\max}&F(A)\\\
\text{s.t.}&A\in\mathscr{F},\end{array}$ (1)

where $F:2^{\mathscr{A}}\rightarrow\mathbf{R}$ is a set function,
$\mathscr{A}$ is a set, and $\mathscr{F}\subseteq 2^{\mathscr{A}}$ is the
collection of feasible sets defining the constraint. Optimization problems
over set functions—such as (1)—are commonplace in machine learning (e.g.
Krause, 2010). For a recent review on the theory of discrete set function
optimization see Wu et al. (2019). Most state of the art developments have
been made on discrete or combinatorial set function optimization, especially
on the class of submodular functions (e.g. Goldengorin, 2009). Our paper
provides one of the first steps towards continuous set function optimization
as here the aim will be to solve (1) when $\mathscr{F}$ is a family of Borel
subsets of a compact set $T\subset\mathbf{R}^{d}$.

As it will be seen in Section 2, the objective set function of interest in our
case will be a measure of proximity between marginal features of the pair of
stochastic processes of interest, whereas the collection of feasible sets
introduces the constraint on the ‘size’ of the feasible regions—i.e. periods
of time or space—over which the comparison is made.

###  1.2 Our Contributions

Our main contributions are:

  1. 1.

We pioneer the study, formulation, and analysis of the learning problem of
tracking down regions of maximum discrimination as described in Section
1.1—and formally defined in Sections 2–3.

  2. 2.

In its most standard version, the proposed learning problem is shown to be
equivalent to a continuous set function optimization on a monotone modular
function, under a Lebesgue measure constraint. Hence, as a byproduct, the
paper contributes to the literature on set function optimization which is
mostly focused on a discrete and combinatorial framework, with a particular
emphasis on monotone submodular functions (e.g. Nemhauser et al., 1978;
Calinescu et al., 2011; Goldengorin, 2009; Buchbinder et al., 2017; Buchbinder
and Feldman, 2018), under cardinality or matroid constraints. Little is known
on continuous set function optimization, and thus the tools, concepts, and
strategies devised herein can be of further interest elsewhere.

  3. 3.

Our approach is fully general in the sense that it applies to most random
process data (including functional data, time series, and point processes).
The functional parameters of the processes of interest (say mean functions,
volatility functions, or intensity functions) are modelled by composing an
inverse link function with a basis function representation, and Bayesian
modeling is conducted via latent Gaussian models and inference is conducted
via INLA (Integrated Nested Laplace Approximations) (Rue et al., 2009, 2017).

  4. 4.

An extension of the proposed method applies also to the context where the
interest is on comparing more than one marginal feature via a multi-objective
version of the set function optimization problem of interest.

  5. 5.

A variant of the proposed approach shows that a well-known performance measure
known as Youden index and the Youden’s J statistic (e.g. Sokolova et al.,
2006; Inácio de Carvalho et al., 2017) as well as the Kolmogorov metric (e.g.,
Gretton et al., 2012) can be regarded as particular cases of the framework
developed herein.

###  1.3 Some Related Prior Work

While the learning problem introduced in Section 1.1 is new—and while there
are novel contributions that will arise from our solution to it—there are some
recent approaches in the context of functional data analysis (Ramsay and
Silverman, 2002, 2006; Ferraty and Vieu, 2006; Horváth and Kokoszka, 2012)
that are tangentially related to it, and that are briefly reviewed below.

Pini and Vantini (2016, 2017) propose an interval testing procedure for
functional data that points out specific differences between functional
populations. Also in the context of functional data, Berrendero et al. (2016)
propose a discretization method consisting on learning to choose a finite
collection of points in the domain of a set of functions in order to improve
the performance of a functional data classifier. Martos and de Carvalho (2018)
propose a Mann–Whitney type of statistic for functional data so to learn about
the regions at which two processes differ the most on aspects related with
symmetry. Finally, Dette and Kokot (2021) develop hypothesis tests for the
equivalence of functional parameters in a two sample functional data setup.

Our approach differs from the ones mentioned above in a number of important
ways. Perhaps the most important ones are that: i) here the goal is not to
test hypothesis but rather to learn about regions with a given volume where
two processes differ the most; ii) our approach applies to random processes in
general whereas the methodologies reviewed above have mainly been designed
with the context of functional data analysis in mind.

###  1.4 Outline of the Paper

The rest of the paper unfolds as follows. In Section 2 we introduce sets of
maximum dissimilarity, present examples, and introduce the inference methods.
In Section 3 we comment on extensions of the main concepts, methods, and ideas
covered in Section 2. In Section 4 we report the main findings of a Monte
Carlo numerical study on artificial data. In Section 5 we showcase the
proposed methods on real data applications. Closing remarks are given in
Section 6. Appendix A includes a selection of auxiliary facts, and the proofs
of main results are included in Appendix B. Table 1 lists symbols and notation
used throughout the article. We write $X(t)$ instead of $X_{t}$ when
typographically convenient.

Symbol | Description
---|---
$X$, $Y$ | stochastic processes under comparison, that is, $X=\\{X_{t}\\}$ and $Y=\\{Y_{t}\\}$
$T$ | ground, or index, set over which processes $X$ and $Y$ are defined
$\theta_{X}$, $\theta_{Y}$ | functional parameters for $X$ and $Y$
$\mathcal{D}_{X},\mathcal{D}_{Y}$ | data on proceses $X$ and $Y$
$\|\cdot\|^{(A)}_{p}$ |  $L^{p}$ sub-norm over $A$
$F(\cdot)$ | set objective function $\|\theta_{X}-\theta_{Y}\|^{(\cdot)}_{p}$
$f(t,r)$ | set objective function evaluated at closed ball $B(t,r)$, that is, $F\\{B(t,r)\\}$
$F_{i}$ | set objective functions in multi-objective context
$\mathcal{A}$ | set of compact and convex subsets of the ground set $T$
$|A|$ | volume functional, i.e., Lebesgue measure of $A$ in $\mathbf{R}^{d}$
$\mathscr{F}_{c}$ | collection of feasible sets
$A^{*}_{c}$ | region of maximum or multi-maximum dissimilarity
$D^{*}_{c}$ | dissimilarity index for region of maximum dissimilarity
$B^{*}_{c}$ | ball of maximum or multi-maximum dissimilarity
$\mathscr{D}^{*}_{c}$ | dissimilarity index for ball of maximum dissimilarity
$t_{c}^{*},r^{*}_{c}$ | center and radius of ball of maximum dissimilarity
$\mathscr{B}_{p}$ | set of all closed balls in $L^{p}$
$d_{\mathrm{H}}(A,B)$ | Hausdorf distance between sets $A$ and $B$
$d(x,B)$ | minimum Euclidean distance between a point $x\in A$ and set $B$
$g|_{A}$ | restriction of function $g$ to set $A$
$\partial A$ | boundary of set $A$
$\mathchoice{{\vbox{\hbox{$\textstyle-$ }}\kern-7.83337pt}}{{\vbox{\hbox{$\scriptstyle-$ }}\kern-5.90005pt}}{{\vbox{\hbox{$\scriptscriptstyle-$ }}\kern-4.75003pt}}{{\vbox{\hbox{$\scriptscriptstyle-$ }}\kern-4.25003pt}}\\!\int$ | average integral symbol
$\twoheadrightarrow$ | notation for defining correspondences (i.e. set-valued functions)
Table 1:  Symbols and notation used throughout the article

##  2 Learning about Sets of Maximum Dissimilarity

###  2.1 Sets of Maximum Dissimilarity

#### A taster—some simple instances

Prior to introducing sets of maximum dissimilarity in a formal fashion, we
introduce some simple instances of the concept based on different parameter
functions $\theta_{X}$ and $\theta_{Y}$. The respective sets of maximum
dissimilarity and parameter functions for Examples 1–3 to be presented below
are depicted in Fig. 1.

|  |
---|---|---

(a)

(b)

(c)

Figure 1: Sets of maximum dissimilarity for Examples 1–3. The sets are
depicted using brackets in the time axis in the case of mean functions and
volatility functions (a, c) and on the top of the box with a line in the case
of intensity functions (b). Raw data alongside the true and estimated
functional parameters $\theta_{X}$ and $\theta_{Y}$ are depicted in red and
blue respectively; the proposed inferences for the sets will be discussed in
Section 2.2.

######  Example 1 (Mean Functions)

Let $X_{t}=\theta_{X}(t)+\varepsilon_{X}(t)$ and
$Y_{t}=\theta_{Y}(t)+\varepsilon_{Y}(t)$ be two stochastic processes defined
on the unit interval, $T=[0,1]$, with

$\begin{cases}\theta_{X}(t)=E(X_{t})=b(t)+4\cos(10t)-2(t-0.75)^{2},\\\
\theta_{Y}(t)=E(Y_{t})=b(t)+3\sin(12t),\end{cases}$ (2)

where $b(t)=1/2\exp\\{10(t-0.5)^{2}\\}$ is a baseline curve, and
$\varepsilon_{X}(t)$ and $\varepsilon_{Y}(t)$ are zero mean Gaussian error
functions. Fig. 1(a) shows the connected set of size $c=0.2$ where both mean
functions differ the most in the $L^{1}$ sense.

######  Example 2 (Intensity Functions)

Consider the intensity functions associated to two point processes

$\begin{cases}\theta_{X}(t)=\lambda_{X}(t)=\gamma\exp\\{-(t_{1}^{2}+t_{2}^{2})/2\\},\\\
\theta_{Y}(t)=\lambda_{Y}(t)=\delta\theta_{X}(t),\end{cases}$ (3)

defined over the region $T=[-3,3]^{2}$ with $\gamma,\delta>0$. Fig. 1(b) shows
the true ball of maximum dissimilarity with area $4\pi$ at which the two
intensity functions differ the most in the $L^{2}$ sense in the case where
$\gamma=100$ and $\delta=1/2$.

######  Example 3 (Volatility Functions)

Suppose $X_{t}$ and $Y_{t}$ are the log returns of two stock markets, with
$E(X_{t})=E(Y_{t})=0$, and that the goal is to search for the period of about
a quarter (think of $T=[0,1]$, so that $c=0.25$), where the volatility between
both markets differed the most. Let

$\begin{cases}\theta_{X}(t)=\sigma_{X}(t):=\\{\operatorname{var}(X_{t})\\}^{1/2},\\\
\theta_{Y}(t)=\sigma_{Y}(t):=\\{\operatorname{var}(Y_{t})\\}^{1/2},\\\
\end{cases}$

Fig. 1(c) depicts a simulated example of two artificial stock prices evolving
during a period of 1 year ($T=[0,1]$), along with the 3 month period over
which the volatilities in both markets differed the most in the $L^{2}$ sense.

To allow for visualizations, in the examples above we focused on low-
dimensional instances but the theory to be presented next holds in general for
any compact ground set $T\subset\mathbf{R}^{d}$.

#### Preparations

We start by laying the groundwork and recalling background. Let
$\theta_{X}\equiv\theta_{X}(t)$ and $\theta_{Y}\equiv\theta_{Y}(t)$ be
functional parameters characterizing $q$ different marginal features of the
processes $X$ and $Y$ for $t\in T\subset\mathbf{R}^{d}$. Throughout, we will
assume that the ground set $T$ is compact and that $\theta_{X}$ and
$\theta_{Y}$ live in the Banach space $(L^{p}(T),\|\cdot\|_{p})$, with
$\|f\|_{p}=(\int_{T}|f|^{p}\,\mathrm{d}\mu)^{1/p}$, where $\mu$ is a measure
over the Borel sets on $T$; we will refer to
$\|f\|_{p}^{(A)}=(\int_{A}|f|^{p}\,\mathrm{d}\mu)^{1/p}$ as the $L^{p}$
sub–norm over $A\subseteq T$.

The learning problem from Section 1.1 boils down to searching for a set
$A\subseteq T$, with volume not greater than $c\geq 0$, over which the
difference between $\theta_{X}$ and $\theta_{Y}$ is highest. As it will be
shown below, a solution to such optimization problem can be ensured to exist
provided that we impose additional structure on the search domain, namely
convexity and compactness. We equip the set $\mathcal{A}$ of compact and
convex subsets of the ground set $T$, with the Hausdorff distance

$d_{\text{H}}(A,B)=\max\left\\{\max_{t\in A}\,d(t,B),\max_{t\in
B}\,d(t,A)\right\\},\qquad A,B\in\mathcal{A}.$

Here, $d(t,B)$ is the minimum Euclidean distance from the point $t\in A$ to
the set $B$. Let $A\subseteq T$ be such that $|A|\leq c$, where $c\geq 0$
defines the maximum size of the regions of interest, and with $|A|$ denoting
the Lebesgue measure of $A$.

Figure 2:  Geometry underlying regions of maximum discrimination. The set
function optimization problem in (4) searches for the region $A$ whose $L^{p}$
sub-norm of $\theta_{Z}=\theta_{X}-\theta_{Y}$ is the largest.

#### Regions of maximum dissimilarity

A compact and convex set $A\subseteq T$ is said to be a region of maximum
dissimilarity, with volume bounded by $c$, if it maximizes
$F(A)=\|\theta_{X}-\theta_{Y}\|^{(A)}_{p}$ and if its volume does not exceed
$c$; the formal definition is as follows.

######  Definition 1 ($L^{p}$ Region of Maximum Dissimilarity)

Suppose $\theta_{X},\theta_{Y}\in L^{p}(T)$ and $T\subset\mathbf{R}^{d}$ is
compact. A region of maximum dissimilarity (RMD) is defined as a set
$A^{*}_{c}\subseteq T$ that solves,

$\begin{array}[]{rl}\underset{A\subseteq
T}{\max}&\|\theta_{X}-\theta_{Y}\|^{(A)}_{p}\\\ \textit{s.t.}&|A|\leq c,\\\
&A~{}\text{is compact and convex},\end{array}$ (4)

for a fixed $c\geq 0$. In addition,
$D_{c}^{*}(\theta_{X},\theta_{Y})=\|\theta_{X}-\theta_{Y}\|^{(A^{*}_{c})}$ is
said to be the dissimilarity index.

Some comments on Definition 1 are in order. Fig. 2 depicts the geometry
underlying the set function optimization problem in (4). The optimization
problem in (4) is a continuous set function optimization problem similar to
(1), with

$\begin{cases}F(A)=\|\theta_{X}-\theta_{Y}\|^{(A)}_{p},\\\ \mathscr{A}=T,\\\
\mathscr{F}_{c}=\\{A:|A|\leq c,A~{}\text{is compact and
convex}\\}.\end{cases}$

Clearly, $F(A)$ is an increasing set function, that is if $A\subseteq B$ then
$F(A)\leq F(B)$. Also, $F(A)$ is modular (or additive) in the sense that
$F(A)+F(B)=F(A\cup B)+F(A\cap B)$. While it is evident that $A^{*}_{c}$ exists
for a few specific examples (say, $A^{*}_{|T|}=T$ almost everywhere, if $T$ is
convex), the existence of $A^{*}_{c}$ in general is not straightforward and it
is proved below in Theorem 1. RMDs and dissimilarity sets also have a number
of attributes which we summarize over Theorem 1.

###### Theorem 1

Suppose $\theta_{X},\theta_{Y}\in L^{p}(T)$ and $T\subset\mathbf{R}^{d}$ is
compact. The quantities $A_{c}^{*}$ and $D_{c}^{*}$ obey the following
properties:

  1. 1.

$A_{c}^{*}$ exists, for every $c\geq 0$.

  2. 2.

Suppose $\mu(A)=|A|$. Then, the RMDs of $(K\theta_{X},K\theta_{Y})$ and
$(L\theta_{X},L\theta_{Y})$ are respectively $A^{K}_{c}=A^{*}_{c}$ and
$A^{L}_{c}=A^{*}_{\alpha+\beta c}$, where $K\theta(t)=\alpha+\beta\theta(t)$,
$L\theta(t)=\theta(\alpha+\beta t)$, $\alpha\in\mathbf{R}$, and $\beta\neq 0$.

  3. 3.

$D_{c}^{*}(\theta_{X},\theta_{Y})$ is a distance over
$\Theta_{c}=\\{\theta|_{A}\colon\theta\in L^{p}(T),A\in\mathscr{F}_{c}\\}$,
for $c>0$.

  4. 4.

$D_{c}^{*}(\theta_{X},\theta_{Y})$ is non-decreasing as a function of $c$, for
fixed $\theta_{X},\theta_{Y}\in L^{p}(T)$.

Theorem 1 warrants some comments. The existence of RMDs (Claim 1) follows from
the continuity of the volume functional and the Blaschke selection theorem
(Appendix A); an RMD needs not however to be unique, as can be easily seen by
considering the limiting case $\theta_{X}=\theta_{Y}$ for which every compact
and convex subset of $T$ with measure $c$ is an RMD. Claim 2 shows how RMDs
are impacted by a group of linear transformations acting either on the
functional parameter or over time; the assumptions $\mu(A)=|A|$ and $\beta\neq
0$ are only required for $A^{L}_{c}=A^{*}_{\alpha+\beta c}$, as indeed
$A_{c}^{K}=A_{c}^{*}$ holds more generally. It is proved in Claim 3 that
$D_{c}^{*}(\theta_{X},\theta_{Y})$ defines a distance on the space of
functional parameters in $L^{p}(T)$ restricted to regions in
$\mathscr{F}_{c}$, which implies that dissimilarity indices have a metric
interpretation. Finally, Claim 4 notes that
$D_{c}^{*}=D_{c}^{*}(\theta_{X},\theta_{Y})$ cannot decrease as $c$ increases.

#### Balls of maximum dissimilarity

Let $\mathscr{B}_{p}$ be the family of all closed balls in $L^{p}(T)$, for
$p\geq 1$, that is,

$\mathscr{B}_{p}=\\{B(t,r):r>0,t\in T\\},$

where $B(t,r)=\\{s:\|s-t\|_{p}\leq r\\}$; to ease notation we drop the
dependence of $B(t,r)$ on $p$. A more parsimonious option for modeling is to
put more structure on the shape of RMDs, and this leads us to the following
definition.

######  Definition 2 ($L^{p}$ Ball of Maximum Dissimilarity)

Suppose $\theta_{X},\theta_{Y}\in L^{p}(T)$ and $T\subset\mathbf{R}^{d}$ is
compact. An $L^{p}$ ball of maximum dissimilarity (BMD) is defined as a set
$B_{c}^{*}\subseteq T$ that solves,

$\begin{array}[]{rl}\underset{B\subseteq
T}{\max}&\|\theta_{X}-\theta_{Y}\|^{(B)}_{p}\\\ \textit{s.t.}&|B|\leq c\\\
&B\in\mathscr{B}_{p},\end{array}$ (5)

for a fixed $c\geq 0$.

As in Definition 5 we refer to
$\mathscr{D}^{*}_{c}=\mathscr{D}^{*}_{c}(\theta_{X},\theta_{Y})=F(B^{*}_{c})$
as the dissimilarity index for BMDs. Since the volume of an $L^{p}(T)$ ball of
radius $r$ in $\mathbf{R}^{d}$ is

$|{B(t,r)}|=\frac{\\{2r\Gamma(1/p+1)\\}^{d}}{\Gamma(d/p+1)},$

for all $t\in T$, where
$\Gamma(z)=\int^{\infty}x^{z-1}\text{e}^{-x}\,\mathrm{d}x$ is the Gamma
function, it follows that the volume constraint on BMDs, $|B|\leq c$, can be
rewritten as a function of the radius, that is

$r\leq\frac{\\{c\Gamma(d/p+1)\\}^{1/d}}{2\Gamma(1/p+1)}\equiv R_{c}.$ (6)

Given the similarities between the definitions of RMDs and BMDs it is not
surprising that have identical properties (say, for $c>0$,
$\mathscr{D}^{*}_{c}(\theta_{X},\theta_{Y})$ is a also distance and
$\mathscr{D}_{c}^{*}$ is non-decreasing). The line of attack for proving
existence of BMDs is analogous to the one from Theorem 1 for RMDs but at this
time Tikhonov’s theorem (Appendix A) implies compactness of the search domain.
And interestingly the smoothness of $R_{c}$ is inherited by
$\mathscr{D}_{c}^{*}$.

###### Theorem 2

Suppose $\theta_{X},\theta_{Y}\in L^{p}(T)$ and $T\subset\mathbf{R}^{d}$ is
compact. The quantities $B_{c}^{*}$ and $\mathscr{D}_{c}^{*}$ obey the
following properties:

  1. 1.

$B^{*}_{c}$ exists, for every $c\geq 0$.

  2. 2.

$\mathscr{D}_{c}^{*}$ is continuous, and the “argmax” correspondence of
center–radius, $\alpha_{c}:[0,\infty)\twoheadrightarrow T\times[0,R_{c}]$,
defined as $\alpha_{c}=\\{(t,r)\in
T\times[0,R_{c}]:f(t,r)=\mathscr{D}_{c}^{*}\\}$, is upper hemicontinuous for
every $c\geq 0$, where $f(t,r)=F\\{B(t,r)\\}$.

#### A modicum on computing and numerical optimization

An important consequence of (6) is that the set function constrained
optimization problem (5) that leads to BMDs, can actually be written as a
standard continuous optimization problem over
$T\times[0,R_{c}]\subset\mathbf{R}^{d+1}$. Indeed, it follows from (6) that
(5) is equivalent to computing

${\max}\\{f(t,r):(t,r)\in T\times[0,R_{c}]\\}.\\\ $ (7)

where $f(t,r)=F\\{B(t,r)\\}$. That is, the BMD is ${B(t_{c}^{*},r^{*}_{c})}$
with $(t_{c}^{*},r^{*}_{c})$ maximizing (7), and hence in practice BMDs can be
computed via derivative free optimization algorithms such as conjugate search,
implicit filtering, pattern search, or Nelder–Mead (Nocedal and Wright, 2006,
Ch. 9). When the center of the optimal BMD ($t^{*}_{c}$) is ‘sufficiently far’
from the boundary of $T$, then the optimal radius ($r^{*}_{c}$) is $R_{c}$ in
(6). That is, when $d(t^{*}_{c},\partial T)>R_{c}$ the optimal radius is
$r^{*}_{c}=R_{c}$ (as $f(t,r)$ is a non-decreasing function of $r$) in which
case the optimization problem in (7) resumes to searching for
$t^{*}_{c}\in\mathbf{R}^{d}$. This also implies that often in practice the
marginal posterior of $r^{*}_{c}$ is essentially degenerated.

###  2.2 Learning from Data

#### Latent Gaussian model specification

To model BMDs in applications we consider a version of the latent Gaussian
model specification in Rue et al. (2009) adapted to our setup; to ease
notation, below we only refer to $X(t)$, and denote its functional parameter
by $\theta(t)\equiv\theta_{X}(t)$, but all comments apply to $Y(t)$ as well. A
latent Gaussian model is essentially a Bayesian generalized additive model
that assigns Gaussian priors to parameters and a possibly non-Gaussian prior
to its hyperparameters. Specifically, suppose that $Z(t)=h\\{X(t)\\}$ is in
the exponential family, with its mean function coinciding with the functional
parameter, and that

$\theta(t)=g\bigg{(}\beta_{0}+\sum_{i=1}^{B}\beta_{i}\phi_{i}(t)\bigg{)}.$ (8)

Here $\\{\phi_{i}\equiv\phi_{i}(t)\\}_{i=1}^{B}$ is a set of basis functions
in $L^{p}(T)$, $\beta=(\beta_{0},\dots,\beta_{B})^{\mathrm{\scriptscriptstyle
T}}$ is a parameter, and $g$ is an inverse link function. Following Rue et al.
we assign a multivariate Normal prior with a sparse precision matrix ($Q$) to
$\beta$, which induces a Gaussian process prior on $g^{-1}(\theta(t))$ with a
conditional independence property. Many functional parameters can be modeled
in this way including those from Examples 1–3 and all numerical instances from
Sections 4–5. The theoretical developments from Section 2.1 apply however more
generally beyond the modeling assumptions made over this section.

#### Inla-based inference for balls of maximum dissimilarity

We now discuss how to conduct inference for balls of maximum dissimilarity. It
is well-known that the latent Gaussian model described above can be fitted
with an Integrated Nested Laplace Approximation (INLA) (Rue et al., 2009); the
method is effective even when the dimension of the precision matrix $Q$ is
large, and is particularly tailored for the case where the number of
hyperparameters, $\alpha$, is moderate (say, 6–12). INLA is a deterministic
method for approximating the marginal posterior of each parameter that is
based on the Laplace approximation; loosely speaking, the Laplace
approximation is an approximation for integrals of the type $\int
e^{-nf(y)}\,\mathrm{d}y$ for large $n$, that approximates the integrand
($e^{-nf(y)}$) with a Gaussian density centered at its mode and sets the
covariance matrix as the inverse of the curvature (around the mode); see, for
instance, Young and Smith (2005, Section 9.7). Below, we sketch some brief
details on INLA; further details can be found elsewhere (Rue et al., 2009;
Blangiardo and Cameletti, 2015; Rue et al., 2017; Wang et al., 2018; Krainski
et al., 2018; Gómez-Rubio, 2020). The first step of INLA approximates the
marginal posterior of $\alpha$ via the Laplace approximation, that is,

$\begin{split}p(\alpha\mid\text{data})&=\frac{p(\alpha,\beta\mid\text{data})}{p(\beta\mid\alpha,\text{data})}\\\
&\approx\frac{p(\alpha,\beta\mid\text{data})}{\widetilde{p}(\beta\mid\alpha,\text{data})}\bigg{|}_{\beta=\beta_{\alpha}^{*}},\end{split}$
(9)

where $\widetilde{p}(\beta\mid\alpha,\text{data})$ is the Gaussian
approximation based on the mode—and the curvature around the mode—of the full
conditional of $\beta$, and where $\beta_{\alpha}$ is the mode of this
approximated full conditional for a given $\alpha$. Next, INLA approximates
the marginal posterior of each component of $\beta$. Let $\beta_{-i}$ be the
elements of $\beta$, except $\beta_{i}$. Similarly to (9) it follows that

$p(\beta_{i}\mid\alpha,\text{data})\propto\frac{p(\alpha,\beta\mid\text{data})}{p(\beta_{-i}\mid\alpha,\beta_{i},\text{data})},$
(10)

which can be approximated using a Laplace approximation for
$p(\beta_{-i}\mid\alpha,\beta_{i},\text{data})$; faster approximations are
also available from Rue et al. (2009). Finally, the marginal posterior density
$p(\beta_{i}\mid\text{data})$ is obtained by numerically integrating out
$\alpha$. Independent samples from the full posterior of $\bm{\beta}$ can then
be generated following Seppä et al. (2019, Section 2.5.2), which can then be
used for estimating functionals of the parameters of interest.

Estimation and inference for balls of maximum dissimilarity can be conducted
with an algorithm that combines the deterministic nature of INLA along with
sampling, according to the steps below. Data from processes $X$ and $Y$ are
respectively denoted by $\mathcal{D}_{X}$ and $\mathcal{D}_{X}$.

Algorithm 1 INLA-Based Posterior Inference for BMDs

Fit the marginal posterior densities, $p(\beta_{X,i}\,|\,\mathcal{D}_{X})$ and
$p(\beta_{Y,i}\,|\,\mathcal{D}_{Y})$, using the Integrated Nested Laplace
Approximation, for $i=0,\dots,B_{X}$ and $j=0,\dots,B_{Y}$.

Sample $m$ posterior draws from the full posterior of
$\beta_{X}=(\beta_{X,0},\dots,\beta_{X,B_{X}})$ and
$\beta_{Y}=(\beta_{Y,0},\dots,\beta_{Y,B_{Y}})$, given by
$\beta_{X}^{(1)},\dots,\beta_{X}^{(m)}\,\operatorname{\overset{iid}{\sim}}\,p(\beta_{X}\,|\,\mathcal{D}_{X}),\qquad\beta_{Y}^{(1)},\dots,\beta_{Y}^{(m)}\,\operatorname{\overset{iid}{\sim}}\,p(\beta_{Y}\,|\,\mathcal{D}_{Y}),$
so to generate $m$ posterior trajectories from the functional parameters using
(8); that is, for $k=1,\dots,m$, do
$\theta_{X}^{(k)}(t)=g\bigg{(}\beta_{X,0}^{(k)}+\sum_{i=1}^{B_{X}}\beta_{X,i}^{(k)}\phi_{i}(t)\bigg{)},\qquad\theta_{Y}^{(k)}(t)=g\bigg{(}\beta_{Y,0}^{(k)}+\sum_{j=1}^{B_{Y}}\beta_{Y,j}^{(k)}\phi_{j}(t)\bigg{)}.$
(11)

Use the posterior trajectories from Step 2 to obtain a sequence of posterior
BMD draws for
$\\{B_{c}^{(k)*}\equiv{B(t_{c}^{(k)*},r_{c}^{(k)*})}\\}_{k=1}^{m},$ with
$(t_{c}^{(k)*},r_{c}^{(k)*})$ solving the optimization problem in (7), for
$k=1,\dots,m$.

Algorithm 1 warrants some comments. Step 1 is deterministic, it follows by
numerically integrating out the hyperparameters in (10), and to facilitate its
implementation we recommend using the R-INLA package (Martins et al., 2013;
Lindgren et al., 2015) from R (R Development Core Team, 2022), which is also
equipped with routines that facilitate the implementation of Step 2 (e.g.
inla.posterior.sample). Step 3 boils down to solving the optimization problem
in (7) using the pair of posterior trajectories in (11). Except where
mentioned otherwise, in all experiments reported below we draw $m=1\,000$
times from the posterior distribution of BMDs using Algorithm 1.

##  3 Variants, Consequences, and Extensions

###  3.1 Multi-Objective RMDs

We now consider the context where the interest is on learning about regions
over which a set of marginal features of two processes most differ. Let

$\begin{cases}\theta_{X}=(\theta_{X,1},\dots,\theta_{X,q})=(\theta_{X,1}(t),\dots,\theta_{X,q}(t)),\\\
\theta_{Y}=(\theta_{Y,1},\dots,\theta_{Y,q})=(\theta_{Y,1}(t),\dots,\theta_{Y,q}(t)),\end{cases}$

be functional parameters characterizing marginal features of these processes.
The learning problem consists of searching for a set $A$, with volume $|A|\leq
c$, over which the overall difference between $\theta_{X}$ and $\theta_{Y}$ is
highest. Formally, the concept below entails the combination of two fields of
optimization on whose interface little is known, namely: Multi-objective
optimization (e.g. Pardalos et al., 2017) and set function optimization. In
this section

$F_{i}(A)=\|\theta_{X,i}-\theta_{Y,i}\|^{(A)},$ (12)

denotes the set objective functions of interest, for $i=1,\dots,q$.

######  Definition 3 (Region of Multi-Maximum Dissimilarity)

Suppose $\theta_{X},\theta_{Y}\in L^{p}(T)$ and $T\subset\mathbf{R}^{d}$ is
compact. The region of multi-maximum dissimilarity (RMMD) is defined as a set
$A^{*}_{c}$ that solves,

$\begin{array}[]{rl}\arg\underset{A\subseteq
T}{\max}&(F_{1}(A),\dots,F_{q}(A))\\\ \textit{s.t.}&|A|\leq c,\\\
&A~{}\text{is compact and convex},\end{array}$ (13)

for a fixed $c\geq 0$, where $F_{i}(A)$ is defined as in (12). In addition,
$D_{c}^{*}=(F_{1}(A^{*}_{c}),\dots,F_{q}(A^{*}_{c}))$ is said to be the multi-
dissimilarity index.

Ideally, we would aim to simultaneously maximize all targets
$F_{i}(A)=\|\theta_{X,i}-\theta_{Y,i}\|^{(A)}$, for all $i$. Yet, in practice
targets may conflict each other, that is, if one is increased some other may
be decreased; we illustrate this situation in Section 5.2. To give a more
concrete meaning to Definition 3, an ordering concept over the collection of
feasible sets $\mathscr{F}$ is required, so to rank the corresponding
objective function values. The following Pareto optimal concept induces that
order, and it defines optimality via a compromise across all objective set
functions, in the sense that improvements on one target cannot be made at the
cost of deteriorating another target.

######  Definition 4 (Pareto Optimal Region of Multi-Maximum Dissimilarity)

Let $F_{i}(A)$ be defined as in (12). The set $A^{*}_{c}\in\mathscr{F}_{c}$ is
a Pareto optimal region of multi-maximum dissimilarity if there exists no
other set $A\in\mathscr{F}_{c}$ such that $F_{i}(A)\geq F_{i}(A^{*}_{c})$, for
all $i\in\\{1,\dots,q\\}$, and $F_{i}(A)>F_{j}(A^{*}_{c})$, for at least one
$j\in\\{1,\dots,q\\}$.

In common with standard theory for multi-objective function, in our set
function context the set of Pareto optimal RMMDs can be analytically obtained
only in very specific cases, and thus we need to resort to scalarization
(Pardalos et al., 2017, Ch. 2). The aim of scalarization is to reduce a multi-
objective problem into a single-objective problem. Here, we define and
characterize the following linear scalarization method for our set function
context

$\mathcal{F}_{w}(A)=\sum_{i=1}^{q}w_{i}F_{i}(A),$ (14)

$w=(w_{1},\dots,w_{q})\in(0,\infty)^{q}$, and we will refer to

$A_{w,c}^{*}=\arg\underset{A\in\mathscr{F}_{c}}{\max}\,\mathcal{F}_{w}(A),$
(15)

as the solution to the set function linear scalarization problem with weight
$w$.

###### Theorem 3

Every solution to the linear scalarization problem is a Pareto optimal RMMD.

Section 5.2 illustrates how Theorem 3 can be applied in practice.

###  3.2 RMDs-Based on Averaging

This section shows that a simple modification of the notion of RMDs based on
‘averaging’ leads to links with other well-known concepts. To streamline the
presentation of ideas we will focus in $L^{1}(T)$, and thus exceptionally over
this section we will write $B(t,r)$ to denote $\\{t:\|t-t\|_{1}\leq r\\}$, and
we will further set $\mu(A)=|A|$. Averaging is here understood in the usual
sense of the well-known average integral symbol,
$\mathchoice{{\vbox{\hbox{$\textstyle-$
}}\kern-7.83337pt}}{{\vbox{\hbox{$\scriptstyle-$
}}\kern-5.90005pt}}{{\vbox{\hbox{$\scriptscriptstyle-$
}}\kern-4.75003pt}}{{\vbox{\hbox{$\scriptscriptstyle-$
}}\kern-4.25003pt}}\\!\int$, which is defined as

$\mathchoice{{\vbox{\hbox{$\textstyle-$
}}\kern-7.83337pt}}{{\vbox{\hbox{$\scriptstyle-$
}}\kern-5.90005pt}}{{\vbox{\hbox{$\scriptscriptstyle-$
}}\kern-4.75003pt}}{{\vbox{\hbox{$\scriptscriptstyle-$
}}\kern-4.25003pt}}\\!\int_{B(t,r)}f(u)\,\mathrm{d}u=\frac{1}{|B(t,r)|}\int_{B(t,r)}f(u)\,\mathrm{d}u,$
(16)

for an absolutely integrable $f$ and for $t\in T$ and $r\geq 0$. For $r=0$,
the expression in (16) should be understood as the limit,
$\mathchoice{{\vbox{\hbox{$\textstyle-$
}}\kern-7.83337pt}}{{\vbox{\hbox{$\scriptstyle-$
}}\kern-5.90005pt}}{{\vbox{\hbox{$\scriptscriptstyle-$
}}\kern-4.75003pt}}{{\vbox{\hbox{$\scriptscriptstyle-$
}}\kern-4.25003pt}}\\!\int_{B(t,r)}f(u)\,\mathrm{d}u\to f$, as $r\to 0$. There
are two motifs in this section:

  * •

Averaging the target: In the optimization problem that defines BMDs
(Definition 2 and Equation (7)), we now replace $\int$ with
$\mathchoice{{\vbox{\hbox{$\textstyle-$
}}\kern-7.83337pt}}{{\vbox{\hbox{$\scriptstyle-$
}}\kern-5.90005pt}}{{\vbox{\hbox{$\scriptscriptstyle-$
}}\kern-4.75003pt}}{{\vbox{\hbox{$\scriptscriptstyle-$
}}\kern-4.25003pt}}\\!\int$ in the objective function.

  * •

Moving beyond random processes: While the focus in the previous sections has
been on modeling RMDs and BMDs for time and spatial processes, it would seem
natural to employ these concepts beyond the context of stochastic processes.
Thus, in this section $\theta_{X}(t)$ will $\theta_{Y}(t)$ will be allowed to
be general parameters, possibly unrelated with stochastic processes, such as
distribution functions associated with random variables $X$ and $Y$.

We define the averaged or Hardy–Littlewood BMD as a ball,
$\mathbb{B}_{c}^{*}=B(t_{c}^{*},r^{*}_{c})$, that maximizes,

$\mathbb{D}_{c}^{*}(\theta_{X},\theta_{Y})=\max\bigg{\\{}\mathchoice{{\vbox{\hbox{$\textstyle-$
}}\kern-7.83337pt}}{{\vbox{\hbox{$\scriptstyle-$
}}\kern-5.90005pt}}{{\vbox{\hbox{$\scriptscriptstyle-$
}}\kern-4.75003pt}}{{\vbox{\hbox{$\scriptscriptstyle-$
}}\kern-4.25003pt}}\\!\int_{B(t,r)}|\theta_{X}(u)-\theta_{Y}(u)|\,\mathrm{d}u:(t,r)\in
T\times[0,R_{c}]\bigg{\\}},$ (17)

where $\theta_{X},\theta_{Y}$ are absolutely integrable and $c\geq 0$. As can
be seen by comparing (17) with (7) the only difference between
Hardy–Littlewood BMDs and standard BMDs in (7) is the use of the average
integral symbol in (17) rather than the standard integral $\int$.111We will
refer to these BMDs as “Hardy–Littlewood”, given that
$\mathbb{D}_{c}^{*}(\theta_{X},\theta_{Y})$ in (17) has links with the so-
called Hardy–Littlewood maximal function (e.g. Tao, 2011, Section 1.6),
$(Mg)(t)=\sup\bigg{\\{}\mathchoice{{\vbox{\hbox{$\textstyle-$
}}\kern-7.83337pt}}{{\vbox{\hbox{$\scriptstyle-$
}}\kern-5.90005pt}}{{\vbox{\hbox{$\scriptscriptstyle-$
}}\kern-4.75003pt}}{{\vbox{\hbox{$\scriptscriptstyle-$
}}\kern-4.25003pt}}\\!\int_{B(t,r)}|g(u)|\,\mathrm{d}u:r>0\bigg{\\}},$ where
$g$ is absolutely integrable. The fact that $\mathbb{B}_{c}^{*}$ exists
follows by noticing that $f(t,r)=\mathchoice{{\vbox{\hbox{$\textstyle-$
}}\kern-7.83337pt}}{{\vbox{\hbox{$\scriptstyle-$
}}\kern-5.90005pt}}{{\vbox{\hbox{$\scriptscriptstyle-$
}}\kern-4.75003pt}}{{\vbox{\hbox{$\scriptscriptstyle-$
}}\kern-4.25003pt}}\\!\int_{B(t,r)}|\theta_{X}(u)-\theta_{Y}(u)|\,\mathrm{d}u$
is continuous and that Tikhonov’s theorem (Appendix A) implies that the search
domain, $T\times[0,R_{c}]$, is compact for every $c\geq 0$, as both $T$ and
$[0,R_{c}]$ are themselves compact. In addition, the following proposition
holds.

###### Proposition 1

Suppose $\theta_{X},\theta_{Y}\in L^{1}(T)$ are continuous and
$T\subset\mathbf{R}^{d}$ is compact. Then, for every $c>0$ it holds that
$\mathbb{D}_{c}^{*}(\theta_{X},\theta_{Y})\leq\max_{t\in
T}|\theta_{X}(t)-\theta_{Y}(t)|$, and as $c\to 0$,

$\mathbb{D}_{c}^{*}(\theta_{X},\theta_{Y})\to\underset{t\in
T}{\max}\,|\theta_{X}(t)-\theta_{Y}(t)|.$

The next example uses the lenses of Proposition 1 to show that a well-known
performance measure known as Youden index (e.g. Sokolova et al., 2006; Inácio
de Carvalho et al., 2017) has links with the averaged BMDs defined in (17).

######  Example 4 (Youden index and Kolmogorov metric)

If $\theta_{X}(t)=F_{X}(t)$ and $\theta_{Y}(t)=F_{Y}(t)$ are distribution
functions supported on $T\subset\mathbf{R}$, and if we set $c\to 0$, then
maximization takes place over singletons and the Lebesgue differentiation
theorem (Appendix A) yields that

$\lim_{r\to 0}\,\mathchoice{{\vbox{\hbox{$\textstyle-$
}}\kern-7.83337pt}}{{\vbox{\hbox{$\scriptstyle-$
}}\kern-5.90005pt}}{{\vbox{\hbox{$\scriptscriptstyle-$
}}\kern-4.75003pt}}{{\vbox{\hbox{$\scriptscriptstyle-$
}}\kern-4.25003pt}}\\!\int_{B(t,r)}|F_{X}(u)-F_{Y}(u)|\,\mathrm{d}u=|F_{X}(t)-F_{Y}(t)|.$

Thus, when $c\to 0$ the set function optimization problem in (17), becomes a
standard continuous optimization problem that can be written as

$\max_{t\in T}|F_{X}(t)-F_{Y}(t)|,$ (18)

which is the well-known Youden index. Thus, the Youden index is an
Hardy–Littlewood dissimilarity index, and the popular Youden’s J statistic
(i.e. the maximizer of (18)) is a singleton averaged BMD, as in (17) with
$c\to 0$. Another consequence of (18) is that the limiting Hardy–Littlewood
distance in (17) when $c\to 0$ is the well-known Kolmogorov metric (e.g.,
Gretton et al., 2012).

##  4 Numerical Study with Artificial Data

In this section we assess the performance of the proposed tools via a Monte
Carlo study.

#### Artificial data generating processes and simulation settings

Examples 1–2 from Section 2.1 will form the basis of this numerical workout.
Namely, we consider the following scenarios.

Scenario 1 (Mean Functions) BMDs between mean functions as in Example 1. Here,
$X_{t}$ and $Y_{t}$ are Gaussian processes with mean functions given in (2)
and where the same Matérn covariance function is assumed for both processes.
Specifically
$\text{cov}(X_{s},X_{t})=\text{cov}(Y_{s},Y_{t})=C_{\nu}(\|s-t\|_{2},\sigma,\ell)$,
where

$\begin{split}C_{\nu}(d,\sigma,\ell)=\sigma^{2}\frac{2^{1-\nu}}{\Gamma(\nu)}\bigg{(}\sqrt{2\nu}\frac{d}{\ell}\bigg{)}^{\nu}K_{\nu}\bigg{(}\sqrt{2\nu}\frac{d}{\ell}\bigg{)},\end{split}$

for $(s,t)\in[0,1]^{2}$, where $K_{\nu}$ is the modified Bessel function
(Abramowitz and Stegun, 1964, Section 9.6), and where $\sigma,\nu$, and $\ell$
are positive parameters, here set as $\sigma=\nu=\ell=1$. The simulated data
are then a discretized version of $n$ simulated Gaussian processes evaluated
over a grid on the unit interval,

$\mathcal{D}_{X}=\\{X_{t,i}:t\in\\{0/J,\dots,(J-1)/J\\}\\}_{i=1}^{n},$

with $n\in\\{10,50,100,200\\}$ and $J\in\\{10,20\\}$; the same comments apply
to $\mathcal{D}_{Y}$.

Scenario 2 (Intensity Functions) BMDs between intensity functions as in
Example 2. Here, points drawn from non-homogeneous bivariate Poisson process
with mean measures,

$\begin{cases}E\\{N_{X}(A)\\}=\int_{A}\gamma\exp\\{-(t_{1}^{2}+t_{2}^{2})/2\\}\,\mathrm{d}t,\\\
E\\{N_{Y}(A)\\}=\delta E\\{N_{X}(A)\\},\end{cases}$ (19)

for $A\subseteq T=[-3,3]^{2}$. While the sample sizes in this scenario are
random quantities, given by $N_{X}=N_{X}(T)$ and $N_{Y}=N_{Y}(T)$, the mean
number of simulated points over $T$ is
$(E(N_{X}),E(N_{Y}))\approx(2\pi\gamma,2\pi\gamma\delta)$, a simple yet
accurate approximation that follows immediately from multiple Gaussian
integrals. The simulated data are given by the following collection of points

$\mathcal{D}_{X}=\\{(X_{1,1},X_{2,1}),\dots,(X_{1,N_{X}},X_{2,N_{X}})\\},$

with $\gamma\in\\{25,50\\}$ and $\delta\in\\{2,4,8,16\\}$, and where
$\mathcal{D}_{Y}$ is analogously defined.

#### Modeling, prior specification, and posterior inference

Inferences for the BMDs were carried out by sampling $m=1\,000$ and 500 times
using Algorithm 1 for Scenarios 1 and 2, respectively. As can be seen from
Algorithm 1, inferences for BMDs are constructed from the functional
parameters, and thus we now comment on what versions of (8) have been used for
fitting the latter. For Scenario 1 the identity link function was used in (8)
along with B-spline basis, and the number of basis functions was selected
using the DIC (Deviance Information Criterion; Spiegelhalter et al. 2002,
2014). The default uninformative priors of R-INLA have been used, which
consist of diffuse priors for the $\beta$’s—i.e. $\beta_{0}\sim N(0,\infty)$
and $\beta_{i}\sim N(0,1000)$—and a long-tailed prior for the variance of the
error term—i.e. a log gamma distribution, where the gamma distribution has
mean $a/b$ and variance $a/b^{2}$, with $a=1$ and $b=10^{-5}$; see Wang et al.
(2018, Section 5.2.1) for further details. For Scenario 2 we follow Simpson et
al. (2016) and specify a log-Gaussian Cox process using (8) by setting a log
link function, that links the intensity function with a Matérn random field
using piecewise linear basis functions over a mesh, and where the $\beta$’s
are Gaussian-distributed. For the parameters of the Matérn covariance function
we use the PC prior approach of Fuglstad et al. (2019) setting
$P(\sigma>1)=0.001$ and $P(\ell<0.05)=0.001$.

Before we move to the Monte Carlo study, we first illustrate the methods on a
single run experiment for some instances of Scenarios 1 and 2.

Scenario 1

Scenario 2

(a)

(b)

(c)

(d)

Figure 3: One shot experiments for Scenarios 1 and 2 (a) The fitted BMDs
([black]) based on the simulated Gaussian process data are depicted in the
time axis and compared against the true ( [gray]); the chart also depicts the
fitted functional parameters (dashed) along with $95\%$ credible bands, and
the true parameters (solid). (b) Marginal posterior density for the center
$t^{*}_{c}$ plotted against the true center (rug). (c) and (d) are identical
to (a) and (b), respectively, but for the simulated point process data from
Scenario 2.

#### One shot experiments

Let’s start with Scenario 1 and consider for illustration BMDs with length
$c=0.1$ and $c=0.2$. In Fig. 3(a) we show the fitted BMDs, along with the
corresponding mean functions, on a one shot experiment with $n=10$ and $J=10$.
As can can be seen from Fig. 3(a), the fitted BMDs accurately recover the true
$B^{*}_{0.1}=[0.165,0.265]$ and $B^{*}_{0.2}=[0.115,0.315]$. In Fig. 3(b) we
also display the marginal posterior density for the optimal center $t^{*}_{c}$
which quantifies the uncertainty surrounding the true. The marginal posterior
for the radius is essentially degenerated, as predicted earlier in the
comments surrounding (7), and hence not shown.

Let’s now move to Scenario 2 and consider for illustration a BMD with area
$c=2\pi$. In Fig. 3(c) we depict the fitted BMD, on a one shot experiment with
$\gamma=25$ and $\delta=2$, and display the fitted intensity functions. As it
is evident from Fig. 3(c), the fitted BMD nicely uncovers the true, and indeed
it completely overlaps it to the point that the true (depicted in gray) is
barely visible.

#### Monte carlo evidence

We now report the main findings of the Monte Carlo simulation study. We redo
the previous one shot analyses $M=1\,000$ times, considering different samples
sizes, and relying on the GHE (Posterior Mean Global Hausdorff Error),

$\text{GHE}=E\bigg{\\{}\int^{|T|}D_{\mathrm{H}}(B^{*}_{c},\widehat{B}^{*}_{c})\,\mathrm{d}c\,\bigg{|}\,\mathcal{D}_{X},\mathcal{D}_{Y}\bigg{\\}}$
(20)

so to quantify how accurate on average are the estimated BMDs,
$\widehat{B}^{*}_{c}$, over $0\leq c\leq|T|$.

Scenario 1 Scenario 2
---
Figure 4: Side-by-side boxplots of GHE (Posterior Mean Global Hausdorff Error)
for Monte Carlo simulation study.

Some comments on the computation of (20) are in order. In Scenario 1 we use
$D_{\mathrm{H}}(B^{*}_{c},\widehat{B}^{*}_{c})=\max\\{|(t_{0,c}^{*}-r^{*}_{c})-(\widehat{t}_{0,c}^{*}-\widehat{r}^{*}_{c})|,|(t_{0,c}^{*}+r^{*}_{c})-(\widehat{t}_{0,c}^{*}+\widehat{r}^{*}_{c})|\\}$,
while in Scenario 2 we use a numerical approximation of
$D_{\mathrm{H}}(B^{*}_{c},\widehat{B}^{*}_{c})$ implemented using Borchers
(2021). Finally, the GHE for each simulated dataset is computed as

$\text{GHE}_{j}=\frac{1}{m}\sum_{i=1}^{m}\int^{|T|}D_{\mathrm{H}}(B^{*}_{c},\widehat{B}_{c}^{*[i,j]})\,\mathrm{d}c,$

where $\widehat{B}_{c}^{*[i,j]}$ is the $i$th posterior sampled BMD based on
the $j$th simulated sample, for $i=1,\dots,m$ and $j=1,\dots,M$. As can be
seen from Fig. 4, GHE tends to decrease as $n,J,\gamma$, and $\delta$
increases. Such behavior confirms the expected frequentist behavior of the
methods, as $n$ and $J$ dictates the amount of simulated data for Scenario 1,
and $\gamma$ and $\delta$ does the same for Scenario 2. To put it differently,
since larger values of these parameters imply larger sample sizes, the
observed reduction in GHE as a function of the latter parameters suggests a
sensible asymptotic performance of the proposed Bayesian inferences.

##  5 Empirical Section

###  5.1 Thefts in Buenos Aires

Buenos Aires is the most dense metropolis in Argentina and its crime rates are
significantly higher in comparison to the rest of the country. In this section
we illustrate how the proposed method can reveal regions of the city where
nonviolent crimes—such as burglary, pickpocketing or nonviolent thefts—have
changed the most, comparing the years 2019 (pre COVID-19) and 2020 (when the
COVID lockdown took place during several months). The data are publicly
available online in the city hall web page, and consist of point process data
on the latitude and longitude where thefts occurred during 2019
($\mathcal{D}_{2019}$) and 2020 ($\mathcal{D}_{2020}$). Here, the functional
parameters of interest are the intensity functions

$\theta_{2019}(\texttt{latitude},\texttt{longitude}),\quad\theta_{2020}(\texttt{latitude},\texttt{longitude}),$

and its BMD will represent region of the city, of a given size, where the most
noteworthy changes in thefts took place. The fitted BMDs were modeled
according to (8) using again a log-Gaussian Cox process and a similar PC prior
specification as in Section 4.

In Fig. 5(a) we depict the estimated BMD corresponding to an area of $8$km2
along with an heat map of the differences in the estimated posterior intensity
functions between consecutive years; the value of $8$km2 was chosen for
illustration as it corresponds to about twice the size of the largest
neighborhood, which is Palermo. We also depict in Fig. 5(b) an heat map of the
posterior density corresponding to the center of the BMD which shows that
these are substantially concentrated, thus suggesting low uncertainty on the
fitted BMD.

To clarify the applied meaning of such BMDs we provide some additional
background on the social context surrounding the empirical analysis. In Fig.
5(c–d) we depict the fitted intensity functions corresponding to both years.
As can be seen from Fig. 5, in 2019 and 2020 thefts were more likely to occur
in APRV (Almagro, Palermo, Recoleta, and Villa Crespo) which are some of the
neighborhoods where several commercial and touristic activities took place.
Yet, important differences on the estimated intensity functions are perceived
between both years. During the first half of year 2020, local authorities took
strong social distancing measures such as the limitation to the access the
public transportation system, restrictions on business and commerce during the
day, limitations on gatherings and tourism activities, restriction to the
capacity in bars and restaurants, among others. The difference on the
estimated intensity functions between consecutive years evident from Fig.
5—and the implied reduction of thefts over 2020—is in line with the findings
of Mohler et al. (2020) that report similar evidence on the effect of COVID-19
lockdown and social distance policies in nonviolent crime.

(a)

(b)

---

(c)

(d)

Figure 5: (a) Fitted BMD corresponding to $c=8$km2 along with a contour of
pointwise differences between fitted posterior intensity functions over
consecutive years. (b) Posterior density estimation corresponding to the
centers of of the fitted BMD. (c) and (d) shows the fitted intensity functions
for 2019 and 2020.

The BMD in Fig. 5(a) suggests that APRV (Almagro, Palermo, Recoleta, and Villa
Crespo) are the neighborhoods where there was a most impactful effect of
COVID-19 lockdown. To put it differently, while nonviolent crime has decreased
during lockdown over the entire city, what the fitted BMD highlights is that
such reduction was relatively much higher in the APRV neighborhoods.

###  5.2 Volatility in Stock Markets

Our second illustration will shed light on the multi-objective approach from
Section 3.1. Data were gathered from Yahoo Finance and consist of monthly
values of the NYSE and NASDAQ composite indices from the New York Stock
Exchange ($\\{X_{t}\\}$) and the NASDAQ Stock Exchange ($\\{Y_{t}\\}$),
respectively. The data ranges from January 1990 to May 2021, thus covering a
variety of episodes of financial turbulence such as the dotcom tech bubble
that peaked around 2000, the subprime crisis that started around 2007, and the
recent COVID–19 global pandemic. In the multi-objective RMD analysis to be
conducted here, we will consider two functional parameters of interest: The
mean values of the indices over time and the volatility of their log returns,
that is,

$\begin{cases}m_{X}(t)=E(X_{t}),\\\
m_{Y}(t)=E(Y_{t}),\end{cases}\qquad\begin{cases}\sigma_{X}(t)=[E\\{\log(X_{t}/X_{t-1})^{2}\\}]^{1/2},\\\
\sigma_{Y}(t)=[E\\{\log(Y_{t}/Y_{t-1})^{2}\\}]^{1/2}.\end{cases}$

These functional parameters were modeled according to (8) exactly as in
Section 4, that is, using an identity link function and B-spline basis
functions, choosing the number of basis functions using the DIC, and using the
same uninformative prior.

We consider intervals of six months, one year, and two years (corresponding to
$c=6,12$, and 24) and we aim to evaluate on what periods of such length these
two stock markets differed the most—in terms of both average returns as well
as volatility. We thus seek for the interval of time
$B^{*}_{c}=[t_{c}^{*}-r_{c}^{*},t_{c}^{*}+r_{c}^{*}]$ that as in (15)
maximizes the following scalarized set function optimization problem,

$\begin{array}[]{rl}{\max}\\{\mathcal{F}_{w}\\{B(t,r)\\}:(t,r)\in
T\times[0,R_{c}]\\},\end{array}$ (21)

where $w$ is the scalarization parameter and

$\mathcal{F}_{w}\\{B(t,r)\\}=w\int_{t-r}^{t+r}|m_{X}(u)-m_{Y}(u)|\,\mathrm{d}u+(1-w)\int_{t-r}^{t+r}|\sigma_{X}(u)-\sigma_{Y}(u)|\,\mathrm{d}u.$

It follows from Theorem 3 that every solution to the linear scalarization
problem (21) is a Pareto optimal BMMD (ball of multi-maximum dissimilarity),
and hence Pareto optimal BMMDs obtained by linear scalarization have the nice
feature of allowing one to put more emphasis on the mean values or on
volatilities according to how we set $w$. That is, by setting $w=0$ or $w=1$,
we only consider volatilities or mean values respectively and the analysis
corresponds to a standard BMD, while for $w\in(0,1)$ absolute values in
differences between mean functions are more important than those in
volatilities as $w$ increases. In terms of computing we adapt Algorithm 1 to
the multi-objective context. That is, inference about the BMMDs is conducted
by sampling $m=1\,000$ times from the posteriors for means and
volatilities—and rather than solving (7) as in Algorithm 1—we now solve the
scalarized set function optimization problem in (21).

BMD for mean $(w=1)$

BMD for volatility $(w=0)$

Multi-objective BMD for mean–volatility $(w=2/3)$

(a)

(b)

(c)

(d)

(e)

(f)

Figure 6: (a, c, e) The brackets on the time axis correspond to the fitted
multi-objective BMDs for six monts, one year, and two years for $w=1$, $w=0$,
and $w=2/3$; the charts also show the realized (solid) and estimated (dotted)
volatilities corresponding to NYSE (red) and NASDAQ (blue). (b, d, f) Marginal
posterior density for the centre of the BMD.

In Fig. 6 we compare the results of the single BMD analysis [panels (a) to
(d)] versus the multi-objective analysis [panels (e) and (f)] considering a
value of $c$ corresponding to a period of six months, one year, and two years.
As can be seen from Fig. 6(a–b), the BMD associated with the mean ($w=1$)
concentrates around the subprime crisis, indicating that July 2006 to July
2008 is the period over which mean levels of NYSE and NASDAQ differed the
most. In Fig. 6(c–d), we see that the BMD associated with volatility ($w=0$)
ranges from October 1999 to October 2001—which corresponds to the dotcom
bubble burst. Finally, the multi-objective approach with $w=2/3$ is depicted
in Fig. 6 (e–f) and it suggests that volatility has a greater control on the
objective function in (21); that is, even when we set $w=2/3$—that is, even
when we set more emphasis on the differences in means rather than the
differences in volatility—we still get a similar result as setting $w=0$, as
we end up recovering the period of the dotcom bubble burst as can be seen from
Fig. 6(e–f).

###  5.3 Electrocardiogram Data (ECG200)

For our final illustration we use the ECG200 dataset contributed by Olszewski
(2001). The data is the result of monitoring electrical activity recorded
during one heartbeat and it consists of 200 ECG signals sampled at 96 time
instants, corresponding to 133 normal heartbeats ($\mathcal{D}_{X}$) and 67
myocardial infarction signals ($\mathcal{D}_{Y}$); the data are publicly
available from the UCR Time Series Classification and Clustering website.

In this illustration, the functional parameters of interest are the mean
functions of ECG signals for both classes (normal heartbeat
$\theta_{X}(t)=\operatorname{E}(X_{t})$, and myocardial infarction
$\theta_{Y}(t)=\operatorname{E}(Y_{t})$) and one of the goals of the analysis
is to track down periods, of a given length, over a cardiac cycle where the
differences between the two classes is most pronounced. To model these
functional parameters the Gaussian process prior specification in (8) was once
more applied using an identity link, B-spline basis functions, the DIC to
select the number of basis functions, and a Matérn covariance function with
the PC prior of Fuglstad et al. (2019) setting $P(\sigma>1)=0.001$ and
$P(\ell<0.05)=0.001$.

|
---|---

(a)

(b)

Figure 7: ECG200 dataset: (a) The brackets on the time axis correspond to the
fitted BMDs,
$\widehat{B}^{*}_{10}\subset\widehat{B}^{*}_{20}\subset\widehat{B}^{*}_{30}$;
the chart also shows the raw data (dots) and fitted mean functions (dashed)
along with $95\%$ credible bands, where red and blue respectively correspond
to normal heartbeats and myocardial infarction signals. (b) Marginal posterior
density for the centre of the set of maximum dissimilarity $B^{*}_{10}$.

In Fig. 7(a), we plot the fitted posterior estimates for a sequence of BMDs
($\widehat{B}^{*}_{10}$, $\widehat{B}^{*}_{20}$, and $\widehat{B}^{*}_{30}$)
using brackets in the time axis, along with the raw data, and the fitted mean
functions with 95% credible bands. The obtained BMDs uncover periods of a
given length where we observe largest differences between the estimated mean
ECG functions. Informally, we may think of such BMDs as corresponding to
intervals of about 10–30 deciseconds, since the 96 time instants cover about a
cardiac cycle for all subjects and thus are not expected to last longer than 1
second. The fitted BMD centers are localized around the time instants 22–55.
All in all, the analysis suggests that while normal heartbeats and myocardial
infarction signals have similar ‘peaks’ at the beginning of the sample period
(i.e. they have similar Q waves, in ECG signal analysis terminology),
immediately in the period right after (i.e. over their ST segments) they
greatly differ.

We close the analysis with two final remarks. First, in this illustration the
fitted BMDs verify the chained inclusions
$\widehat{B}^{*}_{10}\subset\widehat{B}^{*}_{20}\subset\widehat{B}^{*}_{03}$,
but this property does not hold in general—nor for the fitted BMDs, nor for
the true ones; counterexamples can be constructed either numerically or
analytically. Second, although BMDs are unrelated to classification, since the
ECG200 dataset is a popular benchmark for new classifiers it may be sensible
to ask whether the accuracy of some classifiers at discriminating outcome
classes (diseased–nondiseased) can be improved by focusing on BMDs rather than
treating the entire time horizon equally; we leave such open problem for
future analysis.

##  6 Closing Remarks

Regions of maximum dissimilarity and their variants are here proposed as a
tool for acquiring knowledge on the region with a given size where two
stochastic processes differ the most. The proposed learning problem is shown
to be equivalent to a continuous set function optimization on a monotone
modular function, under a Lebesgue measure constraint. As a byproduct, the
paper contributes to the literature on set function optimization which thus
far is focused on a discrete and combinatorial framework, and which has been
focused on monotone submodular functions (e.g. Nemhauser et al., 1978;
Calinescu et al., 2011; Goldengorin, 2009; Buchbinder et al., 2017; Buchbinder
and Feldman, 2018), typically under cardinality or matroid constraints.

The existence of the proposed regions of maximum dissimilarity is nontrivial
but we prove their existence, and illustrate with artificial and real data
that it only requires a moderate computational investment to learn them from
data. The proposed methods are developed in full generality for the setting
where the data of interest are themselves stochastic processes, and thus the
proposed toolbox can be used for unveiling the regions of maximum
dissimilarity with a given volume, for a variety of random process data. All
modeling was framed within a latent Gaussian framework, with inference being
conducted using the Integrated Nested Laplace Approximation; clearly, other
computational approaches could have been employed as well, including, for
example, variational inference (Blei et al., 2017).

A multi-objective version of the proposed framework is also devised so to
learn about multi-objective RMDs, where several functional parameters are
considered—each characterizing a specific feature of the processes being
compared. In addition, another variant of the method to which we refer to as
Hardy–Littlewood BMDs showcases that the current framework includes the Youden
index and the Youden’s J statistic as well as the Kolmogorov metric as
particular cases.

While the theoretical developments from Section 2.1 establish the existence of
general compact and convex sets of maximum dissimilarity, BMDs turn out be a
much convenient simplification for a variety of reasons. First, numerical
optimization is much more challenging with general RMDs, whereas for BMDs it
is relatively simple as can be seen from (7). Second, the inference for
general RMDs would entail averaging posterior simulated RMDs—that is,
averaging sets—whereas with BMDs we just need to compute posterior mean of the
$(d+1)$-dimensional centre-radius pair.

Acknowledgments. We thank, without implicating, Vanda Inácio de Carvalho for
comments and feedback and Finn Lindgren for discussions on INLA. MdC was
partially supported by FCT (Fundação para a Ciência e a Tecnologia, Portugal)
through the projects PTDC/MAT-STA/28649/2017 and UID/MAT/00006/2019.

##  Appendix A Technical Details and Auxiliary Lemmata

In this section we state some auxiliary facts that will be used to prove the
main results of this paper. Beyond the auxiliary lemmata to be stated below we
use some basic facts from measure, topology, and convex analysis, such as, for
example, Lebesgue differentiation theorem (e.g. Tao, 2011, Theorem 1.6.19),
Tikhonov’s theorem (Waldmann, 2014, Theorem 5.3.1), and the well-known fact
that the volume functional is continuous in the space of convex bodies
(Schneider, 2014, Theorem 1.8.16), under the Hausdorff metric.

Recall that Tikhonov’s theorem implies that from the Cartesian product of two
compact sets results a compact set. In addition, recall that Lebesgue
differentiation theorem implies that if $f:\mathbf{R}\to\mathbf{C}$ is an
absolutely integrable function, then for almost every $x\in\mathbf{R}^{d}$,

$\lim_{r\to 0}\frac{1}{|B(x,r)|}\int_{B(x,r)}f(y)\,\mathrm{d}y=f(x),$

where $B(x,r)=\\{y:\|x-y\|_{1}\leq r\\}$.

We now present the auxiliary lemmata. Blaschke selection theorem is a
classical result on convex analysis; the version stated below can be found,
for instance, in Benyamini (1998). In addition, we also recall below two key
results on optimization of correspondences (i.e. set-valued functions), namely
Berge’s maximum theorem and the product of correspondences theorem; the
versions below can be found in Aliprantis and Border (2006, Theorems 17.31 and
17.28).

######  Lemma 1 (Blaschke selection theorem)

The set of all compact convex subsets of a fixed compact subset of
$\mathbf{R}^{d}$ is compact under the Hausdorff metric.

######  Lemma 2 (Berge’s maximum theorem)

Let $\varphi:X\twoheadrightarrow Y$ be a continuous correspondence between
topological spaces with nonempty compact values, and suppose that $g:\text{Gr
}\varphi\to\mathbb{R}$ is continuous. Define the “value function”
$v:X\to\mathbf{R}$ by

$v(x)=\max\\{g(x,y):y\in\varphi(x)\\}$

and the correspondence $\alpha:X\twoheadrightarrow Y$ of maximizers by

$\alpha(x)=\\{y\in\varphi(x):g(x,y)=v(x)\\}.$

Then:

  1. 1.

The value function $v$ is continuous.

  2. 2.

The “argmax” correspondence $\alpha$ has nonempty compact values.

  3. 3.

If $Y$ is Hausdorff, then the “argmax” correspondence $\alpha$ is upper
hemicontinuous.

######  Lemma 3 (Product of correspondences theorem)

The product of correspondences obeys the following properties:

  1. 1.

The product of a family of upper hemicontinuous correspondences with compact
values is upper hemicontinuous with compact values.

  2. 2.

The product of a finite family of lower hemicontinuous correspondences is
lower hemicontinuous.

##  Appendix B Proofs of Main Results

###  B.1 Proof of Theorem 1

Claim 1. We start by showing that

$\mathscr{F}_{c}=\\{A\in\mathcal{A}:|A|\leq c\\}$

is compact, for every $c\geq 0$ where $\mathcal{A}$ is the family of compact
and convex subsets of the ground set $T$. Recall that by the Blaschke
selection theorem (Lemma 1), $\mathcal{A}$ is compact under the Hausdorff
metric. Further, since the volume functional $|\cdot|$ is continuous
(Schneider, 2014, Theorem 1.8.16), and given that $\mathscr{F}_{c}$ is the
preimage of the closed set $[0,c]$, it follows that $\mathscr{F}_{c}$ is a
closed subset of $\mathcal{A}$ and hence it is compact.

Observe next that maximizing $F(A)=\|\theta_{X}-\theta_{Y}\|^{(A)}_{p}$ is
equivalent to maximizing $F^{p}(A)$, for $p>0$, and as we show next $F^{p}(A)$
is upper semicontinuous under the Hausdorff metric, for every
$A\in\mathscr{F}_{c}$. Let $A_{n}\to A$ in $(\mathscr{F}_{c},d_{H})$, for a
fixed $c\geq 0$. It can be easily shown that (e.g. Schneider and Weil, 2008,
Theorem 12.3.6)

$1_{A}(t)\geq\lim\underset{n}{\sup}\,1_{A_{n}}(t),\quad t\in
T\subset\mathbf{R}^{d}.$ (22)

Combining (22) with Fatou’s lemma yields that

$\begin{split}F^{p}(A)&=\int_{T}1_{A}(t)|\theta_{X}(t)-\theta_{Y}(t)|^{p}\,\mu(\mathrm{d}t)\\\
&\geq\int_{T}\lim\underset{n}{\sup}\,1_{A_{n}}(t)|\theta_{X}(t)-\theta_{Y}(t)|^{p}\,\mu(\mathrm{d}t)\\\
&\geq\lim\underset{n}{\sup}\int_{T}1_{A_{n}}(t)|\theta_{X}(t)-\theta_{Y}(t)|^{p}\,\mu(\mathrm{d}t)\\\
&=\lim\underset{n}{\sup}\,F^{p}(A_{n}),\end{split}$

thus showing that $F^{p}(A)$ is upper semicontinuous under the Hausdorff
metric, for every $A\in\mathscr{F}_{c}$. The final step of the proof is
tantamount to a standard argument used for proving Weierstrass theorem. Let
$u_{c}=\sup\\{F^{p}(A):A\in\mathscr{F}_{c}\\}\cup\\{\infty\\}$. By definition,
for every $c$ there exists a maximizing sequence $A_{n}\in\mathscr{F}_{c}$
such that $F^{p}(A_{n})\to u_{c}$. By compactness, we can assume that
$A_{n}\to A^{*}_{c}$. Upper semicontinuity of $F^{p}(A)$ implies that
$u_{c}=\lim{\sup}_{n}\,F^{p}(A_{n})\leq F^{p}(A^{*}_{c}),$ and on the other
hand we have $u_{c}\geq F(A^{*}_{c})$ since $u_{c}$ is the supremum. This
proves that $u_{c}=F^{p}(A^{*}_{c})$ is the maximum of $F^{p}(A)$, subject to
$A\in\mathscr{F}_{c}$, and hence $A^{*}_{c}$ exists that solves (4).

Claim 2. First, observe that

$\begin{split}\|K\theta_{X}-K\theta_{Y}\|_{p}^{(A)}&=\|\alpha+\beta\theta_{X}-(\alpha+\beta\theta_{Y})\|_{p}^{(A)}\\\
&=|\beta|\times\|\theta_{X}-\theta_{Y}\|_{p}^{(A)},\end{split}$

from where it follows that a set $A$ that maximizes
$\|K\theta_{X}-K\theta_{Y}\|_{p}^{(A)}$ also maximizes
$\|\theta_{X}-\theta_{Y}\|_{p}^{(A)}$; that is, $A^{K}_{c}=A^{*}_{c}$ for all
$c\in[0,\infty)$. Second, it follows from the change of variables formula that

$\begin{split}\|L\theta_{X}-L\theta_{Y}\|_{p}^{(A^{L})}&=\bigg{[}\int_{A^{L}}\\{\theta_{X}(\alpha+\beta
t)-\theta_{Y}(\alpha+\beta t)\\}^{p}\mu(\mathrm{d}t)\bigg{]}^{1/p}\\\
&=\frac{1}{\beta^{1/p}}\bigg{[}\int_{A}\\{\theta_{X}(u)-\theta_{X}(u)\\}^{p}\mu(\mathrm{d}u)\bigg{]}^{1/p},\end{split}$

where $A=\\{\alpha+\beta t:t\in A^{L}\\}$. Thus,

$\arg\max\\{\|L\theta_{X}-L\theta_{Y}\|_{p}^{(A^{L})}:A^{L}\in\mathscr{F}_{c}\\}=\arg\max\\{\|\theta_{X}-\theta_{Y}\|_{p}^{(A)}:A\in\mathscr{F}_{\alpha+\beta
c}\\},$

and hence, $A^{L}_{c}=A^{*}_{\alpha+\beta c}$.

Claim 3. Note first that $D_{c}^{*}\geq 0$, for all $0<c<\infty$. Also, it
holds that $D_{c}^{*}=0$ if and only if $\theta_{X}=\theta_{Y}$ in
$\Theta_{c}$; indeed, $D_{c}^{*}=0$ implies that
$0=\|\theta_{X}-\theta_{Y}\|^{(A^{*}_{c})}_{p}\geq\|\theta_{X}-\theta_{Y}\|^{(A)}_{p}$,
for every $A\in\mathscr{F}^{c}$, which is only possible if
$\theta_{X}=\theta_{Y}$ in $\Theta_{c}$, as $c>0$. Finally, the triangle
inequality for $\|\cdot\|_{p}^{(A)}$ and $\|F\|_{\infty}=\max_{A}|F(A)|$
yields that

$\begin{split}D^{*}_{c}(\theta_{X},\theta_{Z})&=\max\\{\|\theta_{X}-\theta_{Z}\|^{(A)}_{p}:A\in\mathscr{F}_{c}\\}\\\
&=\max\\{\|(\theta_{X}-\theta_{Y})+(\theta_{Y}-\theta_{Z})\|^{(A)}_{p}:A\in\mathscr{F}_{c}\\}\\\
&\leq\max\\{\|\theta_{X}-\theta_{Y}\|^{(A)}_{p}+\|\theta_{Y}-\theta_{Z})\|^{(A)}_{p}:A\in\mathscr{F}_{c}\\}\\\
&\leq\max\\{\|\theta_{X}-\theta_{Y}\|^{(A)}_{p}:A\in\mathscr{F}_{c}\\}+\max\\{\|\theta_{Y}-\theta_{Z})\|^{(A)}_{p}:A\in\mathscr{F}_{c}\\}\\\
&=D^{*}_{c}(\theta_{X},\theta_{Y})+D^{*}_{c}(\theta_{Y},\theta_{Z}),\end{split}$

hence concluding the proof.

Claim 4. First, we note that an increase in $c$ represents augmenting the
search domain as

$\mathscr{F}_{a}\subset\mathscr{F}_{b},\quad\text{for }b>a.$

This combined with the fact that the set objective function
$F(A)=\|\theta_{X}-\theta_{Y}\|_{p}^{(A)}$ is non-decreasing ($A\subseteq B$
implies $F(A)\leq F(B)$) yields that $D^{*}_{b}\geq D^{*}_{a}$, for $b>a$.

###  B.2 Proof of Theorem 2

Claim 1. As a consequence of Tikhonov’s theorem (Appendix A), the search
domain $T\times[0,R_{c}]$ is compact for every $c\geq 0$. Next, it is a
routine exercise to prove that $f(t,r)=F\\{B(t,r)\\}$, is continuous for all
$(t,r)\in T\times[0,R_{c}]$ as both $T$ and $[0,R_{c}]$ are compact. The final
result then follows from Weierstrass theorem.

Claim 2. Let $\theta_{X},\theta_{Y}\in L^{p}(T)$ be fixed and set
$\mathscr{D}^{*}_{c}=\mathscr{D}^{*}_{c}(\theta_{X},\theta_{Y}).$ The proof
uses Berge’s maximum theorem (Lemma 2) which, as can be seen from Appendix A,
claims that a value function is continuous provided that both the objective
function and the constraint correspondence are continuous. In our setup, the
value function is

$\begin{split}\mathscr{D}^{*}_{c}&=\max\\{g(c,t,r):(t,r)\in\varphi(c)\\}\\\
&=\max\\{f(t,r):(t,r)\in\varphi(c)\\},\\\ \end{split}$

where $g(c,t,r)=f(t,r)=F\\{B(t,r)\\}$, and the constraint correspondence is
$\varphi:[0,\infty)\twoheadrightarrow\mathbf{R}^{d+1}$, defined by

$\varphi(c)=T\times[0,R_{c}],\qquad c\geq 0.$ (23)

Thus, we just have to prove that $\varphi(c)$ in (23) is continuous, given
that the objective function $g(c,t,r)=f(t,r)=F\\{B(t,r)\\}$ is trivially
continuous for all $(t,r)\in T\times[0,R_{c}]$. Continuity of $\varphi(c)$
follows immediately from the product of correspondences theorem (Lemma 3),
which yields that $\varphi(c)=\varphi_{1}(c)\times\varphi_{2}(c)$ is
continuous as both $\varphi_{1}(c)=T$ and $\varphi_{2}(c)=[0,R_{c}]$ are
compacted-valued, and $R_{c}$ is a continuous function for every $c\geq 0$.
Finally, upper semicontinuity of the argmax correspondence,

$\begin{split}\alpha_{c}&=\\{(t,r)\in
T\times[0,R_{c}]:g(c,t,r)=\mathscr{D}_{c}^{*}\\}\\\ &=\\{(t,r)\in
T\times[0,R_{c}]:f(t,r)=\mathscr{D}_{c}^{*}\\},\\\ \end{split}$

follows also from Berge’s maximum theorem as $\mathbf{R}^{d+1}$ is Hausdorff.

###  B.3 Proof of Theorem 3

Suppose by contradiction that $A^{\prime}_{w,c}$ is the solution to the set
function linear scalarization problem (15), for a fixed $w\in(0,\infty)^{q}$,
but that $A^{\prime}_{w,c}$ was not a Pareto optimal RMMD; then, there would
exist a Pareto improvement $A\in\mathscr{F}_{c}$ , with $A\neq
A^{\prime}_{w,c}$, so that $F_{i}(A)\geq F_{i}(A^{\prime}_{w,c})$, for all
$i$, and $F_{i}(A)>F_{i}(A^{\prime}_{w,c})$, for at least an $i$. But then,

$\sum_{i=1}^{q}w_{i}F_{i}(A)>\sum_{i=1}^{q}w_{i}F_{i}(A^{\prime}_{w,c}),$

which is a contradiction as $A^{\prime}_{w,c}$ solves the set function linear
scalarization problem.

###  B.4 Proof of Proposition 1

The proof of the first claim is straightforward. For every $c>0$,

$\begin{split}\mathbb{D}_{c}^{*}(\theta_{X},\theta_{Y})&=\max\bigg{\\{}\mathchoice{{\vbox{\hbox{$\textstyle-$
}}\kern-7.83337pt}}{{\vbox{\hbox{$\scriptstyle-$
}}\kern-5.90005pt}}{{\vbox{\hbox{$\scriptscriptstyle-$
}}\kern-4.75003pt}}{{\vbox{\hbox{$\scriptscriptstyle-$
}}\kern-4.25003pt}}\\!\int_{B(t,r)}|\theta_{X}(u)-\theta_{Y}(u)|\,\mathrm{d}u:(t,r)\in
T\times[0,R_{c}]\bigg{\\}}\\\
&=\max\bigg{\\{}\frac{1}{|B(t,r)|}\int_{B(t,r)}|\theta_{X}(u)-\theta_{Y}(u)|\,\mathrm{d}u:(t,r)\in
T\times[0,R_{c}]\bigg{\\}}\\\
&\leq\max\bigg{\\{}\frac{1}{\cancel{|B(t,r)|}}\cancel{|B(t,r)|}\max_{u\in
T}|\theta_{X}(u)-\theta_{Y}(u)|:(t,r)\in T\times[0,R_{c}]\bigg{\\}}\\\
&=\max_{t\in T}|\theta_{X}(t)-\theta_{Y}(t)|.\end{split}$

The second claim follows directly from the Lebesgue differentiation theorem
(Appendix A).

## References

  * Abramowitz and Stegun (1964) Milton Abramowitz and Irene A Stegun.  _Handbook of Mathematical Functions_.  Dover, New York, 1964.
  * Aliprantis and Border (2006) CD Aliprantis and KC Border.  _Infinite Dimensional Analysis_.  Springer, New York, 2006.
  * Benyamini (1998) Yoav Benyamini.  Applications of the universal surjectivity of the cantor set.  _The American Mathematical Monthly_ , 105(9):832–839, 1998.
  * Berrendero et al. (2016) José R Berrendero, Antonio Cuevas, and José L Torrecilla.  Variable selection in functional data classification: A maxima-hunting proposal.  _Statistica Sinica_ , pages 619–638, 2016.
  * Berrendero et al. (2020) José R Berrendero, Beatriz Bueno-Larraz, and Antonio Cuevas.  On Mahalanobis distance in functional settings.  _Journal of Machine Learning Research_ , 21(9):1–33, 2020.
  * Blangiardo and Cameletti (2015) Marta Blangiardo and Michela Cameletti.  _Spatial and Spatio-temporal Bayesian Models with R-INLA_.  Wiley, New York, 2015.
  * Blei et al. (2017) David M Blei, Alp Kucukelbir, and Jon D McAuliffe.  Variational inference: A review for statisticians.  _Journal of the American statistical Association_ , 112(518):859–877, 2017.
  * Borchers (2021) Hans W. Borchers.  _pracma: Practical Numerical Math Functions_ , 2021.  URL https://cran.r-project.org/web/packages/pracma/pracma.pdf.  R package version 2.3.6.
  * Buchbinder and Feldman (2018) Niv Buchbinder and Moran Feldman.  Submodular functions maximization problems.  In _Handbook of Approximation Algorithms and Metaheuristics, 2nd ed_ , pages 753–788. Chapman and Hall/CRC, Boca Raton, FL, 2018.
  * Buchbinder et al. (2017) Niv Buchbinder, Moran Feldman, and Roy Schwartz.  Comparing apples and oranges: Query trade-off in submodular maximization.  _Mathematics of Operations Research_ , 42(2):308–329, 2017.
  * Calinescu et al. (2011) Gruia Calinescu, Chandra Chekuri, Martin Pal, and Jan Vondrák.  Maximizing a monotone submodular function subject to a matroid constraint.  _SIAM Journal on Computing_ , 40(6):1740–1766, 2011.
  * Dette and Kokot (2021) Holger Dette and Kevin Kokot.  Bio-equivalence tests in functional data by maximum deviation.  _Biometrika_ , 108(4):895–913, 2021.  ISSN 0006-3444.  doi: 10.1093/biomet/asaa096.  URL https://doi.org/10.1093/biomet/asaa096.
  * Faouzi and Janati (2020) Johann Faouzi and Hicham Janati.  pyts: A python package for time series classification.  _Journal of Machine Learning Research_ , 21(46):1–6, 2020.  URL http://jmlr.org/papers/v21/19-763.html.
  * Ferraty and Vieu (2006) Frédéric Ferraty and Philippe Vieu.  _Nonparametric Functional Data Analysis: Theory and Practice_.  Springer, New York, 2006.
  * Fuglstad et al. (2019) Geir-Arne Fuglstad, Daniel Simpson, Finn Lindgren, and Håvard Rue.  Constructing priors that penalize the complexity of gaussian random fields.  _Journal of the American Statistical Association_ , 114(525):445–452, 2019.
  * Goldengorin (2009) Boris Goldengorin.  Maximization of submodular functions: Theory and enumeration algorithms.  _European Journal of Operational Research_ , 198(1):102–112, 2009.
  * Gómez-Rubio (2020) Virgilio Gómez-Rubio.  _Bayesian Inference with INLA_.  Chapman and Hall/CRC, address=Boca Raton, FL, 2020.
  * Gretton et al. (2012) Arthur Gretton, Karsten M Borgwardt, Malte J Rasch, Bernhard Schölkopf, and Alexander Smola.  A kernel two-sample test.  _Journal of Machine Learning Research_ , 13(1):723–773, 2012.
  * Horváth and Kokoszka (2012) Lajos Horváth and Piotr Kokoszka.  _Inference for Functional Data with Applications_.  Springer, New York, 2012.
  * Inácio de Carvalho et al. (2017) V. Inácio de Carvalho, M. de Carvalho, and A. J. Branscum.  Nonparametric Bayesian covariate-adjusted estimation of the Youden index.  _Biometrics_ , 73(4):1279–1288, 2017.
  * Krainski et al. (2018) Elias Krainski, Virgilio Gómez-Rubio, Haakon Bakka, Amanda Lenzi, Daniela Castro-Camilo, Daniel Simpson, Finn Lindgren, and Håvard Rue.  _Advanced Spatial Modeling with Stochastic Partial Differential Equations using R and INLA_.  Chapman and Hall/CRC, Boca Raton, FL, 2018.
  * Krause (2010) Andreas Krause.  SFO: A toolbox for submodular function optimization.  _Journal of Machine Learning Research_ , 11(38):1141–1144, 2010.  URL http://jmlr.org/papers/v11/krause10a.html.
  * Lindgren et al. (2015) Finn Lindgren, Håvard Rue, et al.  Bayesian spatial modelling with R-INLA.  _Journal of Statistical Software_ , 63(19):1–25, 2015.
  * Martins et al. (2013) Thiago G Martins, Daniel Simpson, Finn Lindgren, and Håvard Rue.  Bayesian computing with INLA: New features.  _Computational Statistics & Data Analysis_, 67:68–83, 2013.
  * Martos and de Carvalho (2018) G. Martos and M. de Carvalho.  Discrimination surfaces with application to region-specific brain asymmetry analysis.  _Statistics in Medicine_ , 11(37):1859–1873, 2018\.
  * Mohler et al. (2020) George Mohler, Andrea L Bertozzi, Jeremy Carter, Martin B Short, Daniel Sledge, George E Tita, Craig D Uchida, and P Jeffrey Brantingham.  Impact of social distancing during COVID-19 pandemic on crime in Los Angeles and Indianapolis.  _Journal of Criminal Justice_ , 68:101692, 2020.
  * Nemhauser et al. (1978) George L Nemhauser, Laurence A Wolsey, and Marshall L Fisher.  An analysis of approximations for maximizing submodular set functions—I.  _Mathematical Programming_ , 14(1):265–294, 1978\.
  * Nocedal and Wright (2006) Jorge Nocedal and Stephen Wright.  _Numerical Optimization_.  Springer, New York, 2006.
  * Olszewski (2001) Robert T Olszewski.  Generalized feature extraction for structural pattern recognition in time-series data.  Technical report, Carnegie-Mellon University, School of Computer Science, 2001.
  * Pardalos et al. (2017) Panos M Pardalos, Antanas Žilinskas, and Julius Žilinskas.  _Non-convex Multi-objective Optimization_.  Springer, New York, 2017.
  * Pini and Vantini (2016) Alessia Pini and Simone Vantini.  The interval testing procedure: A general framework for inference in functional data analysis.  _Biometrics_ , 72(3):835–845, 2016.
  * Pini and Vantini (2017) Alessia Pini and Simone Vantini.  Interval-wise testing for functional data.  _Journal of Nonparametric Statistics_ , 29(2):407–424, 2017.
  * R Development Core Team (2022) R Development Core Team.  _R: A Language and Environment for Statistical Computing_.  R Foundation for Statistical Computing, Vienna, Austria, 2022.
  * Ramsay and Silverman (2006) James O. Ramsay and B. W. Silverman.  _Functional Data Analysis_.  Springer, New York, 2006.
  * Ramsay and Silverman (2002) James O Ramsay and Bernard W Silverman.  _Applied Functional Data Analysis: Methods and Case Studies_.  Springer, New York, 2002.
  * Rue et al. (2009) Håvard Rue, Sara Martino, and Nicolas Chopin.  Approximate Bayesian inference for latent Gaussian models by using integrated nested Laplace approximations.  _Journal of the Royal Statistical Society: Series B (Statistical Methodology)_ , 71(2):319–392, 2009.
  * Rue et al. (2017) Håvard Rue, Andrea Riebler, Sigrunn H Sørbye, Janine B Illian, Daniel P Simpson, and Finn K Lindgren.  Bayesian computing with inla: A review.  _Annual Review of Statistics and Its Application_ , 4:395–421, 2017.
  * Schneider (2014) Rolf Schneider.  _Convex Bodies: The Brunn–Minkowski Theory_.  Cambridge University Press, Cambridge, 2014.
  * Schneider and Weil (2008) Rolf Schneider and Wolfgang Weil.  _Stochastic and Integral Geometry_.  Springer, New York, 2008.
  * Seppä et al. (2019) Karri Seppä, Håvard Rue, Timo Hakulinen, Esa Läärä, Mikko J Sillanpää, and Janne Pitkäniemi.  Estimating multilevel regional variation in excess mortality of cancer patients using integrated nested Laplace approximation.  _Statistics in Medicine_ , 38(5):778–791, 2019\.
  * Simpson et al. (2016) Daniel Simpson, Janine Baerbel Illian, Finn Lindgren, Sigrunn H Sørbye, and Havard Rue.  Going off grid: Computationally efficient inference for log-Gaussian Cox processes.  _Biometrika_ , 103(1):49–70, 2016.
  * Sokolova et al. (2006) Marina Sokolova, Nathalie Japkowicz, and Stan Szpakowicz.  Beyond accuracy, F-score and ROC: A family of discriminant measures for performance evaluation.  In _Australasian Joint Conference on Artificial Intelligence_ , pages 1015–1021. Springer, 2006.
  * Spiegelhalter et al. (2002) David J Spiegelhalter, Nicola G Best, Bradley P Carlin, and Angelika Van Der Linde.  Bayesian measures of model complexity and fit.  _Journal of the Royal Statistical Society: Series B (Statistical Methodology)_ , 64(4):583–639, 2002.
  * Spiegelhalter et al. (2014) David J Spiegelhalter, Nicola G Best, Bradley P Carlin, and Angelika Van der Linde.  The deviance information criterion: 12 years on.  _Journal of the Royal Statistical Society: Series B (Statistical Methodology)_ , 76(3):485–493, 2014.
  * Tao (2011) Terence Tao.  _An Introduction to Measure Theory_.  American Mathematical Society, Providence, RI, 2011.
  * Waldmann (2014) Stefan Waldmann.  _Topology: An Introduction_.  Springer, New York, 2014.
  * Wang et al. (2018) Xiaofeng Wang, Yuryan Yue, and Julian J Faraway.  _Bayesian regression modeling with INLA_.  Chapman and Hall/CRC, Boca Raton, FL, 2018.
  * Wu et al. (2019) Wei-Li Wu, Zhao Zhang, and Ding-Zhu Du.  Set function optimization.  _Journal of the Operations Research Society of China_ , 7(2):183–193, 2019.
  * Xu et al. (2020) Ganggang Xu, Ming Wang, Jiangze Bian, Hui Huang, Timothy R. Burch, Sandro C. Andrade, Jingfei Zhang, and Yongtao Guan.  Semi-parametric learning of structured temporal point processes.  _Journal of Machine Learning Research_ , 21(192):1–39, 2020.
  * Young and Smith (2005) G. A Young and Richard L Smith.  _Essentials of Statistical Inference_.  Cambridge University Press, Cambridge, UK, 2005.