quadtree during the search and did not arrive at the level of $\sigma_{s}$, we
add its four children to the former leaf. We then step down into the newly
created leaf that contains $\sigma_{s}$ and repeat the process until we added
$\sigma_{s}$.

Note that it is possible that even the root of the quadtree containing
$\sigma_{s}$ does not exist. In this case, we insert the root cell into the
balanced binary search tree and create a root node before processing the way
down as described before. As the insertion into the search tree takes
$O\mathopen{}\mathclose{{}\left(\log n}\right)$ time just as in the look-up
step, this does not increase the running time. Additionally, the operation of
searching for $\sigma_{s}$ or accordingly building the quadtree downwards to
the level of $\sigma_{s}$ requires at most
$O\mathopen{}\mathclose{{}\left(\lfloor\log\Psi\rfloor}\right)$ time which is
the maximum height of the quadtrees in the forest.

To complete the insertion of $s$ in the forest component, we also have to
introduce the cells of the neighborhood $N_{15\times 15}(\sigma_{s})$ if not
present yet. Here, the look-up operation as well as the insertion, if
required, are executed in a similar manner as for $\sigma_{s}$. Since the
neighborhood is of constant size, the overall processing time of the quadtree
forest during the insertion of $s$ is in $O\mathopen{}\mathclose{{}\left(\log
n+\log\Psi}\right)$.

When introducing cells into the forest, we also need their associated regions
to be present in the proxy graph $H$. Hence, we insert the regions of the new
cells as isolated vertices to $\mathcal{H}$ and initialize the corresponding
AWNN of the regions for $S_{1}(A)$ and $\overline{S_{2}(A)}$. Subsequently, we
have to insert $s$ itself and its edges as follows: first, we add $s$ as an
isolated vertex to $\mathcal{H}$. Then, we insert $s$ to the sets containing
it and if applicable the associated AWNN. Since we implicitly store $S_{1}(A)$
and $S_{2}(A)$, we add the corresponding edges between the regions and $s$ to
$\mathcal{H}$ along the way.

Recall, that by Lemma 6.1, $s$ can only be in a set $S_{1}(A)$ defined by a
cell in $N(s)$. Therefore, we iterate through the cells of $N(s)$ and add $s$
with weight $-r_{s}$ to the AWNN of the sets $S_{1}(A)$ it belongs to. We also
introduce the respective edges $(s,A)$ to the proxy graph. As $N(s)$ is of
constant size this step takes $O(1)$ amortized time for the insertion of edges
in $\mathcal{H}$ and $O(\lambda_{6}(\log n)\log^{5}n)$ expected time for
insertions into the AWNN.

Furthermore, for a region $A$ changes in a set $S_{1}(A)$ may affect the set
$S_{2}(A)$. Hence, we have to check if the insertion of $s$ into $S_{1}(A)$
yields a transfer of a site $t$ from $\overline{S_{2}(A)}$ to $S_{2}(A)$
because $D_{s}$ intersects $D_{t}$. We identify these sites by iteratively
doing a weighted nearest neighbor search with $s$ on the AWNN on
$\overline{S_{2}(A)}$ until we do not find a site intersecting $s$ anymore.
When finding a site $t$ intersecting $s$, we delete $t$ from the AWNN on
$\overline{S_{2}(A)}$ and introduce the edge $(t,A)$ into $\mathcal{H}$. The
amortized expected running time of $O(\lambda_{6}(\log n)\log^{9}n)$ to delete
a site from $\overline{S_{2}(A)}$ dominates this step. We do not know how many
deletions are required in this step, but a site $t$ can be present in the most
$O(\log\Psi)$ sets $\overline{S_{2}(A)}$, associated to the cells containing
$t$ on the path from the root to $\sigma_{t}$. Furthermore, a site $t$ is only
inserted into $\overline{S_{2}(A)}$, when $t$ is inserted. Thus, these
deletions take overall $O(n\log\Psi\lambda_{6}(\log n)\log^{9}n)$ expected
time, which yields an amortized expected running time of
$O(\log\Psi\lambda_{6}(\log n)\log^{9}n)$ per insertion.

Now, we consider the time needed to insert of $s$ either into $S_{2}(A)$ or
$\overline{S_{2}(A)}$. By definition, we know that $s$ has to be added either
to $S_{2}(A)$ or to $\overline{S_{2}(A)}$ of the regions of the cells along
the path from the root to $\sigma_{s}$ in the quadtree. Let $A$ be a region
defined by a cell along this path. To determine whether to add $s$ to
$S_{2}(A)$ or to $\overline{S_{2}(A)}$, we perform a nearest neighbor search
on the AWNN of $S_{1}(A)$. If we find a site in $S_{1}(A)$ intersecting $s$,
we know that $s$ belongs to $S_{2}(A)$. Then, we insert the edge $(s,A)$ to
$\mathcal{H}$, if it does not exist yet. This takes $O(1)$ amortized time. In
the other case, $s$ is inserted to the AWNN of $\overline{S_{2}(A)}$ in
$O(\lambda_{6}(\log n)\log^{5}n)$ amortized expected time. This last step
concludes the insertion. Since each cell defines a constant number of regions,
the step requires a constant number of look-up and insertion operations on
each level. As the time for insertions into the AWNN of $\overline{S_{2}(A)}$
dominates, the last insertion step yields an amortized expected running time
of $O(\log\Psi\lambda_{6}(\log n)\log^{5}n)$

By the summing over all insertion steps, we achieve an amortized expected
running time of $O(\log\Psi\lambda_{6}(\log n)\log^{9}n)$ per insertion
because of the dominating running time of the deletions from
$\overline{S_{2}(A)}$.

It remains to consider the running time for the query which directly follows
from the query time in $\mathcal{H}$. This concludes our proof.

##  7 Arbitrary Radius Ratio

In this section, we extend the approach from Section 6 to give a data
structure that allows update times independent of $\Psi$. The cost for
dropping the dependence on $\Psi$ is exchanging the additive
$O\mathopen{}\mathclose{{}\left(n\log\Psi\log^{4}n}\right)$ term in the
running time in Theorem 6.13 with an additional $O(\log n)$ factor in the
first term. The additive
$O\mathopen{}\mathclose{{}\left(n\log\Psi\log^{4}n}\right)$ part in Section 6
came from the size of the sets $S_{2}(A)$ and thus from the height of each
quadtree in the forest of quadtrees. We can get rid of this dependency by
using a compressed quadtree $\mathcal{Q}$. Recall from Section 2 that the
height and size of the compressed quadtree do not depend on the diameter of
the point set, but only on $n$. However the height of $\mathcal{Q}$ can still
be $O(n)$ which is not favorable for our application. In order to reduce the
number of edge in our proxy graph to
$O\mathopen{}\mathclose{{}\left(\log^{2}n}\right)$, we consider a heavy path
decomposition on $\mathcal{Q}$ in combination with binary search trees for
each heavy path. The new proxy graph is described in Section 7.1 and the
decremental connectivity data structure based on this graph can be found in
Section 7.2.

###  7.1 The Proxy Graph

The main structure of the proxy graph is the same as that described in Section
6.1, so we will often refer back to it. We still have a bipartite graph with
the sites from $S$ on one side and a set of region on the other side. The
regions will again be used to define sets $S_{1}(A)$ and $S_{2}(A)$ that will
in turn define the edges. However we adapt the definition of the regions and
define them based on certain paths of the compressed quadtree instead of
single cells. Furthermore, we relax the condition on the radii in the
definition of the set $S_{1}(A)$.

Similar to the non-compressed approach, we consider an augmented version of
the quadtree. Again let $N(s)$ be the $15\times 15$ neighborhood of the cell
$\sigma$ containing $s$ such that $|\sigma|\leq r_{s}<2|\sigma|$. Instead of
working with the compressed quadtree on the sites, we work on the tree
$\mathcal{Q}$ that also contains all cells which are part of at least one
neighborhood.

Now let $\mathcal{R}$ be the heavy path decomposition of $\mathcal{Q}$. For
each heavy path $R\in\mathcal{R}$, we find a set $\mathcal{P}_{R}$ of
_canonical paths_ such that every subpath of $P$ can be written as the
disjoint union of $O(\log n)$ canonical paths. To be precise, for each
$R\in\mathcal{R}$ we build a balanced binary search tree with the cells of the
path in the leafs, sorted by increasing diameter. We associate each vertex $v$
in the binary search tree with the path induced by all cells in the subtree
rooted at $v$ and add this path to $\mathcal{P}_{R}$. Using this definition,
we can now write every path in $\mathcal{Q}$ that starts at the root as the
disjoint union of canonical path, as shown in the following lemma:

######  Lemma 7.1.

Let $\sigma$ be a vertex of $\mathcal{Q}$ and let $\pi$ be the path from the
root of $\mathcal{Q}$ to $\sigma$. Then there exists a set $\mathcal{P}_{\pi}$
of canonical paths such that:

  1. 1.

$|\mathcal{P}_{\pi}|=O\mathopen{}\mathclose{{}\left(\log^{2}n}\right)$; and

  2. 2.

$\pi$ is the disjoint union of the canonical paths in $\mathcal{P}_{\pi}$.

######  Proof 7.2.

Figure 19: Illustration of Lemma 7.1. Left we see the decomposition of $R$
into $R_{1},\dots R_{k}$. On the right the vertices defining
$\mathcal{P}_{\pi}$ are depicted in green.

Consider the heavy paths $R_{1},\dots,R_{k}$ encountered along $\pi$. By Lemma
2.4 we have that $k\in O(\log n)$ and that $\pi$ is the disjoint union of the
intersections $\pi\cap R_{i}$. Each of these intersections constitutes a
subpath of $R_{i}$ whose largest cell is also the largest cell of $R_{i}$. Let
$\sigma_{i}$ be the smallest cell of $\pi\cap R_{i}$, then the subpath of
$R_{i}$ consists of all cells in $R_{i}$ having diameter at least
$|\sigma_{i}|$. This subpath can be composed as the disjoint union of all
canonical paths defined by the right children along the search path for
$\sigma_{i}$ in the binary search tree associated with $R_{i}$, together with
the path only consisting of $\sigma_{i}$; see Figure 19. As the height of a
balanced binary search tree is $O(\log n)$, the overall number of canonical
paths needed to represent $\pi$ is $O(\log^{2}n)$.

The vertex set of the proxy graph $H$ again consists of the set of sites $S$
and a set of regions $\mathcal{A}$. We define $O(1)$ regions for each
canonical path in a similar way as we defined the regions in Section 6.1. Let
$\sigma$ be the smallest cell and $\tau$ be the largest cell of a canonical
path. The _inner_ and _middle regions_ are defined as in the bounded case,
using the smallest cell of the path to define the region, instead of the
single cell in Section 6.1. That is the inner region is the disk with center
$a(\sigma)$ and radius $|\sigma|$. The middle regions are the $d_{2}$ regions
which are defined as the intersection of the cones in $\mathcal{C}_{d_{2}}$
with the annulus with inner radius $|\sigma|$ and outer radius
$\frac{5}{2}{|\sigma|}$. For the _outer regions_ , we extend the outer radius
of the annulus. They are defined as the intersections of the cones in
$\mathcal{C}_{d_{1}}$ with the annulus of inner radius $\frac{5}{2}{|\sigma|}$
and outer radius $\frac{5}{2}|\sigma|+2|\tau|$, again centered at $a(\sigma)$.
The set $\mathcal{A}$ now contains the regions defined this way for all
canonical paths.

Given a region $A\in\mathcal{A}$ with smallest cell $\sigma$ and largest cell
$\tau$, we can now define the sets $S_{1}(A)$ and $S_{2}(A)$. The set
$S_{1}(A)$ is defined similarly to the set in Section 6.1, again using
$\sigma$ in the role of the single cell for most parts. The difference is that
the radius range for $t$ is larger, as it’s upper bound depends on the largest
cell of the path. The set $S_{1}(A)$ contains all sites $t$ such that

  1. 1.

$t\in A$

  2. 2.

$|\sigma|\leq r_{t}\leq 2|\tau|$; and

  3. 3.

$\|a(\sigma)t\|\leq r_{t}+\frac{5}{2}|\sigma|$.

The last condition is only relevant, if $A$ is an outer region, as it is
trivially true for middle and inner regions.

The definition for $S_{2}(A)$ is also similar to that in Section 6.1, however
the difference in its definition are larger than for the set $S_{1}(A)$. Let
$\sigma_{s}$ be the the cell such that $s\in\sigma_{s}$ and $|\sigma_{s}|\leq
r_{s}<2|\sigma_{s}|$, let $\pi_{s}$ be the path in $\mathcal{Q}$ from the root
to $\sigma_{s}$ and let $\mathcal{P}_{\pi_{s}}$ be the decomposition of
$\pi_{s}$ into canonical paths as defined in Lemma 7.1. Let $A$ be a region,
defined by a canonical path $P$, then $s\in S_{2}(A)$ if

  1. 1.

$P\in\mathcal{P}_{\pi_{s}}$; and

  2. 2.

$s$ is adjacent in $\mathcal{D}(S)$ to at least one site in $S_{1}(A)$.

If $\sigma$ is the smallest cell in a canonical path defining a region $A$,
then each site $s$ in $S_{2}(A)$ lies in $\sigma$, has a radius less than
$2|\sigma|$ and intersects at least one cell in $S_{2}(A)$. These are
basically the conditions we had in Section 6.1. However, as the definition is
restricted to those canonical paths in $\mathcal{P}_{\pi_{s}}$, not all sites
satisfying these conditions are considered. As we will see below, this
suffices to make sure that the proxy graph represents the connectivity, while
also ensuring that each site $s$ lies in few sets $S_{2}(A)$.

The graph $H$ is now again defined by connecting a site $s\in S_{1}(A)\cup
S_{2}(A)$ to the region $A$. To show that $H$ accurately represents the
connectivity in $D(S)$, we need the following corollary from Lemma 6.3.

######  Corollary 7.3.

Suppose that $d_{1}\geq 23$ and $d_{2}\geq 8$. Then, for any region
$A\in\mathcal{A}$, the associated sites in $S_{1}(A)$ form a clique in
$\mathcal{D}(S)$.

######  Proof 7.4.

Recall that the center of the annuli and disks defining $A$ is the smallest
cell of the associated canonical path. A close inspection of the proof for
Lemma 6.3 shows that we only use the lower bound on the radii of the sites in
$S_{1}(A)$. As this lower bound is unchanged, all arguments carry over for
sites with larger radii.

######  Lemma 7.5.

Two sites $s$ and $t$ are connected in $\mathcal{D}(S)$ if and only if they
are connected in $H$.

######  Proof 7.6.

If $s$ and $t$ are connected in $H$, the same argument as in the proof of
Lemma 6.5 with Corollary 7.3 instead of Lemma 6.3 holds. The more challenging
part is to show that if two sites are connected in $\mathcal{D}(S)$ they are
also connected in $H$.

Figure 20: The cell $\sigma_{s}$ is the smallest cell such that $r_{s}\leq
2|\sigma_{s}|$. The canonical path $P$ contains $\rho$.

If suffices to show that if $s$ and $t$ are connected by an edge in
$\mathcal{D}(S)$, they are connected to the same region $A\in\mathcal{A}$.
Refer to Figure 20 for a depiction of the following argument. Assume without
loss of generality that $r_{s}\leq r_{t}$. Consider the neighborhood $N(t)$ of
$t$. Then by Lemma 6.1 there is a cell $\rho\in N(t)$ that contains $s$.
Consider the path $\pi_{s}$ in $\mathcal{Q}$ from the root to the cell
$\sigma_{s}$ such that $s\in\sigma_{s}$ and $|\sigma_{s}|\leq
r_{s}<2|\sigma_{s}|$. The cell $\sigma_{s}$ is in $\mathcal{Q}$ as it is part
of $N(s)$. Then $\rho$ lies on this path $\pi_{s}$. Let
$\mathcal{P}_{\pi_{s}}$ be the decomposition of $\pi_{s}$ into canonical paths
as defined in 7.1 and let $P$ be the path containing $\rho$. Let $\sigma$ and
$\tau$ be the smallest and largest cell on $P$ respectively. By the definition
of $P$, we have $\sigma_{s}\subseteq\sigma\subseteq\rho\subseteq\tau$. As
$\\{s,t\\}$ is an edge in $D(S)$, we have $\|st\|\leq r_{s}+r_{t}\leq
2|\sigma|+2|\tau|$ and thus $\|a(\sigma)t\|\leq\frac{5}{2}|\sigma|+2|\tau|$.
This implies that $t$ lies in a region $A$ defined by $P$ and thus $t\in
S_{1}(A)$ for this region. As $s$ intersects $t$, it intersects at least one
site in $S_{1}(A)$ for the region $A$. Furthermore, $P$ is a canonical path in
$\mathcal{P}_{\pi_{s}}$ and thus $s\in S_{2}(A)$ by definition. As $s,t\in
S_{1}(A)\cup S_{2}(A)$, they are connected in $H$.

######  Lemma 7.7.

The graph $H$ has $O(n)$ vertices and
$O\mathopen{}\mathclose{{}\left(n\log^{2}n}\right)$ edges.

######  Proof 7.8.

As discussed in Section 2 the compressed quadtree consists of $O(n)$ cells.
Each cell is part of exactly one heavy path, so the total size of the binary
search trees that define the canonical paths is $O(n)$. A balanced binary
search tree with $n$ leafs has $O(n)$ inner vertices, and thus there is a
total of $O(n)$ canonical paths. As each canonical path defines $O(1)$
regions, the number of regions and therefore also the number of vertices in
$H$ follows.

To bound the number of edges, we again count the number of sets $S_{1}(A)$ and
$S_{2}(A)$ a single sites can be contained in. Again let $\sigma$ and $\tau$
be the smallest and largest cell of a canonical path respectively. A site $t$
can only be in a set $S_{1}(A)$ if $|\sigma|\leq r_{t}\leq 2|\tau|$ and
$\sigma$ is contained in a cell of $N(t)$. To see the second claim, let $P$ be
a canonical path such that $t\in S_{1}(A)$ for a region $A$ defined by $P$ and
let $\sigma$ be the smallest cell of $P$. If $|\sigma|\leq r_{t}<2|\sigma|$,
the statement holds by the second property given in Lemma 6.1. In the other
case, let $\rho$ be the cell of the hierarchical grid with
$|\sigma|<|\rho|\leq r_{t}<2|\rho|$ and $\sigma\subseteq\rho$. As $t\in
S_{1}(A)$, we have

$\displaystyle\|ta(\sigma)\|$ $\displaystyle\leq r_{t}+\frac{5}{2}|\sigma|$
$\displaystyle\leq 2|\rho|+\frac{5}{4}|\sigma|$ $\displaystyle\text{since
}r_{t}\leq 2|\rho|\text{ and }|\sigma|\leq 2|\rho|$
$\displaystyle\leq\frac{13}{4}|\rho|$ $\displaystyle\|ta(\rho)\|$
$\displaystyle\leq\frac{15}{4}$ $\displaystyle\text{by triangle inequality as
}\sigma\subseteq\rho$

By the same considerations as in the proof of Lemma 6.1, it follows that
$\rho\in N(t)$ and thus $\rho$ is in $\mathcal{Q}$ and also in $P$.

Let $\rho$ be a cell of $N(t)$, and let $R$ be the heavy path containing
$\rho$. Then $t$ is considered for the set $S_{1}(A)$ for all regions which
are defined by a canonical path containing $\rho$. These are the $O(\log n)$
paths along the search path for $\rho$ in the binary search tree on the cells
of $R$ and thus each site can be part of at most $O(\log n)$ sets $S_{1}(A)$.

The number of sets $S_{2}(A)$, a fixed site $s$ can be part of, is at most the
number of canonical paths decomposing the path from the root to the smallest
cell containing $s$. By Lemma 7.1 there are
$O\mathopen{}\mathclose{{}\left(\log^{2}n}\right)$ such paths. This yields
$\sum_{A\in\mathcal{A}}|S_{1}(A)|+|S_{2}(A)|=O\mathopen{}\mathclose{{}\left(n\log^{2}n}\right)$
as an upper bound for the number of edges.

###  7.2 Decremental Data Structure

The structure of the decremental data structure is the same as in Section 6.2,
recall Figure 17. We again store the sets $S_{1}(A)$ and $S_{2}(A)$ together
with a RDS for each region $A$. The set $B$ for the RDS is again $S_{1}(A)$
and the set $P$ is $S_{2}(A)$. Furthermore, we store the graph $H$ defined by
$S_{1}(A)$ and $S_{2}(A)$ in a Holm et al. data structure $\mathcal{H}$.

Both queries and deletions work exactly as in Section 6.2, however we repeat
them here for completeness. Queries are performed on $\mathcal{H}$. On the
deletion of a site $s$, we first remove all edges incident to $s$ from
$\mathcal{H}$. Then the site is removed from all sets $S_{1}(A)$ containing
it, as well as the associated RDS. The sites $U$ reported as revealed by the
RDS are then removed from the matching $S_{2}(A)$ and the edges $\\{u,A\\}$
for $u\in U\setminus S_{1}(A)$ are removed from $\mathcal{H}$. Finally the
site $s$ is removed from all sets $S_{2}(A)$ and the matching RDS.

######  Lemma 7.9.

The data structure described above can be preprocessed in
$O(n\log^{6}n\lambda_{6}(\log n))$ time.

######  Proof 7.10.

To find the regions $\mathcal{A}$ we first compute the extended compressed
quadtree $\mathcal{Q}$. This can be done by adding $O(n)$ virtual sites to our
site set, similar to a construction of Har-Peled [16]. For each site $s\in S$
and each cell $\sigma\in N(s)$ we add two virtual sites. The virtual sites are
added at the center of two of the cells one level below $\sigma$ that are
contained in $\sigma$, see Figure 21 for an illustration. Now, all cells in
$N(s)$ have at least two children in the non-compressed quadtree and thus the
cells are also present in the compressed quadtree. When having $\mathcal{Q}$
at hand, we can, with additional $O(n)$, time find the heavy paths by Item 3.
of Lemma 2.4, and the binary search trees on the heavy paths by standard
techniques. This gives us the set of regions $\mathcal{A}$.

Figure 21: Four cells from $N_{15\times 15}(\sigma)$ with the virtual sites

To find the set $S_{1}(A)$, recall that a site $t$ can only lie in the set
$S_{1}(A)$ for a region defined by a canonical path that contains a cell in
$N(t)$. The sets $S_{1}(S)$ we can now be found as follows. For each site
$t\in S$, compute the cells in $N(t)$ and for each such cell $\rho$ find the
heavy path $R$ containing it. In the binary search tree defined on $R$, follow
the search path for $\rho$ and for each canonical path defined along this
search path explicitly find the region containing $t$ and check if the
distance condition holds.

When implementing this step naively it takes
$O\mathopen{}\mathclose{{}\left(n\log^{2}n}\right)$ time which is fast enough
for our purposes. We could however reduce this time to $O(n\log n)$ if we
added an additional preprocessing step on $\mathcal{Q}$. As in Lemma 6.11, we
construct a static additively weighted nearest neighbor data structure on each
set $S_{1}(A)$, again assigning the weight $-r_{s}$ to each site $s\in S$.
Recall that the time needed to construct a single nearest neighbor data
structure is $O(|S_{1}(A)|\log n)$ (Lemma 2.6). As
$\sum_{A\in\mathcal{A}}|S_{1}(A)|=O(n\log n)$, the construction of all data
structures takes $O\mathopen{}\mathclose{{}\left(n\log^{2}(n)}\right)$ time.
The query time however remains $O(\log n)$ in each data structure, as each
data structure contains at most $O(n)$ sites.

Recall from the definition of the sets $S_{2}(A)$, that $\pi_{s}$ is the path
in $\mathcal{Q}$ from the root to the cell $\sigma_{s}$ that contains $s$ and
has diameter $|\sigma|\leq r_{s}<2|\sigma|$. To find all sets $S_{2}(A)$
containing $s$, we simply follow the decomposition of $\pi_{s}$ into canonical
paths and query the nearest neighbor data structure with $s$ for all regions
defined by these paths. As there are
$O\mathopen{}\mathclose{{}\left(\log^{2}n}\right)$ canonical path in the
decomposition, this takes an additional
$O\mathopen{}\mathclose{{}\left(n\log^{3}n}\right)$ time for all sites.
Inserting the $O\mathopen{}\mathclose{{}\left(n\log^{2}n}\right)$ edges into
$\mathcal{H}$ takes $O\mathopen{}\mathclose{{}\left(\log^{2}n}\right)$
amortized time each, for a total of
$O\mathopen{}\mathclose{{}\left(n\log^{4}n}\right)$.

Again, the step dominating the preprocessing time is the construction of the
RDS. For a single region this time is
$O\mathopen{}\mathclose{{}\left(|S_{1}(A)|\log^{5}n\lambda_{6}(\log
n)+|S_{2}(A)|\log^{3}(n)}\right)$ by Corollary 5.12. We have
$\sum_{A\in\mathcal{A}}|S_{1}(A)|=O(n\log n)$ and
$\sum_{A\in\mathcal{A}}|S_{2}(A)|=O\mathopen{}\mathclose{{}\left(n\log^{2}n}\right)$
and thus the claimed preprocessing time follows.

######  Theorem 7.11.

The data structure described above correctly answers connectivity queries in
amortized time $O(\frac{\log n}{\log\log n})$ with
$O((n\log^{6}n+m\log^{10}n)\lambda_{6}(\log n))$ overall expected update time
for $m$ deletions.

######  Proof 7.12.

As the only difference in the definition of the data structure is the
definition of the sets $S_{1}(A)$ and $S_{2}(A)$, the correctness follows from
Theorem 6.13. Also the $O(\frac{\log n}{\log\log n})$ bound for the queries in
$\mathcal{H}$ again carries over. The preprocessing of the data structure
takes $O\mathopen{}\mathclose{{}\left(n\log^{6}\lambda_{6}(\log n)}\right)$
time by Lemma 7.9. Deletions from $\mathcal{H}$ take
$O\mathopen{}\mathclose{{}\left(\log^{2}n}\right)$ amortized time for each of
the $O\mathopen{}\mathclose{{}\left(n\log^{2}n}\right)$ edges for an overall
time of $O\mathopen{}\mathclose{{}\left(n\log^{4}n}\right)$. The time for the
sequence of deletions is again dominated by the time needed for the updates in
the RDS. Let $m_{A}$ be the number of sites deleted from $S_{1}(A)$ for a
region $A$. Then the time needed by the RDS associated with $A$ is
$O\mathopen{}\mathclose{{}\left(\mathopen{}\mathclose{{}\left(|S_{1}(A)|\log^{5}n+m_{A}\log^{9}n}\right)\lambda_{6}(\log
n)+|S_{2}(A)|\log^{4}n}\right)$ by Corollary 5.12. We have
$\sum_{A\in\mathcal{A}}m_{A}=O(m\log
n),\sum_{A\in\mathcal{A}}|S_{1}(A)|=O(n\log n)$ and
$\sum_{A\in\mathcal{A}}|S_{2}(A)|=O(n\log^{2}n)$. Summing up the time needed
for the RDS over all $A\in\mathcal{A}$, we get a running time of
$O\mathopen{}\mathclose{{}\left(\mathopen{}\mathclose{{}\left(n\log^{6}n+m\log^{10}n}\right)\lambda_{6}(\log
n)}\right)$ as claimed.

##  8 Conclusion

We considered several problems related to dynamic connectivity data structures
for disk graphs. First of all, we significantly improve the state of the art
for unit disk graphs, by considering data structures tailored to this case.
Furthermore, in the general bounded radius ratio case, we were able to improve
the dependency on $\Psi$ for updates. We then considered the incremental and
decremental setting. For the incremental setting with bounded radius ratio, we
gave an a data structure with an amortized update time, that is logarithmic in
$\Psi$ and polylogarithmic in $n$ and near constant query time.

In order to be able to give a similarly efficient data structure in the
decremental setting, we first considered problems related to the lower
envelopes of hyperplanes and more general functions. Using the these, we were
able to describe a dynamic reveal data structure that is elemental for our
decremental data structure and might be of independent interest. Using the
RDS, we were able to give data structures with $O(\frac{\log n}{\log\log n})$
query time. In the bounded setting, the update time is again logarithmic in
$\Psi$ and polylogarithmic in $n$, while for the setting of unbounded radius
ratio, we managed to achieve a data structure whose update time depends only
on $n$.

For the semi-dynamic data structures described in Sections 6 and 7, this
significantly improves the previously best time bounds that can be derived
from the fully-dynamic data structure of Chan et al. [10].

While giving significant improvements in some regards, there are still some
open questions. First of all, our result in the incremental setting only deals
with the setting of bounded radius ratio. We are currently working on
extending it to the case of general disk graphs. Furthermore, as the
incremental and decremental data structures we developed are significantly
faster than the fully dynamic data structures, an interesting question would
be if similar bounds can be achieved in the fully dynamic setting.

## References

  * [1] P. K. Agarwal and J. Matousek.  On range searching with semialgebraic sets.  Discrete Comput. Geom., 11(4):393–418, 1994.  doi:10.1007/BF02574015.
  * [2] Pankaj K. Agarwal, Ravid Cohen, Dan Halperin, and Wolfgang Mulzer.  Dynamic Maintenance of the Lower Envelope of Pseudo-Lines.  arXiv:1902.09565 [cs], February 2019.  arXiv:1902.09565.
  * [3] Pankaj K. Agarwal, Herbert Edelsbrunner, Otfried Schwarzkopf, and Emo Welzl.  Euclidean minimum spanning trees and bichromatic closest pairs.  Discrete $\&$ Computational Geometry, 6(3):407–422, September 1991.  doi:10.1007/BF02574698.
  * [4] Jon Louis Bentley and James B Saxe.  Decomposable searching problems 1. static-to-dynamic transformation.  J. Algorithms, 1:58, 1980.
  * [5] Jack Bresenham.  A linear algorithm for incremental digital display of circular arcs.  Commun. ACM, 20(2):100–106, 1977.
  * [6] G.S. Brodal and R. Jacob.  Dynamic planar convex hull.  In The 43rd Annual IEEE Symposium on Foundations of Computer Science, 2002. Proceedings., pages 617–626, November 2002.  doi:10.1109/SFCS.2002.1181985.
  * [7] Kevin Buchin, Maarten Löffler, Pat Morin, and Wolfgang Mulzer.  Preprocessing imprecise points for Delaunay triangulation: Simplified and extended.  Algorithmica, 61(3):674–693, 2011.  doi:10.1007/s00453-010-9430-0.
  * [8] Timothy M. Chan.  Dynamic planar convex hull operations in near-logarithmic amortized time.  Journal of the ACM, 48(1):1–12, January 2001.  doi:10.1145/363647.363652.
  * [9] Timothy M. Chan.  A dynamic data structure for 3-D convex hulls and 2-D nearest neighbor queries.  Journal of the ACM, 57(3):16:1–16:15, March 2010.  doi:10.1145/1706591.1706596.
  * [10] Timothy M. Chan, Mihai Pătraşcu, and Liam Roditty.  Dynamic Connectivity: Connecting to Networks and Geometry.  SIAM Journal on Computing, 40(2):333–349, January 2011.  doi:10.1137/090751670.
  * [11] Timothy M. Chan and Konstantinos Tsakalidis.  Optimal deterministic algorithms for 2-d and 3-d shallow cuttings.  Discrete Comput. Geom., 56(4):866–881, 2016.  doi:10.1007/s00454-016-9784-4.
  * [12] Mark de Berg, Otfried Cheong, Marc van Kreveld, and Mark Overmars.  Computational Geometry: Algorithms and Applications.  Springer-Verlag, Berlin, 3rd edition, 2008.  doi:10.1007/978-3-540-77974-2.
  * [13] D. Eppstein.  Dynamic Euclidean minimum spanning trees and extrema of binary functions.  Discrete $\&$ Computational Geometry, 13(1):111–122, January 1995\.  doi:10.1007/BF02574030.
  * [14] David Eppstein, Giuseppe F Italiano, Roberto Tamassia, Robert E Tarjan, Jeffery Westbrook, and Moti Yung.  Maintenance of a minimum spanning forest in a dynamic plane graph.  Journal of Algorithms, 13(1):33–54, March 1992.  doi:10.1016/0196-6774(92)90004-V.
  * [15] Steven Fortune.  A sweepline algorithm for Voronoi diagrams.  Algorithmica, 2(1):153, November 1987.  doi:10.1007/BF01840357.
  * [16] Sariel Har-Peled.  Geometric Approximation Algorithms, volume 173 of Mathematical Surveys and Monographs.  American Mathematical Society, Providence, RI, 2011.  doi:10.1090/surv/173.
  * [17] Sariel Har-Peled.  Quadtrees–Hierarchical grids.  In Geometric Approximation Algorithms, volume 173 of Mathematical Surveys and Monographs, chapter 2. American Mathematical Society, Providence, Rhode Island, 2011.  doi:10.1090/surv/173.
  * [18] Christian Hofer.  Fast Dynamic Planar Convex Set Maintenance Using Finger Trees.  Bachelor Thesis, Freie Universität Berlin, December 2017.
  * [19] Jacob Holm, Kristian de Lichtenberg, and Mikkel Thorup.  Poly-logarithmic deterministic fully-dynamic algorithms for connectivity, minimum spanning tree, 2-Edge, and biconnectivity.  Journal of the ACM, 48(4):723–760, 2001.  doi:10.1145/502090.502095.
  * [20] Haim Kaplan, Eyal Molad, and Robert E. Tarjan.  Dynamic rectangular intersection with priorities.  In Proc. 35th Annu. ACM Sympos. Theory Comput. (STOC), pages 639–648, 2003.
  * [21] Haim Kaplan, Wolfgang Mulzer, Liam Roditty, Paul Seiferth, and Micha Sharir.  Dynamic planar voronoi diagrams for general distance functions and their algorithmic applications.  64(3):838–904.  doi:10.1007/s00454-020-00243-7.
  * [22] Haim Kaplan, Robert E. Tarjan, and Kostas Tsioutsiouliklis.  Faster kinetic heaps and their use in broadcast scheduling.  In Proceedings of the Twelfth Annual ACM-SIAM Symposium on Discrete Algorithms, SODA ’01, pages 836–844, USA, January 2001\. Society for Industrial and Applied Mathematics.
  * [23] Mark H. Overmars and Jan van Leeuwen.  Maintenance of configurations in the plane.  Journal of Computer and System Sciences, 23(2):166–204, October 1981.  doi:10.1016/0022-0000(81)90012-X.
  * [24] Raimund Seidel and Micha Sharir.  Top-Down Analysis of Path Compression.  SIAM Journal on Computing, 34(3):515–525, January 2005.  doi:10.1137/S0097539703439088.
  * [25] Paul Seiferth.  Disk Intersection Graphs: Models, Data Structures, and Algorithms.  PhD Thesis, Freie Universität Berlin, 2016.
  * [26] Micha Sharir.  Intersection and Closest-Pair Problems for a Set of Planar Discs.  SIAM Journal on Computing, 14(2):448–468, May 1985.  doi:10.1137/0214034.
  * [27] Micha Sharir and Pankaj K. Agarwal.  Davenport-Schinzel sequences and their geometric applications.  Cambridge University Press, 1995.
  * [28] Daniel D. Sleator and Robert Endre Tarjan.  A data structure for dynamic trees.  Journal of Computer and System Sciences, 26(3):362–391, 1983.  doi:10.1016/0022-0000(83)90006-5.
  * [29] Mikkel Thorup.  Near-optimal fully-dynamic graph connectivity.  In Proc. 32nd Annu. ACM Sympos. Theory Comput. (STOC), pages 343–350, 2000.
  * [30] Andrew C. Yao and F. Frances Yao.  A general approach to D-dimensional geometric queries.  In Proc. 17th Annu. ACM Sympos. Theory Comput. (STOC), pages 163–168, 1985.  doi:10.1145/22145.22163.


# Performance of Non-Parametric Reconstruction Techniques in the Late-Time
Universe

Celia Escamilla-Rivera  Jackson Levi Said  and Jurgen Mifsud

###### Abstract

In the context of a Hubble tension problem that is growing in its statistical
significance, we reconsider the effectiveness of non-parametric reconstruction
techniques which are independent of prescriptive cosmological models. By
taking cosmic chronometers, Type Ia Supernovae and baryonic acoustic
oscillation data, we compare and contrast two important reconstruction
approaches, namely Gaussian processes (GP) and the Locally weighted
Scatterplot Smoothing together with Simulation and extrapolation method
(LOESS-Simex or LS). In the context of these methods, besides not requiring a
cosmological model, they also do not require physical parameters in their
approach to their reconstruction of data (but they do depend on statistical
hyperparameters). We firstly show how both GP and LOESS-Simex can be used to
successively reconstruct various data sets to a high level of precision. We
then directly compare both approaches in a quantitative manner by considering
several factors, such as how well the reconstructions approximate the data
sets themselves to how their respective uncertainties evolve. In light of the
puzzling Hubble tension, it is important to consider how the uncertain regions
evolve over redshift and the methods compare for estimating cosmological
parameters at current times. For cosmic chronometers and baryonic acoustic
oscillation compiled data sets, we find that GP generically produce smaller
variances for the reconstructed data with a minimum value of $\sigma_{\rm GP-
min}=1.1$, while the situation for LS is totally different with a minimum of
$\sigma_{\rm LS-min}=50.8$. Moreover, some of these characteristics can be
alleviate at low $z$, where LS presents less underestimation in comparison to
GP.

##  1 Introduction

The precise value of late time cosmic acceleration remains a great mystery
that is causing serious tensions with predictions from the early Universe
assuming vanilla $\Lambda$CDM [1]. Since the discovery of the accelerating
expansion of the Universe in the late 90s [2, 3], flat $\Lambda$CDM has
emerged as the concordance model despite several theoretical challenges in the
model such as fine-tuning issues, the horizon and coincidence problems, the
cosmological constant and dark matter problems, among many more [4, 5]. In the
background of these challenges, there have been several compelling potential
solution proposals ranging from beyond $\Lambda$CDM such as dynamical dark
energy [6], extended gravity [7] to beyond general relativity (GR) [8], and
others.

The increasing separation in $H_{0}$ values between early Universe
$\Lambda$CDM-based predictions and local cosmology-independent measurements
has become a serious threat that may soon pass over to other cosmological
parameters such as the value of $f\sigma_{8}$ [9, 10], calling into question
the predictive power of $\Lambda$CDM. In this context, it is extremely
important to have a range of reliable techniques to reconstruct combined data
sets that can accurately determine late time cosmological parameter values
[11, 12]. Given a cosmological model, the well known Markov chain Monte Carlo
(MCMC) approach readily determines best-fit model parameters for cosmological
data [13]. In recent years, interest has grown in non-parametric models which
do not assume a cosmological model at all [14]. These methods have largely
relied on supervised learning in which the relations between the points in a
data set are used to either reduce noise in the data or to simulate other
points not in the observations.

In this work, we develop comparison techniques for two popular approaches that
have gained traction in recent years, namely Gaussian process (GP) regression
[15] , and Locally weighted Scatterplot Smoothing (LOESS or LS) [16]. On the
one hand, GP offers an interesting avenue to constrain a covariance function
which can then be used to simulate new points as well as to smoothen places
where data points exist. The effect of having many points in a data set helps
reduce noise in observations. The drawback of this method is that if the
behavior or density of points in one region of a data set is drastically
different to another then one of these areas may be poorly reconstructed. GP
has been used widely in cosmology from smoothening Hubble data [17, 18, 19,
20, 21, 22, 21], $f\sigma_{8}$ data [23], and gravitational wave analysis [24,
25, 26]. A recent development that has garnered momentum in recent years is
that of using GP to reconstruct various possible extended gravity scenarios in
which a precise Lagrangian function can be left relatively arbitrary. In Ref.
[27] viable Horndeski gravity was shown to adequately reproduce late time
Hubble expansion, while in Refs. [28, 29] the inverse problem in $f(T)$
gravity was approached using Hubble data, this was then extended in Ref. [30]
to growth data. Finally, in Ref. [31] some possible interactions between dark
energy and dark matter were explored.

On the other hand, the LOESS reconstruction method is a robust and
computationally non-parametric regression approach to reconstructing data. The
main purpose of this technique is to avoid assuming any prior or cosmological
model. In this line of thought, LOESS recovers the global trend of a data set
by studying the suitable number of neighbouring data points around the pivot
point. Furthermore, in order to compute the measurement errors of the data we
use Simulation and extrapolation method (Simex). The merging of both
techniques is the so-called LOESS-Simex reconstruction, which has been used in
order to reconstruct the cosmic expansion [32], the Om diagnostic up to first
derivative on the Hubble parameter [33] and to reconstruct rotation curves for
dark matter profiles [34]. Both the LOESS-Simex and GP approaches to
reconstruction are model-independent in that they do not depend on a physical
model of cosmology in order to produce reconstructions of cosmological
parameters. However, they are more than this in that they do not require
physical parameters beyond the one they are reconstructing, which is a
significant property for these methods. On the other hand, they do require the
use of statistical hyperparameters in order to function. However, these
hyperparameters are used for statistical purposes and do not have a direct
relationship to the physical parameters of a particular cosmological model.

Nonparametric and model independent reconstructions approaches
---
Method |  Assumption of prior/models |  Binned |  Low efficiency at high $z$ |  Underestimation of errors $z$ |  High computational cost
Principal Components Analysis (PCA) [35] |  ✓ |  ✗ |  ✓ |  ✓ |  ✓
Nonlinear Inverse Approach (NIA) [36] |  ✓ |  ✗ |  ✓ |  ✗ |  ✓
Dipole of the Luminosity Distance method (DLD) [37] |  ✗ |  ✗ |  ✓ |  ✗ |  ✗
Nodal Reconstruction (NR) [38] |  ✓ |  ✗ |  ✗ |  ✗ |  ✗
Genetic Algorithms (GA) [39] |  ✗ |  ✗ |  ✗ |  ✓ |  ✗
Reconstructions of the Expansion History (MIR-I,II,III) [40, 41, 42] |  ✗ |  ✗ |  ✗ |  ✓ |  ✓
Gaussian Processes (GP) [43] |  ✓ |  ✗ |  ✗ |  ✓ |  ✓
Loess-Simex [44] |  ✗ |  ✗ |  ✗ |  ✗ |  ✗

Table 1: A table of nonparametric and model independent reconstructions
approaches. First column: Method and its reference. Second column: Indicates
when a model requires an assumption of a prior ir a fiducial cosmological
model, usually this leads to biased results. Third column: Indicates when the
method share the same data with the bins selected. This causes fluctuations.
Fourth column: This can leads to numerical instability. Fifth column: Methods
that require implementation of tools to propagate errors. Sixth column: High
computational cost mostly due the increase of the data sample.

The literature contains several other approaches to reconstructing
observational data using model-independent approaches [35, 36, 37, 38, 39, 40,
43, 44] although many of them are aimed at continuing data in regions of
scarcity of data. To this end, we list many of the popular ones in the
literature in Table 1 where we emphasize how they cope with error estimation,
their computational cost, as well as whether they make prior assumptions on
the behaviour of the data or whether they use binning to analyse the data.
There is also a growing tendency to use a combination of these approaches to
better reconstruct data sets such as using genetic algorithms with GP to
better motivate the kernel selection problem in this method [45]. There are
also novel approaches being developed such as using artificial neural networks
to reconstruct cosmological data. In this work, our aim is to compare two of
the most promising and used methods in the literature in the cosmological
context, and to show how they compare against each other.

In the context of cosmology and specifically in terms of Hubble data sets,
these methods have not been directly compared against each other. However,
more generally, reconstruction methods are not generally compared against each
other which leaves open the question of which of the vast array of non-
parametric approaches reproduce the data sets best? In this work, we consider
both GP and LOESS-Simex reconstruction techniques in Sec. 2 which we then
apply to Hubble data in Sec. 3. Here, we also discuss how we use the
observational data sets in the context of priors. In Sec. 4, we explore
several ways of directly comparing both reconstruction approaches against each
other. Finally, we discuss and summarise our core results in Sec. 5.

##  2 Model independent approaches to reconstructing Hubble observations

Non-parametric approaches to reconstructing observational data have been
exhaustively explored in light of the possible tensions in observations and
deviations in model predictions with recent precision tests, such as the
$H_{0}$ tension [46, 1] as well as the growing $\sigma_{8,0}$ tension [9, 47].
We introduce here the two methods that we will then use to compare and
contrast expansion rate reconstructions that follow in the ensuing sections.

###  2.1 LOESS-Simex reconstruction for Cosmology

LOESS is a model independent non-parametric regression method. For instance,
with this method it is not necessarily to assume a prior nor do we assume
anything about the cosmological model. The fact that it is local is because we
use the local neighbourhood of each pivot point in turn to infer the global
trend of data set at hand.

Our goal is to reconstruct the Hubble function which can explain the behaviour
of the data set using LOESS. A recipe to generate the reconstructed values
with this method is as follows [32]:

  1. 1.

We choose a subset with $n$-points out of the $N$ data points in the
neighbourhood of our pivot point $z_{i,0}$. The difference between the pivot
point $z_{i,0}$ and the afar point of the neighbourhood
($\text{max}|z_{i,0}-z_{j}|$) where($j=1,2,\ldots n$), is known as the span or
bandwidth $h$ of this subset. As a first step we have to select the
neighbouring points around the pivot point. This can be fixed by deciding the
value of the the smoothing parameter $\alpha$, which relates the optimal
number of neighbourhood points ($n$) to the total number of data points ($N$),
i.e. $n=\alpha N$, where $\alpha$ ranges from $0$ to $1$. The convenient path
to select the optimum value for the smoothing parameter is via a cross
validation technique.

The reconstructed value of $\eta^{obs}_{i}$ is obtained by removing the
$i^{th}$ observation for a given value of $\alpha$. The reconstructed value of
$\eta^{obs}$ is denoted by $\widehat{\eta}_{-i}$ and this process is repeated
for the whole data set. The cross validation function (CV) is given by

$\text{CV}(\alpha)=\frac{1}{N}\sum_{i=1}^{N}(\eta^{\text{obs}}_{i}-{\widehat{\eta}_{-i}})^{2}\,,$
(2.1)

where CV($\alpha$) denotes a mean squared error for different values of
$\alpha$. We then compute CV$(\alpha)$ for several values of $\alpha$. The
value of $\alpha$ for which CV$(\alpha)$ is a minimum is chosen to be the
optimal choice of $\alpha$ and used as the smoothing parameter for the
subsequent calculations.

  2. 2.

Next, we consider a weight function $w_{ij}$ such that it gives more weight to
the points nearest to the pivot point and less weight to points far from the
pivot point. At this stage we take into account a kernel function of the form
$w_{ij}=F[(z_{j}-z_{i,0})/h)]$. The choice of weight function can be performed
by guessing a behaviour where the points near to each other may be more
correlated in comparison to the points which are far away. In this line of
thought, we need to consider a higher weight to observations that are near to
the pivot point $z_{i,0}$. For this reason we choose a tricube weight function
of the form:

$w_{ij}=\left\\{\begin{array}[]{rcl}(1-u_{ij}^{3})^{3}&\mbox{for}&|u_{ij}|<1\\\
0&\mbox{for}&\text{ otherwise}\,,\\\ \end{array}\right.$ (2.2)

where $u_{ij}$ is defined by

$u_{ij}\equiv\frac{z_{i,0}-z_{j}}{h}\,.$ (2.3)

Here $h=\text{max}|z_{i,0}-z_{j}|$ is the maximum distance between the point
of interest and the $j^{\text{th}}$ element of its window.

  3. 3.

We perform a fit of the subset of the data to a local polynomial
$f(\theta_{i})$, with $\theta_{i}$ being the free parameters of the polynomial
model, up to first or second order using the weight function $w_{ij}$. The
main idea under the LOESS assumption is to fit the small subset of data using
low degree polynomials. In another case, a higher degree polynomial may
increase the computational cost without giving any significant improvement on
the result. Therefore, the polynomial fit for each subset of the data is
usually of first or second degree i.e. linear or quadratic. In this work we
consider a fit neighbourhood subset of each pivot point with a linear
polynomial fit as: $f_{i}(a,b)=a+bz_{i}$, which allows to fit an extrapolant
function to the averaged and error-contaminated estimates decrease. Also, we
define a $\chi^{2}$-statistic over this subset of data as

$\chi^{2}_{i}=\sum_{j=1}^{n}w_{ij}(\eta^{\text{obs}}_{j}-f_{j}(\theta_{i}))^{2}\,.$
(2.4)

  4. 4.

As it is standard, we minimise Eq. (2.4) to obtain the best fit value(s)
$\theta_{i}$. Then, the reconstructed value of $\eta^{obs}_{i}$ at pivot point
$z_{i,0}$ is given by

$\widehat{\eta}_{i}=f_{i}(\theta_{0})\,,$ (2.5)

where $\theta_{0}$, denotes best value at the pivot point. We repeat these
four steps for each point in the data set until we obtain the corresponding
reconstructed value at each selected point. Furthermore, in LOESS we require
to fix three numerical parameters: (i) the bandwidth or smoothing parameter,
(ii) the weight function and (iii) the degree of the local polynomial. To
select viable values of each one we need to consider the following:

Once the LOESS process is done, we proceed with the Simex method. In cosmology
most of the observed quantities come with some degree of noise or measurement
errors. However, in the LOESS method we do not use the observed measurement
errors $\sigma_{i}$ while reconstructing the response parameter. The effect of
observational errors can be accommodated in the LOESS process by using a
statistical technique called Simex. This method is based on a two step
resampling approach. In the simulation step, some additional error is
introduced by hand in the data with some controlling parameter $\xi$. Next, by
using regression analysis on this new data set, we try to trace the effect of
the measured error in our original. To compute this we consider three steps:

  1. (a)

For the simulation step, we introduce a fixed amount of measurement error in
each observation data point and define a new variable as

$\eta^{\text{rec}}_{i}(\eta_{k})=\eta^{\text{obs}}_{i}+\sqrt{\eta_{k}}\,\,\,\sigma_{i}\,,$
(2.6)

where $\sigma_{i}$ is the measurement error associated with the observed data
$\eta^{\text{obs}}_{i}$ (we assume that the data points are independent of
each other). Here the parameter $\eta_{k}$ acts as the controlling parameter
for the variance of the measurement error. It is a vector of length $K$ such
that $\eta_{k}>0$. We thus form a matrix of order $K\times N$. The elements of
this matrix will be the values of $\eta^{\text{rec}}_{i}(\eta_{k})$ where
$i=1,2\ldots N$ and $k=1,2\ldots K$. Notice that we assume a standard normal
distribution of the errors.

  2. (b)

After, we reconstruct each data point given in each row of the matrix by
applying the LOESS technique and by selecting each column of the reconstructed
values of $\eta^{\text{rec}}$, we can apply the simple regression technique to
find the best fit value. At this stage, a quadratic polynomial is a better
choice for our purposes.

Finally, we obtain a row having $N$ elements. Each
$\widehat{\eta}^{\text{rec}}_{i}(\eta_{k})$ can be written as a function of
$\eta_{k}$, i.e.

$\widehat{\eta}^{\text{rec}}_{i}(\eta_{k})=k_{1}+k_{2}\eta_{k}+k_{3}\eta_{k}^{2}\,,$
(2.7)

where $k_{1}$, $k_{2}$ and $k_{3}$ are normalised constants. If a normal
distribution of the errors is assumed, then the error variance associated with
the simulated data points $\widehat{\eta}^{\text{rec}}_{i}$ will be
$(1+\eta_{k})\sigma_{i}^{2}$. On substituting $\eta_{k}=-1$, we obtain error-
free smoothed data points.

  3. (c)

To construct the 1-$\sigma$ and 2-$\sigma$ C.L. around the nonparametric
regression curve of $\widehat{\eta}$ , we assume that the errors are
distributed normally. This can be obtained by using the limiting values:
$\widehat{\eta}_{i}\pm\gamma\sqrt{V(\widehat{\eta}_{i})}$, with
$\gamma=1,2,\ldots$ and

$V(\widehat{\eta}_{i})=\frac{1}{df_{\text{res}}}\sum_{i}^{N}d_{i}^{2}\sum_{j}^{N}w_{ij}^{2}\,,$
(2.8)

where $d_{i}=\eta_{i}-\widehat{\eta}_{i}$ and
$df_{\text{res}}=N-df_{\text{mod}}$, such that $df_{\text{mod}}$ is the number
of effective degrees of freedom or the effective number of parameters used in
this regression. We calculate it by using the normalised smoothing matrix S,
such that $df_{\text{mod}}=\text{Tr}($SST). The smoothing matrix S, which is a
$N\times N$ square matrix of $w_{ij}$ elements, is directly calculated from
the weight function.

LOESS-Simex reconstruction has been used to reconstruct the cosmic expansion
in Ref. [32] along with the transition redshift reconstruction of the cosmic
accelerated expansion cosmokinematics analysed in Ref. [48]. Furthermore, its
implementation to null tests of dark energy has been done in order to perform
an Om diagnostic where a static domain wall network dark energy scenario
results viable [33]. Also, in Ref. [34] it was shown that LOESS-Simex is
practical in astrophysical configurations to reconstruct rotation curves for
dark matter profiles .

###  2.2 Gaussian Processes reconstruction for cosmology

GP regression is a non-parametric approach to reconstructing observational
data which are incorporated as a finite collection of normally distributed
points. This type of supervised learning is implemented for stochastic data
whose data points are related together enough such that a covariance function
can be fit in such a way as to prescribe these relations in a specific way
[49, 11]. In practice a GP is defined in terms of the mean function $\mu(z)$
together with its associated two–point covariant function
$\mathcal{C}(z,\tilde{z})$ which produces the continuous realization

$\xi(z)\sim\mathcal{GP}\left(\mu(z),\mathcal{C}(z,\tilde{z})\right)\,,$ (2.9)

where the uncertainty $\Delta\xi(z)$ is also part of the reconstruction which
produces a realization region $\xi(z)\pm\Delta\xi(z)$. By and large the mean
is set to zero without loss of generality so that the covariance function
becomes the defining feature for each reconstructed data set. For redshift
points $z^{*}$ which represent points where observations do not occur, we can
define a kernel function for the covariance function such that
$\mathcal{K}\left(z^{*},z^{*^{\prime}}\right)=\mathcal{C}\left(z^{*},z^{*^{\prime}}\right)$
where the trained covariance function is utilised. Thus, the kernel will
embody all the information about the strength and amplitude of the
correlations between the redshift data points [17]. The only strong
requirement on the kernel is that it is a symmetric function.

For observational point $\tilde{z}$, we also have available to us information
on the associated uncertainty region covariance matrix
$\mathcal{D}\left(\tilde{z},\tilde{z}^{\prime}\right)$ between data points.
Thus, the covariance function for these points can be written as
$\mathcal{C}\left(\tilde{z},\tilde{z}^{\prime}\right)=\mathcal{K}\left(\tilde{z},\tilde{z}^{\prime}\right)+\mathcal{D}\left(\tilde{z},\tilde{z}^{\prime}\right)$
which will provide information to fit the kernel. Naturally, observational and
reconstructed points can be interrelated by the kernel alone through
$\mathcal{C}\left(z^{*},\tilde{z}^{\prime}\right)=\mathcal{K}\left(z^{*},\tilde{z}^{\prime}\right)$
[50].

Therefore, for a Gaussian distribution, the posterior distribution of a
reconstructed function can be expressed via the joint Gaussian distribution of
different data points in which the kernel expresses the interrelated mean
value and uncertainties of each point in the resulting distribution [21, 17].
In Appendix A we develop a discussion of this GP technique and the GaPP
code111http://ascl.net/1303.027 in order to check how sensitive the chosen
horizontal axis scale is related to thedata points. In the present work, we
use the squared–exponential kernel but there exist a plethora of such kernel
options [11]. Despite each producing a slightly different non-parametric
reconstruction of a distribution, they largely agree to well within $1\sigma$
confidence regions. The general purpose squared–exponential kernel is
explicitly defined as

$\mathcal{K}\left(z,\tilde{z}\right)=\sigma_{f}^{2}\exp\left[-\frac{\left(z-\tilde{z}\right)^{2}}{2l_{f}^{2}}\right]\,,$
(2.10)

where $\sigma_{f}$ and $l_{f}$ are the kernel hyperparameters, which are
parameters of the kernel but do not parameterise the function being
reconstructed, as a model would. The hyperparameters characterise the variance
in the data through $\sigma_{f}$ and the length–scale through $l_{f}$. Thus,
the hyperparameters define the smoothness and reach for fluctuations in a
signal. As a result, larger values of $l_{f}$ lead to smoother GP functions
while higher values of $\sigma_{f}$ express a lower signal–to–noise ratio.

To determine the best suited hyperparameters, their values are derived from
the maximisation of the probability of the GP to generate the data set under
consideration which is implemented via the minimisation of the GP marginal
likelihood similar to a Bayesian approach [51, 52, 21, 53]. GP regression has
now been extensively studied as a tool for reconstructions in the cosmological
context which have largely focused on the late-time behaviour of the expansion
rate [54, 21, 31, 55, 56, 57, 58, 17, 59, 60, 20, 61, 62, 19]. In Refs. [28,
29] the use of GP reconstructions was extended from kinetic parameters to
theories beyond $\Lambda$CDM. As an example in this line of thought, $f(T)$
gravity was then extended into the perturbative section in Ref. [30], where an
$f\sigma_{8}$ data set was used to constrain the form of potential $f(T)$
cosmology models that are in agreement with earlier GP reconstructions of
growth data in Ref. [23]. Other recent work also includes a possible
resolution to the Hubble tension problem within $\Lambda$CDM in Refs. [63, 18]
where the model independent nature of GP is brought into question due to its
dependence on the kernel choice in the regression procedure of the GP
reconstruction algorithm. Furthermore, in quintessence and k-essence models
from Horndeski theories of gravity [27] that represent an extension to
quintessence can reproduce reconstructions of the late expansion of the
Universe within 2$\sigma$.

##  3 Reconstructions applied to cosmic data

With the recipes explained above, in this section we develop the
reconstruction for late-time observations based on luminosity distances and
Hubble flow reconstructions. The LOESS-Simex and GP reconstructions are based
on current data sets of Type Ia supernovae (SNeIa), observational Hubble data
(CC) and baryon acoustic oscillations (BAO).

###  3.1 Cosmic chronometers and baryonic acoustic oscillations

For the observational Hubble data we consider a sample of $l=51$ measurements
in the range $0.07<z<2.36$. A total sample size of 31 data points from passive
galaxies, or cosmic chronometers (CC) [64, 65, 66, 67, 68], and 20 data points
which are estimated from BAO data points [69, 70, 71, 72, 73, 74, 75, 76, 77,
78] under a $\Lambda$CDM prior. Following Ref. [79], we assume that these
points are not correlated. To remove the model dependence in some of the data
points, we use a sound horizon $r_{s}$ computed from Planck 2018 [10]
(henceforth denoted by PL18). To perform the model independent reconstruction,
we set the construction of a $\chi_{H}^{2}$ as

$\chi_{H}^{2}=\sum_{i=1}^{l}\frac{\left[H_{\text{recons}}\left(z_{i},\mathbf{x}\right)-H_{\text{obs}}(z_{i})\right]^{2}}{\sigma^{2}_{H}(z_{i})}\,,$
(3.1)

where $\text{H}_{\textit{obs}}(z_{i})$ is the observed value at $z_{i}$,
$\sigma_{H}(z_{i})$ are the observational errors, and
$\text{H}\left(z_{i},\mathbf{x}\right)$ is the reconstructed H for the same
$z_{i}$ with the specific parameter vector $\mathbf{x}$.

For the BAO sample we consider 15 transversal measurements in a range of
$z=[0.11,\,2.225]$ , obtained in a quasi–model independent approach. As it is
standard, this can be computed by considering the 2-point angular correlation
function tracers with $D_{A}(z;r_{\text{drag}})$ [80]. This data can be
compared with the sound horizon today to the sound horizon at the time of
recombination inferred from the CMB anisotropy measurements. Defining the
fundamental quantities for this observation, we have the BAO distances
$d_{z}\equiv\frac{r_{s}(z_{d})}{D_{V}(z)}$, with
$r_{s}(z_{d})=\frac{c}{H_{0}}\int_{z_{d}}^{\infty}\frac{c_{s}(z)}{E(z)}\,\mathrm{d}z$
being the comoving sound horizon at the baryon dragging epoch, $c$ the speed
of light, $z_{d}$ is the drag epoch redshift, and
$c^{2}_{s}=c^{2}/3[1+(3\Omega_{b0}/4\Omega_{\gamma 0})(1+z)^{-1}]$ the sound
speed with $\Omega_{b0}$ and $\Omega_{\gamma 0}$ being the current values of
baryon and photon density parameters, respectively. We can also define the
dilation scale by

$D_{V}(z,\Omega_{m0};\Theta)=\left[\frac{c\,z(1+z)^{2}D_{A}^{2}}{H(z,\Omega_{m0};\Theta)}\right]^{1/3}\,,$
(3.2)

where $D_{A}$ is the angular diameter distance given by

$D_{A}(z,\Omega_{m0};\Theta)=\frac{1}{1+z}\int_{0}^{z}\frac{c\,\mathrm{d}\tilde{z}}{H(\tilde{z},\Omega_{m0};\Theta)}\,,$
(3.3)

where $\Theta$ is the vector that contains the free cosmological parameters to
be reconstructed (in such cases fitted). Using the comoving sound horizon, the
distance ratio $d_{z}$ can be related to the normalised Hubble parameter $h$
($H\doteq 100h$) and the critical densities $\Omega_{m0}$ and $\Omega_{b0}$.
To connect the BAO measurements with SNeIa to CMB measurements (PL18), we take
into account the Alcock-Paczynski distortion parameter given by

$F(z,\Theta)=(1+z)\frac{D_{A}(z,\Theta)H(z,\Theta)}{c}\,.$ (3.4)

Here we notice that by calibrating the $D_{A}$ from BAO, we can define

$\chi^{2}_{\text{BAO}}=(\Delta\mathcal{F}_{\text{BAO}})^{T}\cdotp
C_{\text{BAO}}^{-1}\cdotp\Delta\mathcal{F}_{\text{BAO}}\,,$ (3.5)

where $\Delta\mathcal{F}_{\text{BAO}}$ is the difference between the data and
the resulting value for $\Theta$, and $C_{\text{BAO}}^{-1}$ is the inverse of
the covariance matrix of this sample. In our analysis the off-diagonal entries
of the covariance matrix are taken to be zero. We do this since the LS method
cannot tackle covariance matrices. However, we did perform the following
analyses for GP both with and without the covariance matrix for this data set
which resulted in very similar results. Thus, for consistency, we are able to
take the GP analysis without a covariance matrix.

Figure 1: Analysis reconstruction using CC+BAO data set (red color points).
Left column: Cross validation (Eq. (2.1)) is shown for $\alpha$ values ranging
between $0$ and $1$. Middle column: The LOESS-Simex reconstructed Hubble
diagram is shown (blue color contours) for each prior against a $\Lambda$CDM
(PL18) reference (dashed line). Right column: The GP reconstructed Hubble
diagram is shown (blue color contours) for each prior against a $\Lambda$CDM
(PL18) reference. In both reconstructions we show the 1$\sigma$ and 2$\sigma$
confidence regions of the reconstructed data. From top to bottom we denote the
$H_{0}$ prior cases as: TRGB, HW and R19, respectively.

###  3.2 Pantheon data set

The Pantheon SNeIa compilation is one of the latest Type Ia Supernovae
compilations [81] and it contains 1048 SNeIa at redshift $0.01<z<2.26$. The
constraining power of this kind of supernovae is due to the fact that this
observations can be used as standarisable candles. This can be implemented
through the use of the distance modulus

$\mathcal{F}(z,\Theta)_{\text{theo}}=5\log_{10}\left[D_{L}(z,\Theta)\right]+\mu_{0}\,,$
(3.6)

where $D_{L}$ is the luminosity distance given by

$D_{L}(\Theta)=(1+z)\int_{0}^{z}{\frac{c\,dz^{\prime}}{H_{0}E(z^{\prime},\Theta)}}\,,$
(3.7)

and $\Theta$ is the vector with the free cosmological parameters to be fitted.
We notice that the factor $c/H_{0}$ can be absorbed in $\mu_{0}$. Furthermore,
we can write
$\Delta\mathcal{F}(\Theta)=\mathcal{F}_{\text{theo}}-\mathcal{F}_{\text{obs}}$,
using for this purpose the distance modulus $\mathcal{F}_{\text{obs}}$
associated with the observed magnitude. At this point, it may be thought that
a possible $\chi_{SN}^{2}$ is given by

$\chi_{SN}^{2}(\Theta)=\left(\Delta\mathcal{F}(\Theta)\right)^{T}\cdotp
C_{SN}^{-1}\cdotp\Delta\mathcal{F}(\Theta)\,,$ (3.8)

where $C_{SN}$ is the total covariance matrix. As in the case of the BAO data
set, we tested GP both with and without the covariance matrix with very
similar results being returned. For this reason, and due to the LS method not
being able to handle covariance matrices, we used the GP method without the
covariance matrix approach for consistency. This equation can be used to
contain the nuisance parameter $\mu_{0}$, which in turn is a function of the
Hubble constant, the speed of light $c$ and the SNeIa absolute magnitude. To
circumvent this issue, $\chi_{SN}^{2}$ is marginalised analytically with
respect to $\mu_{0}$ and we can obtain a new $\chi_{SN}^{2}$ estimator, given
by (using Eq.(6,7) of Ref. [81])

$\chi_{SN}^{2}(\Theta)=\left(\Delta\mathcal{F}(\Theta)\right)^{T}\cdotp
C_{SN}^{-1}\cdotp\Delta\mathcal{F}(\Theta)+\ln{\frac{S}{2\pi}}-\frac{k^{2}(\Theta)}{S}\,,$
(3.9)

where $S$ is the sum of all entries of $C_{SN}^{-1}$. This equation gives an
estimation of the precision of these data points independently of $\Theta$,
and $k$ is $\Delta\mathcal{F}(\Omega_{m},\Omega_{r},\Omega_{\Lambda})$ but
weighed by a covariance matrix as follows:

$k(\Theta)={\left(\Delta\mathcal{F}(\Theta)\right)^{T}\cdotp C_{SN}^{-1}}\,.$
(3.10)

Also, for this sampler we are taking the nuisance parameter $M$ inside the
sample. Consequently, we perform three calibrations taking into account three
different priors of $H_{0}$:

  * •

TRGB [82]: $H_{0}^{\text{TRGB}}=69.8\pm 1.9\,\text{km/s/Mpc}$.

  * •

H0LiCOW (HW) [83]: $H_{0}^{\text{HW}}=73.3^{+1.7}_{-1.8}\,\text{km/s/Mpc}$.

  * •

Riess et al (R19) [84] : $H_{0}^{\text{R}}=73.24\pm 1.74\,\text{km/s/Mpc}$.

For this case we based our calibration on Ref. [85] by extrapolating the
expression $M_{1}=M_{2}+5\log 10(H_{1}/H_{2})$, where the indices $1$ and $2$
denote our different values for $H_{0}$. Finally, our result is $\mu(z)=m-M$.
A compilation of these measurements can be found in [86, 87, 88, 89].

##  4 Comparative analysis of reconstruction methods

In this section, we will compare and contrast the GP and LOESS-Simex
reconstruction methods in terms of their performance against each other at
reproducing the observations along with their relative uncertainties. In
particular, we aim to quantify the ability of both methods to reproduce the
expansion data together with the uncertainties in those reconstructions of the
data.

Figure 2: Analysis reconstruction using Pantheon SNeIa data set (orange color
points). Left column: Cross validation (Eq. (2.1)) is shown for $\alpha$
values ranging between $0$ and $1$. Middle column: The LOESS-Simex
reconstructed modulus distance $\mu(z)$ diagram is shown (blue color contours)
for each prior against a $\Lambda$CDM (PL18) reference (dashed line). Right
column: The GP reconstructed modulus distance $\mu(z)$ diagram is shown (blue
color contours) for each prior against a $\Lambda$CDM (PL18) reference. In
both reconstructions we depict the 1$\sigma$ and 2$\sigma$ confidence regions
of the reconstructed data. From Top to Bottom we denote the $H_{0}$ prior
cases as: TRGB, HW and R19, respectively.

To this end, the reconstructions for the CC+BAO data set combination shown in
Fig. 1 as well as the Pantheon data set reconstructions shown in Fig. 2 show
good agreement between the mean reconstructed functions and their respective
observational data for each method and prior combination. However, the
respective uncertainties between both approaches is the most drastic
difference between the methods. In this context, Fig. 3 shows the discrepancy
between the variances for reconstructions using both the CC+BAO and SNeIa data
sets in the background of all three priors. Here, we observe that the
selection of priors has little to no effect on the different variances that
both reconstructions result in.

Another important property to notice in the CC+BAO combination is that LOESS-
Simex retains a variance that largely remains at the same order of magnitude
while the GP reconstruction has a dependence on the density of actual
observational data points in any particular redshift region. This is to be
expected since LOESS-Simex reconstructs data points with an equal contribution
from each point in the data set while GP is impacted by the density of points
between different regions on the redshift reconstruction space. We also notice
that GP predicts much lower uncertainties as compared with LOESS-Simex which
is a result of the way in which it assumes the data points to be sourced,
namely that each point is a peak in a background of a normally distributed set
of data points. Thus, GP makes stronger assumptions on the underlying data
being probed. Finally, despite the impact of the prior being minor, higher
$H_{0}$ priors do lead to slightly lower variances across the redshift range
being considered with the biggest effect being sourced by the R19 prior since
this has the highest prior value. Precisely the same overall behaviour is
observed for the Pantheon reconstructed variance shown in Fig. 3 in which the
fluctuations in the GP variances is less pronounced.

Figure 3: Top: Comparisons of the variance are shown for the three priors for
the CC+BAO data set. Bottom: Comparisons of the variance are shown for the
three priors for the Pantheon data set. From Left to Right we denote the
$H_{0}$ prior cases as: TRGB, HW and R19, respectively.

Figure 4: Top: The distance between the two reconstructions is shown in
$\sigma$ units for the CC+BAO data set. Bottom: The distance between the two
reconstructions is shown in $\sigma$ units for the Pantheon data set. From
Left to Right we denote the $H_{0}$ prior cases as: TRGB, HW and R19,
respectively.

Both reconstruction approaches can be compared directly by determining the
distance between their predicted reconstruction values (in units of $\sigma$),
which can be defined as

$\displaystyle D(z)$ $\displaystyle=$ $\displaystyle\frac{Q_{\rm GP}(z)-Q_{\rm
LS}(z)}{\sqrt{\sigma_{\rm GP}^{2}(z)+\sigma_{\rm LS}^{2}(z)}}\,,$ (4.1)

where $Q$ denotes the Hubble parameter $H(z)$ for the CC+BAO combination, and
$\mu(z)$ for the Pantheon SNeIa data set. We show the evolution of the
distance parameter $D(z)$ for the CC+BAO and Pantheon data sets in Fig. 4. In
both cases the reconstructions never cross the 1$\sigma$ line and in the
Pantheon reconstruction, this limit is much stricter which is a result of the
much larger quantity of data points. As in Fig. 4, the effect of the priors
has little impact on the general behavior of $D(z)$. On the other hand, there
is a weak correlation between higher distances at higher redshifts for high
priors such as R19. For the CC+BAO data set combination, GP predicts larger
values for the Hubble parameter except for a small neighbourhood at about
$z=1.2$ (mostly prior independent). In the Pantheon data set case, the LOESS-
Simex reconstruction surpasses the GP predicted values only in the initial
neighbourhood of $z=0$ (also prior independent). However, in this second case,
the predicted values are much closer together and, generally, much closer to
the data points themselves.

Furthermore, the value measured by the following $d(z)$ indicator quantifies
the displacement of both reconstructions against the original data. However,
to quantify the quality of the entire reconstructions, we sum its square
through

$d(z)=\frac{Q_{\rm reconstructed}(z)-Q_{\rm observation}(z)}{\sqrt{\sigma_{\rm
reconstructed}^{2}(z)+\sigma_{\rm observation}^{2}(z)}}\,,$ (4.2)

which gives another important measure of the approximation of either approach
to producing the original data (in units of $\sigma$). We denote the sum as

$\Upsilon=\sum_{z}d^{2}(z)\,,$ (4.3)

which is used to produce Table 2.

Prior | Reconstruction method | $\Upsilon_{\rm CC+BAO}$ | $\Upsilon_{\rm SN}$
---|---|---|---
TRGB | GP | 16.53 | 826.41
LS | 16.68 | 16.60
HW | GP | 15.72 | 917.93
LS | 16.84 | 16.35
R19 | GP | 15.66 | 976.11
LS | 16.86 | 16.22
Table 2: Values of $\Upsilon$ are shown for all combinations of reconstruction
methods and priors for the CC+BAO and SNeIa data sets.

In Table 2, the value of $\Upsilon$ is shown for each combination of prior
$H(z)$ value together with both reconstructions. In the background of this
setting, both reconstructions perform at a similar order of magnitude with the
GP reconstruction producing values that systematically have smaller variances.
Table 2 also shows the over-fitting by the GP reconstruction, which in this
instance produces very small variances which thus produce large values of
$\Upsilon$.

##  5 Conclusions

In the current era of precision cosmology, the amount of observational data is
significantly increasing. Moreover, the robustness of cosmological data is
always being enhanced due to state-of-the-art cosmological probes of the
Universe. Therefore it is natural to adopt machine learning techniques to
cosmological data in order to infer model independent results which are purely
driven by observational data. In this analysis we focused on the model
independent reconstruction of functions via two distinct machine learning
techniques, namely LOESS-Simex and GP which are respectively introduced in
Sec. 2.1 and Sec. 2.2.

As described in Sec. 3.1 and Sec. 3.2, we here considered the model
independent reconstruction of $H(z)$ by using the CC and BAO data, as well as
the reconstruction of $\mu(z)$ by adopting the Pantheon SNeIa data. We further
considered the latter data sets in the presence of the TRGB, HW and R19
$H_{0}$ priors. In the Pantheon SNeIa data set, these prior values were
introduced via the calibration of the original SNeIa data set, as described in
Sec. 3.2. In our study, we considered observational data directly to avoid
convoluting our results with that of calibration, for instance. A further
study may better incorporate this together with a fuller understanding of the
tension in the value of the $H_{0}$ constant, among other effects such as the
tendency of BAO data to prefer a lower value of this parameter.

Further to our GP and LOESS-Simex reconstructions of $H(z)$ and $\mu(z)$, we
statistically compared the two methods in Sec. 4. In general, both machine
learning techniques were able to reconstruct $H(z)$ and $\mu(z)$ data very
well. Moreover, the $\Lambda$CDM model was found to be in very good agreement
with all the inferred GP and LOESS-Simex reconstructions. We observed that the
variance in GP reconstructions is smaller than that of LOESS-Simex, although
in the case of LOESS-Simex the variance does not fluctuate as much as in the
GP scenarios. Generically, the LOESS-Simex technique is characterised by a
broader and more conservative uncertainty region than that of GP. We also
showed that these distinguishing features of the considered data driven
reconstruction techniques are independent from the adopted $H_{0}$ prior
values.

Undoubtedly, machine learning techniques are already giving new insights to
the theoretical cosmological framework via a more elaborate exploration and
interpretation of observational data. Therefore, further development and
implementation of machine learning techniques to cosmological data would
facilitate our search for specific cosmological signatures, which would shed
light on the set of viable cosmological models.

## Acknowledgments

CE-R is supported by the Royal Astronomical Society as FRAS 10147 and by
DGAPA-PAPIIT-UNAM Project IA100220. This article is based upon work from
CANTATA COST (European Cooperation in Science and Technology) action CA15117,
EU Framework Programme Horizon 2020. The authors would like to acknowledge
networking support by the COST Action CA18108 and funding support from
Cosmology@MALTA which is supported by the University of Malta. We are grateful
to Reginald Christian Bernardo for important feedback on new methods in GP
regression. This research has been carried out using computational facilities
procured through the European Regional Development Fund, Project No. ERDF-080
“A supercomputing laboratory for the University of Malta”.

##  Appendix A GP subtleties

We would like to note that the GP technique and the GaPP code are sensitive to
the chosen horizontal axis scale of the data points, and hence to the data set
that one aims to reconstruct. Indeed, in Fig. 2 we have adopted a logarithmic
horizontal scale for the data points, which led to very tight constraints on
the reconstructed function of $\mu(z)$. However, when we considered a linear
horizontal scale for the Pantheon SNeIa data set, we inferred a more
conservative and less robust reconstruction, as illustrated in Fig. 5. In the
case of a linear redshift scale, the GP hyperparameters are dominated by the
low redshift data points. Consequently, the derived hyperparameter values do
not fit the high redshift data points very well in this case. Indeed, as
illustrated in Fig. 5, at high redshifts, the fitted hyperparameters do not
lead to a continuously and monotonically smooth reconstruction band. On the
other hand, when a logarithmic redshift scale is considered, the
hyperparameters fit the entire data set very well, as illustrated in Fig. 2.
This observation could be explained via the inferred GP kernel hyperparameter
values of $\sigma_{f}$ and $l_{f}$, as reported in Table 3. It is clear that
the optimal value of $l_{f}$ is significantly altered by the chosen redshift
scale of the data points, which is expected since the length separation from
one data point to another is altered when rescaling from a linear to a
logarithmic axis.

Figure 5: GP reconstructions of the Pantheon SNeIa data set using a linear redshift scale with the TRGB (top left), R19 (top right) and HW (bottom) $H_{0}$ prior values. Prior | Scale | $\sigma_{f}$ | $l_{f}$
---|---|---|---
TRGB | Linear | 42.89 | 0.18
Logarithmic | 61.32 | 7.36
HW | Linear | 44.61 | 0.18
Logarithmic | 59.37 | 6.98
R19 | Linear | 46.55 | 0.18
Logarithmic | 58.76 | 6.82
Table 3: A comparison of the derived hyperparameter values of $\sigma_{f}$ and
$l_{f}$ with the TRGB, HW and R19 Hubble constant priors when adopting a
linear or a logarithmic redshift scale in the corresponding GP
reconstructions.

## References

  * [1] J. L. Bernal, L. Verde and A. G. Riess, _The trouble with $H_{0}$_, _JCAP_ 10 (2016) 019 [1607.05617].
  * [2] Supernova Cosmology Project collaboration, _Measurements of $\Omega$ and $\Lambda$ from 42 high redshift supernovae_, _Astrophys.J._ 517 (1999) 565 [astro-ph/9812133].
  * [3] Supernova Search Team collaboration, _Observational evidence from supernovae for an accelerating universe and a cosmological constant_ , _Astron.J._ 116 (1998) 1009 [astro-ph/9805201].
  * [4] S. Weinberg, _The cosmological constant problem_ , _Rev. Mod. Phys._ 61 (1989) 1.
  * [5] P. Bull et al., _Beyond $\Lambda$CDM: Problems, solutions, and the road ahead_, _Phys. Dark Univ._ 12 (2016) 56 [1512.05356].
  * [6] E. J. Copeland, M. Sami and S. Tsujikawa, _Dynamics of dark energy_ , _Int. J. Mod. Phys. D_ 15 (2006) 1753 [hep-th/0603057].
  * [7] S. Capozziello and M. De Laurentis, _Extended Theories of Gravity_ , _Phys. Rept._ 509 (2011) 167 [1108.6266].
  * [8] T. Clifton, P. G. Ferreira, A. Padilla and C. Skordis, _Modified Gravity and Cosmology_ , _Phys. Rept._ 513 (2012) 1 [1106.2476].
  * [9] E. Di Valentino et al., _Cosmology Intertwined III: $f\sigma_{8}$ and $S_{8}$_, 2008.11285.
  * [10] Planck collaboration, _Planck 2018 results. VI. Cosmological parameters_ , _Astron. Astrophys._ 641 (2020) A6 [1807.06209].
  * [11] C. E. Rasmussen and C. K. I. Williams, _Gaussian Processes for Machine Learning (Adaptive Computation and Machine Learning)_. The MIT Press, 2005.
  * [12] J. Surma, _Hacking Machine Learning: Towards The Comprehensive Taxonomy of Attacks Against Machine Learning Systems_ , ICIAI 2020. Association for Computing Machinery, New York, NY, USA, 2020, 10.1145/3390557.3394126.
  * [13] J. Møller, _Spatial Statistics and Computational Methods_ , Lecture Notes in Statistics, 173. Springer New York, New York, NY, 1st ed. 2003. ed., 2003, 10.1007/978-0-387-21811-3.
  * [14] J. Kiusalaas, _Numerical Methods in Engineering with Python 3_. Cambridge University Press, 3 ed., 2013, 10.1017/CBO9781139523899.
  * [15] B. C. S. Banerjee and A. Gelfand, _Hierarchical modeling and analysis for spatial data_ , .
  * [16] W. S. Cleveland, _Robust locally weighted regression and smoothing scatterplots_ , _Journal of the American Statistical Association_ 74 (1979) 829.
  * [17] A. Gómez-Valent and L. Amendola, _$H_{0}$ from cosmic chronometers and Type Ia supernovae, with Gaussian Processes and the novel Weighted Polynomial Regression method_, _JCAP_ 04 (2018) 051 [1802.01505].
  * [18] E. Colgáin and M. M. Sheikh-Jabbari, _Elucidating cosmological model dependence with $H_{0}$_, 2101.08565.
  * [19] M. K. Yennapureddy and F. Melia, _Reconstruction of the HII Galaxy Hubble Diagram using Gaussian Processes_ , _JCAP_ 11 (2017) 029 [1711.03454].
  * [20] E.-K. Li, M. Du, Z.-H. Zhou, H. Zhang and L. Xu, _Test the effects of $H_{0}$ on $f\sigma_{8}$ tension with Gaussian Process method_, 1911.12076.
  * [21] M. Seikel and C. Clarkson, _Optimising Gaussian processes for reconstructing dark energy dynamics from supernovae_ , 1311.6678.
  * [22] M. Seikel, C. Clarkson and M. Smith, _Reconstruction of dark energy and expansion dynamics using gaussian processes_ , _Journal of Cosmology and Astroparticle Physics_ 2012 (2012) 036–036.
  * [23] D. Benisty, _Quantifying the $S_{8}$ tension with the Redshift Space Distortion data set_, _Phys. Dark Univ._ 31 (2021) 100766 [2005.03751].
  * [24] E. Belgacem, S. Foffa, M. Maggiore and T. Yang, _Gaussian processes reconstruction of modified gravitational wave propagation_ , _Phys. Rev. D_ 101 (2020) 063505 [1911.11497].
  * [25] C. J. Moore, C. P. L. Berry, A. J. K. Chua and J. R. Gair, _Improving gravitational-wave parameter estimation using Gaussian process regression_ , _Phys. Rev. D_ 93 (2016) 064001 [1509.04066].
  * [26] G. Cañas Herrera, O. Contigiani and V. Vardanyan, _Learning how to surf: Reconstructing the propagation and origin of gravitational waves with Gaussian Processes_ , 2105.04262.
  * [27] M. Reyes and C. Escamilla-Rivera, _Improving data-driven model-independent reconstructions and new constraints in Horndeski cosmology_ , 2104.04484.
  * [28] R. Briffa, S. Capozziello, J. Levi Said, J. Mifsud and E. N. Saridakis, _Constraining teleparallel gravity through Gaussian processes_ , _Class. Quant. Grav._ 38 (2020) 055007 [2009.14582].
  * [29] Y.-F. Cai, M. Khurshudyan and E. N. Saridakis, _Model-independent reconstruction of $f(T)$ gravity from Gaussian Processes_, _Astrophys. J._ 888 (2020) 62 [1907.10813].
  * [30] J. Levi Said, J. Mifsud, J. Sultana and K. Z. Adami, _Reconstructing teleparallel gravity with cosmic structure growth and expansion rate data_ , 2103.05021.
  * [31] T. Yang, Z.-K. Guo and R.-G. Cai, _Reconstructing the interaction between dark energy and dark matter using Gaussian Processes_ , _Phys. Rev. D_ 91 (2015) 123533 [1505.04443].
  * [32] A. Montiel, R. Lazkoz, I. Sendra, C. Escamilla-Rivera and V. Salzano, _Nonparametric reconstruction of the cosmic expansion with local regression smoothing and simulation extrapolation_ , _Phys. Rev. D_ 89 (2014) 043007 [1401.4188].
  * [33] C. Escamilla-Rivera and J. C. Fabris, _Nonparametric reconstruction of the O m diagnostic to test $\Lambda$CDM_, _Galaxies_ 4 (2016) 76 [1511.07066].
  * [34] L. M. Fernández-Hernández, A. Montiel and M. A. Rodríguez-Meza, _Galaxy rotation curves using a non-parametric regression method: core, cuspy and fuzzy scalar field dark matter models_ , _Mon. Not. Roy. Astron. Soc._ 488 (2019) 5127 [1809.06875].
  * [35] D. Huterer and G. Starkman, _Parameterization of dark-energy properties: A Principal-component approach_ , _Phys. Rev. Lett._ 90 (2003) 031301 [astro-ph/0207517].
  * [36] C. Espana-Bonet and P. Ruiz-Lapuente, _Dark energy as an inverse problem_ , hep-ph/0503210.
  * [37] C. Bonvin, R. Durrer and M. Kunz, _The dipole of the luminosity distance: a direct measure of H(z)_ , _Phys. Rev. Lett._ 96 (2006) 191302 [astro-ph/0603240].
  * [38] J. Alberto Vazquez, M. Bridges, M. P. Hobson and A. N. Lasenby, _Reconstruction of the Dark Energy equation of state_ , _JCAP_ 09 (2012) 020 [1205.0847].
  * [39] C. Bogdanos and S. Nesseris, _Genetic algorithms and supernovae type ia analysis_ , _Journal of Cosmology and Astroparticle Physics_ 2009 (2009) 006–006.
  * [40] R. A. Daly and S. G. Djorgovski, _A model‐independent determination of the expansion and acceleration rates of the universe as a function of redshift and constraints on dark energy_ , _The Astrophysical Journal_ 597 (2003) 9–20.
  * [41] S. Fay and R. Tavakol, _Model-independent dark energy reconstruction scheme using the geometrical form of the luminosity-distance relation_ , _Physical Review D_ 74 (2006) .
  * [42] S. Benitez-Herrera, F. Röpke, W. Hillebrandt, C. Mignone, M. Bartelmann and J. Weller, _Model-independent reconstruction of the expansion history of the universe from type ia supernovae_ , _Monthly Notices of the Royal Astronomical Society_ 419 (2011) 513–521.
  * [43] T. Holsclaw, U. Alam, B. Sansó, H. Lee, K. Heitmann, S. Habib et al., _Nonparametric dark energy reconstruction from supernova data_ , _Physical Review Letters_ 105 (2010) .
  * [44] A. Montiel, R. Lazkoz, I. Sendra, C. Escamilla-Rivera and V. Salzano, _Nonparametric reconstruction of the cosmic expansion with local regression smoothing and simulation extrapolation_ , _Physical Review D_ 89 (2014) .
  * [45] R. C. Bernardo and J. Levi Said, _Towards a model-independent reconstruction approach for late-time Hubble data_ , 2106.08688.
  * [46] R.-Y. Guo, J.-F. Zhang and X. Zhang, _Can the $H_{0}$ tension be resolved in extensions to $\Lambda$CDM cosmology?_, _JCAP_ 02 (2019) 054 [1809.02340].
  * [47] G. Lambiase, S. Mohanty, A. Narang and P. Parashari, _Testing dark energy models in the light of $\sigma_{8}$ tension_, _Eur. Phys. J. C_ 79 (2019) 141 [1804.07154].
  * [48] N. Rani, D. Jain, S. Mahajan, A. Mukherjee and N. Pires, _Transition Redshift: New constraints from parametric and nonparametric methods_ , _JCAP_ 12 (2015) 045 [1503.08543].
  * [49] D. J. C. MacKay, _Information Theory, Inference & Learning Algorithms_. Cambridge University Press, USA, 2002.
  * [50] V. C. Busti, C. Clarkson and M. Seikel, _The Value of $H_{0}$ from Gaussian Processes_, _IAU Symp._ 306 (2014) 25 [1407.5227].
  * [51] V. C. Busti, C. Clarkson and M. Seikel, _Evidence for a Lower Value for $H_{0}$ from Cosmic Chronometers Data?_, _Mon. Not. Roy. Astron. Soc._ 441 (2014) 11 [1402.5429].
  * [52] L. Verde, P. Protopapas and R. Jimenez, _The expansion rate of the intermediate Universe in light of Planck_ , _Phys. Dark Univ._ 5-6 (2014) 307 [1403.2181].
  * [53] Z. Li, J. E. Gonzalez, H. Yu, Z.-H. Zhu and J. S. Alcaniz, _Constructing a cosmological model-independent Hubble diagram of type Ia supernovae with cosmic chronometers_ , _Phys. Rev. D_ 93 (2016) 043014 [1504.03269].
  * [54] A. Shafieloo, A. G. Kim and E. V. Linder, _Gaussian Process Cosmography_ , _Phys. Rev. D_ 85 (2012) 123530 [1204.2272].
  * [55] R.-G. Cai, Z.-K. Guo and T. Yang, _Null test of the cosmic curvature using $H(z)$ and supernovae data_, _Phys. Rev. D_ 93 (2016) 043517 [1509.06283].
  * [56] D. Wang and X.-H. Meng, _Improved constraints on the dark energy equation of state using Gaussian processes_ , _Phys. Rev. D_ 95 (2017) 023508 [1708.07750].
  * [57] L. Zhou, X. Fu, Z. Peng and J. Chen, _Probing the Cosmic Opacity from Future Gravitational Wave Standard Sirens_ , _Phys. Rev. D_ 100 (2019) 123539 [1912.02327].
  * [58] P. Mukherjee and N. Banerjee, _Revisiting a non-parametric reconstruction of the deceleration parameter from observational data_ , 2007.15941.
  * [59] M.-J. Zhang and H. Li, _Gaussian processes reconstruction of dark energy from observational data_ , _Eur. Phys. J. C_ 78 (2018) 460 [1806.02981].
  * [60] M. Aljaf, D. Gregoris and M. Khurshudyan, _Constraints on interacting dark energy models through cosmic chronometers and Gaussian process_ , 2005.01891.
  * [61] K. Liao, A. Shafieloo, R. E. Keeley and E. V. Linder, _A model-independent determination of the Hubble constant from lensed quasars and supernovae using Gaussian process regression_ , _Astrophys. J. Lett._ 886 (2019) L23 [1908.04967].
  * [62] H. Yu, B. Ratra and F.-Y. Wang, _Hubble Parameter and Baryon Acoustic Oscillation Measurement Constraints on the Hubble Constant, the Deviation from the Spatially Flat $\Lambda$CDM Model, the Deceleration–Acceleration Transition Redshift, and Spatial Curvature_, _Astrophys. J._ 856 (2018) 3 [1711.03437].
  * [63] C. Krishnan, E. Colgáin, M. M. Sheikh-Jabbari and T. Yang, _Running Hubble Tension and a H0 Diagnostic_ , 2011.02858.
  * [64] M. Moresco, _Raising the bar: new constraints on the Hubble parameter with cosmic chronometers at z $\sim$ 2_, _Mon. Not. Roy. Astron. Soc._ 450 (2015) L16 [1503.01116].
  * [65] M. Moresco et al., _A 6% measurement of the Hubble parameter at $z\sim 0.45$: direct evidence of the epoch of cosmic re-acceleration_, _JCAP_ 05 (2016) 014 [1601.01701].
  * [66] D. Stern, R. Jimenez, L. Verde, M. Kamionkowski and S. A. Stanford, _Cosmic chronometers: constraining the equation of state of dark energy. I: H(z) measurements_ , _JCAP_ 2010 (2010) 008 [0907.3149].
  * [67] M. Moresco, A. Cimatti, R. Jimenez, L. Pozzetti, G. Zamorani, M. Bolzonella et al., _Improved constraints on the expansion rate of the Universe up to z ~1.1 from the spectroscopic evolution of cosmic chronometers_, _JCAP_ 2012 (2012) 006 [1201.3609].
  * [68] C. Zhang, H. Zhang, S. Yuan, S. Liu, T.-J. Zhang and Y.-C. Sun, _Four new observational H(z) data from luminous red galaxies in the Sloan Digital Sky Survey data release seven_ , _Research in Astronomy and Astrophysics_ 14 (2014) 1221 [1207.4541].
  * [69] BOSS collaboration, _Quasar-Lyman $\alpha$ Forest Cross-Correlation from BOSS DR11 : Baryon Acoustic Oscillations_, _JCAP_ 05 (2014) 027 [1311.1767].
  * [70] BOSS collaboration, _Baryon acoustic oscillations in the Ly $\alpha$ forest of BOSS DR11 quasars_, _Astron. Astrophys._ 574 (2015) A59 [1404.1801].
  * [71] J. E. Bautista et al., _Measurement of baryon acoustic oscillation correlations at $z=2.3$ with SDSS DR12 Ly$\alpha$-Forests_, _Astron. Astrophys._ 603 (2017) A12 [1702.00176].
  * [72] BOSS collaboration, _The clustering of galaxies in the SDSS-III Baryon Oscillation Spectroscopic Survey: baryon acoustic oscillations in the Data Releases 10 and 11 Galaxy samples_ , _Mon. Not. Roy. Astron. Soc._ 441 (2014) 24 [1312.4877].
  * [73] C. Blake, S. Brough, M. Colless, C. Contreras, W. Couch, S. Croom et al., _The WiggleZ Dark Energy Survey: joint measurements of the expansion and growth history at z ¡ 1_ , _Mon. Not. Roy. Astron. Soc._ 425 (2012) 405 [1204.3674].
  * [74] S. Alam, M. Ata, S. Bailey, F. Beutler, D. Bizyaev, J. A. Blazek et al., _The clustering of galaxies in the completed SDSS-III Baryon Oscillation Spectroscopic Survey: cosmological analysis of the DR12 galaxy sample_ , _Mon. Not. Roy. Astron. Soc._ 470 (2017) 2617 [1607.03155].
  * [75] C.-H. Chuang and Y. Wang, _Modelling the anisotropic two-point galaxy correlation function on small scales and single-probe measurements of H(z), DA(z) and $f(z)\sigma_{8}(z)$ from the Sloan Digital Sky Survey DR7 luminous red galaxies_, _Mon. Not. Roy. Astron. Soc._ 435 (2013) 255 [https://academic.oup.com/mnras/article-pdf/435/1/255/3865066/stt1290.pdf].
  * [76] BOSS collaboration, _The clustering of galaxies in the completed SDSS-III Baryon Oscillation Spectroscopic Survey: tomographic BAO analysis of DR12 combined sample in configuration space_ , _Mon. Not. Roy. Astron. Soc._ 469 (2017) 3762 [1607.03154].
  * [77] A. Oka, S. Saito, T. Nishimichi, A. Taruya and K. Yamamoto, _Simultaneous constraints on the growth of structure and cosmic expansion from the multipole power spectra of the SDSS DR7 LRG sample_ , _Mon. Not. Roy. Astron. Soc._ 439 (2014) 2515 [1310.2820].
  * [78] E. Gaztanaga, A. Cabre and L. Hui, _Clustering of Luminous Red Galaxies IV: Baryon Acoustic Peak in the Line-of-Sight Direction and a Direct Measurement of H(z)_ , _Mon. Not. Roy. Astron. Soc._ 399 (2009) 1663 [0807.3551].
  * [79] J. Magana, M. H. Amante, M. A. Garcia-Aspeitia and V. Motta, _The Cardassian expansion revisited: constraints from updated Hubble parameter measurements and type Ia supernova data_ , _Mon. Not. Roy. Astron. Soc._ 476 (2018) 1036 [1706.09848].
  * [80] R. C. Nunes, S. K. Yadav, J. F. Jesus and A. Bernui, _Cosmological parameter analyses using transversal BAO data_ , _Mon. Not. Roy. Astron. Soc._ 497 (2020) 2133 [2002.09293].
  * [81] D. M. Scolnic et al., _The Complete Light-curve Sample of Spectroscopically Confirmed SNe Ia from Pan-STARRS1 and Cosmological Constraints from the Combined Pantheon Sample_ , _Astrophys. J._ 859 (2018) 101 [1710.00845].
  * [82] W. L. Freedman, B. F. Madore, D. Hatt, T. J. Hoyt, I. S. Jang, R. L. Beaton et al., _The Carnegie-Chicago Hubble Program. VIII. An independent determination of the Hubble constant based on the tip of the red giant branch_ , _The Astrophysical Journal_ 882 (2019) 34.
  * [83] K. C. Wong, S. H. Suyu, G. C. Chen, C. E. Rusu, M. Millon, D. Sluse et al., _H0LiCOW–XIII. A 2.4 percent measurement of H0 from lensed quasars: 5.3 $\sigma$ tension between early-and late-Universe probes_, _Monthly Notices of the Royal Astronomical Society_ 498 (2020) 1420.
  * [84] A. G. Riess et al., _A 2.4% Determination of the Local Value of the Hubble Constant_ , _Astrophys. J._ 826 (2016) 56 [1604.01424].
  * [85] H. W. Chiang, A. E. Romano, F. Nugier and P. Chen, _Probing homogeneity with standard candles_ , _JCAP_ 11 (2019) 016 [1706.09734].
  * [86] E. Di Valentino, O. Mena, S. Pan, L. Visinelli, W. Yang, A. Melchiorri et al., _In the realm of the hubble tension—a review of solutions *_ , _Classical and Quantum Gravity_ 38 (2021) 153001.
  * [87] E. Di Valentino et al., _Snowmass2021 - Letter of interest cosmology intertwined II: The hubble constant tension_ , _Astropart. Phys._ 131 (2021) 102605 [2008.11284].
  * [88] W. L. Freedman, _Measurements of the Hubble Constant: Tensions in Perspective_ , 2106.15656.
  * [89] G. Alestas and L. Perivolaropoulos, _Late-time approaches to the Hubble tension deforming H(z), worsen the growth tension_ , _Mon. Not. Roy. Astron. Soc._ 504 (2021) 3956 [2103.04045].