A Practical Algorithm for (Geometry-Aware) Interleavings Between Merge Trees

Consider two merge trees $(T_1,f_1)$ and $(T_2,f_2)$. The interleaving distance between $T_1$ and $T_2$ was originally defined in terms of two maps, one from $T_1$ to $T_2$, and one from $T_2$ to $T_1$ [11]. We use an equivalent definition [7] in terms of a single map from $T_1$ to $T_2$ (see Figure 2):

Touli and Wang [17] describe a decision procedure to determine if a $\delta$-good map³³3We note that Definition 2 slightly differs from the one given in [17]. A proof that these definitions are equivalent can be found in [7, Appendix A]. exists for a given value $\delta$. We sketch this procedure below; see the original paper for a full description. To compute the interleaving distance, they search on a finite set of candidate values. In particular, it has been shown [1, 17] that $d(T_1,T_2)$ is always present in the set $\mathrm{\Delta }(T_1,T_2)$ that contains (1) all height differences between a vertex in $T_1$ and a vertex in $T_2$, and (2) all height differences, divided by 2, between vertex pairs in the same tree.

Consider a pair $(X,y)$, where $X$ is a set of points of $T_1$, all at the same height $h$, and $y$ is a single point of $T_2$ at height $h+\delta$. A continuous map $\alpha$ from $T_1[X]$ to $T_2[y]$ is a partial $\delta$-good map if P1–P3 hold for all points of $T_1[X]$ and $T_2[y]$. Note that for P2, we still determine the lowest common ancestor in the complete tree $T_1$. We say $(X,y)$ is feasible if a partial $\delta$-good map from $T_1[X]$ to $T_2[y]$ exists. Touli and Wang [17, Claim 2] show that a $\delta$-good map between the complete trees $T_1$ and $T_2$ exists if and only if there is a feasible pair $(X,y)$ for which $X$ and $y$ both lie above the top vertices in their respective trees.

We maintain for each point $y$ of $T_2$ on the sweepline, all non-empty sets $X$ of points of $T_1$ on the sweepline such that $(X,y)$ is feasible. Let ${\mathrm{\ell }}_1$ and ${\mathrm{\ell }}_2$ be the number of leaves in $T_1$ and $T_2$, respectively. We assume the leaves in each tree $T_k$ are indexed $1,2,\mathrm{\dots },{\mathrm{\ell }}_k$ from left to right in such a way that the indices of the set of leaves in a subtree always form an integer interval; this can be achieved using an in-order tree walk on the input. We index each point on the sweepline by the lowest-index leaf in its subtree (see Figure 4b); when we discuss points on the sweepline, we write “vertex $i$ in $T_k$” instead of “vertex in $T_k$ with index $i$”. If an internal vertex $i$ is on the sweepline, its two descendants just below the sweepline have indices $i$ and $i^{\prime }$, where . We refer to $i$ and $i^{\prime }$ as the children of $i$.

We maintain the set of feasible pairs in a list $F$ of ${\mathrm{\ell }}_2$ elements. Consider an index $j$ from $T_2$. If no point $y$ of $T_2$ on the sweepline has index $j$, then $F[j]=\text{nil}$. Otherwise, if there is a point $y$, then $F[j]$ represents all non-empty sets $X$ of points of $T_1$ on the sweepline such that $(X,y)$ is feasible. We store every such set $X$ as the corresponding set $I$ of indices (see Figure 4c). When the sweepline hits a vertex, an event happens and we need to update $F$:

1. (a)

   Leaf $j$ in $T_2$: Initialize $F[j]=\mathrm{\varnothing }$.

2. (b)

   Internal vertex $j$ in $T_2$: Let $j,j^{\prime }$ be the children of $j$. Set $G=\mathrm{\varnothing }$. For each $I_1\in F[j]$ and $I_2\in F[j^{\prime }]$, add $I_1\cup I_2$ to $G$ if and only if $I_1$ and $I_2$ (i) have the same $2\delta$-ancestor, and (ii) are disjoint (see Figure 5). Additionally, if the depth of the left (right) subtree is at most $2\delta$, add all elements of $F[j^{\prime }]$ ($F[j]$) to $G$. Finally, set $F[j]=G$ and $F[j^{\prime }]=\text{nil}$.

3. (c)

   Leaf $i$ in $T_1$: For each index $j$ from $T_2$ on the sweepline: (1) for each set $I$ in $F[j]$, add the set $I\cup \{i\}$ to $F[j]$ if and only if $I$ and $i$ have the same $2\delta$-ancestor, and (2) if the depth of the point of $T_2$ on the sweepline with index $j$ is at most $2\delta$, add $\{i\}$ to $F[j]$.

4. (d)

   Internal vertex $i$ in $T_1$: Let $i,i^{\prime }$ be the children of $i$. For each index $j$ from $T_2$ on the sweepline, for each set $I$ in $F[j]$, if $I$ contains exactly one of $i$ and $i^{\prime }$, delete $I$ from $F[j]$; if $I$ contains both, replace $I$ by $I\setminus \{i^{\prime }\}$.

The sweepline terminates after handling the last event, that is, when the top vertices of both $T_1$ and $T_2$ have been handled. Consider the set $F$ after the sweepline terminates. Let $y$ be the top vertex in $T_2$. By construction, $y$ has index $1$. We conclude that a $\delta$-good map exists if and only if (a) the depth of $y$ is at most $2\delta$, or (b) $F[1]$ contains at least one element.

To prove that our sweep is correct, we make two observations. Consider a partial $\delta$-good map $\alpha$ from $T_1[X]$ to $T_2[y]$, for some set of points $X$ of $T_1$ and a point $y$ of $T_2$. Then we can extend $\alpha$ until we hit an internal vertex: let $X^{\prime }$ be the ancestors of points in $X$ at a fixed height, and assume that there are no internal vertices between points in $X$ and their ancestor in $X^{\prime }$. We define the extension of $\alpha$ to $T_1[X^{\prime }]$: for every point $x$ of $T_1[X^{\prime }]$ that is not in $T_1[X]$, let $y^{\prime }$ be the ancestor of $y$ at height $f_1(x)+\delta$. We set $\alpha (x)≔y^{\prime }$. Since there are no new internal vertices in $T_1[X^{\prime }]$, the extension is well-defined for all points of $T_1[X^{\prime }]$. Moreover, $\alpha$ is continuous and satisfies P1–P3.

On the other hand, we can also restrict $\alpha$ such that its image remains a subtree of $T_2$. The resulting map is again a partial $\delta$-good map: P1–P3 follow directly.

We are now ready to prove that the sweepline maintains all non-empty feasible pairs.

For efficient event handling, we store $F$ in a suitable data structure. Specifically, we store each (non-nil) element $F[j]$ in a balanced binary search tree, such as a red–black tree, which allows efficient searches, insertions, and deletions. When storing some $I\in F[j]$ in the tree, we use the lowest-indexed element of $I$ as its key. This way, the sets in $F[j]$ are always ordered by their lowest-indexed element, which allows us to efficiently handle the $2\delta$-ancestor check in (b) events as follows.

Event (b) asks us to try all combinations of $I_1\in F[j]$ and $I_2\in F[j^{\prime }]$ and see if they have the same $2\delta$-ancestor. Instead, we check for each potential $2\delta$-ancestor which sets $I_1$ and $I_2$ descend from it, and only try these combinations. Let $A=[a_1,\mathrm{\dots },a_k]$ be the (ordered) sequence of indices of $T_1$ at height $2\delta$ above the sweepline, that is, $A$ is the sequence of potential $2\delta$-ancestors of $I_1$ and $I_2$. The problem boils down to finding, for a given ancestor $a_i$, the subset of $F[j]$ that have $a_i$ as its $2\delta$-ancestor. This is the subset of all $I_1\in F[j]$ for which the lowest-indexed element of $I_1$ is between $a_i$ (inclusive) and $a_{i+1}$ (exclusive). As $F[j]$ is sorted on the lowest-indexed element, this subset of elements is a contiguous range, which we can delineate by searching in the binary search tree for $\{a_i\}$ and $\{a_{i+1}\}$.

We note that this data structure admits an easy and efficient implementation in practice if we represent each set $I\in F[j]$ as a bit string of which the $i$-th bit is set if and only if $i\in I$. Such a representation not only allows us to compute unions and intersections efficiently, but also sorting the bit strings lexicographically ensures that they are sorted by their lowest-indexed element.

Let $(T,f)$ be a merge tree, and let $\epsilon >0$ such that $\epsilon$ is smaller than the height difference between any two adjacent vertices in $T$. A merge tree $(T_{\epsilon },f_{\epsilon })$ is an $\epsilon$-perturbation of $(T,f)$ if there is a homeomorphism $\eta$ from $T$ to $T_{\epsilon }$ such that $|f(x)-f_{\epsilon }(\eta (x))|\le \epsilon$ for all points $x$ of $T$ (see Figure 6a). It is not hard to show that the interleaving distance between $(T,f)$ and an $\epsilon$-perturbation $(T_{\epsilon },f_{\epsilon })$ is at most $\epsilon$.

Suppose $T$ has a vertex $x$ that has $k>2$ children $x_1,\mathrm{\dots },x_k$. We expand $x$ into $k-1$ distinct vertices ${\widehat{x}}_1,\mathrm{\dots },{\widehat{x}}_{k-1}$ such that each vertex ${\widehat{x}}_i$ has two children: $x_i$ and ${\widehat{x}}_{i+1}$ (see Figure 6b). We define the expansion $\widehat{T}$ of $T$ by expanding all vertices in $T$, together with the height function $\widehat{f}$ that assigns to each new vertex ${\widehat{x}}_i$ the same height as $x$. An extension $(\widehat{T},\widehat{f})$ is not a merge tree; however, for any value $\epsilon >0$ there does exists an $\epsilon$-perturbation of $(\widehat{T},\widehat{f})$ that is a merge tree. We define the interleaving distance between two expansions ${\widehat{T}}_1$ and ${\widehat{T}}_2$ as the limit of the interleaving distance between such $\epsilon$-perturbed merge trees ${\widehat{T}}_1^{\epsilon }$ and ${\widehat{T}}_2^{\epsilon }$:

Lemma 6 then directly implies Corollary 7.

In the algorithm we handle expansions by placing the expanded vertices on exactly the same height, and separating them symbolically. That is, the vertices all lie on the same level, and when we process that level, we handle them one by one.

We compute the $ℛ$-restricted interleaving distance between $T_1$ and $T_2$ using a similar approach as for computing the unrestricted interleaving distance. More precisely, we can still search on the set of candidate values, but instead of searching over all values in $\mathrm{\Delta }(T_1,T_2)$, we search over all values in $\mathrm{\Delta }(T_1,T_2,ℛ)$. Moreover, to decide for a given value $\delta$ whether an $ℛ$-restricted $\delta$-good map exists, we use a modified version of either the DP as described in Section 2, or the sweepline algorithm as described in Section 3.

Recall that ${\mathrm{\ell }}_1$ and ${\mathrm{\ell }}_2$ denote the number of leaves in $T_1$ and $T_2$, respectively. We represent a restriction $ℛ$ using a matrix $R$ of size ${\mathrm{\ell }}_1\times {\mathrm{\ell }}_2$. Consider a leaf $x$ in $T_1$ with index $i$, and a leaf $y$ in $T_2$ with index $j$. The entry $R_{ij}$ stores the height of the lowest ancestor of $y$ that is in the pruning $ℛ[x]$. We refer to $R$ as the restriction matrix of $ℛ$ (see Figure 7b).

Consider a (finite) sequence of scalar fields $(M,f_t)$, that is, a 2D or 3D domain $M$ whose height function changes over time. Let $d_e$ denote the distance between points of $M$. Fix two instances $(M,f_1)$ and $(M,f_2)$, and let $(T_1,f_1)$ and $(T_2,f_2)$ be the corresponding scalar field-based merge trees. Without loss of generality, assume $T_1$ and $T_2$ are based on sublevel sets; their leaves correspond to minima in their respective scalar fields. Intuitively, assuming the two instances are not too far apart and the scalar field evolves smoothly over time, we do not expect big changes in the geometry of the extrema: each minimum in $(M,f_1)$ has a corresponding minimum in $(M,f_2)$ close by. We use this to define a restriction from $T_1$ to $T_2$: a leaf in $T_1$ is allowed to map to a point of $T_2$ if the corresponding minima are not “too far” apart. Formally, we define a restriction ${ℛ}_{\rho }$, where $\rho \ge 0$ is a threshold variable that captures how far is too far. Let $x$ be a leaf in $T_1$. Then $x$ corresponds to a minimum in $(M,f_1)$; let $p_x$ be the point of $M$ where the minimum lies. Similarly, each point $y$ of $T_2$ corresponds to a sublevel set component; let $P_y$ denote the sublevel set component that contains $y$. We define $y\in {ℛ}_{\rho }[x]$ if and only if $d_e(p_x,P_y)\le \rho$. This defines a pruning: if $y\in {ℛ}_{\rho }[x]$, then the sublevel set component $P_{y^{\prime }}$ of any ancestor $y^{\prime }$ of $y$ contains $P_y$, so $y^{\prime }\in {ℛ}_{\rho }$.

Our idea of a geometry-aware heuristic is inspired by a heuristic by Yan et al. [19]. Specifically, they also use the geometry of the underlying scalar fields as a heuristic to compute the interleaving distance. However, their method uses a matching between extrema of the scalar fields to compute a labeling on the merge trees; the resulting interleaving distance is then the labeled interleaving distance for that specific labeling. Our method is more general: each leaf is allowed to map to more than one point. In particular, if we set $\rho =\mathrm{\infty }$, then ${ℛ}_{\rho }[x]=T_2$ for each leaf $x$ in $T_1$ and we recover the unrestricted interleaving distance.

We ran our algorithm on all instances of HeatedFlow and RedSea. We used only the exponential search procedure, as running the linear search procedure was exceedingly slow and did not finish in a reasonable amount of time. We also set a timeout at 300 seconds for each run. The DP procedure reached the timeout 5 times for HeatedFlow, while the sweepline procedure never reached the timeout. We summarize the accumulated results in Table 1; in this table we count timeouts as contributing just 300 seconds to the total time.