# Optimal Decision-Making for Autonomous Agents via Data Composition

Émiland Garrabé, Martina Lamberti and Giovanni Russo∗ ∗ emails:
<EMAIL_ADDRESS><EMAIL_ADDRESS>E. Garrabé and G.
Russo with the DIEM at the University of Salerno, 84084, Salerno, Italy.

###### Abstract

We consider the problem of designing agents able to compute optimal decisions
by composing data from multiple sources to tackle tasks involving: (i)
tracking a desired behavior while minimizing an agent-specific cost; (ii)
satisfying safety constraints. After formulating the control problem, we show
that this is convex under a suitable assumption and find the optimal solution.
The effectiveness of the results, which are turned in an algorithm, is
illustrated on a connected cars application via in-silico and in-vivo
experiments with real vehicles and drivers. All the experiments confirm our
theoretical predictions and the deployment of the algorithm on a real vehicle
shows its suitability for in-car operation.

##  I Introduction

We often make decisions by composing knowledge gathered from others [1] and a
challenge transversal to control and learning is to devise mechanisms allowing
autonomous decision-makers to emulate these abilities. Systems based on
sharing data [2] are examples where agents need to make decisions based on
some form of crowdsourcing [3]. Similar mechanisms can also be useful for the
data-driven control paradigm when e.g., one needs to re-use policies
synthesized on plants for which data are available to solve a control task on
a new plant, for which data are scarcely available [3, 4, 5].

Motivated by this, we consider the problem of designing decision-making
mechanisms that enable autonomous agents to compute optimal decisions by
composing information from third parties to solve tasks that involve: (i)
tracking a desired behavior while minimizing an agent-specific cost; (ii)
satisfying safety constraints. Our results enable computation of the optimal
behavior and are turned into an algorithm. This is experimentally validated on
a connected car application.

#### Related works

we briefly survey a number of works related to the results and methodological
framework of this paper. The design of context-aware switches between multiple
datasets for autonomous agents has been recently considered in [3, 4], where
the design problem, formalized as a data-driven control (DDC) problem, did not
take into account safety requirements. Results in DDC include [6, 7, 8], which
take a behavioral systems perspective, [9], which finds data-driven formulas
towards a model-free theory of geometric control. We also recall e.g., [10,
11, 12] for results inspired from MPC, [5] that considers data-driven control
policies transfer and [13] that tackles he problem of computing data-driven
minimum-energy control for linear systems. In our control problem (see Section
III) we formalize the tracking of a given behavior via Kullback-Leibler (KL)
divergence minimization and we refer to e.g., [14, 15] for examples across
learning and control that involve minimizing this functional. Further, the
study of mechanisms enabling agents to re-use data, also arises in the design
of prediction algorithms from experts [16] and of learning algorithms from
multiple simulators [17]. In a yet broader context, studies in neuroscience
[18] hint that our neocortex might implement a mechanism composing the output
of the cortical columns and this might be the basis of our ability to re-use
knowledge.

_Contributions:_ we consider the problem of designing agents that dynamically
combine data from heterogeneous sources to fulfill tasks that involve tracking
a target behavior while optimizing a cost and satisfying safety requirements
expressed as box constraints. By leveraging a probabilistic framework, we
formulate a data-driven optimal control problem and, for this problem, we: (i)
prove convexity under a suitable condition; (ii) find the optimal solution;
(iii) turn our results into an algorithm, using it to design an intelligent
parking system for connected vehicles. Validations are performed both in-
silico and in-vivo, with real cars. As such, the main purpose of this paper is
twofold: (i) we introduce, and rigorously characterize our algorithm; (ii)
propose a stand-alone implementation of our results, suitable for in-vivo
experiments on real cars. In-vivo validations were performed via an hardware-
in-the-loop platform allowing to embed real cars/drivers in the experiments.
Using the platform, we deploy our algorithm on a real vehicle showing its
suitability for in-car operation. All experiments confirm the effectiveness of
our approach (documented code/data for our simulations at
https://tinyurl.com/3ep4pknh).

While our results are inspired by the ones in [3, 4], our paper extends these
in several ways. First, the results in [3, 4] cannot consider box constraints
and hence cannot tackle the control problem of this paper. Second, even when
there are no box constraints, the results in [3, 4] only solve an approximate
version of the problem considered here. That is, the results from [3, 4] only
find an approximate, non-optimal, solution of the problem considered here. As
we shall see, this means that the solutions from [3, 4] cannot get a better
cost than the one obtained with the results of this paper. Third, the
algorithm obtained from the results in this paper is deployed, and validated,
on a real car and this was not done in [3, 4]. The _in-vivo_ implementation is
novel.

#### Notation

sets are in calligraphic and vectors in bold. Given the measurable space
$(\mathcal{X},\mathcal{F}_{x})$, with $\mathcal{X}\subseteq\mathbb{R}^{d}$
($\mathcal{X}\subseteq\mathbb{Z}^{d}$) and $\mathcal{F}_{x}$ being a
$\sigma$-algebra on $\mathcal{X}$, a random variable on
$(\mathcal{X},\mathcal{F}_{x})$ is denoted by $\mathbf{X}$ and its realization
by $\mathbf{x}$. The probability density (resp. mass) function or pdf (pmf) of
a continuous (discrete) $\mathbf{X}$ is denoted by $p(\mathbf{x})$. The convex
subset of such probability functions (pfs) is $\mathcal{D}$. The expectation
of a function $\mathbf{h}(\cdot)$ of the continuous variable $\mathbf{X}$ is
$\mathbb{E}_{{p}}[\mathbf{h}(\mathbf{X})]:=\int\mathbf{h}(\mathbf{x})p(\mathbf{x})d\mathbf{x}$,
where the integral (in the sense of Lebesgue) is over the support of
$p(\mathbf{x})$, which we denote by $\mathcal{S}(p)$. The joint pf of
$\mathbf{X}_{1}$, $\mathbf{X}_{2}$ is $p(\mathbf{x}_{1},\mathbf{x}_{2})$ and
the conditional pf of $\mathbf{X}_{1}$ given $\mathbf{X}_{2}$ is
$p\left(\mathbf{x}_{1}\mid\mathbf{x}_{2}\right)$. Countable sets are denoted
by $\\{w_{k}\\}_{k_{1}:k_{n}}$, where $w_{k}$ is the generic element of the
set and $k_{1}:k_{n}$ is the closed set of consecutive integers between
$k_{1}$ and $k_{n}$. The KL divergence between $p(\mathbf{x})$ and
$q(\mathbf{x})$, where $p$ is absolutely continuous w.r.t. $q$, is
$\mathcal{D}_{\text{KL}}\left(p\mid\mid
q\right):=\int_{{\color[rgb]{0,0,0}\mathcal{S}(p)}}p\;\ln\left({p}/{q}\right)\,d\mathbf{x}$:
it is non-negative and $0$ if and only if $p(\mathbf{x})=q(\mathbf{x})$. In
the expressions for the expectation and KL divergence, the integral is
replaced by the sum if the variables are discrete. Finally: (i) we let
$\mathds{1}_{\mathcal{A}}(\mathbf{x})$ denote the indicator function being
equal to $1$ if $\mathbf{x}\in\mathcal{A}\subseteq\mathcal{X}$ and $0$
otherwise; (ii) set exclusion is instead denoted by $\setminus$.

##  II The Setup

The agent seeks to craft its behavior by combining a number of sources to
fulfill a task that involves tracking a target/desired behavior while
maximizing an agent-specific reward over the time horizon $\mathcal{T}:=0:T$,
$T>0$. The agent’s state at time step $k\in\mathcal{T}$ is
$\mathbf{x}_{k}\in\mathcal{X}$ and the target behavior that the agent seeks to
track is $p_{0:T}:=p_{0}(\mathbf{x}_{0})\prod_{k\in
1:T}p(\mathbf{x}_{k}\mid\mathbf{x}_{k-1})$. As in [3, 4], we design the
behavior of the agent by designing its joint pf
$\pi_{0:T}:=\pi(\mathbf{x}_{0},\ldots,\mathbf{x}_{T})$ and we have:

$\pi_{0:T}=\pi_{0}(\mathbf{x}_{0})\prod_{k\in
1:T}\pi(\mathbf{x}_{k}\mid\mathbf{x}_{k-1}).$ (1)

That is, the behavior of the agent can be designed via the pfs
$\pi(\mathbf{x}_{k}\mid\mathbf{x}_{k-1})$, i.e., the transition probabilities.
To do so, the agent has access to $S$ sources and we denote by
$\pi^{(i)}(\mathbf{x}_{k}\mid\mathbf{x}_{k-1})$, with support
$\mathcal{S}(\pi)\subseteq\mathcal{X}$, the behavior made available by source
$i$, $i\in\mathcal{S}:=1:S$, at $k-1$. We also let $r_{k}(\mathbf{x}_{k})$ be
the agent’s reward for being in state $\mathbf{x}_{k}$ at $k$.

##  III Formulation of the Control Problem

Let $\alpha_{k}^{(1)},\ldots,\alpha_{k}^{(S)}$ be weights and
$\boldsymbol{\alpha}_{k}$ be their stack. Then, the control problem we
consider can be formalized as:

###### Problem 1

find the sequence ${\left\\{\boldsymbol{\alpha}_{k}^{\ast}\right\\}_{1:T}}$
solving

$\displaystyle\underset{\left\\{\boldsymbol{\alpha}_{k}\right\\}_{1:T}}{\text{min}}$
$\displaystyle\mathcal{D}_{\text{KL}}\left(\pi_{0:T}\mid\mid
p_{0:T}\right)-\sum_{k=1}^{T}\mathbb{E}_{\pi(\mathbf{x}_{k-1})}\left[\tilde{r}_{k}(\mathbf{X}_{k-1})\right]$
s.t. $\displaystyle\
\mathbb{E}_{\pi(\mathbf{x}_{k}\mid\mathbf{x}_{k-1})}\left[\mathds{1}_{\mathcal{X}_{k}}(\mathbf{x}_{k})\right]\geq
1-\epsilon_{k},\ \ \ \forall k,$ $\displaystyle\
\pi(\mathbf{x}_{k}\mid\mathbf{x}_{k-1})=\sum_{i\in\mathcal{S}}\alpha_{k}^{(i)}\pi^{(i)}(\mathbf{x}_{k}\mid\mathbf{x}_{k-1}),\
\ \ \forall k,$ $\displaystyle\sum_{i\in\mathcal{S}}\alpha_{k}^{(i)}=1,\ \
\alpha_{k}^{(i)}\in[0,1],\ \ \ \forall k.$

In Problem 1,
$\tilde{r}_{k}(\mathbf{x}_{k-1}):=\mathbb{E}_{\pi(\mathbf{x}_{k}\mid\mathbf{x}_{k-1})}\left[r_{k}(\mathbf{X}_{k})\right]$
and we note that
$\mathbb{E}_{\pi(\mathbf{x}_{k-1})}\left[\tilde{r}_{k}(\mathbf{X}_{k-1})\right]=\mathbb{E}_{\pi(\mathbf{x}_{k})}\left[r_{k}(\mathbf{X}_{k})\right]$
is the expected reward for the agent when the behavior in (1) is followed. The
problem is a finite-horizon optimal control problem with the
$\boldsymbol{\alpha}_{k}$’s as decision variables. As we shall see, these are
generated as feedback from the agent state (Section IV). We say that
$\left\\{\pi^{\ast}\left(\mathbf{x}_{k}\mid\mathbf{x}_{k-1}\right)\right\\}_{1:T}$,
with
$\pi^{\ast}\left(\mathbf{x}_{k}\mid\mathbf{x}_{k-1}\right)=\sum_{i\in\mathcal{S}}\alpha_{k}^{(i),\ast}\pi^{(i)}(\mathbf{x}_{k}\mid\mathbf{x}_{k-1})$,
is the optimal behavior for the agent, obtained by linearly combining the
sources via the $\boldsymbol{\alpha}_{k}^{\ast}$’s. In the problem, the cost
formalizes the fact that the agent seeks to maximize its reward, while
tracking (in the KL divergence sense) the target behavior. Minimizing the KL
term amounts at minimizing the discrepancy between $\pi_{0:T}$ and $p_{0:T}$.
This term can also be thought as a divergence regularizer and, when $p_{0:T}$
is uniform, it becomes an entropic regularizer. The second and third
constraints formalize the fact that, at each $k$,
$\pi^{\ast}\left(\mathbf{x}_{k}\mid\mathbf{x}_{k-1}\right)\in\mathcal{D}$ and
that this is a convex combination of the pfs from the sources. The first
constraint is a box constraint and models the fact that the probability that
the agent behavior is, at each $k$, inside some (e.g., safety) measurable set
$\mathcal{X}_{k}\subseteq\mathcal{X}$ is greater than some
${\color[rgb]{0,0,0}\epsilon}_{k}\geq 0$. We now make the following

###### Assumption 1

the optimal cost of Problem 1 is bounded.

###### Remark 1

the cost in Problem 1 can be recast as
$\mathcal{D}_{\text{KL}}\left(\pi_{0:T}\mid\mid\tilde{p}_{0:T}\right)$, where
$\tilde{p}_{0:T}\propto
p_{0:T}\exp{\left(\sum_{k=1}^{T}r_{k}(\mathbf{x}_{k})\right)}$. This means
that Assumption 1 is satisfied whenever there exists some $\tilde{\pi}_{0:T}$
that is feasible for Problem 1 and that is absolutely continuous w.r.t.
$\tilde{p}_{0:T}$. See also Remark 3.

Algorithm 1 Pseudo-code

1:Input: time horizon $T$, target behavior $p_{0:T}$, reward $r_{k}(\cdot)$,
sources $\pi^{(i)}(\mathbf{x}_{k}\mid\mathbf{x}_{k-1})$, box constraints
(optional)

2:Output: optimal agent behavior
$\left\\{\pi^{\ast}\left(\mathbf{x}_{k}\mid\mathbf{x}_{k-1}\right)\right\\}_{1:T}$

3:$\hat{r}_{T}(\mathbf{x}_{T})\leftarrow 0$

4:for $k=T:1$ do

5:$\bar{r}_{k}(\mathbf{x}_{k})\leftarrow
r_{k}(\mathbf{x}_{k})-\hat{r}_{k}(\mathbf{x}_{k})$

6:$\boldsymbol{\alpha}^{\ast}_{k}(\mathbf{x}_{k-1})\leftarrow\text{minimizer
of the problem in \eqref{eqn:subproblem_alg}}$;

7:$\pi^{\ast}(\mathbf{x}_{k}\mid\mathbf{x}_{k-1})\leftarrow\sum_{i\in\mathcal{S}}\boldsymbol{\alpha}^{(i),\ast}_{k}(\mathbf{x}_{k-1})\pi^{(i)}(\mathbf{x}_{k}\mid\mathbf{x}_{k-1})$

8:$\hat{r}_{k-1}(\mathbf{x}_{k-1})\leftarrow
c_{k}(\pi^{\ast}(\mathbf{x}_{k}\mid\mathbf{x}_{k-1}))$

9:end for

##  IV Main Results

We propose an algorithm to tackle Problem 1. The algorithm takes as input $T$,
the target behavior, the reward, the behaviors of the sources and the box
constraints of Problem 1 (if any). Given the inputs, it returns the optimal
behavior for the agent. The key steps of the algorithm are given as pseudo-
code in Algorithm 1. An agent that follows Algorithm 1 computes
$\left\\{\boldsymbol{\alpha}_{k}\right\\}_{1:N}$ via backward recursion (lines
$4-9$). At each $k$, the $\boldsymbol{\alpha}_{k}$’s are obtained as the
minimizers of

$\displaystyle\underset{\boldsymbol{\alpha}_{k}}{\min}$ $\displaystyle\
c_{k}(\pi(\mathbf{x}_{k}\mid\mathbf{x}_{k-1}))$ (2) s.t. $\displaystyle\
\mathbb{E}_{\pi(\mathbf{x}_{k}\mid\mathbf{x}_{k-1})}\left[\mathds{1}_{\mathcal{X}_{k}}(\mathbf{x}_{k})\right]\geq
1-\epsilon_{k}$ $\displaystyle\
\pi(\mathbf{x}_{k}\mid\mathbf{x}_{k-1})=\sum_{i\in\mathcal{S}}\alpha_{k}^{(i)}\pi^{(i)}(\mathbf{x}_{k}\mid\mathbf{x}_{k-1})$
$\displaystyle\sum_{i\in\mathcal{S}}\alpha_{k}^{(i)}=1,\ \
\alpha_{k}^{(i)}\in[0,1],$

where

$\begin{split}c_{k}(\pi(\mathbf{x}_{k}\mid\mathbf{x}_{k-1}))&:=\mathcal{D}_{\text{KL}}\left(\pi(\mathbf{x}_{k}\mid\mathbf{x}_{k-1})\mid\mid
p(\mathbf{x}_{k}\mid\mathbf{x}_{k-1})\right)\\\
&-\mathbb{E}_{\pi(\mathbf{x}_{k}\mid\mathbf{x}_{k-1})}\left[\bar{r}_{k}(\mathbf{X}_{k})\right],\end{split}$
(3)

with $\bar{r}_{k}(\cdot)$ iteratively built within the recursion (lines $5$,
$8$). The weights are used (line $7$) to compute
$\pi^{\ast}(\mathbf{x}_{k}\mid\mathbf{x}_{k-1})$.

###### Remark 2

results are stated for continuous variables (proofs for discrete variables
omitted for brevity). Note that integrals/summations in the cost are over
$\mathcal{S}(\pi)$.

###### Remark 3

following Remark 1, the optimal cost of the problem in (2) is bounded if there
exists some feasible $\tilde{\pi}(\mathbf{x}_{k}\mid\mathbf{x}_{k-1})$ that is
absolutely continuous w.r.t.
$\tilde{p}(\mathbf{x}_{k}\mid\mathbf{x}_{k-1})\propto
p(\mathbf{x}_{k}\mid\mathbf{x}_{k-1})\exp\left(\bar{r}_{k}(\mathbf{x}_{k})\right)$.
From the design viewpoint, this can satisfied if it holds for at least one
$\pi^{(i)}(\mathbf{x}_{k}\mid\mathbf{x}_{k-1})$.

Finally, we make the following

###### Assumption 2

$\forall i\in\mathcal{S}$ and $\forall\mathbf{x}_{k-1}\in\mathcal{X}$, there
exist some constants, say $m$ and $M$, with $0<m\leq M<+\infty$, such that
$m\leq\pi^{(i)}(\mathbf{x}_{k}\mid\mathbf{x}_{k-1})\leq M$,
$\forall\mathbf{x}_{k}\in\mathcal{S}(\pi)$.

###### Remark 4

Assumption 1 is satisfied for e.g., Gaussian distributions. As we shall see
(Section V) the assumption can be fulfilled by injecting noise in the data.

###  IV-A Properties of Algorithm 1

Before characterizing convexity of the problems recursively solved in
Algorithm 1 and optimality, we give a condition ensuring feasibility of the
problem in (2).

###### Lemma 1

the problem in (2) is feasible if and only if there exists at least one
source, say $\pi^{(j)}(\mathbf{x}_{k}\mid\mathbf{x}_{k-1})$, such that
$\mathbb{E}_{\pi^{(j)}(\mathbf{x}_{k}\mid\mathbf{x}_{k-1})}\left[\mathds{1}_{\mathcal{X}_{k}}(\mathbf{x}_{k})\right]\geq
1-\epsilon_{k}$.

###### Proof:

the if part clearly holds. For the only if part we prove that if problem (2)
is infeasible, then
$\max_{i}\mathbb{E}_{\pi^{(i)}(\mathbf{x}_{k}\mid\mathbf{x}_{k-1})}\left[\mathds{1}_{\mathcal{X}_{k}}(\mathbf{x}_{k})\right]<1-\epsilon_{k}$.
In fact, if the problem is infeasible, then for all $\boldsymbol{\alpha}_{k}$
such that $\sum_{i\in\mathcal{S}}\alpha_{k}^{(i)}=1$ and
$\alpha_{k}^{(i)}\in[0,1]$ it must hold that

$\int_{\mathcal{X}_{k}}\sum_{i\in\mathcal{S}}\alpha_{k}^{(i)}\pi^{(i)}(\mathbf{x}_{k}\mid\mathbf{x}_{k-1})d\mathbf{x}_{k}<1-\epsilon_{k}.$

Note that this must also hold for all $\boldsymbol{\alpha}_{k}$ such that
$\sum_{i\in\mathcal{S}}\alpha_{k}^{(i)}=1$ and $\alpha_{k}^{(i)}\in\\{0,1\\}$,
as these are contained in the set of real-valued $\boldsymbol{\alpha}_{k}$’s.
We conclude the proof by noticing that, if $\boldsymbol{\alpha}_{k}$ is such
that $\alpha_{k}^{(i)}=0\forall i\neq j$ and $\alpha_{k}^{(j)}=1$, then

$\displaystyle\int_{\mathcal{X}_{k}}\sum_{i\in\mathcal{S}}\alpha_{k}^{(i)}\pi^{(i)}(\mathbf{x}_{k}\mid\mathbf{x}_{k-1})d\mathbf{x}_{k}$
$\displaystyle=\int_{\mathcal{X}_{k}}\pi^{(j)}(\mathbf{x}_{k}\mid\mathbf{x}_{k-1})d\mathbf{x}_{k}$
$\displaystyle=E_{\pi^{(j)}(\mathbf{x}_{k}\mid\mathbf{x}_{k-1})}\left[\mathds{1}_{\mathcal{X}_{k}}(\mathbf{x}_{k})\right].$

It then follows that, $\forall j$,

$E_{\pi^{(j)}(\mathbf{x}_{k}\mid\mathbf{x}_{k-1})}\left[\mathds{1}_{\mathcal{X}_{k}}(\mathbf{x}_{k})\right]<1-\epsilon_{k}.$

∎

#### Convexity

we are now ready to prove the following

###### Proposition 1

let Assumption 2 hold. Then, the problem in (2) is convex.

###### Proof:

Clearly, the second and third constraint in the problem are convex. For the
first constraint, we get

$\displaystyle\mathbb{E}_{\pi(\mathbf{x}_{k}\mid\mathbf{x}_{k-1})}\left[\mathds{1}_{\mathcal{X}_{k}}(\mathbf{x}_{k})\right]$
$\displaystyle=\int\pi(\mathbf{x}_{k}\mid\mathbf{x}_{k-1})\mathds{1}_{\mathcal{X}_{k}}(\mathbf{x}_{k})d\mathbf{x}_{k}$
$\displaystyle=\int_{\mathcal{X}_{k}}\sum_{i\in\mathcal{S}}\alpha_{k}^{(i)}\pi^{(i)}(\mathbf{x}_{k}\mid\mathbf{x}_{k-1})d\mathbf{x}_{k}$
$\displaystyle=\sum_{i\in\mathcal{S}}\alpha_{k}^{(i)}\int_{\mathcal{X}_{k}}\pi^{(i)}(\mathbf{x}_{k}\mid\mathbf{x}_{k-1})d\mathbf{x}_{k},$

which is therefore convex in the decision variables.
Now, we show that the cost is also convex in these variables and we do so by
explicitly computing, for each $\mathbf{x}_{k-1}$, its Hessian, say
$\mathbf{H}(\mathbf{x}_{k-1})$. Specifically, after embedding the second
constraint of the problem in (2) in the cost and differentiating with respect
to the decision variables we get, for each $j\in\mathcal{S}$:

$\begin{split}&\frac{\partial
c_{k}}{\partial\alpha_{k}^{(j)}}:=\frac{\partial}{\partial\alpha_{k}^{(j)}}\int_{{\color[rgb]{0,0,0}\mathcal{S}(\pi)}}\sum_{i\in\mathcal{S}}\alpha_{k}^{(i)}\pi^{(i)}(\mathbf{x}_{k}\mid\mathbf{x}_{k-1})\left(\log\left(\frac{\sum_{i\in\mathcal{S}}\alpha_{k}^{(i)}\pi^{(i)}(\mathbf{x}_{k}\mid\mathbf{x}_{k-1})}{p(\mathbf{x}_{k}\mid\mathbf{x}_{k-1})}\right)-\bar{r}_{k}(\mathbf{x}_{k})\right)d\mathbf{x}_{k}\\\
&=\int_{{\color[rgb]{0,0,0}\mathcal{S}(\pi)}}\frac{\partial}{\partial\alpha_{k}^{(j)}}\sum_{i\in\mathcal{S}}\alpha_{k}^{(i)}\pi^{(i)}(\mathbf{x}_{k}\mid\mathbf{x}_{k-1})\left(\log\left(\frac{\sum_{i\in\mathcal{S}}\alpha_{k}^{(i)}\pi^{(i)}(\mathbf{x}_{k}\mid\mathbf{x}_{k-1})}{p(\mathbf{x}_{k}\mid\mathbf{x}_{k-1})}\right)-\bar{r}_{k}(\mathbf{x}_{k})\right)d\mathbf{x}_{k}\\\
&=\int_{{\color[rgb]{0,0,0}\mathcal{S}(\pi)}}\pi^{(j)}(\mathbf{x}_{k}\mid\mathbf{x}_{k-1})\Bigg{(}\Bigg{.}\log\left(\sum_{i\in\mathcal{S}}\alpha_{k}^{(i)}\pi^{(i)}(\mathbf{x}_{k}\mid\mathbf{x}_{k-1})\right)-\log\left(p(\mathbf{x}_{k}\mid\mathbf{x}_{k-1})\right)-\bar{r}_{k}(\mathbf{x}_{k})+1\Bigg{.}\Bigg{)}d\mathbf{x}_{k}.\\\
\end{split}$

The above chain of identities was obtained by swapping integration and
differentiation, leveraging the fact that the cost is smooth in the decision
variables. Similarly, we get

$\begin{split}\frac{\partial^{2}c_{k}}{\partial\alpha_{k}^{(j)^{2}}}&=\frac{\partial}{\partial\alpha_{k}^{(j)}}\int_{{\color[rgb]{0,0,0}\mathcal{S}(\pi)}}\pi^{(j)}(\mathbf{x}_{k}\mid\mathbf{x}_{k-1})\Bigg{(}\Bigg{.}\log\left(\sum_{i\in\mathcal{S}}\alpha_{k}^{(i)}\pi^{(i)}(\mathbf{x}_{k}\mid\mathbf{x}_{k-1})\right)-\log\left(p(\mathbf{x}_{k}\mid\mathbf{x}_{k-1})\right)-\bar{r}_{k}(\mathbf{x}_{k})+1\Bigg{.}\Bigg{)}d\mathbf{x}_{k}\\\
&={\color[rgb]{0,0,0}\int_{\mathcal{S}(\pi)}\frac{\pi^{(j)}(\mathbf{x}_{k}\mid\mathbf{x}_{k-1})^{2}}{\sum_{i\in\mathcal{S}}\alpha_{k}^{(i)}\pi^{(i)}(\mathbf{x}_{k}\mid\mathbf{x}_{k-1})}d\mathbf{x}_{k}}\\\
&=:h_{{\color[rgb]{0,0,0}jj}}(\mathbf{x}_{k-1}),\end{split}$

and, for each $m\neq j$, $m\in\mathcal{S}$,

$\frac{\partial^{2}c_{k}}{\partial\alpha_{k}^{(j)}\partial\alpha_{k}^{(m)}}=\int_{{\color[rgb]{0,0,0}\mathcal{S}(\pi)}}\frac{\pi^{(j)}(\mathbf{x}_{k}\mid\mathbf{x}_{k-1})\pi^{(m)}(\mathbf{x}_{k}\mid\mathbf{x}_{k-1})}{\sum_{i\in\mathcal{S}}\alpha_{k}^{(i)}\pi^{(i)}(\mathbf{x}_{k}\mid\mathbf{x}_{k-1})}d\mathbf{x}_{k}=:h_{{\color[rgb]{0,0,0}j}m}(\mathbf{x}_{k-1}).$

Also, following Assumption 2, $\forall j,m\in\mathcal{S}$ we have that

$\int_{{\color[rgb]{0,0,0}\mathcal{S}(\pi)}}\left|\frac{\pi^{(j)}(\mathbf{x}_{k}\mid\mathbf{x}_{k-1})\pi^{(m)}(\mathbf{x}_{k}\mid\mathbf{x}_{k-1})}{\sum_{i\in\mathcal{S}}\alpha_{k}^{(i)}\pi^{(i)}(\mathbf{x}_{k}\mid\mathbf{x}_{k-1})}\right|d\mathbf{x}_{k}\leq\frac{M}{m}$

where we used the third constraint. That is, the above integrals are well
defined and thus we can conclude the proof by computing
$\mathbf{v}^{T}\mathbf{H}(\mathbf{x}_{k-1})\mathbf{v}$ for some non-zero
$\mathbf{v}\in\mathbb{R}^{S}$:

$\displaystyle\mathbf{v}^{T}\mathbf{H}(\mathbf{x}_{k-1})\mathbf{v}$
$\displaystyle=\sum_{j,m}v_{j}v_{m}h_{jm}(\mathbf{x}_{k-1})$
$\displaystyle=\sum_{j,m}v_{j}v_{m}\int_{{\color[rgb]{0,0,0}\mathcal{S}(\pi)}}a_{jm}(\mathbf{x}_{k},\mathbf{x}_{k-1})d\mathbf{x}_{k}$
$\displaystyle=\int_{{\color[rgb]{0,0,0}\mathcal{S}(\pi)}}\sum_{j,m}v_{j}v_{m}a_{jm}(\mathbf{x}_{k},\mathbf{x}_{k-1})d\mathbf{x}_{k},$

where the $a_{jm}$’s are the elements of the matrix

$A(\mathbf{x}_{k},\mathbf{x}_{k-1}):=\bar{\pi}(\mathbf{x}_{k},\mathbf{x}_{k-1})\left[\begin{array}[]{*{20}c}\pi^{(1)}(\mathbf{x}_{k}\mid\mathbf{x}_{k-1})\\\
\vdots\\\
\pi^{(S)}(\mathbf{x}_{k}\mid\mathbf{x}_{k-1})\end{array}\right]\cdot\\\
\cdot\left[\begin{array}[]{*{20}c}\pi^{(1)}(\mathbf{x}_{k}\mid\mathbf{x}_{k-1})\ldots\pi^{(S)}(\mathbf{x}_{k}\mid\mathbf{x}_{k-1})\end{array}\right],$

with
$\bar{\pi}(\mathbf{x}_{k},\mathbf{x}_{k-1}):={1}/\left({\sum_{i\in\mathcal{S}}\alpha_{k}^{(i)}\pi^{(i)}(\mathbf{x}_{k}\mid\mathbf{x}_{k-1})}\right)$.
The above expression is indeed positive semi-definite for each
$\mathbf{x}_{k}$, $\mathbf{x}_{k-1}$ and we can draw the desired conclusion. ∎

#### Optimality

we can now prove the following

###### Proposition 2

let Assumption 2 and Assumption 1 hold. Then, Algorithm 1 gives an optimal
solution for Problem 1.

###### Proof:

the chain rule for the KL divergence and the linearity of expectation imply
that the cost can be written as

$\mathcal{D}_{\text{KL}}\left(\pi_{0:T-1}\mid\mid
p_{0:T-1}\right)-\sum_{k=1}^{T-1}\mathbb{E}_{\pi(\mathbf{x}_{k-1})}\left[\tilde{r}_{k}(\mathbf{X}_{k-1})\right]+\mathbb{E}_{\pi(\mathbf{x}_{T-1})}\left[c_{T}(\pi(\mathbf{x}_{T}\mid\mathbf{x}_{T-1}))\right],$
(4)

where $c_{T}(\pi(\mathbf{x}_{T}\mid\mathbf{x}_{T-1}))$ is defined as in (3)
with $\bar{r}_{T}(\mathbf{x}_{T})$ given by Algorithm 1 – see lines $3$ and
$5$ and note that, at time step $T$,
$\bar{r}_{T}(\mathbf{x}_{T})={r}_{T}(\mathbf{x}_{T})$. To obtain the above
expression, the fact that $c_{T}(\pi(\mathbf{x}_{T}\mid\mathbf{x}_{T-1}))$
only depends on $\mathbf{x}_{T-1}$ was also used. Hence, Problem 1 can be
split into the sum of two sub-problems: a first problem over $k\in 0:T-1$ and
the second for $k=T$. For this last time step, the problem can be solved
independently on the others and is given by:

$\displaystyle\underset{\boldsymbol{\alpha}_{T}}{\text{ min }}$
$\displaystyle\
\mathbb{E}_{\pi(\mathbf{x}_{T-1})}[c_{T}(\pi(\mathbf{x}_{T}\mid\mathbf{x}_{T-1}))]$
(5) s.t. $\displaystyle\
\mathbb{E}_{\pi(\mathbf{x}_{T}\mid\mathbf{x}_{T-1})}\left[\mathds{1}_{\mathcal{X}_{T}}(\mathbf{x}_{T})\right]\geq
1-\epsilon_{T},$ $\displaystyle\
\pi(\mathbf{x}_{T}\mid\mathbf{x}_{T-1})=\sum_{i\in\mathcal{S}}\alpha_{T}^{(i)}\pi^{(i)}(\mathbf{x}_{T}\mid\mathbf{x}_{T-1}),$
$\displaystyle\sum_{i\in\mathcal{S}}\alpha_{T}^{(i)}=1,\ \
\alpha_{T}^{(i)}\in[0,1].$

Using linearity of the expectation and the fact that the decision variable is
independent on $\pi(\mathbf{x}_{T-1})$, we have that the minimizer of the
problem in (5) is the same as the problem in (2) with $k=T$. Following
Proposition 1, such a problem is convex and we denote its optimal solution as
$\boldsymbol{\alpha}^{\ast}_{T}(\mathbf{x}_{T-1})$ – see line $6$ of Algorithm
1 – and the optimal cost of the problem, which is bounded by Assumption 1, is
$c_{T}(\pi^{\ast}(\mathbf{x}_{T}\mid\mathbf{x}_{T-1}))$, where:

$\pi^{\ast}(\mathbf{x}_{T}\mid\mathbf{x}_{T-1})=\sum_{i\in\mathcal{S}}\boldsymbol{\alpha}^{(i),\ast}_{T}(\mathbf{x}_{T-1})\pi^{(i)}(\mathbf{x}_{T}\mid\mathbf{x}_{T-1}){\color[rgb]{0,0,0}.}$

This gives $\hat{r}_{T-1}(\mathbf{x}_{T-1})$ in Algorithm 1 (lines $7-8$),
thus yielding the steps for the backward recursion of the Algorithm 1 at time
step $T$. Now, the minimum value of the problem in (5) is given by
$\mathbb{E}_{\pi(\mathbf{x}_{T-1})}\left[\hat{r}_{T-1}(\mathbf{X}_{T-1})\right]$.
Hence, the cost of Problem 1 becomes

$\mathcal{D}_{\text{KL}}\left(\pi_{0:T-1}\mid\mid
p_{0:T-1}\right)-\sum_{k=1}^{T-1}\mathbb{E}_{\pi(\mathbf{x}_{k-1})}\left[\tilde{r}_{k}(\mathbf{X}_{k-1})\right]+\mathbb{E}_{\pi(\mathbf{x}_{T-1})}\left[\hat{r}_{T-1}(\mathbf{X}_{T-1})\right].$

Then, following the same reasoning used to obtain (4) and by noticing that

$\mathbb{E}_{\pi(\mathbf{x}_{T-1})}\left[\hat{r}_{T-1}(\mathbf{X}_{T-1})\right]=\mathbb{E}_{\pi(\mathbf{x}_{T-2})}\left[\mathbb{E}_{\pi(\mathbf{x}_{T-1}\mid\mathbf{x}_{T-2})}\left[\hat{r}_{T-1}(\mathbf{X}_{T-1})\right]\right],$

we get:

$\mathcal{D}_{\text{KL}}\left(\pi_{0:T-2}\mid\mid
p_{0:T-2}\right)-\sum_{k=1}^{T-2}\mathbb{E}_{\pi(\mathbf{x}_{k-1})}\left[\tilde{r}_{k}(\mathbf{X}_{k-1})\right]+\mathbb{E}_{\pi(\mathbf{x}_{T-2})}\left[c_{T-1}(\pi(\mathbf{x}_{T-1}\mid\mathbf{x}_{T-2}))\right],$

where $c_{T-1}(\pi(\mathbf{x}_{T-1}\mid\mathbf{x}_{T-2}))$ is again given in
(3) with $\bar{r}_{T-1}(\mathbf{x}_{T-1})$ again defined as in Algorithm 1. By
iterating the arguments above, we find that at each time step Problem 1 can
always be split as the sum of two sub-problems, where the last sub-problem can
be solved independently on the previous ones. Moreover, the minimizer of this
last sub-problem is always the solution of a problem of the form

$\displaystyle\underset{\boldsymbol{\alpha}_{k}}{\text{min}}$
$\displaystyle\mathcal{D}_{\text{KL}}\left(\pi(\mathbf{x}_{k}\mid\mathbf{x}_{k-1})\mid\mid
p(\mathbf{x}_{k}\mid\mathbf{x}_{k-1})\right)-\mathbb{E}_{\pi(\mathbf{x}_{k}\mid\mathbf{x}_{k-1})}\left[\bar{r}_{k}(\mathbf{X}_{k})\right]$
(6) s.t. $\displaystyle\
\mathbb{E}_{\pi(\mathbf{x}_{k}\mid\mathbf{x}_{k-1})}\left[\mathds{1}_{\mathcal{X}_{k}}(\mathbf{x}_{k})\right]\geq
1-\epsilon_{k},$ $\displaystyle\
\pi(\mathbf{x}_{k}\mid\mathbf{x}_{k-1})=\sum_{i\in\mathcal{S}}\alpha_{k}^{(i)}\pi^{(i)}(\mathbf{x}_{k}\mid\mathbf{x}_{k-1}),$
$\displaystyle\sum_{i\in\mathcal{S}}\alpha_{k}^{(i)}=1,\ \
\alpha_{k}^{(i)}\in[0,1],$

where
$\bar{r}_{k}(\mathbf{x}_{k}):=r_{k}(\mathbf{x}_{k})-\hat{r}_{k}(\mathbf{x}_{k})$,
with
$\hat{r}_{k}(\mathbf{x}_{k}):=c_{k+1}({\pi}^{\ast}(\mathbf{x}_{k+1}\mid\mathbf{x}_{k}))$
and
$\pi^{\ast}(\mathbf{x}_{k}\mid\mathbf{x}_{k-1})=\sum_{i\in\mathcal{S}}\boldsymbol{\alpha}^{(i),\ast}_{k}(\mathbf{x}_{k-1})\pi^{(i)}(\mathbf{x}_{k}\mid\mathbf{x}_{k-1})$.
This yields the desired conclusions. ∎

###### Remark 5

the results from [3, 4] solve an approximate version of Problem 1 when there
are no box constraints. Hence, even in this special case, the results from [3,
4] do not lead to the optimal solution found with Algorithm 1. Specifically,
the approximate solution from [3, 4] corresponds to the optimal solution of a
problem with additional binary constraints on the decision variables. As a
result, the algorithm from [3, 4] searches solutions in a feasibility domain
that is contained in the feasibility domain of Problem 1. Hence, the solutions
found in [3, 4] cannot achieve a better cost than the ones obtained via
Algorithm 1.

##  V Designing an Intelligent Parking System

We now use Algorithm 1 to design an intelligent parking system for connected
cars and validate the results via in-silico and in-vivo experiments. For the
latter set of experiments, Algorithm 1 is deployed on a real car and
validation is performed via an hardware-in-the-loop (HiL) platform inspired
from [19]. Before reporting the results, we describe the validation scenarios
and the experimental set-up. Code, maps and parameters with instructions to
replicate the simulations are given at https://tinyurl.com/3ep4pknh.

Figure 1: campus map. The magnified areas show the obstructed road link (in
blue) and links used within the validations. Colors online.

###  V-A Validation Scenarios and Experimental Set-up

We consider the problem of routing vehicles in a given geographic area to find
parking. In all experiments we consider a morning rush scenario at the
University of Salerno campus (see Figure 1). Specifically, cars arrive to the
campus through a highway exit and, from here, users seek to park in one of the
parking locations: Biblioteca and Terminal.

In this context, vehicles are agents equipped with Algorithm 1. The set of
road links within the campus is $\mathcal{X}$ and time steps are associated to
the time instants when the vehicle changes road link. The state of the agent,
$x_{k}$, is the road link occupied by the vehicle at time step $k$. Given this
set-up, at each $k$ Algorithm 1 outputs the turning probability for the car
given the current car link, $\pi^{\ast}({x}_{k}\mid{x}_{k-1})$. The next
direction for the car is then obtained by sampling from this pf. Agents have
access to a set of sources, each providing different routes. As discussed in
[20], sources might be third parties navigation services, routes collected
from other cars/users participating to some sharing service. Agents wish to
track their target/desired behavior (driving them to the preferred parking –
Terminal in our experiments) and the reward depends on the actual road
conditions within the campus. Links adjacent to a parking lot are assigned a
constant reward of: (i) $3.8$ if the parking has available spaces; (ii) 0 when
it becomes full. Unparked cars already on full parking lots are assigned a
target behavior leading them to another parking. In normal road conditions,
the reward for the other links is $0$ and becomes $-20$ when there is an
obstruction. In our experiments, the reward was selected heuristically so that
it would be: (i) sufficiently penalizing for links affected by obstruction;
(ii) encouraging cars to drive towards parking lots with available spaces. In
the first scenario (Scenario $1$) there are no box constraints: this is done
to benchmark Algorithm 1 with [3, 4]. To this aim, we use the campus map from
[20] in which [3, 4] were thoroughly validated via simulations. Then, we show
that by introducing box constraints Algorithm 1 can effectively regulate
access of vehicles through selected road links. This is Scenario $2$ and we
considered three situations: (A) the road towards the Biblioteca parking lot
is forbidden. To account for this, we set in Algorithm 1
$\mathcal{X}_{k}=\mathcal{X}\setminus{{l_{2}}}$, where the link ${l_{2}}$ is
shown in Figure 1, and $\epsilon_{k}=0.027$; (B) the set $\mathcal{X}_{k}$ as
before but now $\epsilon_{k}=0.5$; (C) the road towards the Terminal parking
lot is forbidden and in this case we set
$\mathcal{X}_{k}=\mathcal{X}\setminus{l_{1}}$, $\epsilon_{k}=0.027$ (see
Figure 1 for link ${l_{1}}$). For this last scenario, Algorithm 1 is validated
both in-silico and in-vivo. Next, we describe the simulation and HiL
experimental set-ups.

#### Simulation set-up

simulations were performed in SUMO [21]; see also [20] for a description of
the pipeline to import maps and traffic demands. In our simulations, each
parking lot can accommodate up to $50$ cars and we generated the traffic
demand so that $100$ cars would arrive on campus at $5$-second intervals. All
the cars seek to park and, by doing so, the parking capacity is saturated.
Given this setting, we simulated a road obstruction on the main link (in blue
in Figure 1) from the highway exit to the campus entrance. This was done by
restricting, in SUMO, the speed of the link to less than one kilometer per
hour. Information on the cars in the simulation are contained in the stand-
alone file agent.npy. Instead, the pfs associated with the sources (described
below) are all stored in behaviors.npy.

#### HiL set-up

the platform embeds a real vehicle into a SUMO simulation. By doing so,
performance of the algorithm deployed on a real car can be assessed under
arbitrary synthetic traffic conditions generated via SUMO. The high-level
architecture of the platform is shown in Figure 2. The platform creates an
avatar of the real car in the SUMO simulation. Namely, as shown in Figure 2,
the position of the real car is embedded in SUMO by using a standard
smartphone to collect its GPS coordinates. These coordinates are then sent via
bluetooth to a computer, also hosted on the car in the current implementation.
The connection is established via an off-the-shelf app, which writes the
coordinates in a rfcomm file. By using the pySerial library, an interface was
developed to read data in Python. Here, a script was designed leveraging
pynmeaGPS to translate the data in the NMEA format for longitude/latitude
coordinates. With this data format, a Python script was created to place the
avatar of the real car in the position given by the coordinates. A GUI is also
included to highlight the trajectory of the real car on the map and an off-
the-shelf text-to-speech routine is used to provide audio feedback to the
driver on the vehicle.

Figure 2: HiL functional architecture.

###  V-B In-car Implementation of the Algorithm

For both the in-silico and in-vivo experiments, Algorithm 1 was implemented in
Python as a stand-alone class so that each car equipped with the algorithm
could function as a stand-alone agent. The class has methods accepting all the
input parameters of Algorithm 1 and providing as output the car behavior
computed by the algorithm. The optimization problems within the algorithm loop
were solved via the Python library scipy.optimize. Additionally, the class
also implements methods to compute the cost and to support receding horizon
implementations of Algorithm 1. In our experiments, we used this receding
horizon implementation: the width of the receding horizon window was $T=5$ and
every time the car entered in a new link/state computations were triggered.
The pfs from the sources in our experiments were such that, at each $k$,
feasibility of the problem was ensured (see Lemma 1). Following [20], we also
implemented an utility function that restricts calculations of the agent only
to the road links that can be reached in $T$ time steps (rather than through
the whole state space/map). With this feature, in our experiments the
algorithm took on average approximately half a second to output a behavior,
less than the typical time taken to drive through a road link.111this duration
was measured between the moment where the GUI shows the car merging on a new
link and the moment where new directions are displayed. The simulation was run
on a modern laptop. Finally, the pfs of the sources were obtained via built-in
routing functions in SUMO and we added noise to the routes so that Assumption
2 would be fulfilled (for each road link, $\mathcal{S}(\pi)$ is the set of
outgoing links). See our github for the details.

###  V-C Experimental Results

#### Simulation results

first, we benchmarked the performance obtained by Algorithm 1 against these
from the algorithm in [3, 4], termed as crowdsourcing algorithm in what
follows. To this aim, we considered Scenario $1$ and performed two sets of
$10$ simulations. In the first set of experiments, Algorithm 1 was used to
determine the behavior of cars on campus (note that Assumption 1 is
fulfilled). In the second set of simulations, the cars instead used the
crodwourcing algorithm. Across the simulations, we recorded the number of cars
that the algorithms were able to park. The results are illustrated in Figure 3
(top panel). The figure shows that the crowdsourcing algorithm was not able to
park all the cars within the simulation. This was instead achieved by
Algorithm 1, which outperformed the algorithm from [3, 4]. To further quantify
the performance, we also computed the average time spent by a car looking for
a parking space after it enters the simulation (ATTP: average time-to-
parking). Across the simulations, the ATTP for the algorithm in [3] was of
$224.74\pm 19.67$, while for Algorithm 1 it was of $151.32\pm 30.59$ (first
quantities are means, second quantities are standard deviations). That is,
Algorithm 1 yielded an average improvement of $32.7\%$ in the ATTP. Then, we
simulated the three cases of Scenario $2$ to verify that Algorithm 1 can
effectively regulate access through specific links. The constraints for the
three cases of Scenario $2$ were given as an input to the algorithm and in
Figure 3 (bottom panel) the optimal solution $\pi^{\ast}({x}_{k}\mid
x_{k-1}=l_{r})$ is shown. The figure shows that the optimal solution indeed
fulfills the constraints (the road link $l_{r}$ is in Figure 1).

Figure 3: Top: unparked cars over time for crowdsourcing and Algorithm 1.
Solid lines are means across the simulations, shaded areas are confidence
intervals (one standard deviation). Bottom: $\pi^{\ast}({x}_{k}\mid
x_{k-1}=l_{r})$ for the three cases of Scenario $2$. The pfs satisfy the
constraints. Link definitions in Figure 1.

#### HiL results

we deployed Algorithm 1 on a real vehicle using the HiL platform and validated
its effectiveness in Scenario $2$ (C): the target behavior of the agent would
make the real car reach the Terminal parking but this route is forbidden. What
we observed in the experiment was that, once the car entered in the campus,
this was re-routed towards the Biblioteca parking. The re-routing was an
effect of Algorithm 1 computing the rightmost pf in Figure 3 (bottom panel). A
video of the HiL experiment is available on our github. The video shows that
the algorithm is suitable for real car operation: it would run smoothly during
the drive, providing feedback to the driver on the vehicle. Figure 4 shows the
car’s route recorded during the experiment.

Figure 4: Route of the real vehicle. The continuous line shows the GPS
position during the HiL experiment (map from OpenStreetMaps).

##  VI Conclusions

We considered the problem of designing agents able to compute optimal
decisions by re-using data from multiple sources to solve tasks involving: (i)
tracking a desired behavior while minimizing an agent-specific cost; (ii)
satisfying certain safety constraints. After formulating the control problem,
we showed that this is convex under a mild condition and computed the optimal
solution. We turned the results in an algorithm and used it to design an
intelligent parking system. We evaluated the algorithm via in-silico and in-
vivo experiments with real vehicles/drivers. All experiments confirmed the
effectiveness of the algorithm and its suitability for in-car operation.
Besides considering the use of other divergences in the cost and deploying our
results in more complex urban scenarios that would use data from pedestrians
and sensors on-board vehicles, our future research will involve devising
mechanisms for the composition of policies for the tasks with actuation
constraints in [22].

## References

  * [1] B. M. Lake _et al._ , “Human-level concept learning through probabilistic program induction,” _Science_ , vol. 350, no. 6266, pp. 1332–1338, 2015\.
  * [2] E. Crisostomi _et al._ , _Analytics for the sharing economy: Mathematics, Engineering and Business perspectives_. Springer, 2020.
  * [3] G. Russo, “On the crowdsourcing of behaviors for autonomous agents,” _IEEE Cont. Sys. Lett._ , vol. 5, pp. 1321–1326, 2020.
  * [4] É. Garrabé and G. Russo, “On the design of autonomous agents from multiple data sources,” _IEEE Cont. Sys. Lett._ , vol. 6, pp. 698–703, 2021\.
  * [5] L. Li, C. De Persis, P. Tesi, and N. Monshizadeh, “Data-based transfer stabilization in linear systems,” 2022. [Online]. Available: https://arxiv.org/abs/2211.05536
  * [6] J. Coulson, J. Lygeros, and F. Dörfler, “Data-enabled predictive control: In the shallows of the DeePC,” in _European Control Conference_ , 2019, pp. 307–312.
  * [7] H. J. van Waarde, J. Eising, H. L. Trentelman, and M. K. Camlibel, “Data informativity: a new perspective on data-driven analysis and control,” _IEEE Trans. Automatic Control_ , vol. 65, pp. 4753–4768, 2020.
  * [8] C. De Persis and P. Tesi, “Formulas for data-driven control: Stabilization, optimality, and robustness,” _IEEE Trans. Automatic Control_ , vol. 65, pp. 909–924, 2020.
  * [9] F. Celi and F. Pasqualetti, “Data-driven meets geometric control: Zero dynamics, subspace stabilization, and malicious attacks,” _IEEE Cont. Sys. Lett._ , vol. 6, pp. 2569–2574, 2022.
  * [10] U. Rosolia and F. Borrelli, “Learning model predictive control for iterative tasks. A data-driven control framework,” _IEEE Trans. Automatic Control_ , vol. 63, pp. 1883–1896, 2018.
  * [11] K. P. Wabersich and M. N. Zeilinger, “Bayesian model predictive control: Efficient model exploration and regret bounds using posterior sampling,” ser. Proc. of ML Research, vol. 120, 2020, pp. 455–464.
  * [12] J. Berberich, J. Kohler, M. A. Muller, and F. Allgower, “Data-driven model predictive control with stability and robustness guarantees,” _IEEE Trans. Aut. Contr._ , vol. 66, pp. 1702–1707, 2021.
  * [13] G. Baggio, V. Katewa, and F. Pasqualetti, “Data-driven minimum-energy controls for linear systems,” _IEEE Cont. Sys. Lett._ , vol. 3, pp. 589–594, 2019\.
  * [14] Émiland Garrabé and G. Russo, “Probabilistic design of optimal sequential decision-making algorithms in learning and control,” _Annual Rev, in Control_ , vol. 54, pp. 81–102, 2022.
  * [15] N. Cammardella, A. Busic, Y. Ji, and S. P. Meyn, “Kullback-Leibler-Quadratic optimal control of flexible power demand,” in _IEEE 58th Conference on Decision and Control_ , 2019, pp. 4195–4201.
  * [16] N. Cesa-Bianchi and G. Lugosi, _Prediction, learning, and games._ Cambridge University Press, 2006.
  * [17] M. Cutler, T. J. Walsh, and J. P. How, “Real-world reinforcement learning via multifidelity simulators,” _IEEE Trans. on Robotics_ , vol. 31, pp. 655–671, 2015.
  * [18] V. B. Mountcastle, “The columnar organization of the neocortex.” _Brain_ , vol. 120, pp. 701–722, Apr 1997.
  * [19] W. Griggs _et al._ , _A Vehicle-in-the-Loop Emulation Platform for Demonstrating Intelligent Transportation Systems_. Cham: Springer International Publishing, 2019, pp. 133–154.
  * [20] É. Garrabé and G. Russo, “CRAWLING: a Crowdsourcing Algorithm on Wheels for Smart Parking,” 2022, preprint submitted to Scientific Reports. [Online]. Available: https://arxiv.org/abs/2212.02467
  * [21] Pablo Alvarez Lopez and others, “Microscopic traffic simulation using SUMO,” in _21st IEEE International Conference on Intelligent Transportation Systems_ , 2018, pp. 2575–2582.
  * [22] D. Gagliardi and G. Russo, “On a probabilistic approach to synthesize control policies from example datasets,” _Automatica_ , vol. 137, p. 110121, 2022\.

## Appendix

We report here an investigation of how the time taken by Algorithm 1 changes
as a function of the number of sources, $S$, and time horizon $T$. This time
is a crucial aspect to investigate whether the approach we propose would scale
to more complex urban scenarios than the one presented in Section V, which
would e.g., include more parking locations (note that these are seen as links
by Algorithm 1) and more complex road networks together with more sources that
the agent could use to determine its optimal behavior. The time Algorithm 1
takes to output a behavior depends on the number of available sources, the
time horizon $T$ and the number of links accessible to the car within the time
horizon. We analyze the computation time w.r.t. the
To investigate Algorithm 1’s computation time, we considered the same
implementation as in Section V-B and ran the algorithm by varying the receding
horizon window between $0$ and $5$ time steps, and the number of sources
available to the agent between $1$ and $6$. For this, additional sources were
taken from [20], where more sources were used to implement the algorithm from
[4]. For each pair of these parameters, we measured the time taken by
Algorithm 1 to output a behavior, by running the algorithm over each link in
the network on a standard computer and taking the average of such times. This
gives a fair estimate of the average runtime, as the amount of states
considered depends on the density of the connections in the neighborhood of
each link.

The results of this numerical investigation are shown in Figure 5. The figure
shows that the highest computation time is about one second, which appears to
be suitable for our reference application.

Figure 5: Computation time as a function of time horizon and number of sources


# Intersection properties of finite disk collections

Jesús F. Espinoza<EMAIL_ADDRESS>and  Cynthia G. Esquer-
Pérez<EMAIL_ADDRESS>Departamento de Matemáticas, Universidad de
Sonora, México

###### Abstract.

In this work we study the intersection properties of a finite disk system in
the euclidean space. We accomplish this by utilizing subsets of spheres with
varying dimensions and analyze specific points within them, referred to as
poles. Additionally, we introduce two applications: estimating the common
scale factor for the radii that makes the re-scaled disks intersects in a
single point, this is the Čech scale, and constructing the minimal Axis-
Aligned Bounding Box (AABB) that encloses the intersection of all disks in the
system.

##  1\. Introduction

One of the new techniques developed for the analysis of large clusters of
information, known as Big Data, is Topological Data Analysis (TDA). In TDA,
simplicial complexes associated with the data are constructed. These
structures include the Vietoris-Rips complex, the Čech complex, and the
piecewise linear lower star complex, among others. Of special interest to us
is the generalized Čech complex structure. Although the standard Čech complex
is formed by intersecting a collection of disks with a fixed radius, the
generalized version allows varying radii. This flexibility enables us to
highlight specific data points by assigning or weighting them with larger
and/or more rapidly expanding balls to them, while de-emphasizing others by
using smaller and/or slower growing balls. This approach proves valuable for
handling noisy data sets, offering an alternative to discarding data that may
not meet a specific significance threshold [1].

Understanding the patterns of intersections and the timing of intersections
among a set of disks in $\mathbb{R}^{d}$, each with potentially different
radii, is a fundamental problem. This leads to the exploration of the
generalized Čech complex structure, which captures the intersection
information of these disks, regardless of their radii. Rescaling the radii by
the same factor, we obtain a filtered generalized Čech complex, where the
associated simplicial complexes evolve as the scale parameter varies. In
particular, in [4], algorithms are provided to calculate the generalized Čech
complex in $\mathbb{R}^{2}$, and [3] presents an algorithm to determine the
Čech scale for a collection of disks in the plane.

To establish the necessary foundation for our study, Section 2 introduces
crucial concepts and notation that will be used throughout the article and we
focus on analyzing the intersection of a disk system in $\mathbb{R}^{d}$. We
start by investigating the intersection of two disks in Subsection 2.1 and
then expand our analysis to a system of $m$ disks in Subsection 2.2. By
applying Helly’s Theorem, we prove that it is sufficient to examine the
intersection of all subsystems consisting of $d+1$ disks in order to determine
if the system has a nonempty intersection.

In Section 3, we define Vietoris-Rips systems and Čech systems, together with
presenting results regarding the Rips scale and the Čech scale, as well as
their connections. In Subsection 3.1, we present an algorithm that can
determine whether the intersection of the system is empty or non-empty. This
is achieved by exclusively computing the poles of subsystems of disks (or
spheres). Finally, in Subsection 3.2, we introduce the algorithm that computes
an approximation to the Čech scale using the numerical bisection method.

Additionally, in Section 4 we incorporate the concept of an minimal axis-
aligned bounding box (AABB) into our methodology. An AABB is a rectangular
parallelepiped whose faces are perpendicular to the basis vectors. These
bounding boxes frequently arise in spatial subdivision problems, such as ray
tracing [5] and collision detection [2]. In this paper, we study AABBs to
enclose the intersection of a finite collection of disks. This approach proves
valuable for discerning whether the collection intersects at a singular point
or not. In this section, we also provide an algorithm for constructing the
AABB of a disk system.

##  2\. Intersection properties of sphere systems

Throughout this work, we will refer to a $d$-disk system $M$, or simply a disk
system, as a finite collection of closed disks in $\mathbb{R}^{d}$ with
positive and not necessarily equal radii, i.e.,

$M=\\{D_{i}(c_{i};r_{i})\subset\mathbb{R}^{d}\mid
c_{i}\in\mathbb{R}^{d},r_{i}>0,1\leq i\leq m<\infty\\}.$

Moreover, in order to study the intersection properties of a disk system $M$
with the approach addressed in Sections 4 and 3 of this work, we will conduct
a study in this section of the intersection properties of the spheres
corresponding to the boundaries of each disk in $M$, which we call a sphere
system and denote by $\partial M$,

$\partial M=\\{\partial D_{i}\subset\mathbb{R}^{d}\mid D_{i}\in M\\},$

where $\partial$ denotes the topological boundary operator.

Following the notation in [6], we introduce the following generalization of
the sphere.

######  Definition 1.

An $i$-sphere in $\mathbb{R}^{d}$ is the intersection of a sphere with an
affine subspace of dimension $i$.

Of course, the notions of a sphere (as a $(d-1)$-dimensional surface) and a
$d$-sphere in $\mathbb{R}^{d}$ agree. However, an $i$-sphere in
$\mathbb{R}^{d}$ can also be viewed as the intersection of $d$-spheres. For
instance, the intersection of two spheres typically occurs in a hyperplane,
forming a $(d-1)$-sphere in $\mathbb{R}^{d}$. When another $d$-sphere
intersects this configuration, the result may be a $(d-1)$-sphere, a
$(d-2)$-sphere, a $0$-sphere (a single point), or it might even be empty, all
within the same hyperplane. For a disk system $M=\\{D_{i}(c_{i};r_{i})\\}$
composed of $m$ disks, where $\\{c_{1},\dots,c_{m}\\}$ is a set in general
position in $\mathbb{R}^{d}$, the maximum dimension of the affine subspace
associated with the $i$-sphere, obtained from the intersection of all the
spheres in $\partial M$, is at most $d-m+1$, or equivalently, $i=d-m+1$. This
conclusion is drawn from [6, Theorem 2.1] and the fact that the affine hull of
$\\{c_{1},\dots,c_{m}\\}$ is of dimension $m-1$. Consequently, the following
result is proven.

######  Lemma 2.

Let $M=\\{D_{1}(c_{1};r_{1}),\ldots,D_{m}(c_{m};r_{m})\\}$ be disk system such
that $\\{c_{1},\dots,c_{m}\\}$ is a set in general position in
$\mathbb{R}^{d}$. Then, the possibilities for the set $\cap_{D_{i}\in
M}\partial D_{i}$ are:

  1. (1)

the empty set;

  2. (2)

a single point;

  3. (3)

a $(d-m+1)$-sphere.

Remarkable points in $i$-spheres that will play a key role in the rest of the
article are the poles. Let $\pi_{i}:\mathbb{R}^{d}\longrightarrow\mathbb{R}$
be the canonical projection on the $i$-th factor for $i=1,...,d$, and let
$\\{e_{1},e_{2},...,e_{d}\\}$ be the standard basis of $\mathbb{R}^{d}$.

######  Definition 3.

Let $e_{q}$ be the $q$-th vector of the canonical base of $\mathbb{R}^{d}$. An
$e_{q}$-north (south) pole of an $i$-sphere $S$ in $\mathbb{R}^{d}$ is a point
on $S$ whose projection on the $q$-th coordinate is maximum (minimum). In
other words, $x\in S$ is the $e_{q}$-north pole if $\pi_{q}(y)\leq\pi_{q}(x)$
for all $y\in S-\\{x\\}$, where $\pi_{q}$ represents the projection onto the
$q$-th coordinate.

We denote the $e_{q}$-north pole of $S$ by $s_{q}^{+}$ and the $e_{q}$-south
pole by $s_{q}^{-}$.

An $i$-sphere can have a single $e_{q}$-pole (north or south) or an infinite
number of them, which occurs when a normal vector to the affine space
containing the $i$-sphere is aligned with the vector $e_{q}$. We are
interested in finding the $e_{q}$-poles of $(d-m+1)$-spheres originating from
disk systems $M=\\{D_{1}(c_{1};r_{1}),\ldots,D_{m}(c_{m};r_{m})\\}$, by taking
the intersection $\cap_{j=1}^{m}\partial D_{j}$. Such $(d-m+1)$-spheres will
be denoted by $S_{M}(c;r)$, to emphasize the disk system $M$, as well as its
center and radius.

######  Lemma 4.

Let $M=\\{D_{1},\ldots,D_{m}\\}$ be a $d$-disk system such that
$\bigcap_{j=1}^{m}D_{j}\neq\emptyset$, and let $p$ be a point in
$\bigcap_{j=1}^{m}D_{j}$ such that $\pi_{q}(p)\leq\pi_{q}(x)$ (resp.
$\pi_{q}(p)\geq\pi_{q}(x)$) for every $x$ in $\bigcap_{j=1}^{m}D_{j}$. Then,
there exists an $i$-sphere $S=\partial D_{j_{1}}\cap\cdots\cap\partial
D_{j_{i}}$ such that $p$ is in $S$ and $p$ is the $e_{q}$-south pole (resp.
$e_{q}$-north pole) of $S$.

###### Proof.

Since $\cap_{j=1}^{m}D_{j}\neq\emptyset$, then
$\partial(\cap_{j=1}^{m}D_{j})\neq\emptyset$,
$\partial(\cap_{j=1}^{m}D_{j})\subset\cap_{j=1}^{m}D_{j}$ due to the
closedness of the sets $D_{j}$, for $j=1,\ldots,m$, and
$p\in\partial(\cap_{j=1}^{m}D_{j})$.

On the other hand, since
$\partial(\cap_{j=1}^{m}D_{j})\subset\cup_{i=1}^{m}\partial D_{j}$, there
exist indices $j_{1},\ldots,j_{i}$ such that $p\in\partial D_{j_{r}}$ for any
$r=1,\ldots,i$; let
$\Lambda(p)=\\{j_{1},\ldots,j_{i}\\}\subseteq\\{1,\ldots,m\\}$ be a maximal
subset of indices such that $p\in D_{j}$ if and only if $j\in\Lambda(p)$. We
claim that $p\in\cap_{r=1}^{i}\partial D_{j_{r}}$ is the $e_{q}$-south pole of
$S:=\cap_{r=1}^{i}\partial D_{j_{r}}$.

In effect, let $V_{p}\subset\mathbb{R}^{d}$ be an open neighborhood of $p$
sufficiently small such that:

  1. (1)

Every $x\in V_{p}\cap\partial(\cap_{j=1}^{m}D_{j})$ has as maximal set of
indices a proper subset of $\Lambda(p)$,

  2. (2)

$S\cap V_{p}\subset\partial(\cap_{j=1}^{m}D_{j})$.

The first condition can be guaranteed by the finiteness of the disk system
$M$, and the second condition is a consequence of the maximality of the set
$\Lambda(p)$. Therefore, $\pi_{q}(p)\leq\pi_{q}(x)$ for every $x\in S\cap
V_{p}$, which is equivalent to the fact that $\pi_{q}(p)\leq\pi_{q}(x)$ for
every $x\in S$, in the case of $i$-spheres.

∎

###  2.1. Sphere systems with two spheres

In the following two lemmas we provide the computations to determine the
center, radius and poles for a $(d-1)$-sphere given by the intersection of two
disks in $\mathbb{R}^{d}$.

######  Lemma 5.

Let $M=\\{D_{1}(c_{1};r_{1}),D_{2}(c_{2};r_{2})\\}$ be a disk system with two
$d$-disks such that $\partial D_{1}\cap\partial D_{2}$ is a $(d-1)$-sphere
$S=S_{M}(c;r)$ with center $c$ and radius $r$. Then,

$c=\frac{1}{2}\left(1+\frac{r_{2}^{2}-r_{1}^{2}}{\|c_{2}-c_{1}\|^{2}}\right)c_{1}+\frac{1}{2}\left(1+\frac{r_{1}^{2}-r_{2}^{2}}{\|c_{2}-c_{1}\|^{2}}\right)c_{2}$

and

$r=\frac{2\sqrt{s(s-\|c_{2}-c_{1}\|)(s-r_{1})(s-r_{2})}}{\|c_{2}-c_{1}\|}$

where $s=\frac{1}{2}(\|c_{2}-c_{1}\|+r_{1}+r_{2})$.

###### Proof.

Let $\Pi$ be the hyperplane containing the $(d-1)$-sphere $S$, which is
defined by the equation:

(1)
$\sum_{i=1}^{d}(k_{i}-h_{i})x_{i}-\frac{1}{2}\sum_{i=1}^{d}(k_{i}^{2}-h_{i}^{2})=\frac{r_{1}^{2}-r_{2}^{2}}{2},$

where $c_{1}=(h_{1},\ldots,h_{d})$ and $c_{2}=(k_{1},\ldots,k_{d})$. Then the
normal vector of the hyperplane $\Pi$ is given by
$N:=c_{2}-c_{1}=(k_{1}-h_{1},\ldots,k_{d}-h_{d})$, and the center $c$ of $S$
is determined by the intersection point of the hyperplane $\Pi$ with the
perpendicular line that passes through the center $c_{1}$ of $D_{1}$. This
line can be parameterized as $\gamma:t\mapsto
c_{1}+tN=(x_{1}(t),\ldots,x_{d}(t))$, such that $\gamma(0)=c_{1}$ and
$\gamma(1)=c_{2}$. We can compute the intersection point $c=\gamma(t_{*})$ of
$\Pi$ and $\gamma([0,1])$, for any $t_{*}\in[0,1]$, by substituting it in (1),

$\sum_{i=1}^{d}(k_{i}-h_{i})(h_{i}+t_{*}(k_{i}-h_{i}))=\frac{1}{2}\left(\sum_{i=1}^{d}(k_{i}^{2}-h_{i}^{2})+(r_{1}^{2}-r_{2}^{2})\right).$

And solving for $t_{*}$, we obtain that
$t_{*}=\frac{1}{2}+\frac{r_{1}^{2}-r_{2}^{2}}{2\|c_{2}-c_{1}\|^{2}}$. Hence,
the center of $S$ is given by:

$\displaystyle
c=c_{1}+t_{*}N=\frac{1}{2}\left(1+\frac{r_{2}^{2}-r_{1}^{2}}{\|c_{2}-c_{1}\|^{2}}\right)c_{1}+\frac{1}{2}\left(1+\frac{r_{1}^{2}-r_{2}^{2}}{\|c_{2}-c_{1}\|^{2}}\right)c_{2}$

Next, we will compute the radius $r$ of $S$. This radius can be determined as
the height $r$ of the triangle with base $\|c_{2}-c_{1}\|$ formed by the
points $c_{1}$, $c_{2}$, and a point on $S$. Thus, by the Heron’s formula we
have that

$r=\frac{2\sqrt{s(s-\|c_{2}-c_{1}\|)(s-r_{1})(s-r_{2})}}{\|c_{2}-c_{1}\|},$

where $s=\frac{1}{2}(\|c_{2}-c_{1}\|+r_{1}+r_{2})$ correspond to the semi-
perimeter. ∎

We can proceed now to compute the poles of the $(d-1)$-sphere $\partial
D_{1}\cap\partial D_{2}$.

######  Lemma 6.

Let $D_{1}(c_{1};r_{1})$ and $D_{2}(c_{2};r_{2})$ be two $d$-disks such that
$\partial D_{1}\cap\partial D_{2}$ is a $(d-1)$-sphere $S=S(c;r)$ with center
$c$ and radius $r$. Then, the $e_{q}$-poles of $S$ are
$s_{q}^{\pm}=c\pm\sum_{i}^{d}x_{i}e_{i}$, where

$x_{i}=\begin{cases}\dfrac{r|\pi_{i}(c_{2}-c_{1})\pi_{q}(c_{2}-c_{1})|}{\|c_{2}-c_{1}\|\sqrt{\|c_{2}-c_{1}\|^{2}-\pi_{q}(c_{2}-c_{1})^{2}}},&i\neq
q\\\ \\\
\dfrac{r{\sqrt{\|c_{2}-c_{1}\|^{2}-\pi_{q}(c_{2}-c_{1})^{2}}}}{\|c_{2}-c_{1}\|},&i=q.\end{cases}$

###### Proof.

For simplicity, we translate the hyperplane $\Pi$, which contains the
$(d-1)$-sphere $S$, as well as the sphere itself, to the origin; in such case,
the corresponding equations are given by,

$\displaystyle\sum_{i=1}^{d}(k_{i}-h_{i})x_{i}$ $\displaystyle=0,$
$\displaystyle\sum_{i=1}^{d}x_{i}^{2}$ $\displaystyle=r^{2},$

where $h_{i}:=\pi_{i}(c_{1})$ and $k_{i}:=\pi_{i}(c_{2})$ for $i=1,\ldots,d$.
In the case that $k_{q}-h_{q}=\pi_{q}(c_{2}-c_{1})=0$, the normal vector
$N=c_{2}-c_{1}$ of the hyperplane $\Pi$ is orthogonal to the basis vector
$e_{q}$. Therefore, the $e_{q}$-poles of $S$ are $c\pm re_{q}$, which agree
with the formulae of the lemma.

On the other hand, suppose that $k_{q}-h_{q}\neq 0$. To find the $e_{q}$-poles
of $S$, we will use the Lagrange multiplier method. Consider the following
function:

(2) $x_{q}=f(x_{1},x_{2},...,\widehat{x_{q}},...,x_{d})=\frac{-\sum_{j\neq
q}^{d}(k_{j}-h_{j})x_{j}}{k_{q}-h_{q}},$

subject to the restriction:

$g(x_{1},x_{2},...,\widehat{x_{q}},...,x_{d})=\sum_{j\neq
q}^{d}x_{j}^{2}+\left(\frac{-\sum_{j\neq
q}^{d}(k_{j}-h_{j})x_{j}}{k_{q}-h_{q}}\right)^{2}-r^{2}=0$

Let $\lambda$ be the Lagrange multiplier, we define

$h(x_{1},...,\widehat{x_{q}},...,x_{d},\lambda)=f(x_{1},...,\widehat{x_{q}},...,x_{d})+\lambda
g(x_{1},...,\widehat{x_{q}},...,x_{d})$

For any $i\neq q$, consider the following system of equations:

$\frac{\partial h}{\partial x_{i}}=-\frac{k_{i}-h_{i}}{k_{q}-h_{q}}+2\lambda
x_{i}+2\lambda\left(\frac{-\sum_{j\neq
q}^{d}(k_{j}-h_{j})x_{j}}{k_{q}-h_{q}}\right)\left(-\frac{k_{i}-h_{i}}{k_{q}-h_{q}}\right)=0.$

Then

$\displaystyle-\frac{k_{i}-h_{i}}{k_{q}-h_{q}}+2\lambda\left(x_{i}+\frac{k_{i}-h_{i}}{(k_{q}-h_{q})^{2}}{\sum_{j\neq
q}^{d}(k_{j}-h_{j})x_{j}}\right)$ $\displaystyle=0$
$\displaystyle-(k_{i}-h_{i})(k_{q}-h_{q})+2\lambda\left((k_{q}-h_{q})^{2}x_{i}+(k_{i}-h_{i}){\sum_{j\neq
q}^{d}(k_{j}-h_{j})x_{j}}\right)$ $\displaystyle=0$

Solving this system of equations for $\lambda$, we obtain that,

$\lambda=\frac{(k_{i}-h_{i})(k_{q}-h_{q})}{2\left((k_{q}-h_{q})^{2}x_{i}+(k_{i}-h_{i}){\sum_{j\neq
q}^{d}(k_{j}-h_{j})x_{j}}\right)}$

Comparing the last expression for two indices $i\neq\tilde{i}$, we have that,

$\displaystyle x_{i}^{2}$
$\displaystyle=\frac{r^{2}(k_{i}-h_{i})^{2}}{(k_{\tilde{i}}-h_{\tilde{i}})^{2}+\sum_{j\neq\tilde{i},q}{(k_{j}-h_{j})}^{2}+\frac{1}{(k_{q}-h_{q})^{2}}\left(\sum_{j\neq
q}(k_{j}-h_{j})^{2}\right)^{2}}$
$\displaystyle=\frac{r^{2}(k_{i}-h_{i})^{2}}{\sum_{j\neq
q}{(k_{j}-h_{j})}^{2}+\frac{1}{(k_{q}-h_{q})^{2}}\left(\sum_{j\neq
q}(k_{j}-h_{j})^{2}\right)^{2}}$
$\displaystyle=\frac{r^{2}(k_{i}-h_{i})^{2}(k_{q}-h_{q})^{2}}{\left((k_{q}-h_{q})^{2}+\sum_{j\neq
q}(k_{j}-h_{j})^{2}\right)\left(\sum_{j\neq q}{(k_{j}-h_{j})}^{2}\right)}$
$\displaystyle=\frac{r^{2}\pi_{i}(c_{2}-c_{1})^{2}\pi_{q}(c_{2}-c_{1})^{2}}{\|c_{2}-c_{1}\|^{2}\left(\|c_{2}-c_{1}\|^{2}-\pi_{q}(c_{2}-c_{1})^{2}\right)}$

Finally, for $i=q$, we can use the last expression to substitute it in (2) and
obtain the desired result. ∎

###  2.2. Sphere systems with more than two spheres

Now, let us proceed with the explicit calculation of the coefficients for the
center $c$ of the $(d-m+1)$-sphere $S=\cap_{j=1}^{m}\partial D_{j}$. We can
achieve this by considering the disk system translated to $c_{m}$, denoted as
$\\{D_{j}(c_{j}-c_{m};r_{j})\\}_{j=1}^{m}$, and by defining the
$(d-m+1)$-sphere $S-\\{c_{m}\\}=\cap_{j=1}^{m}\partial
D_{j}(c_{j}-c_{m};r_{j})$. This sphere is positioned at the intersection of
hyperplanes (for more details, refer to [6]).

(3) $(c_{k}-c_{m})^{T}x=\frac{1}{2}(r_{m}^{2}+\|c_{k}-c_{m}\|^{2}-r_{k}^{2})$

for all $k=1,...,m-1$. Utilizing the information that the center of
$S-\\{c_{m}\\}$ can be expressed as a combination of the centers $c_{k}-c_{m}$
and substituting it into (3), we obtain a linear system of equations with
dimensions $(m-1)\times(m-1)$:

$\sum_{j=1}^{m-1}\lambda_{j}(c_{k}-c_{m})\cdot(c_{j}-c_{m})=\frac{1}{2}(r_{m}^{2}+\|c_{k}-c_{m}\|^{2}-r_{k}^{2})$

for $k=1,...,m-1$. Solving the system of equations for
$\lambda=(\lambda_{1},...,\lambda_{m-1})$, we find the center of $S$ as
follows:

$c=\lambda_{1}(c_{1}-c_{m})+...+\lambda_{m-1}(c_{m-1}-c_{m})+c_{m}$

The radius of the sphere $S$ can be computed using the equation:

$r^{2}=r_{k}^{2}-\|c-c_{k}\|^{2}$

for any $k\in\\{1,...,m-1\\}$.

Now that we have determined the center and radius of $S$, as well as the
affine space that contains it, we can proceed to compute its $e_{q}$-poles for
each $q\in\\{1,2,...,d\\}$. These poles reside in the affine space that
contains $S$ and within a set that we define below.

Let $S$ be an $i$-sphere in $\mathbb{R}^{d}$, and let $n_{1},...,n_{d-i}$ be
orthogonal vectors to the affine space $L$ that contains $S$. Consider the
space $M$ generated by these vectors together with the vector $e_{q}$ from the
canonical basis of $\mathbb{R}^{d}$. Let us denote
$n_{j}=\left(n_{l}^{(j)}\right)_{l=1}^{d}$ for each $j=1,...,d-i$. Then, we
can define $L_{0}$, the set $L$ translated to the origin, as follows:

$\displaystyle L_{0}$
$\displaystyle=\left<\\{n_{j}\\}_{j=1}^{d-i}\right>^{\perp}$
$\displaystyle=\left\\{x\in\mathbb{R}^{d}\mid x\cdot n_{j}=0,\hskip
5.69046pt\forall j=1,...,d-i\right\\}$

The set $M$ is defined as:

$\displaystyle M$
$\displaystyle=\left<\\{n_{j}\\}_{j=1}^{d-i}\cup\\{e_{q}\\}\right>$
$\displaystyle=\left\\{x=\sum_{j=1}^{d-i}\lambda_{j}n_{j}+\lambda_{d-i+1}e_{q}\in\mathbb{R}^{d}\mid\lambda_{j}\in\mathbb{R},\hskip
5.69046pt\forall j=1,...,d-i+1\right\\}$
$\displaystyle=\left\\{\left(\sum_{j=1}^{d-i}\lambda_{j}n_{1}^{(j)},...,\sum_{j=1}^{d-i}\lambda_{j}n_{q}^{(j)}+\lambda_{d-i+1},...,\sum_{j=1}^{d-i}\lambda_{j}n_{d}^{(j)}\right)\mid\lambda_{j}\in\mathbb{R}\right\\}$

Refer to Figure 1 for a visual representation of the subspaces $L$ and
$M+\\{c\\}$.

Figure 1. Visualization of the subspaces $L$ and $M+\\{c\\}$.

As mentioned above, the $e_{q}$-poles of $S$ lie at the intersection of $L$
and $M+\\{c\\}$, where $c$ is the center of $S$. To simplify the calculations,
we will utilize $L_{0}$ and $M$, and then translate them into $c$. The
intersection of $M$ and $L_{0}$ can be expressed as follows:

$\displaystyle M\cap L_{0}$ $\displaystyle=\left\\{x\in M\mid x\cdot
n_{k}=0,\hskip 5.69046pt\forall k=1,...,d-i\right\\}$
$\displaystyle=\left\\{\sum_{j=1}^{d-i}\lambda_{j}n_{j}+\lambda_{d-i+1}e_{q}\in\mathbb{R}^{d}\left|\left(\sum_{j=1}^{d-i}\lambda_{j}n_{j}+\lambda_{d-i+1}e_{q}\right)\cdot
n_{k}=0,\hskip 5.69046pt\lambda_{j}\in\mathbb{R}\hskip 5.69046pt\forall
k=1,...,d-i\right.\right\\}$
$\displaystyle=\left\\{\sum_{j=1}^{d-i}\lambda_{j}n_{j}+\lambda_{d-i+1}e_{q}\in\mathbb{R}^{d}\left|\sum_{j=1}^{d-i}\lambda_{j}n_{j}\cdot
n_{k}+\lambda_{d-i+1}n_{q}^{k}=0,\hskip
5.69046pt\lambda_{j}\in\mathbb{R}\hskip 5.69046pt\forall
k=1,...,d-i\right.\right\\}$

Let us consider a disk system in $\mathbb{R}^{d}$, denoted
$\\{D_{j}(c_{j};r_{j})\\}_{j=1}^{m}$, where $m<d$. The intersection of their
boundaries forms a $(d-m+1)$-sphere $S$. In this case, the subspace $M$ has
dimension $m$, or $dim(M)=m-1$ if $e_{q}\in M$. We choose the normal vectors
for the affine space containing $S$ as $n_{j}=c_{j}-c_{m}$, where
$j=1,...,m-1$. Then

$M=\left\\{\left(\sum_{j=1}^{m-1}\lambda_{j}\left(c_{1}^{(j)}-c_{1}^{(m)}\right),...,\sum_{j=1}^{m-1}\lambda_{j}\left(c_{q}^{(j)}-c_{q}^{(m)}\right)+\lambda_{m},...,\sum_{j=1}^{m-1}\lambda_{j}\left(c_{d}^{(j)}-c_{d}^{(m)}\right)\right)\left|\lambda_{j}\in\mathbb{R}\right.\right\\}$

By rewriting, we have

$\displaystyle M\cap L_{0}$
$\displaystyle=\left\\{\sum_{j=1}^{m-1}\lambda_{j}n_{j}+\lambda_{m}e_{q}\in\mathbb{R}^{d}\left|\sum_{j=1}^{m-1}\lambda_{j}n_{j}\cdot
n_{k}+\lambda_{m}n_{q}^{k}=0,\hskip 5.69046pt\lambda_{j}\in\mathbb{R},\hskip
5.69046pt\forall k=1,...,m-1\right.\right\\}$

If $S(c;r)=\cap_{j=1}^{m}\partial D_{j}$ is the $(d-m+1)$-sphere with center
in $c$ and radius $r$, then the $e_{q}$-poles of $S$ are the $e_{q}$-poles of
$S-\\{c_{m}\\}$ but translated by $c$. The poles of $S-\\{c_{m}\\}$ are
located in $M\cap L_{0}$. If $p$ is an $e_{q}$-pole of $S-\\{c_{m}\\}$, then
it can be expressed as

$p=\sum_{j=1}^{m-1}\lambda_{j}n_{j}+e_{q}\lambda_{m}$

for some $\lambda_{j}\in\mathbb{R}$, $j=1,...,m$ and the following conditions
holds:

$\sum_{j=1}^{m-1}\lambda_{j}n_{j}\cdot n_{k}+\lambda_{m}n_{q}^{k}=0$

for each $k=1,...,m-1$ and

$\|p\|^{2}=r^{2}$

Thus, if $p=\sum_{j=1}^{m-1}\lambda_{j}n_{j}+e_{q}\lambda_{m}$ is an
$e_{q}$-pole of $S-\\{c_{m}\\}$, the following equations are satisfied for
$\lambda_{1},\lambda_{2},...,\lambda_{m}\in\mathbb{R}$:

(4) $\sum_{j=1}^{m-1}\lambda_{j}n_{j}\cdot n_{k}+\lambda_{m}n_{q}^{k}=0$

(5)
$\sum_{i=1}^{d}\left(\sum_{j=1}^{m-1}\lambda_{j}n_{i}^{j}\right)^{2}+2\lambda_{m}\sum_{j=1}^{m-1}\lambda_{j}n_{q}^{j}+\lambda_{m}^{2}=r^{2}$

for all $k=1,...,m-1$, with $r$ the radius of the $(d-m+1)$-sphere $S$. From
(4) we have the system

$\begin{pmatrix}n_{1}\cdot n_{1}&n_{1}\cdot n_{2}&....&n_{1}\cdot n_{m-1}\\\
n_{2}\cdot n_{1}&n_{2}\cdot n_{2}&....&n_{2}\cdot n_{m-1}\\\ ...\\\
n_{m-1}\cdot n_{1}&n_{m-1}\cdot n_{2}&....&n_{m-1}\cdot n_{m-1}\\\
\end{pmatrix}\begin{pmatrix}\lambda_{1}\\\ \lambda_{2}\\\ ...\\\
\lambda_{m-1}\end{pmatrix}+\lambda_{m}\begin{pmatrix}n_{q}^{1}\\\ n_{q}^{2}\\\
...\\\ n_{q}^{m-1}\end{pmatrix}=0.$

Let us denote $A$ as the matrix $(n_{i}\cdot n_{j})_{i,j}$ and $B$ as the
vector $(-n_{q}^{j})_{j=1}^{m-1}$. Then, we have $A\lambda=\lambda_{m}B$,
where $\lambda=(\lambda_{1},...,\lambda_{m-1})$. Solving for $\lambda$, we
obtain $\lambda_{j}=\lambda_{m}(A^{-1}B)[j]$ for each $j=1,...,m-1$ (where
$(A^{-1}B)[j]$ denotes the entry $j$ of the $(m-1)\times 1$ vector
$(A^{-1}B)$). By substituting the value of $\lambda_{j}$ into (5), we obtain
the quadratic equation:

$\lambda_{m}^{2}\left[\sum_{i=1}^{d}\left(\sum_{j=1}^{m-1}(A^{-1}B)[j]n_{i}^{j}\right)^{2}+2\sum_{j=1}^{m-1}(A^{-1}B)[j]n_{q}^{j}+1\right]-r^{2}=0$

Let us define

$\Gamma_{i}=\sum_{j=1}^{m-1}(A^{-1}B)[j]n_{i}^{j}$

for all $i=1,...,d$. Solving this equation, we find:

$\lambda_{m}=\frac{\pm r}{\sqrt{\sum_{i=1}^{d}\Gamma_{i}^{2}+2\Gamma_{q}+1}}$
$\lambda_{j}=\frac{\pm
r(A^{-1}B)[j]}{\sqrt{\sum_{i=1}^{d}\Gamma_{i}^{2}+2\Gamma_{q}+1}}$

for $j=1,...,m-1$. Therefore, the $e_{q}$-poles of $S$, for $q=1,...,d$, are:

$p=\sum_{j=1}^{m-1}\lambda_{j}n_{j}+e_{q}\lambda_{m}+c=\frac{\pm
r}{\sqrt{\sum_{i=1}^{d}\Gamma_{i}^{2}+2\Gamma_{q}+1}}\left(\sum_{j=1}^{m-1}(A^{-1}B)[j]n_{j}+e_{q}\right)+c$

##  3\. Vietoris-Rips and Čech systems

Our goal in this section is to provide a comprehensive understanding of the
disk system, the Vietoris-Rips system, and the Čech system. Additionally, we
introduce some results that establish a certain connection between both disk
systems. Investigating the features and qualities of data and spaces can
provide us with useful knowledge about their geometric and topological
characteristics.

Before we look into the definitions of Vietoris-Rips and Čech systems, let us
give a brief overview. These systems are essential in the field of topological
data analysis for recognizing and comprehending the geometric structure of
point cloud data. Vietoris-Rips complex and Čech complex share the goal of
capturing the topology of the underlying metric space, both provide different
ways of recognizing connections and associations among data points. The
Vietoris-Rips complex tends to be more efficient and scalable for large
datasets, while the Čech complex can be more accurate but computationally more
expensive. The choice between the two depends on the nature of the dataset and
the specific goals of the topological analysis. Now, let us move on to
defining these fundamental concepts.

######  Definition 7.

Let $M=\\{D_{1},D_{2},\ldots,D_{m}\\}$ be a $d$-disk system. We say $M$ is a
Vietoris-Rips system if $D_{i}\cap D_{j}\neq\emptyset$ for each pair
$i,j\in\\{1,2,\ldots,m\\}$. Furthermore, if the $d$-disk system $M$ has the
nonempty intersection property $\bigcap_{D_{i}\in M}D_{i}\neq\emptyset$, then
$M$ is called a Čech system.

For each $\lambda\geq 0$, we define a collection of $d$-disks $M_{\lambda}$
with the same centers as those in the $d$-disk system $M$, but with radii
rescaled by $\lambda$. When $\lambda>0$, $M_{\lambda}$ is a $d$-disk system
again. $M_{1}$ is equal to $M$, and $M_{0}$ is the set of the centers of the
$d$-disks in $M$.

In the field of topological data analysis, the Rips scale and the Čech scale
are essential parameters for determining the closeness and connectivity
between data points. These two scales offer different perspectives on how we
measure and comprehend geometric relationships within point-cloud data. To
understand their importance in capturing the underlying topological structure,
let us look at their definitions. The Vietoris-Rips scale $\nu_{M}$, of a
$d$-disk system $M$ is the smallest $\lambda\in\mathbb{R}$ such that
$M_{\lambda}$ is a Vietoris-Rips system. Similarly, the Čech scale $\mu_{M}$,
of $M$ is the smallest $\lambda\in\mathbb{R}$ such that $M_{\lambda}$ is a
Čech system. This is

$\displaystyle\nu_{M}$ $\displaystyle=\inf\\{\lambda\in\mathbb{R}\mid
M_{\lambda}\mbox{ is a Vietoris-Rips system}\\}$ $\displaystyle\mu_{M}$
$\displaystyle=\inf\\{\lambda\in\mathbb{R}\mid M_{\lambda}\mbox{ is a \v{C}ech
system}\\}$

Next, we present some easily observable properties for both scales. It can be
easily seen that $M$ is a Vietoris-Rips system if and only if $\nu_{M}\leq 1$
(in particular $\nu_{M_{\nu_{M}}}=1$); similarly, $M$ is a Čech system if and
only if $\mu_{M}\leq 1$.

Note that for a given $d$-disk system $M=\\{D_{1},D_{2},\ldots,D_{m}\\}$ the
Vietoris-Rips scale is

$\nu_{M}=\max_{i<j}\\{\|c_{i}-c_{j}\|/(r_{i}+r_{j})\\}$

where $c_{i}$ and $r_{i}$ are the center and radii of $D_{i}$. An additional
observation is that, in cases where the disk system consists of either one or
two disks, the Vietoris-Rips scale coincides with the Čech scale. It is
evident that every Čech system is also a Vietoris-Rips system; however, the
reverse assertion, in general, is not true.

Conversely, if the $d$-disk system contains at least three disks, determining
the Čech scale becomes more complex. In the context of Čech scale, the
following remark is important and play a key role in implementation (see [3]
for details).

######  Remark 8.

If $\mu_{M}$ is the Čech scale for $M$, then the $\mu_{M}$-rescaled system
$M_{\mu_{M}}$, has only one point in the intersection $\bigcap_{D_{i}\in
M}D_{i}(c_{i};\mu_{M}r_{i})$.

As we have mentioned, a Čech system is also a Vietoris-Rips system, but the
converse is not true. What we can affirm is that if a system is a Vietoris-
Rips system, then the system rescaled by the factor $\sqrt{2d/(d+1)}$ is also
a Čech system. This is established by the following lemma, the proof of which
can be found in [3].

######  Lemma 9.

Let $M=\\{D_{i}(c_{i};r_{i})\\}$ be a $d$-disk system in euclidean space
$\mathbb{R}^{d}$. If $D_{i}(c_{i};r_{i})\cap D_{j}(c_{j};r_{j})\neq\emptyset$
for every pair of disks in $M$, then

$\bigcap_{D_{i}\in M}D_{i}(c_{i};\sqrt{2d/(d+1)}\,r_{i})\neq\emptyset.$

One of the implications of the previous result is that, for any given disk
system $M$, we can bound the Čech scale using the Vietoris-Rips scale
$\nu_{M}$. This is stated by the following corollary.

######  Corollary 10.

If $M$ is an arbitrary $d$-disk system and $\nu_{M}$ is its Vietoris-Rips
scale, then its Čech scale satisfies
$\nu_{M}\leq\mu_{M}\leq\sqrt{2d/(d+1)}\,\nu_{M}$. Therefore, for every
$d$-disk system $M$, the rescaled disk system $M_{\sqrt{2d/(d+1)}\,\nu_{M}}$
is always a Čech system. In particular, if $\sqrt{2d/(d+1)}\,\nu_{M}\leq 1$
then $M_{\nu_{M}}$ is a Čech system.

###  3.1. Algorithm for determining Čech system.

In the previous section, we have determined the $e_{q}$-poles for the
intersection of any number of disks in $\mathbb{R}^{d}$. If any of these poles
is in all disks of the disk system, it indicates that the system conforms to
the criteria of a Čech system. It is important to recognize that this result
streamlines our calculation process, focusing on specific points to establish
whether the system exhibits a non-empty intersection.

Given a system of $m$ disks in $\mathbb{R}^{d}$ where $m>d$, it is enough to
verify if every subsystem of $d+1$ disks qualifies as a Čech system to
conclude that the entire system of disks has a non-empty intersection. This
assertion is supported by the Helly’s Theorem.

Now, we introduce an algorithm that determines whether a disk system qualifies
as a Čech system. In simpler terms, if the disk system exhibits a non-empty
intersection, the algorithm outputs ”TRUE”; otherwise, it outputs ”FALSE”. The
algorithm operates by seeking poles within the intersections of the disk
boundaries, which, as we have observed, correspond to $i$-spheres. It
initiates the search for poles within individual disks and then progresses to
the pairwise intersections of the disk boundaries ($(d-2)-sphere$), continuing
the process iteratively. If a pole is found within the remaining disks, the
system is classified as a Čech system.

1

Input : A $d$-disk system $M=\\{D_{j}\\}_{j=1}^{m}$

Output : A logical TRUE/FALSE to indicate if $M$ is a Čech system

2 Initialize: $\texttt{Is\\_Cech\\_System}\leftarrow\texttt{FALSE}$

3 for  _$k\leftarrow 1$ to $m$_ do

4 Let $\mathcal{S}$ be the set of $(d-k+1)$-spheres of $\partial M$

5 for  _$S$ in $\mathcal{S}$_ do

6 for  _$q\leftarrow 1$ to $d$_ do

7 Compute the set $\\{s_{q}^{\pm}\\}$ of $e_{q}$-poles of $S$

8 for  _$s\leftarrow\\{s_{q}^{\pm}\\}$_ do

9 if  _$s\in\cap_{j=1}^{d}D_{j}$_ then

10 $\texttt{Is\\_Cech\\_System}\leftarrow\texttt{TRUE}$

11 Go to line 16

12 end if

13

14 end for

15

16 end for

17

18 end for

19

20 end for

return (Is_Cech_System)

Algorithm 1 Cech.system

######  Theorem 11.

Let $M=\\{D_{1},\ldots,D_{m}\\}$ be a $d$-disk system. Then, $M$ is a Čech
system if and only if Cech.system$(M)=$ TRUE.

###### Proof.

If Cech.system$(M)=$ TRUE, then the Cech.system algorithm (Algorithm 1) found
a pole contained in the intersection $\cap_{j=1}^{m}D_{j}$, therefore
$\cap_{j=1}^{m}D_{j}\neq\emptyset$ and it follows that $M$ is a Čech system.

On the other hand, if $\bigcap_{j=1}^{m}D_{j}\neq\emptyset$, let $p$ be a
point in $\bigcap_{j=1}^{m}D_{j}$ satisfying $\pi_{1}(p)\leq\pi_{1}(x)$ for
every $x$ in $\bigcap_{j=1}^{m}D_{j}$. By Lemma 4, it follows that $p$ belongs
to an $i$-sphere and must be an $e_{1}$-south pole. Therefore, by the
exhaustive search of Algorithm 1 across all poles, its output is
Cech.system$(M)=$ TRUE. ∎

###  3.2. Algorithm to compute the Čech scale.

Finding the minimum parameter for which the rescaled system of disks has a
non-empty intersection is significant because it helps identify a critical
threshold at which the disks come into contact. This parameter, known as the
Čech scale, provides valuable information about the proximity or overlap of
the disks, which can be crucial in various applications such as collision
detection in computer graphics, spatial packing problems, and modeling
physical phenomena. In this section, we introduce an algorithm to compute an
approximation of the Čech scale for a system of $m$ disks in $\mathbb{R}^{d}$.

1

Input : A $d$-disk system $M$ in $\mathbb{R}^{d}$ and precision parameter
$\eta>0$

Output : $\mu_{M}$, a Čech scale approximation

2 Compute the Vietoris-Rips scale of $M$: $\nu_{M}$

3 if  _Cech.system $(M_{\nu_{M}})=$ TRUE_ then

4 $\mu_{M}\leftarrow\nu_{M}$

5 Go to line 16

6

7else

8 Initialize: $\mu_{M}^{*}\leftarrow\nu_{M}$,
$\mu_{M}\leftarrow\sqrt{2d/(d+1)}\nu_{M}$

9 while  _$\mu_{M}-\mu_{M}^{*} >\eta$_ do

10 Compute: $\lambda\leftarrow\dfrac{\mu_{M}^{*}+\mu_{M}}{2}$

11 if  _Cech.system $(M_{\lambda})=$ TRUE_ then

12 Update: $\mu_{M}\leftarrow\lambda$

13 else

14 Update: $\mu_{M}^{*}\leftarrow\lambda$

15 end if

16

17 end while

18

19 end if

return ($\mu_{M}$)

Algorithm 2 Cech.scale

The given code presents an algorithm to compute an approximation of the Čech
scale of a disk system in Euclidean space using Algorithm 1 and a precision
parameter $\eta>0$. It initializes the scale factor $\lambda$ to the Rips
scale $\nu_{M}$. If this scale satisfies Cech.system$(M_{\nu_{M}})=$ TRUE, it
indicates that the Čech scale has been found. Otherwise, we initiate a cycle
in which we compute Cech.system of the system rescaled by a factor $\lambda$.
The Čech scale is known to fall between the Rips scale and the value
$\sqrt{{2d}/{(d+1)}}\nu_{M}$ (Generalized Vietoris-Rips, Corollary 10). To
approximate the Čech scale, we employ the bisection method as long as the
interval enclosing the Čech scale has a length greater than $\eta$. Finally,
the algorithm returns an approximation of the Čech scale.

Utilizing the previously described algorithm, we can construct the filtered
generalized Čech complex for a disk system $M$. Let $\mathscr{C}(M)$ denote
the set of Čech subsystems, and $\mathscr{C}_{\lambda}(M)$ the set of Čech
subsystems for the rescaled disk system $M_{\lambda}$. The Čech filtration of
the $M$ system forms a maximal chain of Čech complexes

$\mathscr{C}_{*}(M):\mathscr{C}_{0}(M)\subsetneq\mathscr{C}_{\lambda_{1}}(M)\subsetneq\mathscr{C}_{\lambda_{2}}(M)\subsetneq...\subsetneq\mathscr{C}_{\mu_{M}}(M),$

where each $\lambda_{i}$ represents the Čech scale of the system
$M_{\lambda_{i}}$. Since the Čech scale of a disk system indicates the factor
by which we must rescale the system to make it Čech, defining a level of the
filtration, $\mathscr{C}_{\lambda}(M)$, simply requires determining the Čech
scale of the system $M_{\lambda}$.

##  4\. Minimal Axis-Aligned Bounding Box.

In this section, we introduce the concept of the minimal axis-aligned bounding
box (AABB) for the intersection of $d$-disks and present methods for its
computation. The AABB provides a simplified representation of the disk
intersection, making it easier to obtain valuable information about the disks.
This information could be useful for computing the Čech scale of a disk
system.

######  Definition 12.

Let $M$ be a disk system in $\mathbb{R}^{d}$. The minimal axis-aligned
bounding box of $M$, denoted as $AABB(M)$, is defined as the smallest axis-
aligned bounding box that contains the intersection $D=\cap_{D_{i}\in M}D_{i}$
, given by

$AABB(M):=\bigcap_{D\subset\tilde{B}}\tilde{B}$

where $\tilde{B}$ ranges over all axis-aligned bounding boxes that contain
$D$.

Note that the AABB can be expressed as

$AABB(M)=\prod_{k=1}^{d}[\inf\pi_{k}(\partial D),\sup\pi_{k}(\partial D)]$

where $\pi_{k}:\mathbb{R}^{d}\longrightarrow\mathbb{R}$ is the canonical
projection onto the $k$-th factor, and $\partial D$ denotes the boundary of
$D$. In other words, the AABB of a disk system $M$ is given by the Cartesian
product of intervals, where each interval is determined by the minimum and
maximum values of the corresponding projection of the disk boundaries.

###  4.1. Minimal axis-aligned bounding box for two disks

Let’s consider the situation when the disk system $M$ is composed of two disks
$D_{1}$ and $D_{2}$ in $\mathbb{R}^{d}$. If $D_{1}\cap D_{2}\neq\emptyset$ and
$D_{1}\neq D_{2}$, the subset $\partial D_{1}\cap\partial D_{2}$ can take one
of three forms: an empty set, a single common point (when the disks are
tangent), or a $(d-2)$-dimensional sphere. In the last case, we denote the
$(d-2)$-dimensional sphere or the $(d-1)$-sphere $\partial D_{1}\cap\partial
D_{2}$ by $S_{1,2}$. To calculate $\inf\pi_{i}(\partial(D_{1}\cap D_{2}))$ and
$\sup\pi_{i}(\partial(D_{1}\cap D_{2}))$ for each $i\in\\{1,2,...,d\\}$, we
can use:

(6) $\displaystyle\inf\pi_{i}(\partial(D_{1}\cap
D_{2}))=\begin{cases}\pi_{i}(c_{1}-r_{1}e_{i})&\text{if }c_{1}-r_{1}e_{i}\in
D_{2},\\\ \pi_{i}(c_{2}-r_{2}e_{i})&\text{if }c_{2}-r_{2}e_{i}\in D_{1},\\\
\inf(\pi_{i}(\partial D_{1}\cap\partial D_{2}))&\text{otherwise},\end{cases}$

$\displaystyle\sup\pi_{i}(\partial(D_{1}\cap
D_{2}))=\begin{cases}\pi_{i}(c_{1}+r_{1}e_{i})&\text{if }c_{1}+r_{1}e_{i}\in
D_{2},\\\ \pi_{i}(c_{2}+r_{2}e_{i})&\text{if }c_{2}+r_{2}e_{i}\in D_{1},\\\
\sup(\pi_{i}(\partial D_{1}\cap\partial D_{2}))&\text{otherwise}.\end{cases}$

Indeed, by Lemma 4, the extremes of the AABB are the projections of certain
poles, either from the $(d-1)$-sphere or some $d$-sphere. The $d$-spheres
represent the boundaries of each disk, with poles given by $c_{j}\pm
r_{j}e_{i}$, and the $(d-1)$-sphere is $S_{1,2}=\partial D_{1}\cap\partial
D_{2}$, whose poles are computed using Lemma 6. It is worth noting that there
are no further options for $(d-m+1)$-spheres in the case of a two-disk system.

To simplify the notation, we will use $B_{i,j}$ to denote the AABB of the
intersection of disks $D_{i}$ and $D_{j}$. Figure 2 illustrates the AABB of
the intersection of two disks in the plane.

Figure 2. AABB of two disks.

According to (6) and Lemma 5, we have a method to calculate the axis-aligned
boundary box (AABB) for systems of two disks.

Knowing how to compute the AABB of two disks is not sufficient to determine
the AABB for a disk system with more than two disks in $\mathbb{R}^{d}$. In
the following examples, we demonstrate that the AABB of a disk system is not
simply the intersection of all AABBs of two disks in $\mathbb{R}^{3}$.

######  Example 13.

Let

$M=\\{D_{1}((4,1,0);\sqrt{2}),D_{2}((4,-1,0);\sqrt{2}),D_{3}((0,0,0);3)\\}\subset\mathbb{R}^{3}$

be a Vietoris-Rips system in $\mathbb{R}^{3}$ with the following projection
onto the $xy$-plane:

Figure 3. Disk system $M$ projected onto the $xy$-plane.

By computing the boxes $B_{i,j}$ for all $i<j$, we obtain:

$\displaystyle B_{1,2}$
$\displaystyle=[3,5]\times[-\sqrt{2},\sqrt{2}]\times[-1,1]$ $\displaystyle
B_{1,3}$ $\displaystyle=[4-\sqrt{2},3]\times[0,1.41]\times[-0.72,0.72]$
$\displaystyle B_{2,3}$
$\displaystyle=[4-\sqrt{2},3]\times[-1.4,0]\times[-0.72,0.72]$

Therefore, $\cap_{i<j}B_{i,j}=[3,3]\times[0,0]\times[-0.72,0.72]$. However,
the disk system intersects at the point $P=(3,0,0)$, which means
$AABB(M)=\\{P\\}$. In other words, $AABB(M)$ is not equal to
$\cap_{i<j}B_{i,j}$.

We know that if $D\neq\emptyset$, then $\cap_{i<j}B_{i,j}\neq\emptyset$.
However, the converse is not always true. The following example illustrates
this fact.

######  Example 14.

Let
$N=\\{D_{1},D_{2},D_{3}((0,1,0);\sqrt{10}),D_{4}((3,0,1);0.9)\\}\subset\mathbb{R}^{3}$
be a Vietoris-Rips system. The projection of disks $D_{1},D_{2},$ and $D_{3}$
onto the $xy$-plane is illustrated in Figure 4.

Figure 4. Disk system $N$ projected onto the $xy$-plane.

We will now compute the intersections $B_{i,j}$ for different pairs of disks:

$\displaystyle B_{1,2}$
$\displaystyle=[3,5]\times[-\sqrt{2},\sqrt{2}]\times[-1,1]$ $\displaystyle
B_{1,3}$ $\displaystyle=[4-\sqrt{2},\sqrt{10}]\times[0,2]\times[-1,1]$
$\displaystyle B_{1,4}$
$\displaystyle=[2.6,3.9]\times[-0.4,0.9]\times[0.100007,1.3]$ $\displaystyle
B_{2,3}$ $\displaystyle=[2.6,3]\times[-0.8,0]\times[-0.44,4.47]$
$\displaystyle B_{2,4}$
$\displaystyle=[2.6,3.9]\times[-0.9,0.4]\times[0.100007,1.3]$ $\displaystyle
B_{3,4}$ $\displaystyle=[2.1,3.11]\times[-0.73,0.9]\times[0.1,1.7]$

The intersection of all pairwise intersections, $\cap_{i<j}B_{i,j}$, is given
by $[3,3]\times[0,0]\times[0.1,1.7]$. However, the intersection of disks
$D_{1},D_{2},$ and $D_{3}$ is a single point $P$, which is not contained in
$D_{4}$ (by construction). Therefore, $D$ is an empty set, but
$\cap_{i<j}B_{i,j}$ is not.

These examples clearly illustrate that when dealing with the AABB of three
disks or more, knowing the AABB for pairs of disks is insufficient. In Example
13, we observe that the intersection of $AABB(\\{D_{i},D_{j}\\})$ is not equal
to the AABB of the disk system $M$. Similarly, in Example 14, we find that the
intersection of three disks is empty, yet the intersection of
$AABB(\\{D_{i},D_{j}\\})$ contains points.

Therefore, the next crucial step is to determine how to calculate the AABB of
a disk system consisting of more than two disks in $\mathbb{R}^{d}$.

###  4.2. Minimal axis-aligned bounding box for more than two disks

Given a system $M$ consisting of $m$ disks in $\mathbb{R}^{d}$, we can compute
the $e_{q}$-poles for any subcollection of disks. Using these $e_{q}$-poles,
we can determine the axis-aligned bounding box (AABB) of $M$.

If $m\leq d$, we calculate the $e_{q}$-poles of the $(d-m+1)$-sphere
$\cap\partial D_{i}$, and with these poles, we define the AABB of $M$ by
taking $\inf\pi_{q}\partial D=\pi_{q}(p)$, where $p$ is the $e_{q}$-south pole
of the $(d-m+1)$-sphere (similarly for $\sup\pi_{q}\partial D$).

In the case where $m=d+1$, we consider $AABB(M(1)),...,AABB(M(m))$ as a
collection of minimal axis-aligned boxes for the disk system
$M(i)=M-\\{D_{i}\\}$ in $\mathbb{R}^{d}$. If the intersection of any $d+1$ of
these sets is nonempty, then the intersection of the entire collection gives
us the minimal axis-aligned box of the disk system $M$. This can be expressed
as $AABB(M)=\cap AABB(M(i))$. The next theorem confirms this finding.

######  Theorem 15 (Helly’s theorem for minimal axis-aligned boxes).

Let $M$ be a Rips system with $d+1$ disks in $\mathbb{R}^{d}$. Then, the
minimal axis-aligned bounding box for the intersection set
$D=\cap_{j=1}^{d+1}D_{j}$ satisfies:

$AABB(M)=\bigcap_{j=1}^{d+1}AABB(M(j))$

where $M(j)=M-\\{D_{j}\\}$.

###### Proof.

We know that $AABB(M)\subseteq\bigcap_{j=1}^{d+1}AABB(M(j))$. Now, our
objective is to establish the reverse inclusion, that is,
$\bigcap_{j=1}^{d+1}AABB(M(j))\subseteq AABB(M)$. In order to derive a
contradiction, suppose that the reverse inclusion is not true. By definition,
we have:

$AABB(M)=\prod_{i=1}^{d}\left[\inf\pi_{i}(\partial D),\sup\pi_{i}(\partial
D)\right]$

and

$\bigcap_{j=1}^{d+1}AABB(M(j))=\prod_{i=1}^{d}\left[\max_{k}\\{\inf\pi_{i}(\partial\cap_{j\neq
k}D_{j})\\},\min_{k}\\{\sup\pi_{i}(\partial\cap_{j\neq k}D_{j})\\}\right]$

Without loss of generality, let’s assume that:

$\pi_{1}(p)>\max_{k}\left\\{\inf\pi_{1}(\cap_{j\neq k}\partial
D_{j})\right\\}_{k=1}^{d+1}$

where $p\in D$ satisfies $\pi_{1}(p)=\inf\pi_{1}(\partial D)$.

Let $q_{j}$ be the point in $\cap_{k\neq j}D_{k}$ such that
$\pi_{1}(q_{j})=\inf\pi_{1}(\cap_{k\neq j}D_{k})$ for each $j=1,2,...,d+1$.
Note that $q_{j}\notin D_{j}$ because $q_{j}$ is not in $D$. Now, let
$\gamma_{j}$ be the line segment that connects $q_{j}$ and $p$.

Choose $\epsilon>0$ small enough such that the hyperplane
$P:x_{1}=\pi_{1}(p)-\epsilon$ does not contain any $q_{j}$, and
$\pi_{1}(q_{j})<\pi_{1}(p)-\epsilon$ for all $j=1,2,...,d+1$ (such hyperplane
$P$ exists because $\pi_{1}(p)>\pi_{1}(q_{j})$ for all $j=1,...,d+1$). Since
$q_{j}$ is in $\cap_{k\neq j}D_{k}$, the hyperplane $P$ intersects every disk
in $M(j)$ for each $j$, and $\gamma_{j}\subset\cap_{k\neq j}D_{k}$ intersects
$P$ at a point that is in $\cap_{k\neq j}D_{k}$ (see Figure 5). Furthermore,
$D_{k}\cap P$ is a $(d-1)$-dimensional disk. Therefore, we have a collection
$\mathcal{D}=\\{D_{j}\cap P|j=1,2,...,d+1\\}$ of $d+1$ disks, each of
dimension $d-1$, such that every subset $A$ of $\mathcal{D}$ consisting of $d$
disks has the non-empty intersection property. By Helly’s Theorem, the
intersection of all $(d-1)$-disks in $\mathcal{D}$ is not empty. Therefore,
there exists a point $q\in D$ with $\pi_{1}(q)<\pi_{1}(p)$. However, this
contradicts the fact that $p\in D$ is such that
$\pi_{1}(p)=\inf\pi_{1}(\partial D)$.

Therefore, we conclude that $\bigcap_{j=1}^{d+1}AABB(M(j))=AABB(M)$.

Figure 5.

∎

######  Lemma 16.

Let $M$ be a Vietoris-Rips system of $d+1$ disks in $\mathbb{R}^{d}$. If
$D=\cap_{i=1}^{d+1}D_{i}=\emptyset$ and $AABB(M(j))\neq\emptyset$ for each
$j=1,...,d+1$, then $\cap_{i=1}^{d+1}AABB(M(i))$ consists only of intervals of
the form $[a,b]$ with $a>b$ (inverted intervals).

###### Proof.

Suppose that $\cap_{i=1}^{d+1}AABB(M(i))$ contains an interval. Without loss
of generality, let’s assume it is the interval

$[\inf\pi_{1}(\partial D(k)),\sup\pi_{1}(\partial D(l))]$

for some $k,l\in\\{1,2,...,d+1\\}$ where $D(k)=\cap_{i\neq k}D_{i}$. By
definition of $\cap_{i=1}^{d+1}AABB(M(i))$, we have $\inf\pi_{1}(\partial
D(k))\geq\inf\pi_{1}(\partial D(i))$ for all $i\neq k$ and
$\sup\pi_{1}(\partial D(l))\leq\sup\pi_{1}(\partial D(j))$ for all $j\neq l$.
Now, consider a hyperplane $P:x_{1}=p$ with $\inf\pi_{1}(\partial D(k))\leq
p\leq\sup\pi_{1}(\partial D(k))$. The hyperplane $P$ intersects $D(j)$ for
every $j=1,...,d+1$ (since $P$ cuts through all the boxes $AABB(M(j))$). We
also know that $P\cap D_{i}$ is a $d-1$-dimensional disk for each $i$.

Therefore, we have a collection of $d$ disks in a $d-1$ dimensional space, and
by Helly’s Theorem, this collection must have a non-empty intersection. This
implies that $D\neq\emptyset$, which contradicts our assumption that $D$ is
empty. Hence, all intervals in $\cap_{i=1}^{d+1}AABB(M(i))$ must be inverted
intervals of the form $[a,b]$ with $a>b$. ∎

As an example, we provide the lower bounds of the AABB for a system of three
disks in $\mathbb{R}^{d}$, this is, we compute $\inf\pi_{i}(\partial D)$ for
each $i=1,...,d$ (analogous $\sup\pi_{i}\partial D$) with
$D=\cap_{j=1}^{3}D_{j}$. Let $M=\\{D_{1},D_{2},D_{3}\\}$ be a Rips system in
$\mathbb{R}^{d}$, then, for each $i=1,...,d$, the computation of the AABB for
the disk system $M$ is given by:

$\displaystyle\inf\pi_{i}(\partial
D)=\begin{cases}\pi_{i}(c_{k}-r_{k}e_{i})&\text{if $c_{k}-r_{k}e_{i}\in D$}\\\
\max_{j<k}\\{\inf\pi_{i}(\partial D_{j}\cap\partial D_{k})\\}&\text{if (*)}\\\
\inf\pi_{i}(\cap_{j=1}^{3}\partial D_{j})&\text{otherwise}\end{cases}$

where (*) denotes the case where $q\in D$ with $q\in\partial D_{j}\cap\partial
D_{k}$ such that

$\pi_{i}(q)=\max_{j<k}\\{\inf\pi_{i}(\partial D_{j}\cap\partial D_{k})\\}.$

The preceding calculations are a result of Lemma 4. It is known that the
extremes of the AABB lie in the projections of specific poles, either from a
$d$-sphere, $(d-1)$-sphere, or the $(d-2)$-sphere.

Given a disk system $M$ in $\mathbb{R}^{d}$, we can compute the minimal axis-
aligned bounding box (AABB) for the intersection of all disks in $M$. If the
AABB is a point, then it represents the intersection of the disks. This
property allows us to identify when the AABB is a point. If the AABB of $M$ is
not a point, we can rescale $M$ by a scale factor $\lambda$ such that the AABB
of $M_{\lambda}$ becomes a point. The value of $\lambda$ is referred to as the
Čech scale of the system $M$.

###  4.3. Minimal axis-aligned bounding box Algorithm.

Now, we present an algorithm to calculate the minimal axis-aligned bounding
box (AABB) of a system of $m$ disks in $\mathbb{R}^{d}$. As previously
explained, the AABB’s extremes are defined by projecting the poles of specific
$i$-spheres. Thus, in computing the AABB, we will determine the poles for each
$i$-sphere within the disk system.

1

Input : A $d$-disk system $M=\\{D_{j}\\}_{j=1}^{m}$

Output : The minimal AABB for the system $M$

2 Initialize: $P=\emptyset$

3 for  _$k\leftarrow 1$ to $m$_ do

4 Let $\mathcal{S}$ be the set of $(d-k+1)$-spheres of $\partial M$

5 for  _$S$ in $\mathcal{S}$_ do

6 for  _$q\leftarrow 1$ to $d$_ do

7 Compute the set $\\{s_{q}^{\pm}\\}$ of $e_{q}$-poles of $S$

8 for  _$s\leftarrow\\{s_{q}^{\pm}\\}$_ do

9 if  _$s\in\cap_{j=1}^{d}D_{j}$_ then

10 Add: $P\leftarrow s$

11

12 end if

13

14 end for

15

16 end for

17

18 end for

19

20 end for

21if  _$P\neq\emptyset$_ then

22 for  _$q\leftarrow 1$ to $d$_ do

23 $a_{q}=min_{P}\\{\pi_{q}(s_{q}^{-})\\}$

24 $b_{q}=max_{P}\\{\pi_{q}(s_{q}^{+})\\}$

25

26 end for

27 return ($\Pi_{q=1}^{d}[a_{q},b_{q}]$)

28else

29 return (The disk system does not intersect)

30 end if

Algorithm 3 AABB.minimal

The algorithm starts by initializing the set of poles, $P$, as an empty set.
In each iteration of the first loop, we determine the spheres formed by the
intersection of the boundaries of the disk subcollections in the system $M$
(for this, we require the center and radius, which are computed in Subsection
2.2). Next, we identify the poles of each sphere and if any of them are
present in all the disks of $M$, we add them to the set $P$. If $P$ is not
empty, we proceed to calculate the extremes for each dimension of the AABB
using the set of north poles for the upper bounds and the set of south poles
for the lower bounds. In this case, the output is the product of the intervals
defined by the computed extremes. If $P$ is empty, it indicates that there is
no intersection in the disk system.

##  5\. Acknowledgements

C.G.E.P acknowledges CONACYT for the financial support provided through a
National Fellowship (CVU-638165).

## References

  * [1] G. Bell, A. Lawson, J. Martin, J. Rudzinski, and C. Smyth, Weighted persistent homology, Involve, a Journal of Mathematics, 12 (2019), pp. 823–837.
  * [2] P. Cai, Y. Cai, I. Chandrasekaran, and J. Zheng, Collision detection using axis aligned bounding boxes, Simulation, Serious Games and Their Applications, (2013), pp. 1–14.
  * [3] J. F. Espinoza, R. Hernández-Amador, H. A. Hernández-Hernández, and B. Ramonetti-Valencia, A numerical approach for the filtered generalized čech complex, Algorithms, 13 (2019), p. 11.
  * [4] N. K. Le, P. Martins, L. Decreusefond, and A. Vergne, Construction of the generalized cech complex, 2015.
  * [5] J. Mahovsky and B. Wyvill, Fast ray-axis aligned bounding box overlap tests with plucker coordinates, Journal of Graphics Tools, 9 (2004), pp. 35–46.
  * [6] D. S. Maioli, C. Lavor, and D. S. Gonçalves, A note on computing the intersection of spheres in $\mathbb{R}^{n}$, The ANZIAM Journal, 59 (2017), p. 271–279.