# On scaling limits of random Halin-like maps with large faces

Daniel Amankwah  Mathematics Division, Science Institute, University of
Iceland, Dunhaga 5, 107 Reykjavík, Iceland<EMAIL_ADDRESS>and  Sigurdur Örn
Stefánsson<EMAIL_ADDRESS>
###### Abstract.

We consider maps which are constructed from plane trees by connecting each
pair of consecutive leaves on their contour by a single edge. When every non-
leaf vertex of the underlying tree has exactly one child which is a leaf,
these maps are bijectively related to plane trees with marked corners. We
study the scaling limits of such maps where Boltzmann weights are assigned to
each face. The main result is that when the degree of a typical face is in the
domain of attraction of an $\alpha$-stable law, with $\alpha\in(1,2)$, the
maps, properly rescaled, converge to the $\alpha$-stable looptrees of Curien
and Kortchemski.

##  1\. Introduction

A planar map is an embedding of a plane graph onto a 2-dimensional sphere up
to continuous deformations. We refer to a planar map as rooted if it has a
distinguished and oriented edge. Rooting a planar map makes it easier to
handle combinatorially since it eliminates symmetry issues. Also, the root
gives a starting point for recursive decomposition. The connected components
of the complement of a planar map are called faces. If $f$ is a face of a
planar map, we will denote its degree (the number of edge sides that surround
it, counted with multiplicity) by $\mathsf{deg}(f)$. A rooted map with a
single face is called a plane tree.

Research into asymptotic behavior of random planar maps has been active in the
last two decades. This has been motivated in part by the connection with two-
dimensional Liouville quantum gravity(see [16, 18] together with references
for detailed account). Much of the advances in this field have been made
possible by the use of bijections between specific classes of planar maps and
certain plane trees. See for example [9, 19, 7]. This is due to the fact that
plane trees are easier to analyse since they are coded by functions on which
one can apply the standard tools from analysis and probability theory. In [2,
3, 4], Aldous gives an exposition on the scaling limit of various classes of
discrete plane trees conditioned to be large. For a Galton-Watson tree
conditioned to be large, whose offspring distribution has mean $1$ and a
finite variance, he showed that the scaling limit is a continuous random tree
called the Brownian continuum random tree (CRT). The CRT admits a universality
property in the sense that it is the limit of several classes of trees and
maps. This includes random dissections of polygons [10] which is of interest
in this work. Generalizations to Galton-Watson trees whose offspring
disribution is in the domain of attracion of a stable law with index
$\alpha\in(1,2)$ were given by Duquesne and Le-Gall (see [13, 12]). They
showed that these trees, with a properly rescaled graph distance admit a weak
limit which is referred to as the $\alpha$-stable tree,
$\mathcal{T}_{\alpha}$. The $\alpha$-stable trees are qualitatively different
from the Brownian CRT, and have the feature that they possess a countably
infinite number of vertices of infinite degree.

Curien and Kortchemski introduced the study of a family of random compact
metric spaces $(\mathcal{L}_{\alpha})_{1<\alpha<2}$ called $\alpha$-stable
looptrees (see [11]). On the discrete level, looptrees are graphs constructed
on the vertex set of a plane tree by connecting each vertex to its leftmost
child by an edge, connecting adjacent siblings by edges, and finally
connecting the rightmost child to its parent by an edge (see Fig. 1). The
space $\mathcal{L}_{\alpha}$ may be informally described in the same spirit by
replacing each vertex of large degree in the $\alpha$-stable tree
$\mathcal{T}_{\alpha}$, by a circle of length proportional to this degree.

Figure 1. A looptree associated to a plane tree.

Curien and Kortchemski gave sufficient conditions for which a random sequence
of discrete looptrees (see figure 1) converges in distribution to the stable
looptree and as an application showed that the random dissection of a polygon
sampled according to a Boltzmann distribution on the face degrees, converges
to the stable looptree for the Gromov-Hausdorff topology. Furthermore,
Stefánsson and Stufler used this invariance principle to show that Boltzmann
weighted outerplanar maps (maps whose blocks are dissections of a polygon) may
also converge to stable looptrees [20].

These examples have the common feature that the objects involved are tree-like
and when their faces are forced to have large degrees, in the limit, they
appear as circles arranged along a tree structure. We are interested in
understanding how generic this behaviour is and provide in this paper an
additional example of a family of maps which shares this feature of being
tree-like and converges to the stable looptrees.

The maps we consider in this work are defined as follows: They are two-
connected, each bounded face in the map shares exactly one edge with the
unbounded face and each boundary vertex has degree $3$. The _root edge_ is
chosen as as an oriented edge which does not belong to the unbounded face. The
tail of the root edge is referred to as the _root vertex_. The internal face
to the left of the root edge is called the _root face_ and denoted by $f_{0}$.
Furthermore, we assume for notational convenience, that a half edge is
attached to the root vertex and drawn inside $f_{0}$, see Fig. 2. For
comparison, an outerplanar map has the property that each vertex lies on the
boundary of the unbounded face and its two-connected components are
dissections of a polygon.

There is a natural bijection between our object of interest and a rooted plane
tree which is described by deleting all edges on the boundary of the map to
obtain a plane tree and inversely by connecting consecutive leaves of the
plane tree in a cyclic order to obtain the map.111This description is in the
spirit of the construction of graph from trees with a given surplus of edges
(see [1]), where in the current case the number of leaves in the tree are
exactly the same as the surplus (the number of edges needed to be removed to
obtain a tree) Figure 2 illustrates this bijection. We will refer to this tree
as the _underlying tree_. When the underlying tree has at least four edges and
no vertex of degree two, the corresponding map has been referred to as a Halin
graph in the literature. We will use the term Halin map in the following,
keeping in mind this minor distinction. A generic Halin map will be denoted by
$H$ and its underlying tree by $H^{\circ}$.

We make one further assumption on the maps which is best described as a
condition on the underlying tree:

($H_{\star}$) The underlying tree $H^{\circ}$ of the Halin map $H$ consist of
vertices that are either leaves or have exactly one child which is a leaf.

For the Halin map, condition $(H_{\star})$ means that every vertex of the map
either lies on the boundary of the unbounded face or has exactly one neighbour
on the boundary of the unbounded face.

Another relevant bijection is obtained by considering the weak dual of the
Halin map $H$, namely the usual geometric dual but with the vertex
corresponding to the external face removed. This weak dual is always a
dissection of a polygon $D$. If $v_{0}$ is the vertex in $D$ dual to $f_{0}$,
then $D$ is rooted at the leftmost internal edge emanating from $v_{0}$.
Furthermore, for this duality to be a bijection, the corner of $v_{0}$ which
lies in the face of $D$ dual to the root vertex in $H$, receives a mark (see
Fig. 2, middle). The condition ($H_{\star}$) makes our analysis simpler due to
properties of this dissection. Under ($H_{\star}$), the map obtained by
removing the edges of $D$ is a tree $T$, whose vertices correspond to faces of
the Halin map. This tree has been studied before in the literature under the
name _non-crossing tree_ , see e.g. [15] in the context of scaling limits. We
note that our rooting convention is different from that of [15] which changes
the combinatorics slightly.

Figure 2. Left: A rooted Halin map $H$. The underlying tree $H^{\circ}$ is
indicated by solid lines. A half edge is attached to the root vertex. Middle:
The dual marked dissection $D$ in black. Its boundary is dotted and the solid
lines indicate the non-crossing tree $T$. Right: The non-crossing tree $T$.
Each corner receives a mark which indicates how to reconstruct the marked
dissection from the tree.

The non-crossing tree $T$ is bijectively reletated to the original Halin map
if one assigns marks to the corners of the tree which allow to reconstruct the
marked dissection, see Fig. 2 and precise construction in Section 2.2.

######  Remark 1.1.

We do not attempt to relax condition $(H_{\star})$ in the current work. We
note however, that by assuming that each non-leaf vertex has _at least_ one
child which is a leave, gives a dual dissection with the property that if its
boundary is removed one gets a forest of non-crossing trees (this follows from
the construction in Section 2.2). We plan to consider this case in more detail
in future work.

In the most general case where non-leaf vertices could have no children which
are leaves, the dual dissection becomes even more complicated and we do not
currently know how to solve this case.

###  1.1. The main result

For every integer $n\geq 1$, denote by $\mathbb{H}_{n}$ the set of Halin maps
satisfying $(H_{\star})$ and having $2n$ vertices. Let $H\in\mathbb{H}_{n}$
and fix a sequence $w=(w(k):k\geq 4)$ of non-negative _face weights_. We
define the weight of the Halin map $H$ by

$W(H)=\prod_{\begin{subarray}{c}f\;\text{face
of}\;H\end{subarray}}w(\mathsf{deg}(f)).$

When we assign the weight to the root face $f_{0}$, we count the half-edge as
having contribution one to the degree of $f_{0}$. This is only for notational
convenience later on and this shift does not significantly affect the results.

We define the normalization

$Z_{n}=\sum_{H^{\prime}\in\mathbb{H}_{n}}W(H^{\prime})$

and if $Z_{n}>0$, we define a random element $\mathcal{H}_{n}^{w}$ in
$\mathbb{H}_{n}$ by

$\mathbb{P}(\mathcal{H}_{n}^{w}=H)=\frac{1}{Z_{n}}W(H)\hskip
28.45274pt\text{for all}\;H\in\mathbb{H}_{n}.$ (1)

We view $\mathcal{H}_{n}^{w}$ as a random metric space equipped with the graph
distance as its metric (where every edge of $\mathcal{H}_{n}^{w}$ has unit
length).

We will assume throughout this paper that the generating function of the
sequence $w$ has a finite radius of convergence. This allows us to define a
probability measure, $\mu$ on $\mathbb{Z}_{+}=\\{0,1,2,\ldots\\}$ by
$\mu(k)=ab^{k}(k+1)w(k+4)$ for appropriate choice of $a,b>0$. In order to
guarantee that $Z_{n}>0$ for sufficiently large $n$ we assume that $\mu(k)$ is
aperiodic. The general case may also be considered but the results are
essentially the same and we aim to keep the notation simple.

The probability measure $\mu$ is said to be in the domain of attraction of an
$\alpha$-stable law, with $\alpha\in(1,2)$, if

$\mu([j,\infty))=j^{-\alpha}L(j),$

where $L:\mathbb{R}_{+}\rightarrow\mathbb{R}_{+}$ is a slowly varying
function, i.e. $L(x)>0$ and $\lim_{x\rightarrow\infty}\frac{L(tx)}{L(x)}=1$
for all $t>0$ (see [6] for details). In this case, the variance of $\mu$ is
infinite. Common examples of slowly varying functions are constant functions,
powers of logarithms, and iterated logarithms.

The main result of the current work is the following scaling limit theorem for
$\mathcal{H}_{n}^{w}$.

######  Theorem 1.2.

Fix $\alpha\in(1,2)$ and let $\mu$ be a probability measure on
$\mathbb{Z}_{+}$ which has mean $1$ and belongs to the domain of attraction of
an $\alpha$-stable law. Then there is a sequence of positive numbers
$b_{n}\to\infty$ such that the following convergence holds in distribution for
the Gromov-Hausdorff topology

$b_{n}^{-1}\cdot\mathcal{H}_{n}^{w}\xrightarrow[n\rightarrow\infty]{\enskip(d)\enskip}\mathcal{L}_{\alpha}.$

The notation $s\cdot\mathcal{H}_{n}^{w}$ represents the metric space obtained
from $\mathcal{H}_{n}^{w}$ by multiplying all distances by $s>0$. By [14,
Theorem 1.10], it holds that

$\displaystyle b_{n}=n^{1/\alpha}\Lambda(n)$ (2)

where $\Lambda$ is a slowly varying function. An explicit example of a measure
$\mu$ which satisfies the conditions of the theorem is

$\displaystyle\mu(k)=ck^{-1-\alpha},\qquad c>0$

in which case

$\displaystyle
b_{n}=\frac{n^{\frac{1}{\alpha}}}{(c|\Gamma(-\alpha)|)^{\frac{1}{\alpha}}}.$

We refer the reader to e.g. [11, 14] for more information on random variables
in the domain of attraction of stable laws, in a relevant context.

######  Remark 1.3.

When $\mu$ is critical and has finite variance, we expect that
$n^{-1/2}\mathcal{H}_{n}^{w}$ converges weakly towards a multiple of the
Brownian CRT. This includes for example the case when $\mathcal{H}_{n}^{w}$ is
uniformly distributed on $\mathbb{H}_{n}$, in which case $w_{k}=1$ for all
$k\geq 0$. However, in order to prove this, it seems one needs to control very
precisely the geodesics in $\mathcal{H}_{n}^{w}$, in a similar way as is done
in the case of dissections in [10]. We do not pursue this further in the
current work and leave it as an open question.

##  2\. Halin maps encoded by size conditioned Galton-Watson trees

In this section we define plane trees and Halin maps and describe a useful
bijection between these objecs. The bijection allows us to related the random
Halin map $\mathcal{H}_{n}^{w}$ to a size conditioned Galton-Watson tree.

###  2.1. Rooted plane trees and Halin maps

We briefly recall the formalism of rooted plane trees (or rooted ordered
trees). Let $\mathbb{N}=\\{1,2,...\\}$ and define the set of labels

$\mathcal{U}=\bigcup_{n=0}^{\infty}\mathbb{N}^{n},$

where by convention $\mathbb{N}^{0}=\\{\emptyset\\}$. An element of
$\mathcal{U}$ is thus a sequence $u=(u^{1},...,u^{n})$ of elements of
$\mathbb{N}$ and we set $|u|=n$, so that $|u|$ represents the generation of
$u$. If $u=(u^{1},...,u^{m})$ and $v=(v^{1},...,v^{n})$ belong to
$\mathcal{U}$, we write the concatenation of $u$ and $v$ as
$uv=(u^{1},...,u^{m},v^{1},...,v^{n})$.

A plane tree $T$ is a finite or infinite subset of $\mathcal{U}$ such that,

  1. (1)

$\emptyset\in T$.

  2. (2)

if $v\in T$ and $v=uj$ for some $j\in\mathbb{N}$ then $u\in T$.

  3. (3)

for every $u\in T$, there exists an integer $k_{u}(T)\geq 0$ ( the number of
children of $u$) such that for every $j\in\mathbb{N}$, $uj\in T$ if and only
if $1\leq j\leq k_{u}(T)$.

The vertex, $\emptyset$ is called the root vertex and the edge
$\\{\emptyset,(1)\\}$ is called the root edge. Let us denote the set of all
plane trees by T and the set of plane trees with $n$ vertices by
$\textbf{T}_{n}$ . For $T\in\textbf{T}$, we will view each of its vertex as an
individual of a population whose $T$ is the geneological tree. If
$T\in\textbf{T}$ and $u\in T$, we define the shift of $T$ at $u$ by
$\sigma_{u}T=\\{v\in\mathcal{U}:uv\in T\\}$ which is itself a tree. The total
progeny of $T$(total number of vertices of $T$) will be denoted by $\zeta(T)$.
The number of leaves of the tree $T$ is denoted by $\lambda(T)$. The maximal
generation of the tree $T$ is called its height and we denote it by
$\mathsf{Height}(T)$.

######  Definition 2.1.

Let $\rho$ be a probability measure on $\mathbb{Z}_{+}$ such that $\rho(1)<1$.
We call $\rho$, the offspring distribution. The Galton-Watson measure with
offspring distribution $\rho$ is the probability measure, $GW_{\rho}$ on T
such that:

  1. (1)

$GW_{\rho}(k_{\emptyset}=j)=\rho(j)$ for $j\geq 0$,

  2. (2)

If $\tau$ is distrubted by $GW_{\rho}$, then for every $j\geq 1$ with
$\rho(j)>0$, conditionally on $\\{k_{\emptyset}=j\\}$, the subtrees
$\sigma_{1}\tau...\sigma_{j}\tau$ are i.i.d. with distribution $GW_{\rho}$.

If $\rho$ has mean less than or equal to $1$ then a $\rho$-Galton-Watson tree,
$\tau$ is almost surely finite and

$GW_{\rho}(\tau=T)=\prod_{u\in T}\rho(k_{u}(\tau)).$ (3)

(See [17], Proposition 1.4).

######  Definition 2.2.

Let $T$ be a rooted plane tree and let $\ell_{1},\cdots,\ell_{\lambda(T)}$ be
its leaves in a lexicographical order. A Halin map $H$, is a map constructed
by adding to $T$ the edges $(\ell_{i},\ell_{i+1})$, $1\leq i\leq\lambda(T)-1$,
and the edge $(\ell_{\lambda(T)},\ell_{1})$. We refer to $T$ as the
_underlying tree_ and use the notation $H^{\circ}:=T$. The root edge $H$ is
the same as in $T$ and a half edge is attached to its tail. The Halin maps
considered in this work must satisfy condition $(H_{\star})$ in the
introduction. The set of Halin maps with $2n$ vertices, satisfying
$(H_{\star})$, is denoted by $\mathbb{H}_{n}$.

###  2.2. The bijection

We now describe a bijection between the set of Halin maps and the set of plane
trees with marked corners. This allows us to relate random Halin maps to size-
conditioned Galton-Watson trees which gives us the necessary tools to prove
Theorem 1.2. The marked plane trees are defined as followed Let
$T\in\textbf{T}_{n}$. Consider its vertices $\emptyset=v(0)<v(1)<...<v(n-1)$
in lexicographical order and set

$\mathbb{M}(T)\textasciidieresis=\\{(l_{0},l_{1},l_{2},...,l_{n-1}):0\leq
l_{j}\leq k_{v(j)}(T)\;\text{for every}\;0\leq j\leq n-1\\}$

and

$\textbf{T}_{n}^{\ast}=\\{(T,m):T\in\textbf{T}_{n}\,\text{and}\,m\in\mathbb{M}(T)\\}.$

We call elements of $\textbf{T}_{n}^{\ast}$, _marked trees_. Given a marked
tree $T^{\ast}=(T,m)$ we call $T$ its _shape_ and denote it by
$\mathcal{S}(T^{\ast})=T$. Note that leaves always receive the mark 0 and we
will therefore not indicate their mark in figures.

We now construct a function $\phi_{n}:\mathbb{H}_{n}\to\mathbf{T}_{n}^{\ast}$
in several steps as follows. Let $H\in\mathbb{H}_{n}$:

$\bullet$ A weak dual $D$ of $H$ is obtained by placing a vertex inside each
bounded face of $H$ and joining two vertices by an edge if their corresponding
faces share an edge which is not incident to the unbounded face of the Halin
map.

$\bullet$ The week dual $D$ is a dissection of a polygon: Every vertex in $D$
lies on the boundary of its unbounded face since every face of $H$ shares an
edge with its unbounded face. $D$ is connected since $H$ is two-connected, and
$D$ is two-connected since each edge in $H$ shares exactly one edge with its
unbounded face.

$\bullet$ When the edges of $D$ which lie on the unbounded face are removed,
one obtaines a map $T$ which is plane tree:

First of all, $T$ is connected. Assume the opposite, then there exists a face
$F$ in $D$ with a boundary of the form $e_{1}P_{1}e_{2}P_{2}$ where $e_{1}$
and $e_{2}$ are edges inside $D$ (one of which is possible empty) and $P_{1}$
and $P_{2}$ are paths of length $\geq 1$ on the boundary of $D$. This violates
$(H_{\star})$ since the vertex $v$ dual to $F$ in $H^{\circ}$ would have at
least two children which are leaves (the vertices connected to $v$ through the
edges in $P_{1}$ and $P_{2}$). See Fig. 3 (left).

Secondly, $T$ has no cycles. If $T$ had a cycle it would be the boundary of a
face $F$ with perimeter of length at least 2 which shared no edges with the
boundary of $D$. Then, its dual vertex $v$ in $H^{\circ}$ would be a non-leaf
and none of its children would be leaves, violating $(H_{\star})$, see Fig. 3
(middle).

$\bullet$ A root is assigned to $T$ as follows: The root vertex $\varnothing$
in $T$ is chosen to be the vertex dual to the face $f_{0}$. The directed root
edge is chosen as the leftmost edge eminating from $\varnothing$ which lies
inside $D$.

$\bullet$ Marks are assigned to the corners of $T$ as follows, yielding a
marked tree $T^{\ast}=(T,m)$: Firstly, we mark a corner around the root. If no
edge in $T$ is dual to the root edge in $H$ the mark is taken to be $0$.
Otherwise, if $(\varnothing,(i))$ is the edge in $T$ which is dual to the root
edge in $H$, then the mark is taken to be $i$. For other vertices, the marked
corner is chosen as the corner which is adjacent to the unbounded face of $D$,
see Fig. 3 (right). Note that the mark on the root depends on the location of
the root of the Halin map $H$ whereas the marks on other vertices are
completely determined by the dissection $D$.

It is easily shown by induction that if $H$ contains $2n$ vertices, then $T$
contains $n$ vertices. Finally, we let $\phi_{n}(H)=(T,m)$.

Figure 3. Illustration of the fact that the dissection $D$, with its boundary
removed, is a tree $T$ and an explanation of how the corners of $T$ receive
marks. The boundary of the dissection is in grey and the edges in the tree $T$
are in black. The edges in the corresponding Halin map are dotted and its
vertices are white. Left: $T$ is connected. Middle: $T$ has no cycles. Right:
Marks assigned to the corners of $T$.

Since every step in the of $\phi_{n}$ can easily be seen to be reversible we
have shown that

######  Proposition 2.3.

The function $\phi_{n}:\mathbb{H}_{n}\to\textbf{T}_{n}^{\ast}$ is a bijection.

The bijection enjoys the useful property that the faces in the Halin map $H$
correspond to vertices in the marked tree $T^{\ast}$ which follows from the
fact that $D$ is the weak dual of $H$. Furthermore, the degree of a vertex $v$
in $D$ is equal to the degree of the corresponding face $f$ in $H$, minus 1
since we consider the weak dual. The degree of $v$ in $T^{\ast}$ is 2 smaller
than its degree in $D$ since it is adjacent to two boundary edges of $D$. We
conclude that $k_{v}(T^{\ast})=\mathsf{deg}(f)-4$. This is also true for the
root vertex in $T^{\ast}$, since we added a half-edge to the root vertex in
$H$ which contributes 1 to the degree of the root face (in fact, this is the
reason why we add the half-edge).

######  Proposition 2.4.

Suppose $\mathcal{T}_{n}^{\ast,w}$ is the random marked tree which corresponds
to the Halin map $\mathcal{H}_{n}^{w}$ via the bijection $\phi_{n}$. Then,
$\mathcal{T}_{n}^{w}=\mathcal{S}(\mathcal{T}_{n}^{\ast,w})$, is a $GW_{\mu}$
tree conditioned on having $n$ vertices, where

$\displaystyle\mu(k)=ab^{k}(k+1)w(k+4)$

with $a,b>0$ chosen such that $\mu$ is a probability distribution.

###### Proof.

Let $H\in\mathbb{H}_{n}$, $T^{\ast}=\phi_{n}(H)$ and
$T=\mathcal{S}(T^{\ast})$. Then

$\displaystyle\mathbb{P}(\mathcal{T}_{n}^{w}=T)$
$\displaystyle=\sum_{\tau^{\ast}:~{}\mathcal{S}(\tau^{\ast})=T}\mathbb{P}(\mathcal{T}_{n}^{\ast,w}=\tau^{\ast})$
$\displaystyle=\mathbb{P}(\mathcal{T}_{n}^{\ast,w}=T^{\ast})\prod_{v\;\text{vertex
in}\;T^{\ast}}(k_{v}(T^{\ast})+1)$
$\displaystyle=\mathbb{P}(\mathcal{H}_{n}^{w}=H)\prod_{v\;\text{vertex
in}\;T^{\ast}}(k_{v}(T^{\ast})+1)$
$\displaystyle=Z_{n}^{-1}\prod_{\begin{subarray}{c}f\;\text{face
of}\;H\end{subarray}}w(\mathsf{deg}(f))\prod_{v\;\text{vertex
in}\;T^{\ast}}(k_{v}(T^{\ast})+1)$
$\displaystyle=Z_{n}^{-1}a^{-n}b^{-n+1}\prod_{\begin{subarray}{c}v\;\text{vertex
in}\;T^{\ast}\end{subarray}}ab^{k_{v}(T^{\ast})}(k_{v}(T^{\ast})+1)w(k_{v}(T^{\ast})+4)$
$\displaystyle=Z_{n}^{-1}a^{-n}b^{-n+1}\prod_{\begin{subarray}{c}v\;\text{vertex
in}\;T^{\ast}\end{subarray}}\mu(k)=GW_{\mu}(\tau=T~{}|~{}\zeta(\tau)=n).$

The product in the second line counts the number of marked trees which have
the shape $T$ and follows from the fact that the marks are assigned uniformly
to $T$. ∎

##  3\. Invariance principle for discrete looptrees

Our main result is an application of the Invariance principles for discrete
looptrees given by Curien and Kortchemski in [11]. In the following subsection
we give the description of the required material in order to state their
theorem.

###  3.1. Lukasiewicz Path

Let $T$ be a plane tree whose vertices are listed in depth first-search-
order(or lexicographical order), $\emptyset=v(0)<v(1)<...<v(\zeta(T)-1)$. The
Lukasiewicz path $W(T)=(W_{n}(T),0\leq n\leq\zeta(T))$ is defined by
$W_{0}(T):=0$ and

$W_{n}(T):=\sum_{j=0}^{n-1}(k_{v(j)}-1)$

for $1\leq n\leq\zeta(T)$. Figure 4 gives an example of a tree encoded by a
Lukasiewicz path. $W$ is viewed as a a càdlàg functions (right continuous with
left limits) by assuming it is constant in between successive integers. It may
be seen that $W_{i}(T)\geq 0$ for $0\leq i\leq\zeta(T)-1$, but
$W_{\zeta(T)}(T)=-1$.

Figure 4. A tree encoded by a Lukasiewicz path

Consider a sequence ($\tau_{n}$) distributed by
$GW\mu(\cdot~{}|~{}\zeta(\tau)=n)$. If $\mu$ is in the domain of attraction of
a stable law with index $\alpha\in(1,2)$, the properly rescaled Lukasiewicz
path $(W_{\lfloor tn\rfloor}(\tau_{n}))_{t\in[0,1]}$, converges weakly
torwards the excursion of the $\alpha$-stable spectrally positive Lévy process
denoted by $X^{exc,\alpha}$(see [12, 13]). Background information on Lévy
process can be found in [5] and for more details, see [11]. It is important to
note that this process lives on the Skorokhod space $D([0,1],\mathbb{R})$ of
real valued càdlàg functions.

###  3.2. Invariance principles for discrete looptrees

A discrete looptree is constructed on the vertex set of a plane tree, $T$ as
follows: For each vertex $v$ of outdegree $k\geq 1$, add the edges $(v,v1)$
and $(v,vk)$ (with repetition if $k=1$), and add the edges $(vi,v(i+1))$, for
$i=1,\ldots,k-1$. See Fig. 1 for an illustration. We denote the associated
looptree by $\mathsf{Loop}(T)$. $\mathsf{Loop}(T)$ will be viewed as a metric
space equipped with the graph distance as its metric (each edge with length
one). The following invariance theorem gives a sufficient condition on a
sequence of random plane trees $(\tau_{n})_{n\geq 1}$ that ensures that the
associated looptrees, $(\mathsf{Loop}(\tau_{n}))_{n\geq 1}$ appropriately
rescaled, converge towards the $\alpha$-stable looptree
$\mathcal{L}_{\alpha}$. Informally, one can view the $\alpha$-stable looptree
$\mathcal{L}_{\alpha}$ as ”$\mathsf{Loop}(\mathcal{T}_{\alpha})$”, where
$\mathcal{T}_{\alpha}$ is the $\alpha$-stable tree mentioned in the
introduction. The precise definition requires more work and we refer to [11,
Sec. 2] for the details of the construction.

######  Theorem 3.1 (Curien & Kortchemski, [11]).

Let $(\tau_{n})_{n\geq 1}$ be a sequence of random trees such that there
exists a sequence $(b_{n})_{n\geq 0}$ of positive real numbers satisfying

  1. (i)

$\left(\frac{1}{b_{n}}W_{\left\lfloor\zeta(\tau_{n})t\right\rfloor}(\tau_{n});0\leq
t\leq
1\right)\xrightarrow[n\rightarrow\infty]{\enskip(d)\enskip}X^{exc,\alpha}$,

  2. (ii)

$\frac{1}{b_{n}}\mathsf{Height}(\tau_{n})\xrightarrow[n\rightarrow\infty]{\enskip(\mathbb{P})\enskip}0$

where the first convergence holds in distribution for the Skorokhod topology
on $D([0,1],\mathbb{R})$ and the second convergence holds in probability. Then
the convergence

$\frac{1}{b_{n}}\cdot\textsf{Loop}(\tau_{n})\xrightarrow[n\rightarrow\infty]{\enskip(d)\enskip}\mathcal{L}_{\alpha}$

holds in distribution for the Gromov-Hausdorff topology.

######  Remark 3.2.

This result is applicable to a sequence $(\tau_{n})_{n\geq 1}$ of critical
$GW_{\mu}$ trees, whose offspring distribution is in the domain of attraction
of a stable law of index $\alpha\in(1,2)$, which are conditioned to have $n$
vertices. Specifically, (i) follows from [12, Proposition 4.3]. Secondly, (ii)
is due to the fact that $\frac{b_{n}}{n}\cdot\mathsf{Height}(\tau_{n})$
converges in distribution towards a positive real valued random variable as
$n\rightarrow\infty$ [12, Theorem 3.1] and that $b_{n}^{2}/n\to 0$ as
$n\to\infty$ by (2).

###  3.3. Application to scaling limit of random Halin maps

We are now in a position to give an application of the invariance principle
given in the previous subsection by showing that the stable looptrees appear
as Gromov-Hausdorff limit of a sequence of Halin maps satisfying the
conditions in Theorem 1.2.

The Gromov-Hausdorff topology. Let $(E_{1},d_{1})$,$(E_{2},d_{2})$ be two
compact metric spaces. The Gromov-Hausdorff distance between $E_{1}$ and
$E_{2}$ is

$\mathrm{d}_{\mathrm{GH}}(E_{1},E_{2})=\inf_{\begin{subarray}{c}(E,d)\,\mathrm{metric}\,\mathrm{space}\\\
\varphi_{1}:E_{1}\rightarrow E\\\ \varphi_{2}:E_{2}\rightarrow
E\end{subarray}}\mathrm{d}_{\mathrm{Hau}}^{E}(\varphi_{1}(E_{1}),\varphi_{2}(E_{2}))$

where $\varphi_{1}$ and $\varphi_{2}$ are isometries. Recall that the
Hausdorff distance between two compact subsets, $K$ and $K\textprime$ of the
metric space $E$ is

$\mathrm{d}_{\mathrm{Hau}}^{E}(K,K\textprime)=\inf\\{\varepsilon>0,\,K\subseteq
K\textprime_{\varepsilon}\,\mathrm{and}\,K\textprime\subseteq
K_{\varepsilon}\\}$

where $K_{\varepsilon}$ denotes the $\varepsilon$-neighborhood of $K$. An
equivalent definition of the Gromov-Hausdorff distance involves
correspondences. A correspondence between two metric spaces, $(E_{1},d_{1})$
and $(E_{2},d_{2})$ is a subset $\mathcal{R}$ of $E_{1}\times E_{2}$ such that
for every $x_{1}\in E_{1}$, there exists at least one point $x_{2}\in E_{2}$
such that $(x_{1},x_{2})\in\mathcal{R}$ and conversely for every $y_{2}\in
E_{2}$, there exists at least one point $y_{1}\in E_{1}$ such that
$(y_{1},y_{2})\in\mathcal{R}$. The distortion of the correspondence
$\mathcal{R}$ is defined by

$\mathrm{dis}(\mathcal{R})=\sup\\{|d_{1}(x_{1},y_{1})-d_{2}(x_{2},y_{2})|:(x_{1},x_{2}),(y_{1},y_{2})\in\mathcal{R}\\}.$

The Gromov-Hausdorff distance can be written in term of correspondences by the
formula,

$\mathrm{d}_{\mathrm{GH}}(E_{1},E_{2})=\frac{1}{2}\inf\\{\mathrm{dis}(\mathcal{R})\\}.$

where the infimum is over all correspondences $\mathcal{R}$ between $E_{1}$
and $E_{2}$. We refer the reader to [8] for detail about this metric.

The following lemma is deterministic and shows that the Gromov-Hausdorff
distance between a Halin map satisfying $(H_{\star})$ and the looptree of the
shape of its associated marked tree is bounded by a linear function of the the
height of the tree.

######  Lemma 3.3.

Let $H\in\mathbb{H}_{n}$ be a Halin map satisfying $(H_{\star})$ and let
$T=\mathcal{S}(\phi_{n}(H))$ be the shape of its corresponding marked tree.
Then

$d_{GH}(H,\mathrm{Loop}(T))\leq\mathsf{Height}(T)+3/2.$ (4)

###### Proof.

We denote the extension of $\phi_{n}$ to the set of all finite Halin maps by
$\phi$. Let $H$ be a Halin map and define $T^{\ast}=\phi(H)$,
$T=\mathcal{S}(T)$ and $L=\mathsf{Loop}(T)$. To simplify, we divide the
argument into three steps by introducing spaces $\hat{H}$ and $\hat{L}$ which
are close to $H$ and $L$ respectively, in the Gromov-Hausdorff metric. Start
by contracting every non-boundary edge in $H$ which is connected to a leaf in
the underlying tree $H^{\circ}$ and denote the resulting map by $\hat{H}$. The
vertex set of $\hat{H}$ may be viewed as the the vertex set of $H$ apart from
the leaves of $H^{\circ}$. With this in mind, we construct a larger space $E$
from a copy of $H$ and a copy of $\hat{H}$, and connect copies of the same
vertices by edges of length 1. The natural embeddings of $H$ and $\hat{H}$ are
then isometries and their Hausdorff distance is equal to 2 (see Fig. 5 for an
explanation). It is then clear that,

$\displaystyle d_{GH}(H,\hat{H})\leq 2.$ (5)

Figure 5. The isometric embeddings of $H$ and $\hat{H}$ into a common metric
space $E$. A geodesic $\gamma_{H}$ in $H$ is depicted in orange. A priori it
might be possible to find a shorter path $\gamma_{E}$ in $E$ between the
corresponding vertices by making a short-cut into $\hat{H}$. However, the
edges connecting the vertices in $H$ and $\hat{H}$ in $E$ provide exactly the
correct penalty so that one can just as well choose a geodesic which lies
entirely in $H$.

Next, modify the looptree $L$ in the following way. Let
$\delta=k_{\varnothing}(T)$ and denote the subtrees attached to the root of
$T$ by $T^{1},T^{2},\ldots,T^{\delta}$ and the corresponding looptrees by
$L^{1},L^{2},\ldots,L^{\delta}$. Let $m$ be the label of the marked corner
around the root in $T^{\ast}$. Define the map $\hat{L}$ by detaching
$L^{\delta}$ and then identifying its root with the root in $L$ and by
replacing $L_{k}$ by $L_{k-1}$ for each $k=m+2,m+3,\ldots,\delta$, see Fig. 6.
It is clear that

$\displaystyle d_{GH}(L,\hat{L})\leq 1/2.$ (6)

Figure 6. The map $\hat{L}$ is obtained by shifting the subtrees around the
root of $L$, following the mark $m$, as indicated in the figure.

Finally, we compare $\hat{H}$ and $\hat{L}$. The map $\hat{L}$ may be
constructed on the vertex set of $\hat{H}$, yielding a canonical
correspondence $\mathcal{R}$ between $\hat{H}$ and $\hat{L}$, i.e.
$(x,x^{\prime})\in\mathcal{R}$ if and only if $x=x^{\prime}$. Let $f$ be a
bounded face in $\hat{H}$. Let $v$ be the vertex in $T$ dual to $f$. Then
$\mathsf{deg}(f)-2=k_{v}$. Denote the vertices on the boundary of $f$ by
$u_{1},u_{2},\ldots,u_{\mathsf{deg}(f)}$ in clockwise order, where $u_{1}$
denotes the vertex closest to the root in the underlying tree $H^{\circ}$.
Construct for each bounded face $f$, the circle
$(u_{1},u_{2},u_{3},\ldots,u_{\mathsf{deg}(f)},u_{1})$ if $f$ is the root
face, and the circle $(u_{2},u_{3},\ldots,u_{\mathsf{deg}(f)},u_{2})$
otherwise, which in both cases is of length $\mathsf{deg}(f)-1=k_{v}+1$. The
resulting map is equivalent to the looptree $\hat{L}$, see Fig. 7.

Figure 7. The natural correspondence between $\hat{H}$ and $\hat{L}$.

We will show that

$\displaystyle\mathrm{dis}(\mathcal{R})\leq 2\mathsf{Height}(T)$

which implies that

$\displaystyle d_{GH}(\hat{H},\hat{L})\leq\mathsf{Height}(T).$ (7)

We show this by induction over the height of $T$. It is easy to see that the
result holds when the height is 0 since both $\hat{H}$ and $\hat{L}$ are then
single vertices. Assume the bound holds for all $T$ with
$\mathsf{Height}(T)\leq h$. Let $H_{h+1}$ be a Halin map so that
$T^{\ast}_{h+1}=\phi(H_{h+1})$ is of height $h+1$. Define
$L_{h+1}=\mathsf{Loop}(\mathcal{S}(T^{\ast}_{h+1}))$. Now, remove all leaves
in $T^{\ast}_{h+1}$ at height $h+1$ and denote the resulting marked tree of
height $h$ by $T^{\ast}_{h}$. Let $H_{h}=\phi^{-1}(T^{\ast}_{h})$, and
$L_{h}=\mathsf{Loop}(\mathcal{S}(T^{\ast}_{h}))$.

The vertex set of $\hat{H}_{h}$ may be naturally viewed as a subset of the
vertex set of $\hat{H}_{h+1}$ and this relationship may be seen to be an
isometric embedding of $\hat{H}_{h}$ into $\hat{H}_{h+1}$ (adding leaves to
$T_{h}$ does not create shortcuts in $\hat{H}_{h}$). An analogous statement
applies to $\hat{L}_{h+1}$ and $\hat{L}_{h}$. Denote the metrics on
$\hat{H}_{h+1}$ and $\hat{L}_{h+1}$ by $d_{H}$ and $d_{L}$ respectively. Let
$x$ and $y$ be vertices in $\hat{H}_{h+1}$ (and thus also in $\hat{L}_{h+1}$)
and let $\gamma_{H}(x,y)$ be a geodesic between $x$ and $y$ in $\hat{H}_{h+1}$
and $\gamma_{L}(x,y)$ be a geodesic between $x$ and $y$ in $\hat{L}_{h+1}$.

Figure 8. The grey shaded area represents $\hat{H}_{h}$ and the three peaks at
the top represent the part of $\hat{H}_{h+1}$ which is in addition to
$\hat{H}_{h}$. The geodesics between $x$ and $y$ do not visit the vertices
strictly above the peak in the middle.

If no vertices on $\gamma_{H}(x,y)$ and $\gamma_{L}(x,y)$ lie in
$\hat{H}_{h}$, then

$\displaystyle|d_{H}(x,y)-d_{L}(x,y)|\leq 1<2(h+1).$

Otherwise, define $x_{H}$ (respectively $y_{H}$) as the vertex on
$\gamma_{H}(x,y)$ which is closest to $x$ (respectively $y$) in $\hat{H}_{h}$
and define $x_{L}$ (respectively $y_{L}$) as the vertex on $\gamma_{L}(x,y)$
which is closest to $x$ (respectively $y$) in $\hat{H}_{h}$ (it is possible
that $x=x_{H}=x_{L}$ or $y=y_{H}=y_{L}$). Then
$d_{H}(x,x_{H})=d_{L}(x,x_{L})$, $d_{H}(y,y_{H})=d_{H}(y,y_{L})$,
$d_{H}(x_{H},x_{L})\leq 1$ and $d_{H}(y_{H},y_{L})\leq 1$, see Fig. 8. Thus,

$\displaystyle|d_{H}(x,y)-d_{L}(x,y)|$
$\displaystyle=|d_{H}(x_{H},y_{H})-d_{L}(x_{L},y_{L})|$
$\displaystyle\leq|d_{H}(x_{H},y_{H})-d_{H}(x_{L},y_{L})|+|d_{H}(x_{L},y_{L})-d_{L}(x_{L},y_{L})|$
$\displaystyle\leq 2+2h=2(h+1)$

where the last inequality follows from the induction hypothesis and the fact
that $\hat{H}_{h}$ and $\hat{L}_{h}$ are isometrically embedded in
$\hat{H}_{h+1}$ and $\hat{L}_{h+1}$.

The bounds (5), (6) and (7) complete the proof of (4).

∎

###### Proof of Theorem 1.2.

Let $\mu$ be a probability measure on $\mathbb{Z}_{+}$ that satisfies the
assumptions of Theorem 1.2. The random tree
$\mathcal{T}_{n}^{w}=\mathcal{S}(\phi_{n}(\mathcal{H}_{n}^{w}))$ is, by
Proposition 2.4, a $GW_{\mu}$ tree conditioned to have $n$ vertices.
Therefore, by Theorem 3.1 and Remark 3.2,
$\frac{1}{b_{n}}\cdot\textsf{Loop}(\mathcal{T}_{n}^{w})$ converges torwards
$\mathcal{L}_{\alpha}$ in the Gromov-Hausdorff sense. By Lemma 3.3 it holds
that

$\displaystyle\frac{1}{b_{n}}\cdot
d_{GH}(\mathcal{H}_{n}^{w},\text{Loop}(\mathcal{T}_{n}^{w}))\xrightarrow[n\rightarrow\infty]{\enskip(\mathbb{P})\enskip}0$

since

$\displaystyle\frac{1}{b_{n}}\cdot\mathsf{Height}(\mathcal{T}_{n}^{w})\xrightarrow[n\rightarrow\infty]{\enskip(\mathbb{P})\enskip}0$

as is explained in Remark 3.2. This implies the statement of Theorem 1.2 ∎

### Acknowledgement

We thank Thomas Vallier and Thomas Selig for discussions in an early stage of
this work. The research is partially supported by the Icelandic Research Fund,
grant number 185233-051.

## References

  * [1] L. Addario-Berry, N. Broutin and C. Goldschmidt. The continuum limit of critical random graphs. Probab. Theory Related Fields , 152 (2012), no. 3-4, 367–406. MR2892951
  * [2] D. Aldous. The continuum random tree. I, Ann. Probab., 19(1) :1–28, 1991.
  * [3] D. Aldous. The continuum random tree. II. An overview, Stochastic analysis (Durham, 1990), volume 167 of London Math. Soc. Lecture Note Ser., pages 23–70. Cambridge Univ. Press, Cambridge, 1991.
  * [4] D. Aldous. The continuum random tree. III, Ann. Probab., 21(1) :248–289, 1993.
  * [5] J. Bertoin. Lévy Process, vol. 121 of Cambridge Tracts in Mathematics, Cambridge University Press, Cambridge, 1996.
  * [6] N. H. Bingham, C. M. Goldie, and J. L. Teugels. Regular variation, vol. 27 of Encyclopedia of Mathematics and its Applications, Cambridge University Press, Cambridge, 1989. MR1015093.
  * [7] J. Bouttier, P. Di Francesco, and E. Guitter. Planar maps as labeled mobiles, Electron. J.Combin., 11:no. 69, 27 pp. (electronic), 2004.
  * [8] D. Burago, Y. Burago and S. Ivanov.  A Course in Metric Geometry, AMS Graduate Studies in Math, vol. 33. Am. Math. Soc., Providence (2001).
  * [9] R. Cori and B. Vauquelin. Planar maps are well labeled trees, Canad. J. Math., 33(5):1023–1042, 1981.
  * [10] N. Curien, B. Haas, and I. Kortchemski. The CRT is the scaling limit of random dissections, Random Structures Algorithms, 47 (2015), pp. 304–327.
  * [11] N. Curien and I. Kortchemski. Random stable looptrees, Electron. J. Probab., 19 (2014), no. 108, 1-35.
  * [12] T. Duquesne. A limit theorem for the contour process of conditioned Galton-Watson trees, Ann. Probab., 31 (2003), pp. 996–1027. MR-1964956.
  * [13] T. Duquesne and J.-F. Le Gall. Random trees, Lévy processes and spatial branching processes, Astérisque,(2002), pp. vi+147. MR-1954248.
  * [14] I. Kortchemski. Invariance principles for Galton-Watson trees conditioned on the number leaves, Stochastic Process. Appl., 122 (2012), pp. 3126-3172.
  * [15] I. Kortchemski and C. Marzouk. Triangulating stable laminations, Electron. J. Probab. 21(11), 1–31(2016).
  * [16] J.-F. Le Gall. Random geometry on the sphere, In Proceedings of the ICM, 2014.
  * [17] J.-F. Le Gall.  Random trees and applications. Probab. Surv., 322(10):245–311, 1905.
  * [18] G. Miermont.  Aspects of random maps, In Saint-Flour Lecture Notes, 2014.
  * [19] G. Schaeffer. Conjugaison d’arbres et cartes combinatoires aléatoires, PhD thesis, Université Bordeaux I, 1998.
  * [20] S. Stefánsson and B. Stufler. Geometry of large Boltzmann outerplanar maps, Random Structures Algorithm, 55 (2019), no. 3, pp. 742–771. MR-3997486.


#  tachyon: Efficient Shared Memory Parallel Computation of Extremum Graphs

Abhijath Ande0000-0003-0575-9350, Varshini Subhash0000-0003-1889-6821, and
Vijay Natarajan0000-0002-7956-1470
Department of Computer Science and Automation, Indian Institute of Science,
Bangalore, India

###### Abstract

The extremum graph is a succinct representation of the Morse decomposition of
a scalar field. It has increasingly become a useful data structure that
supports topological feature directed visualization of 2D / 3D scalar fields,
and enables dimensionality reduction together with exploratory analysis of
high dimensional scalar fields. Current methods that employ the extremum graph
compute it either using a simple sequential algorithm for computing the Morse
decomposition or by computing the more detailed Morse-Smale complex. Both
approaches are typically limited to two and three dimensional scalar fields.
We describe a GPU-CPU hybrid parallel algorithm for computing the extremum
graph of scalar fields in all dimensions. The proposed shared memory algorithm
utilizes both fine grained parallelism and task parallelism to achieve
efficiency. An open source software library, tachyon, that implements the
algorithm exhibits superior performance and good scaling behavior.

<ccs2012> <concept> <concept_id>10003120.10003145.10003146</concept_id>
<concept_desc>Human-centered computing Visualization techniques</concept_desc>
<concept_significance>300</concept_significance> </concept> <concept>
<concept_id>10010147.10010169.10010170</concept_id> <concept_desc>Computing
methodologies Parallel algorithms</concept_desc>
<concept_significance>500</concept_significance> </concept> </ccs2012>

[300]Human-centered computing Visualization techniques [500]Computing
methodologies Parallel algorithms

##  1 Introduction

Topological analysis has become a key tool for the analysis and visualization
of data from various branches of scientific research [CFST20, HMST21]. In the
context of data represented as scalar fields, topological descriptors such as
contour trees, Reeb graphs, merge trees, Morse-Smale complexes, and extremum
graphs have been widely studied [HLH∗16]. Each abstraction provides a
different perspective to analyze and interact with the scalar field. The
extremum graph is a useful abstraction in situations where the extrema and
gradient flow behavior that determines their connectivity are crucial for
understanding the scientific phenomenon [CLB11, TN13]. It is a subset of the
Morse-Smale complex that is considerably simpler yet intuitive and detailed
enough for various applications. It is a bipartite graph whose nodes are
either extrema or saddles and arcs determine connectivity via a gradient path
between a saddle-extremum pair. The simple structure of the graph makes it
amenable for various applications including segmentation in 2D and 3D scalar
fields, feature tracking in time-varying data, and clustering in higher
dimensional data. While there exist efficient algorithms to compute the more
general Morse-Smale complex, there is a need for a lightweight, fast, and
scalable algorithm that can compute the extremum graph in all dimensions.

###  1.1 Related work

The Morse complex and Morse-Smale complex were introduced as topological
structures to represent the gradient flow behavior of Morse functions [BHEP04,
EHNP03]. While the extremum graph may be considered as a 1-skeleton of the
Morse complex, it was first described in the currently known form by Correa et
al. [CLB11] together with a new 2D visual representation called the
topological spine. The extremum graph preserves the topological and geometric
structure of a scalar field, while the topological spine provides a succinct
and easily digestible 2D representation to the user. It essentially preserves
the relative location of extrema and knowledge of their neighborhoods with
respect to gradient paths connecting them. Useful applications of extremum
graphs include the robust detection of symmetry in noisy scalar fields [TN13],
extrema proximity awareness and a feature-aware comparison of similar scalar
fields using distance measures [NTN15].

The abundance of large-scale scientific data with increasing feature
complexity and precision has directed the attention of the scientific
visualization community towards developing parallel algorithms that
effectively leverage modern compute power and massively parallel
architectures. Specifically, the parallel computation of 3D Morse-Smale
complexes has been studied extensively to this effect. Gyulassy et al.
[GBHP08] developed a memory efficient algorithm for 3D Morse-Smale complex
computation, where they partitioned large data into chunks called parcels that
fit in memory. By computing the Morse-Smale complexes for each parcel and
performing a cancellation based merging of parcels, they were able to obtain
Morse-Smale complexes for large-scale data which originally did not fit in
memory. This approach was adapted to a distributed memory setting by Peterka
et al. [PRG∗11] and Gyulassy et al. [GPPR12] where they leveraged massively
parallel clusters to handle the parcels in parallel. Robins et al. [RWS11] and
Shivashankar and Natarajan [SSN11] introduced novel locally independent
definitions for gradient pairs which facilitated an embarrassingly parallel
gradient pair assignment. Novel traversal algorithms for the extraction of
ascending and descending manifolds of the extrema and saddles further enhanced
the runtime performance of the Morse-Smale complex computation [GRWH12, SN12,
DFRS14].

Additionally, there have been advances in obtaining a higher degree of
geometric accuracy during parallel Morse-Smale complex computation by Gyulassy
et al. [GBP18, GGL∗14, GBP12] and Bhatia et al. [BGL∗18]. Several approaches
dabbled in CPU based shared memory parallelization techniques with the
exception of Shivashankar and Natarajan [SN12], who introduced a hybrid
approach. Their algorithm and the associated software library, pyms3d [SN17],
leveraged GPU parallelism for gradient assignment and extrema traversal and
CPU parallelism for saddle-saddle traversal. Recently, Subhash et al. [SPN20,
SPN22] introduced gmsc, an end-to-end GPU parallel algorithm for 3D Morse-
Smale complex computation, which resulted in substantial speedups over the CPU
and hybrid approaches attempted thus far. gmsc transforms the saddle-saddle
computation into a sequence of matrix operations that are amenable to fast
parallel computation on the GPU, leading to significant improvement in
runtime. We note that the gradient assignment and saddle-extrema traversal
steps in gmsc are identical to those employed by pyms3d. Both pyms3d and gmsc
are restricted to 3D scalar fields.

A discrete Morse theory based approach is employed by several of the above-
mentioned parallel algorithms and has yielded combinatorial and numerically
robust algorithms [GNP∗06, GRWH12, DFRS14, SSN11, SN12]. Fugacci et al.
[FIDF19] also demonstrated an extension to higher dimensions by computing the
Morse complex from simplicial complexes. This discrete Morse complex can be
leveraged to compute homology [RWS11, HMMN14], perform shape analysis and
study scalar fields [DFFIM15].

As mentioned earlier, the extremum graph is a subgraph of the combinatorial
structure or the 1-skeleton of the Morse-Smale complex. The succinct yet
highly informative abstraction and its demonstrated role in many applications
motivates the need for its efficient computation. To the best of our
knowledge, no previous work specifically explores the efficient and scalable
parallel computation of extremum graphs in all dimensions. We attempt to
address this gap here and describe an open source implementation.

###  1.2 Contributions

In this paper, we describe a fast shared memory parallel algorithm for
computing the extremum graph of an $n$-dimensional scalar field defined on a
grid. The hybrid GPU-CPU parallel computation leverages the GPU for fine-
grained parallel computation of nodes of the extremum graph and the CPU for
effective task parallel computation of arcs of the graph. Key contributions of
the paper include:

  * •

A hybrid GPU-CPU algorithm that efficiently uses both computational resources
while enabling the algorithm to be scalable to data that does not fit in the
GPU memory.

  * •

An efficient implicit representation of the neighborhood of a grid vertex that
supports fast parallel critical point classification of all vertices.

  * •

tachyon, a lightweight library that provides a fast implementation of the
hybrid algorithm (bitbucket.org/vgl_iisc/tachyon).

  * •

Multiple strategies for effective simplification of the extremum graph that
makes it amenable for feature identification and applications.

  * •

Experimental studies on multiple datasets to demonstrate superior performance
and scalability.

##  2 Background

In this section, we define and briefly introduce the necessary terms required
to define extremum graphs [EHNP03, CLB11, NTN15]. Given a smooth function
$f:\mathbb{M}\rightarrow\mathbb{R}$ over a manifold $\mathbb{M}$ of dimension
$n$, we say that a point $x\in\mathbb{M}$ is _critical_ iff $\nabla f(x)=0$,
$x$ is called _regular_ otherwise. Also, $f$ is a _Morse function_ if all
critical points have pairwise distinct function values and none of them are
degenerate i.e., the Hessian evaluated at the point is non-singular. For a
critical point $x$ of $f$, its Morse index is defined as the number of
negative eigenvalues of its Hessian matrix evaluated at $x$. Critical points
of index $0$ are named _minima_ , index $n$ are _maxima_ and the others with
index $k$, $1<k<n-1$, are called $k$-_saddles_.

An _integral line_ is a maximal curve in $\mathbb{M}$ whose tangent at every
point is equal to the gradient of $f$ at that point. Naturally, $f$
monotonically increases along the integral line and its two end points (limit
points) are critical points of $f$. The Morse function $f$ determines a
decomposition of $\mathbb{M}$ based on the integral lines. The union of
integral lines that terminate at a critical point is called its _descending
manifold_. The descending manifold of a maximum is an $n$-dimensional
manifold. The collection of descending manifolds of critical points of $f$
partition $\mathbb{M}$ and is called the _Morse decomposition_. The _ascending
manifold_ of a minimum is similarly defined as the union of integral lines
that originate at a minimum. Again, the collection of ascending manifolds of
critical points of $f$ partition $\mathbb{M}$. The _extremum graph_ is a
representation of the Morse decomposition. It is called a maximum / minimum
graph if it represents the decomposition into descending / ascending
manifolds, respectively. Without loss of generality, we restrict the
discussion in this paper to maximum graphs, and refer to them as extremum
graphs.

An $(n-1)$-saddle $s$ of $f$ lies on the boundary of the descending manifold
of a maximum $m$ i.e., an integral line originating at $s$ terminates at $m$.
The extremum graph captures this relationship between maxima and
$(n-1)$-saddles of $f$, and thereby captures the combinatorial structure of
the Morse decomposition. The node set of the extremum graph consists of the
maxima and $(n-1)$-saddles of $f$. An arc $(s,m)$ belongs to the extremum
graph if $s$ lies on the boundary of the descending manifold of $m$. Figure 1
shows the extremum graph for a 3D scalar field, the output of a simulation of
a silicium grid.

Figure 1: Extremum graph for the silicium grid. Maxima (red) and 2-saddles
(green) are connected by arcs. The geometry of an arc is determined by the
integral line (black) that connects the maximum and the 2-saddle. The scalar
field is shown as a volume rendering using a simple banded color map.

##  3 Computation of extremum graphs

There are two popular approaches to compute extremum graphs, namely flood fill
and gradient path tracing. While both methods produce equivalent results, each
method has its own associated advantages and disadvantages. We begin with a
brief discussion about the two approaches and justify why we propose gradient
path tracing towards the development of an efficient and scalable parallel
algorithm.

###  3.1 Flood fill vs. gradient path tracing

The flood fill approach computes the extremum graph in one sweep over the list
of grid vertices by incrementally growing the descending $n$-dimensional
manifolds of all maxima. The vertices are processed in decreasing order of
scalar value. Processing a vertex $v$ includes expanding the descending
manifolds that have reached $v$ in the steepest descent direction, namely the
vertex with the smallest value in the lower link. This approach has been
employed to compute the Morse decomposition [EHNP03]. The method implicitly
recognizes a maximum as a vertex where no descending manifold have reached and
an $(n-1)$-saddle as a vertex where multiple descending manifolds merge. The
adjacency between maxima and saddles is recorded when the saddle vertex is
processed. While this method computes the connectivity of the extremum graph,
it does not explicitly capture the gradient paths between the maxima and
saddles. A subsequent step needs to be invoked to compute these gradient paths
as necessary, a limitation of this approach. Further, this method has a high
memory footprint because it requires the collection of all descending
manifolds to be stored.

The gradient path tracing approach computes the extremum graph in two steps.
The first step locates and classifies the maxima and $(n-1)$-saddles. It
analyzes the link of a vertex to classify the vertex. Importantly, the
classification of a vertex is independent of the classification of other
vertices and depends on a local neighborhood. The second step traces gradient
paths from each $(n-1)$-saddle towards extrema. This tracing requires
gradients to be computed only for a small subset of regular vertices, thereby
resulting in a low computational load. The memory footprint is also low
because the first step requires a simple query within a constant sized local
neighborhood and the second step requires only the end point of the path to be
stored.

There exist multiple parallel algorithms for flood fill. However, these method
do not scale due to the communication required to handle the merge events and
subsequent serialization of the descending manifold growth computation. In
contrast, the first step of the gradient path tracing is amenable to fine
grained parallelism, is simple, and highly scalable. The second step can also
be accelerated via task parallelism. These advantages motivate us to develop a
gradient path tracing algorithm for parallel computation of the extremum
graph. In the following, we describe the two steps of the algorithm with a
focus on how the steps are parallelized, followed by a discussion on effective
simplification of the extremum graph to make it suitable for applications.

###  3.2 Grid tessellation

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

(a) (b)

(c) (d)

(e) (f)

Figure 2: Tessellating 2D and 3D grids. (a) A 2D cell is decomposed into two
triangles. (b) Tessellating a 2D grid. (d) A 3D grid cell is decomposed into
six tetrahedra. (e) Tessellating a 3D grid. The decomposition of the common
faces between two adjacent grid cells is consistent. (c,f) Edges of the
neighborhood graph (blue) of the vertex at the centre. The vertex has 6
neighbors in a 2D grid and 14 neighbors in a 3D grid.

We assume that the scalar field is available as a collection of samples at
vertices of an $n$-dimensional grid. Vertices of the grid that represent the
domain $\mathcal{D}^{n}$ have integral coordinates,

$V(\mathcal{D}^{n})=\\{v\mid
v\in\underbrace{\mathbb{Z}\times\mathbb{Z}\times\cdots\times\mathbb{Z}}_{n}\\}.$

Edges of the grid are between two vertices that differ in exactly one
coordinate by a value of 1. Methods for topological analysis and visualization
of this scalar field typically assume an input piecewise multilinear or
piecewise linear scalar function. The samples at vertices may be extended into
a scalar function by either using a multilinear interpolant (trilinear in 3D)
within each grid cell or by tessellating the grid cells into linear cells
(tetrahedra in 3D) followed by linear interpolation within each linear cell.
In either case, the local maxima of the scalar field are located at grid
vertices.

We tessellate the uniform grid into irregular linear cells (tetrahedra in 3D)
with the aim of simplifying the extremum graph computation and to support
future extensions of the method to unstructured grids. We choose a
tessellation that decomposes all grid cells in the same manner and one that is
consistent on the common face between adjacent cells, called the Freudenthal
subdivision [CMS06]. Figure 2 shows the tessellated grid in 2D and 3D (for
illustration). The resulting edge set of the tessellated grid
$E(\mathcal{D}^{n})$ can be stored implicitly and recovered using a simple
routine, see Algorithm 1. The algorithm essentially dictates that two distinct
vertices are connected by an edge if and only if the difference vector between
the two points consists entirely of non-negative or non-positive values and
the magnitude of the non-zero values is 1.

The edge set of the tessellated grid

$E(\mathcal{D}^{n})=\\{(u,v)\mid u,v\in
V(\mathcal{D}^{n})\wedge\textnormal{{GridAdjacency}}(u,v)\\}.$

Figure 2 also shows the edges incident on a given vertex of the tessellated
grid. It follows from the symmetry of the edge set that the number of edges
incident on a vertex $v$ is equal to twice the number of non-zero difference
vectors, namely $2\times(2^{n}-1)$ for nD grids. So, the number of vertices in
the neighborhood of $v$ in a 3D grid equals 14. Critical points are identified
and classified based on a connected component labeling of this neighborhood.
We discuss this classification next.

Input: $p,q:$ Two vertices. Test if they are adjacent in the tessellated
$n$-dimensional grid.

Result: $True$ if $p,q$ are adjacent and $False$ otherwise

1

2$n\leftarrow len(p)$ $\triangleright$ $p$ is a vertex in an $n$-dimensional
grid

3 $U\leftarrow\phi$

4foreach  _$i\in 1..n$_ do

5 $d\leftarrow p_{i}-q_{i}$ $\triangleright$ difference between $i^{th}$
coordinates

6 $U\leftarrow U\cup\\{d\\}$ $\triangleright$ collect component wise
difference

7

8 end foreach

9if  _$U{\subseteq}\\{0,1\\}$ or $U{\subseteq}\\{0,-1\\}$_ then

10 return _$True$_

11 end if

12else

13 return _$False$_

14 end if

15

Algorithm 1 GridAdjacency

###  3.3 Critical point classification

Critical points of the piecewise linear function defined over a tessellated
grid can be identified and classified based on a local neighborhood [Ban70,
EHNP03]. All critical points are located at vertices of the grid. The _star_
of a vertex $v$ consists of $v$ together with the collection of edges,
triangles, tetrahedra, and higher dimensional cells incident on it. Figure
2(f) shows the edges (blue) of the star of a vertex in a cube grid. The
collection of end points of these edges, excluding $v$, together with the
induced edges and triangles form the _link_ of $v$. The star and link are two
useful notions of neighborhood of a vertex within the tessellated grid.
Assuming that there are no degeneracies, scalar values at vertices in the link
are either lower or higher than at $v$. Vertices of the link with scalar value
lower than at $v$ together with the induced edges and triangles form the
_lower link_ of $v$. Similarly, the _upper link_ of $v$ is defined as the
collection of vertices of the link with scalar values greater than at $v$
together with the induced edges and triangles.

Critical points are classified based on the number of connected components of
the upper and lower link, denoted $\beta_{0}^{+}$, $\beta_{0}^{-}$ and called
the zeroth Betti number [EHNP03]. Figure 3 illustrates the upper (red) and
lower (green) link for various types of critical points and for a regular
point. Since our objective is to compute the extremum graph, we are
specifically interested in identifying the maxima and 2-saddles. The first
step of the algorithm visits each grid vertex and classifies it as critical or
regular as shown in Table 1.

| $\beta_{0}^{+}$ | $\beta_{0}^{-}$
---|---|---
maximum | 0 | 1
$(n-1)$-saddle |  $\geq$2 | $\bullet$
1-saddle | $\bullet$ |  $\geq$2
minimum | 1 | 0
regular | 1 | 1
Table 1: Classifying a vertex based on the topology of the upper and lower
link. $\beta_{0}^{+}$ and $\beta_{0}^{-}$ count the number of connected
components of the upper and lower link, respectively. $\bullet$ indicates that
the value does not affect the classification.

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

(a) Maximum (b) 2-Saddle

(c) Regular point (d) Minimum

(e) 1-Saddle (f) Multi-saddle

Figure 3: Critical points can be identified and classified based on the upper
(red) and lower (green) link connectivity. The upper link of a maximum (a) is
empty and the lower link is equivalent to a sphere, whereas the upper link of
a minimum (d) is equivalent to a sphere and the lower link is empty. The upper
link of a simple 2-saddle (b) consists of two connected components and its
lower link consists of a single component, whereas the upper link of a
1-saddle (e) consists of a single component and the lower link consists of two
connected components. Both the upper and lower links of a regular point (c)
consist of a single component. The multi-saddle (f) is a degenerate structure
that does exist in piecewise linear functions.

Both $\beta_{0}^{+}$ and $\beta_{0}^{-}$ are typically computed by performing
a BFS graph traversal on the upper and lower link of the vertex $v$. The BFS
traversal of the upper and lower link is an embarrassingly parallel task
across all vertices. We utilize GPU parallelism to perform this traversal in a
massively parallel fashion by launching a CUDA thread for each vertex. Our
implementation moves this task automatically to a multicore CPU in the absence
of a CUDA enabled GPU.

In order to make the computation GPU friendly, we utilize a union-find based
connected component tracking instead of the BFS traversal, similar to previous
approaches [CWS∗19]. The union-find data structure is initialized with a
collection of singleton sets, each containing a vertex of the link of $v$.
Next, we test for the existence of an edge between a pair of vertices in the
link using Algorithm 1. If the edge exists and both end points of the edge
belong to the upper link (lower link), we perform a union operation on the
corresponding sets. Degeneracies are handled using a simulated perturbation
[EH09, Section 1.4], which consistently determines if the scalar value at a
link vertex is lower or higher than the value at $v$.

After identifying the connected components in the link, we label each
component as upper or lower link by testing one vertex within the component
and hence compute $\beta_{0}^{+}$ and $\beta_{0}^{-}$. This method is GPU
friendly as all working threads test the same number of edges and require
almost the same time, thereby causing a very low thread divergence. If a
vertex is identified as regular, we additionally compute the gradient of the
scalar field at the vertex. The gradient is approximated as the vector towards
the vertex with the highest scalar value in the upper link.

###  3.4 Path tracing

We opt to use a direct path tracing algorithm to trace the saddle-maximum arcs
that constitute the extremum graph over the flood-fill approach. The primary
reason for this choice is the benefit in terms of scalability of the
computation. In addition, path tracing does not require each vertex of the
grid be visited. Further, a flood-fill based algorithm requires an ordered
list of vertices, sorted on the scalar values.

Input: $s,M$: Source $(n-1)$ saddle and the set of all maxima

Result: Set of gradient paths originating from $s$

1 $P\leftarrow\phi$

2foreach  _$u\in\textit{UpperLinkRep}(s)$_ do

3

4 $p\leftarrow s$ $\triangleright$ Initialize path $p$ with $s$

5

6 while  _$u\not\in M$_ do

7 $p\leftarrow p\,\,||\,\,u$ $\triangleright$ Append $u$ to $p$

8 $u\leftarrow\textit{gradient}(u)$ $\triangleright$ Follow gradient at $u$

9

10 end while

11

12 $P\leftarrow P\,\,\cup\,\,\\{p\\}$

13 end foreach

14

15return _$P$_

Algorithm 2 TraceGradientPaths

The path tracing algorithm requires only the gradient at vertices. The
gradient paths in an extremum graph originate from $(n-1)$-saddles, so the
algorithm visits only a small fraction of vertices of the grid. However, this
approach is not without challenges. The operation of tracing a gradient path
is inherently serial in nature. One approach towards parallelism is to trace
the gradient paths for all $(n-1)$-saddles concurrently.

Algorithm 2 computes the collection of gradient paths that originate at an
$(n-1)$-saddle. The subroutine UpperLinkRep returns the highest vertex in each
upper link component of the saddle. For each such vertex, it iteratively
follows the gradient vectors until termination at a maximum. This algorithm is
executed in parallel for each $(n-1)$-saddle. Unlike the critical point
classification, Algorithm 2 is not amenable to efficient GPU parallelism.
Varying lengths of the gradient paths causes unequal division of work among
CUDA threads, and results in significant thread divergence.

Figure 4: Gradient path tracing for computing arcs of extremum graph.
Ascending 1-dimensional manifolds are computed for each 2-saddle. The gradient
paths that constitute these ascending manifolds represent the arcs of the
extremum graph. Note that they pass through regular points, entering and
exiting its neighborhood via the lower and upper link, respectively.

Figure 4 illustrates the ascending 1-dimensional manifolds for a 2-saddle.
Notice how the gradient paths originate from the 2-saddle via each upper link
component and terminate at maxima. Given two adjacent regular points (orange)
$(p,q)$ on a gradient path, $p$ lies in the lower link of $q$ and $q$ lies in
the upper link of $p$. The neighborhood of one regular point is shown in
detail. The gradient path enters through its lower link and exits via its
upper link.

###  3.5 Extremum graph simplification

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

(a) Fuel (b) Lobster

(c) Foot (d) Boston Teapot

(e) Nucleon (f) Neghip

(g) Hydrogen (h) Engine

(i) Heptane Gas

Figure 5: Simplified extremum graphs for various datasets. Nodes (extrema and
saddles) are not displayed to reduce clutter. The scalar field is mapped to
color as indicated by the legend.

Scientific data is often noisy, which manifests as an increased number of
critical points and hence a larger sized extremum graph. Topological
simplification identifies noisy topological features and removes them in a
controlled manner, often as a sequence of critical point pair cancellation
operations. The simplification is typically directed by the notion of
persistence [ELZ02]. Indeed, the extremum graph computed using the algorithm
that we describe in the previous sections contains a large number of saddles
and extrema for all datasets. Several critical points and their incident
gradient paths correspond to noisy topological features. The presence of these
noisy elements cause occlusion, thereby adversely affecting data visualization
tasks, and hinder feature selection. We propose three simplification
operations to aid the proper identification of key features from the extremum
graph. The persistence directed cancellation and saturated persistence
directed simplification algorithms are implemented in tachyon as they are
useful in removing noise and potentially uninteresting features from the
extremum graph [CLB11]. We introduce the arc-bundling algorithm as a method to
remove multiple connections between a pair of extrema. Any sequence of these
operations may be applied on an extremum graph. We describe the simplification
operations below. Their implementations in tachyon are serial in nature.

Arc bundling. The gradient path tracing algorithm computes individual arcs
between an $(n-1)$-saddle-maximum pair independent of each other. So, it is
likely that a pair of maxima contain more than one common $(n-1)$-saddle in
their neighborhood. We select a single representative saddle based on a
specified criterion and discard the remaining saddles that are shared by the
two maxima within their neighborhood. We term this pruning operation as arc
bundling, as this reduces the number of arcs between an extrema pair. Note
that this operation does not disturb the connectivity between extrema and the
resulting graph continues to represent the Morse decomposition. In our current
implementation, we select the saddle with the highest scalar value as the
representative.

Persistence directed cancellation. This operation cancels an $(n-1)$-saddle-
maximum pair that is connected by an arc in the extremum graph and reconnects
the neighborhood of both nodes. The critical point pairs are ordered based on
their difference in scalar value, motivated by the notion of topological
persistence, and scheduled for cancellation. The simplification begins by
initializing a priority queue with all $(n-1)$-saddles, where the priority is
inversely related to the persistence of the saddle. Computing the topological
persistence of all saddles upfront is unnecessary and computationally
expensive. We use lazy updates where the true persistence of a saddle is
computed when it reaches the top of the queue and is ready for cancellation.
We initialize the cost of canceling a simple $(n-1)$-saddle as the smallest
difference in scalar value between the saddle and its two adjacent maxima in
the extrema graph. This cost is recomputed whenever a saddle is removed from
the top of the queue.

The next step of the simplification iteratively removes the saddle $s$ from
the top of the queue, recomputes its cost, and decides whether to discard,
reinsert, or cancel the saddle. If the updated cost is above a user-specified
persistence threshold, the saddle should not be cancelled. Note that
subsequent cancellations will not decrease the cost of this saddle and so, it
is safely discarded. If the updated cost is greater than the cost of the top
of the queue (but smaller than the persistence threshold) then the saddle is
reinserted into the queue because the queue contains other saddles with lower
costs that need to be canceled. Else, the saddle should be canceled. The
gradient path from $s$ towards its persistence pair $m$ is reversed to cancel
the critical point pair. The above step is repeated until the queue is empty.
This simplification removes all low persistence critical point pairs (arcs)
from the extremum graph, while retaining features of potential interest.

Multi-saddles have greater than two arcs incident on them. The cost of a
multi-saddle is computed as the difference in scalar value with the second
highest maximum adjacent to it. The multi-saddle is canceled by reversing
gradient paths that originate from it towards all maxima, except for the
highest maximum, also called the surviving maximum.

Saturated persistence directed simplification. While persistence directed
cancellation removes saddles whose persistence falls below a given threshold,
it does not remove long and potentially unintuitive arcs that connect
spatially distant maxima via a saddle. These additional arcs often do not
represent structural features. In addition, they cause clutter and occlude
potentially interesting features. We propose the use of saturated persistence
directed extremum graph simplification [CLB11], which prunes such arcs.

We apply all three operations described above to obtain simplified extremum
graphs for various scientific datasets as shown in Figure 5.

##  4 Hybrid GPU-CPU parallel computation

The parallel algorithm described in the previous section works when the
dataset together with the associated data structures fit within the GPU and
main memory. The GPU memory is often much smaller and determines the size of
the data that can be processed. We now describe an extension of both steps of
the gradient path tracing algorithm to handle datasets that are too large to
fit within the GPU memory. We do assume, however, that the data fits within
the main memory. In order to classify critical points in large datasets, we
view the dataset as an array of scalar field values in row major order. Thus,
each vertex of the grid is mapped to an index of this array. We partition the
array into contiguous subarrays called _blocks_. These contiguous blocks
partition the $n$-dimensional dataset along one of its dimensions, for example
along the y-axis for 2D data and along the z-axis for 3D data. A block is
typically adjacent to two blocks, except for the first and last block which
are adjacent to exactly one block. The domain is partitioned such that each
block fits in the available GPU memory during runtime. The size of a block is
computed based on a linear regression model that predicts the number of
vertices that can be accommodated in the GPU. The linear regression model is
constructed from a series of experiments on datasets of different sizes
executed on different GPUs.

Critical point classification of vertices that lie on the boundary of a block
requires access to their neighboring vertices. So, we append each block with a
set of ghost vertices constituted by vertices from the link of the boundary.
These ghost vertices are not classified into regular or critical while
processing the block. Critical point classification of all vertices in a block
is computed on the GPU and the results are transferred to main memory. This
process is repeated for all blocks in sequence. If multiple GPUs are
available, then multiple blocks could be processed concurrently. tachyon
currently does not utilize multiple GPUs.

---
(a) Sequential, GPU followed by CPU execution. (b) Hybrid GPU-CPU mode.

Figure 6: Hybrid GPU-CPU parallel computation for a dataset partitioned into 4
blocks. In the sequential execution model (a), point classification step
(green) processes the blocks one by one on the GPU. Gradient path tracing
(red) for the blocks executes on the CPU after point classification is
complete for all blocks. The CPU cores lie idle, waiting for the completion of
critical point classification for all blocks on the GPU. In the hybrid GPU-CPU
mode (b), gradient path tracing within the first block begins soon after it is
processed by the GPU. The transfer of control from GPU to CPU for a block is
indicated by arrows. The base of the arrow lies at the finishing time for a
block on GPU and the tip lies at the time when that block is taken up for
processing on the CPU.

We could employ a simple approach for parallel gradient path tracing after all
vertices are processed in the first step. However, this leads to a waste of
CPU resources while the GPU is executing the first step. An optimal approach
utilizes the CPU simultaneously for gradient path tracing within blocks that
are already processed in the GPU. However, this path tracing is non-trivial.
The gradient path tracing is indeed executed within each block but the task
needs to be temporarily paused when the path reaches the block boundary and
attempts to exit the block. Figure 6 illustrates how such an approach reduces
the idle CPU time.

Gradients paths that need to exit a block are recorded and their tracing is
resumed once the adjacent block is processed by the GPU. Path tracing in this
hybrid mode is performed by maintaining a list of partially traced paths,
whose tracing could not be completed because it entered a block that is not
yet processed by the GPU to classify the critical points. A partial path is
stored as a triple – the origin saddle, the first vertex on the path, and the
last visited vertex. The first vertex is required to distinguish between all
paths that originate from the saddle. For each partial path, the last visited
vertex is updated when path tracing is paused at a block boundary. Path
tracing is resumed when the adjacent block containing the vertex is processed
by the GPU. The tracing of a specific path continues until a maximum is
reached and the overall tracing process terminates when the list of partial
paths is exhausted. Thus, each vertex on a saddle-maximum path is visited
once, which implies that the worst case running time is determined by the
total number of vertices across all traced paths. The space required for path
tracing is equal to the maximum size of the partial paths list.

The above method can be further optimized by noticing that at any instant the
collection of blocks processed by the GPU thus far can be viewed as a single
but larger contiguous block. This is true if the blocks are processed in
sequence. We fuse the blocks, iteratively including the next block after point
classification is complete. This fuse step reduces the number of boundary
crossings that may occur during gradient path tracing by simply reducing the
number of possible boundaries to one, thus allowing for a majority of threads
to complete tracing. By nature of this scheme, the extremum graph is computed
without the need of any further processing, once gradient path tracing is
completed for the last block.

##  5 The tachyon library

tachyon (bitbucket.org/vgl_iisc/tachyon) is a C++ software library that
enables efficient computation of the extremum graph by offloading major
portions of the computation onto an Nvidia GPU that supports CUDA. We list
below a few salient features of the software library:

  * •

Leverages GPU computation power via CUDA kernels to perform point
classification, an embarrassingly parallel task.

  * •

Written in C++14. Utilizes standard threading libraries and synchronization
utilities to implement parallelism for most tasks that cannot be performed on
GPUs, thereby ensuring ease-of-installation and portability.

  * •

Generalized to work for $n$-dimensional scalar fields. Computes both maximum
and minimum graph, as required.

  * •

Automatically switches to a hybrid GPU-CPU mode, on the fly, for datasets that
do not fit in GPU memory.

  * •

Provides graph simplification and cancellation algorithms for downstream
applications of the extremum graph.

  * •

Utilizes code and memory optimization tricks, such as bit manipulation
operations and caching, for effective use of memory.

  * •

Provides a user friendly command line utility to specify input data format,
choice of graph simplification, cancellation algorithms, etc.

The library is designed in a modular fashion. For example, it is possible to
configure tachyon to perform point classification on the CPU if NVidia GPUs
are not available. On a multi-GPU setup, it allows a user to select which GPU
to utilize for computations. Further, it can be extended to support multi-GPU
computation. The library has only a small set of dependencies, making it a
versatile tool that can be built and used on a wide variety of workstations.

tachyon is bundled with a user friendly interface that enables a non-expert to
use the library for computing the extremum graph. The user interface supports
various desired functionalities such as computation of maximum and/or minimum
graphs based on user requirement, graph cancellation and simplification
operations, support for various data types for storing the scalar field
($8/16/32/64$-bit signed/unsigned integers and single/double precision
floating point values).

##  6 Experimental Results

Experiments were conducted on datasets of sizes varying from $64\times
64\times 64$ to $2048\times 2048\times 2612$. These datasets were obtained
from the Open Scientific Visualization Dataset Repository [Ope22]. All
experiments, unless stated otherwise, were run on a workstation with an
Intel(R) Xeon(R) Gold 6258R CPU @ 2.70GHz 28 Cores/56 Threads, 640GB RAM, and
an Nvidia GeForce RTX 3080 GPU with 10GB GDDR6X RAM. Runtimes were computed as
an average over five samples after dropping the best and worst timings out of
seven runs. GPU temperatures were also kept within optimal operating ranges of
$<85^{o}$ C. We have performed experiments to highlight different properties
of tachyon and to evaluate its performance. Table 2 lists the various datasets
used in the following experiments.

Dataset | Size | #Critical Pts
---|---|---
Nucleon | 41$\times$41$\times$41 | 421
Silicium | 98$\times$34$\times$34 | 548
Neghip | 64$\times$64$\times$64 | 1646
Fuel | 64$\times$64$\times$64 | 296
Hydrogen | 128$\times$128$\times$128 | 11370
Shockwave | 64$\times$64$\times$512 | 991
Lobster | 301$\times$324$\times$56 | 347207
Head Mri Ventricles | 256$\times$256$\times$124 | 1719145
Engine | 256$\times$256$\times$128 | 501799
Statue Leg | 341$\times$341$\times$93 | 454011
Boston Teapot | 256$\times$256$\times$178 | 101967
Skull | 256$\times$256$\times$256 | 1877366
Foot | 256$\times$256$\times$256 | 788482
Aneurism | 256$\times$256$\times$256 | 60767
Bonsai | 256$\times$256$\times$256 | 210570
Mrt Angio | 416$\times$512$\times$112 | 5146370
Heptane Gas | 302$\times$302$\times$302 | 61212
Stent | 512$\times$512$\times$174 | 3083028
Pancreas | 240$\times$512$\times$512 | 7215261
Backpack | 512$\times$512$\times$373 | 7066081
Magnetic Reconnection | 512$\times$512$\times$512 | 31147751
Zeiss | 680$\times$680$\times$680 | 3209058
Marmoset Neurons | 1024$\times$1024$\times$314 | 56384791
Stag Beetle | 832$\times$832$\times$494 | 852383
Pawpawsaurus | 958$\times$646$\times$1088 | 83014508
Spathorhynchus | 1024$\times$1024$\times$750 | 47569170
Kingsnake | 1024$\times$1024$\times$795 | 32483997
Chameleon | 1024$\times$1024$\times$1080 | 49325513
Beechnut | 1024$\times$1024$\times$1546 | 159472969
Richtmyer Meshkov | 2048$\times$2048$\times$1920 | 31585707
Woodbranch | 2048$\times$2048$\times$2048 | 1520839571
3D Neurons | 2048$\times$2048$\times$2384 | 1913765935
Pig Heart | 2048$\times$2048$\times$2612 | 906730898
Schwefel 3D | 128$\times$128$\times$128 to | 1688
1024$\times$1024$\times$1024
Table 2: Datasets used in the computational experiments [Ope22], their size,
and the number of critical points. The Schwefel dataset [Sch81] is resampled
on domains of different sizes.

###  6.1 Internal memory computation

We first evaluate the runtime performance of the algorithm for datasets that
fit within the GPU memory. In this experiment, we compute the wall clock time
for computing the extremum graphs on a large collection of datasets, see Table
2. These computations did not require the hybrid GPU-CPU execution. Figure 7
shows the total runtime and the split up between the two major steps and the
time for data structure updates following the computation. Datasets that are
smaller in size compared to Magnetic Reconnection require under $2.3$s. The
figure shows runtimes only for larger datasets that fit in memory. Our key
observations are:

  * •

Time required for critical point classification on the GPU scales linearly
with the domain size.

  * •

Time required for gradient path tracing on the CPU scales roughly linear with
the number of identified saddles.

Figure 7: Extremum graph computation for datasets that fit in GPU memory. The
time taken increases with size of dataset and naturally depends on the number
of critical points.

We also perform this experiment for higher dimensional data by sampling the
Schwefel function [Sch81] on 3-to-6 dimensional domains of varying sizes.
Results are shown in Table 3. We again observe an increase with size of data,
particularly the time taken for critical point classification.

Dataset | Size | #Critical Pts | Computation time (in secs)
---|---|---|---
Schwefel 4D | $32^{4}$ | 8449 | 0.20
$64^{4}$ | 15072 | 0.71
$128^{4}$ | 15072 | 7.84
Schwefel 5D | $16^{5}$ | 9525 | 0.33
$32^{5}$ | 62167 | 3.87
$64^{5}$ | 128808 | 110.01
Schwefel 6D | $16^{6}$ | 46345 | 9.13
$24^{6}$ | 444241 | 104.46
$32^{6}$ | 444241 | 432.54
Table 3: Extremum graph computation for higher dimensional scalar fields. A
significant fraction of time was taken for the critical point classification
step due to the increased neighborhood size.

###  6.2 Hybrid GPU-CPU hybrid computation

In this experiment, we choose the large datasets that do not fit in GPU memory
and required partitioning into multiple blocks. Figure 8 shows the runtime
results for these datasets. We make the following observations:

  * •

Critical point classification time is nearly identical for all datasets
because they have similar size.

  * •

The gradient path tracing step consumes a large fraction of total runtime.

  * •

An additional factor that influences gradient path tracing time is the number
of paths that cross block boundaries. Such paths carry a computational
overhead.

Figure 8: Performance for large datasets that do not fit in GPU memory.
Runtime broadly depends on data size. Pig Heart has a small number of critical
points, which explains the smaller runtime. Gradient path tracing constitutes
a significant fraction of total runtime in all cases.

###  6.3 GPU scaling - varying number of cores

We evaluate the performance of the critical point classification step across
different GPUs. These experiments are performed on the dataset Zeiss. The aim
is to understand whether the computation scales to the number of available
compute units on the GPU. Total computation time is not compared in this
experiment because it is strongly influenced by the number of CPU cores and
additional hardware features. The results are shown in Figure 9. We observe a
monotonically decreasing trend with a good slope, indicating efficient usage
of the GPU cores.

Figure 9: GPU scaling results.

###  6.4 GPU scaling - varying domain size

This experiment aims to establish how well tachyon utilizes the available GPU
compute resources. We study how performance varies upon increasing the domain
size while fixing the number of critical points. The Schwefel function [Sch81]
is sampled on a 3D domain. The data, originally available over a $500\times
500\times 500$, is sampled to obtain datasets over domains whose sizes range
from $128\times 128\times 128$ to $1024\times 1024\times 1024$. Across all
these sampled functions, the number of critical points and gradient paths do
not change thereby allowing a study of the variation of GPU computation time
with domain size. All GPU cores are used to compute the extremum graph for
this dataset available at different resolutions. The results are shown in
Figure 10. We observe a strong linear dependency on domain size.

Figure 10: GPU scaling study for varying domain size. Experiments on datasets
obtained by sampling the 3D Schwefel function [Sch81] over different domain
sizes reveal a linear relationship between computation time and domain size.

###  6.5 CPU scaling - varying CPU threads

Next, we study scaling behavior by varying the number of CPU threads. The aim
is to study how well tachyon utilizes the available parallelism provided by
the CPU. In the ideal case, we expect the computation time to reduce by a
factor of two on doubling the number of threads. The results of the experiment
are shown in Figure 11. We observe that the performance improvements continue
up to 8 threads and taper off afterwards.

Figure 11: CPU scaling study for varying number of CPU threads. Experiments on
the dataset Zeiss show good scaling up to 8 threads.

###  6.6 Performance comparisons

We now describe results of runtime comparisons against TTK and pyms3d. Both
TTK [TFL∗17] and pyms3d [SN17] support computation of the complete Morse-Smale
complex. The key difference from tachyon is that they employ a discrete Morse
theory based approach. As a first step, they compute a discrete gradient field
on the GPU to locate all critical cells. In contrast, tachyon uses a lower
link based critical point classification to locate nodes of the extremum
graph. They employ a parallel root finding algorithm to compute the entire
descending 3-manifold of all extrema and hence compute saddle-maxima arcs. In
contrast, tachyon employs gradient path tracing from each saddle to compute
the arcs. The gradient field and saddle-extrema arc computation in gmsc
[SPN22] and pyms3d are identical. The only difference between the two
implementations is that gmsc uses CUDA whereas pyms3d uses OpenCL.

In order to ensure a fair comparison, we configure pyms3d and TTK to
selectively compute the ascending 1-manifolds of 2-saddles and hence output
the extremum graph. The results show that our implementation is able to
compute extremum graphs faster. In some instances, both TTK and pyms3d fail
due to system (GPU/main memory) resource limitations. We find that tachyon
performs better, both in terms of running time and utilization of system
resources. Figures 12 and 13 plot the runtimes for increasing dataset sizes.
The running time for tachyon rises slower compared to TTK and pyms3d. Further,
pyms3d ran out of memory and could not complete execution for larger data
sizes. The comparison with pyms3d was performed on a workstation with an
Intel(R) Core(TM) i5-4590 CPU @ 3.30GHz 2 Cores/4 Threads, 16GB RAM, and an
Nvidia GeForce GTX 1060 GPU with 6GB GDDR5 RAM.

Figure 12: Runtime comparison between tachyon and TTK[TFL∗17]. tachyon
performs consistently better, often several orders of magnitude faster. Figure
13: Runtime comparison between tachyon and pyms3d [SN17]. tachyon performs
consistently better, often several orders of magnitude faster.

###  6.7 Simplification

Finally, we study the ability of the simplification operations to remove noise
and clutter. We perform three operations on various datasets with a fixed
parameter set: edge bundling, persistence directed cancellation with $5$%
threshold and saturated persistence directed simplification with high and low
thresholds $p_{lo}=5\%$ and $p_{hi}=95\%$. Figure 5 shows the simplified
extremum graphs. We also study statistics of the size of extremum graph before
and after simplification, see Figure 14. In some noisy datasets, a small
threshold already removes a large fraction of nodes and arcs whereas in other
datasets that contain a smaller number of features, the reduction is not
significant.

Across all results shown in Figure 14, we observe that the number of arcs is a
little over twice the number of saddles, although these absolute numbers are
not depicted in the graph plot. This is expected, because all saddles have at
least two outgoing arcs.

---
Figure 14: Statistics on the fraction of extremum graph nodes and arcs that
survive post simplification. The number of surviving extrema, saddles, and
arcs after applying the three simplification operations in sequence are shown
as a fraction of the corresponding numbers in the unsimplified extremum graph.

##  7 Conclusions and Future Work

We have developed tachyon, a dimension independent, scalable extremum graph
computation library that also supports a hybrid GPU-CPU mode to process large
datasets. This software library will be released in the public domain for use
by the community. We also describe three methods for simplifying the extremum
graph with the aim of identifying and removing noise, resulting in a clutter
free graph that is more amenable for further analysis. Computational
experiments demonstrate good running times and scaling properties.

Future work includes incorporating support for generic non-uniform grids and
an extension to a distributed setting where the dataset does not fit in system
memory. For the latter problem, an approach similar to the hybrid GPU-CPU mode
where different processors work on different blocks of the same dataset may
work but the challenge is to reduce the communication necessary to consolidate
the output.

## Acknowledgments

This work is partially supported by a Swarnajayanti Fellowship from the
Department of Science and Technology, India (DST/SJF/ETA-02/2015-16), an SERB
grant CRG/2021/00527, MoE Govt. of India, and a Mindtree Chair research
grant).

## References

  * [Ban70] Banchoff T. F.:  Critical points and curvature for embedded polyhedral surfaces.  _Am. Math. Monthly 77_ (1970), 475–485.
  * [BGL∗18] Bhatia H., Gyulassy A. G., Lordi V., Pask J. E., Pascucci V., Bremer P.-T.:  Topoms: Comprehensive topological exploration for molecular and condensed-matter systems.  _Journal of computational chemistry 39_ , 16 (2018), 936–952.
  * [BHEP04] Bremer P.-T., Hamann B., Edelsbrunner H., Pascucci V.:  A topological hierarchy for functions on triangulated surfaces.  _IEEE Transactions on Visualization and Computer Graphics 10_ , 4 (2004), 385–396.
  * [CFST20] Carr H., Fujishiro I., Sadlo F., Takahashi S.:  _Topological Methods in Data Analysis and Visualization V_.  Springer, 2020.
  * [CLB11] Correa C., Lindstrom P., Bremer P.-T.:  Topological spines: A structure-preserving visual representation of scalar fields.  _IEEE Transactions on Visualization and Computer Graphics 17_ , 12 (2011), 1842–1851.
  * [CMS06] Carr H., Moller T., Snoeyink J.:  Artifacts caused by simplicial subdivision.  _IEEE Transactions on Visualization and Computer Graphics 12_ , 2 (2006), 231–242.
  * [CWS∗19] Carr H. A., Weber G. H., Sewell C. M., Rübel O., Fasel P., Ahrens J. P.:  Scalable contour tree computation by data parallel peak pruning.  _IEEE transactions on visualization and computer graphics 27_ , 4 (2019), 2437–2454.
  * [DFFIM15] De Floriani L., Fugacci U., Iuricich F., Magillo P.:  Morse complexes for shape segmentation and homological analysis: discrete models and algorithms.  _Computer Graphics Forum 34_ , 2 (2015), 761–785.
  * [DFRS14] Delgado-Friedrichs O., Robins V., Sheppard A.:  Skeletonization and partitioning of digital images using discrete Morse theory.  _IEEE Transactions on Pattern Analysis and Machine Intelligence 37_ , 3 (2014), 654–666.
  * [EH09] Edelsbrunner H., Harer J.:  _Computational Topology: An Introduction_.  Amer. Math. Soc., Providence, Rhode Island, 2009.
  * [EHNP03] Edelsbrunner H., Harer J., Natarajan V., Pascucci V.:  Morse-Smale complexes for piecewise linear 3-manifolds.  In _Proceedings of the nineteenth annual symposium on Computational geometry_ (2003), pp. 361–370.
  * [ELZ02] Edelsbrunner H., Letscher D., Zomorodian. A.:  Topological persistence and simplification.  _Disc. Comput. Geom. 28_ , 4 (2002), 511–533.
  * [FIDF19] Fugacci U., Iuricich F., De Floriani L.:  Computing discrete Morse complexes from simplicial complexes.  _Graphical Models 103_ (2019), 101023.
  * [GBHP08] Gyulassy A., Bremer P.-T., Hamann B., Pascucci V.:  A practical approach to Morse-Smale complex computation: Scalability and generality.  _IEEE Transactions on Visualization and Computer Graphics 14_ , 6 (2008), 1619–1626.  doi:10.1109/TVCG.2008.110.
  * [GBP12] Gyulassy A., Bremer P.-T., Pascucci V.:  Computing Morse-Smale complexes with accurate geometry.  _IEEE Transactions on Visualization and Computer Graphics 18_ , 12 (2012), 2014–2022.
  * [GBP18] Gyulassy A., Bremer P.-T., Pascucci V.:  Shared-memory parallel computation of Morse-Smale complexes with improved accuracy.  _IEEE Transactions on Visualization and Computer Graphics 25_ , 1 (2018), 1183–1192.
  * [GGL∗14] Gyulassy A., Günther D., Levine J. A., Tierny J., Pascucci V.:  Conforming Morse-Smale complexes.  _IEEE Transactions on Visualization and Computer Graphics 20_ , 12 (2014), 2595–2603.
  * [GNP∗06] Gyulassy A., Natarajan V., Pascucci V., Bremer P.-T., Society C., Hamann B.:  A topological approach to simplification of three-dimensional scalar functions.  _IEEE transactions on visualization and computer graphics 12_ (08 2006), 474–84.  doi:10.1109/TVCG.2006.57.
  * [GPPR12] Gyulassy A., Pascucci V., Peterka T., Ross R.:  The parallel computation of Morse-Smale complexes.  In _2012 IEEE 26th International Parallel and Distributed Processing Symposium_ (2012), IEEE, pp. 484–495.
  * [GRWH12] Günther D., Reininghaus J., Wagner H., Hotz I.:  Efficient computation of 3D Morse–Smale complexes and persistent homology using discrete Morse theory.  _The Visual Computer 28_ , 10 (2012), 959–969.
  * [HLH∗16] Heine C., Leitte H., Hlawitschka M., Iuricich F., De Floriani L., Scheuermann G., Hagen H., Garth C.:  A survey of topology-based methods in visualization.  _Computer Graphics Forum 35_ , 3 (2016), 643–667.
  * [HMMN14] Harker S., Mischaikow K., Mrozek M., Nanda V.:  Discrete Morse theoretic algorithms for computing homology of complexes and maps.  _Foundations of Computational Mathematics 14_ , 1 (2014), 151–184.
  * [HMST21] Hotz I., Masood T. B., Sadlo F., Tierny J.:  _Topological Methods in Data Analysis and Visualization VI_.  Springer, 2021.
  * [NTN15] Narayanan V., Thomas D. M., Natarajan V.:  Distance between extremum graphs.  In _2015 IEEE Pacific Visualization Symposium (PacificVis)_ (2015), IEEE, pp. 263–270.
  * [Ope22] Open scientific visualization datasets.  https://klacansky.com/open-scivis-datasets/, 2022.  [Online; accessed 01-July-2022].
  * [PRG∗11] Peterka T., Ross R., Gyulassy A., Pascucci V., Kendall W., Shen H.-W., Lee T.-Y., Chaudhuri A.:  Scalable parallel building blocks for custom data analysis.  In _2011 IEEE Symposium on Large Data Analysis and Visualization_ (2011), IEEE, pp. 105–112.
  * [RWS11] Robins V., Wood P. J., Sheppard A. P.:  Theory and algorithms for constructing discrete Morse complexes from grayscale digital images.  _IEEE Transactions on Pattern Analysis and Machine Intelligence 33_ , 8 (2011), 1646–1658.
  * [Sch81] Schwefel H.-P.:  _Numerical Optimization of Computer Models_.  John Wiley & Sons, Inc., 1981.
  * [SN12] Shivashankar N., Natarajan V.:  Parallel computation of 3D Morse-Smale complexes.  _Computer Graphics Forum 31_ , 3 (2012), 965–974.
  * [SN17] Shivashankar N., Natarajan V.:  Efficient software for programmable visual analysis using Morse-Smale complexes.  In _Topological Methods in Data Analysis and Visualization_ (2017), pp. 317–331.  doi:10.1007/978-3-319-44684-4_19.
  * [SPN20] Subhash V., Pandey K., Natarajan V.:  GPU parallel computation of Morse-Smale complexes.  In _Proc. IEEE Visualization Conference, IEEE VIS 2020 - Short Papers_ (2020), IEEE, pp. 36–40.  doi:10.1109/VIS47514.2020.00014.
  * [SPN22] Subhash V., Pandey K., Natarajan V.:  A gpu parallel algorithm for computing Morse-Smale complexes.  _IEEE Transactions on Visualization and Computer Graphics_ (2022), 1–1.  doi:10.1109/TVCG.2022.3174769.
  * [SSN11] Shivashankar N., Senthilnathan M., Natarajan V.:  Parallel computation of 2d Morse-Smale complexes.  _IEEE Transactions on Visualization and Computer Graphics 18_ , 10 (2011), 1757–1770.
  * [TFL∗17] Tierny J., Favelier G., Levine J. A., Gueunet C., Michaux M.:  The Topology ToolKit.  _IEEE Transactions on Visualization and Computer Graphics_ (2017).  https://topology-tool-kit.github.io/.
  * [TN13] Thomas D. M., Natarajan V.:  Detecting symmetry in scalar fields using augmented extremum graphs.  _IEEE transactions on visualization and computer graphics 19_ , 12 (2013), 2663–2672.