# Computing Persistence Diagram Bundles

Abigail Hickok

###### Abstract.

Persistence diagram (PD) bundles, a generalization of vineyards, were
introduced as a way to study the persistent homology of a set of filtrations
parameterized by a topological space $B$. In this paper, we present an
algorithm for computing piecewise-linear PD bundles, a wide class that
includes many of the PD bundles that one may encounter in practice. Full
details are given for the case in which $B$ is a triangulated surface, and we
outline the generalization to higher dimensions and other cases.

##  1\. Introduction

In persistent homology, one is given a filtration (e.g., a filtered complex),
and one studies how the homology changes as the filtration parameter
increases. Here, we consider the case in which one is instead given a _fibered
filtration function_ , which is a set
$\\{f_{p}:\mathcal{K}^{p}\to\mathbb{R}\\}_{p\in B}$ of filtration functions
that is parameterized by some topological space $B$ (the _base space_), where
$\mathcal{K}^{p}$ is a simplicial complex for each $p\in B$. At each $p\in B$,
the sublevel sets of $f_{p}$ form a filtered complex. The associated
_persistence diagram (PD) bundle_ is the space of persistence diagrams
$\text{PD}(f_{p})$ as they vary with $p\in B$.

PD bundles arise naturally in many circumstances. The prototypical PD bundle
is a _vineyard_ , introduced in [3], which is the case where $B$ is an
interval in $\mathbb{R}$. For example, given a time-varying point cloud, one
obtains a vineyard by constructing a filtration (such as the Vietoris–Rips
filtration) at each time. Figure 1 shows an illustration of a vineyard,
visualized as a continuously-varying “stack of PDs”. More generally, however,
one might have a set $\\{X(p)\\}_{p\in B}$ of point clouds parameterized by an
arbitrary parameter space $B$. One obtains a PD bundle by constructing a
filtration for the point cloud at each $p\in B$. For example, biological
aggregation models (e.g., the D’Orsogna model [2]) produce time-varying point
clouds whose coordinates also depend on various real-valued system parameters
$\mu_{1},\ldots,\mu_{k}$. The parameter space $B$ in this example is a subset
of $\mathbb{R}^{k+1}$. For each $(t,\mu_{1},\ldots,\mu_{k})\in B$, there is a
point cloud from which one can construct a filtration. Another special case is
the persistent homology transform ($B=\mathbb{S}^{d}$), which is used in the
field of shape analysis [15]. Other concrete examples of PD bundles are given
in [7].

Figure 1. An illustration of a vineyard, which consists of a persistence
diagram for each time $t$. (This figure is a slightly modified version of a
figure that appeared originally in [9], which is available under a Creative
Commons license.)

###  1.1. Contributions

I generalize Cohen-Steiner et. al’s algorithm for computing vineyards [3] to
an algorithm for efficiently computing PD bundles. I restrict to the case in
which the PD bundle is _piecewise linear_. This means that $B$ is a simplicial
complex, $\mathcal{K}^{p}\equiv\mathcal{K}$ for all $p\in B$, and for every
simplex $\sigma\in\mathcal{K}$, the function $f_{\sigma}(p):=f_{p}(\sigma)$ is
linear in $p$ on every simplex of $B$. The restriction to piecewise-linear PD
bundles allows us to take advantage of methods in computational geometry, such
as the Bentley–Ottman planesweep algorithm [4] for finding intersections of
lines in a plane and algorithms for solving the point-location problem in a
line arrangement [13]. An analogous piecewise-linear restriction was helpful
for computing vineyards in [3].

The idea of the algorithm is to subdivide the base space $B$ into polyhedrons
and compute a PD “template” for each polyhedron. The subdivision is given by
Proposition 2.7 ([7]). For any $p\in B$, the persistence diagram
$\text{PD}(f_{p})$ can then be computed in $\mathcal{O}(N)$ time from the PD
template for the polyhedron that contains $p$, where $N$ is the number of
simplices in $\mathcal{K}$. By contrast, computing PD$(f_{p})$ from scratch
takes $\mathcal{O}(N^{3})$ time in the worst case.

The piecewise-linear restriction is reasonable for most applications. For
example, suppose that we have a point cloud $X(t,\mu)$ whose coordinates
depend on time $t$ and a parameter $\mu\in\mathbb{R}$. If the data set arises
from either real-world data collection or through numerical simulation, then
we likely only know the coordinates of the point cloud $X(t,\mu)$ at a
discrete set $\\{t_{i}\\}$ of time steps and a discrete set $\\{\mu_{j}\\}$ of
system-parameter values. For every $(t_{i},\mu_{j})$, there is the filtration
function $f_{(t_{i},\mu_{j})}:\mathcal{K}\to\mathbb{R}$ associated with the
Vietoris–Rips filtration (or any other filtration) of $X(t_{i},\mu_{j})$. For
the Vietoris–Rips filtration, $\mathcal{K}$ is the simplicial complex that
contains a simplex for every subset of points in the point cloud. To obtain a
fibered filtration function, we define $B$ to be a triangulation of $[\min
t_{i},\max t_{i}]\times[\min\mu_{j},\max\mu_{j}]$ whose vertices are
$\\{(t_{i},\mu_{j})\\}_{ij}$. We can extend $\\{f_{(t_{i},\mu_{j})}\\}_{ij}$
to a fibered filtration function with base space $B$ by defining the
filtration values of a simplex $\sigma$ via linear interpolation of
$\\{f_{(t_{i},\mu_{j})}(\sigma)\\}_{ij}$. By construction, the resulting PD
bundle is piecewise linear.

Full details are only given for the case in which $\dim(B)$ is a triangulated
surface, but the generalization to higher dimensions is discussed in Section
3.2. When the base space $B$ is a triangulated surface, it is already an
improvement over a vineyard because it allows three parameters in total: a
filtration parameter as well as two parameters that locally parameterize $B$.

###  1.2. Related Work

PD bundles were introduced in [7] as a generalization of vineyards [3]. The
algorithm that I present in this paper for computing PD bundles is a broad
generalization of the algorithm in [3]. The algorithm in this paper is also
related to the Rivet algorithm for computing fibered barcodes of 2D
multiparameter persistence modules [8].

###  1.3. Organization

The paper proceeds as follows. I review the relevant background on persistent
homology, vineyards, and PD bundles in Section 2. I present my algorithm for
computing piecewise-linear PD bundles in Section 3. Finally, I conclude and
discuss possible directions for future research in Section 4.

##  2\. Background

We begin by reviewing persistent homology, vineyards, and PD bundles; for more
details on persistent homology, see [5, 12], for more on vineyards, see [3],
and for an introduction to PD bundles, see [7].

###  2.1. Persistent homology

Let $\mathcal{K}$ be a simplicial complex. A _filtration function_
$f:\mathcal{K}\to\mathbb{R}$ is a real-valued function on $\mathcal{K}$ that
is _monotonic_. That is, $f(\tau)\leq f(\sigma)$ if $\tau$ is a face of
$\sigma$. Monotonicity guarantees that the $r$-sublevel sets
$\mathcal{K}_{r}:=\\{\sigma\in\mathcal{K}\mid f(\sigma)\leq r)\\}$ are
simplicial complexes.

In persistent homology, we study how the homology of $\mathcal{K}_{r}$ changes
as $r$ increases. Let $\\{r_{i}\\}_{i=1}^{k}$ be the image of $f$, ordered
such that $r_{i}<r_{i+1}$. These are the critical values at which
$\mathcal{K}_{r}$ changes; for $r\in[r_{i},r_{i+1})$, we have
$\mathcal{K}_{r}=\mathcal{K}_{r_{i}}$. For every $i\leq j$, the inclusion
$\iota^{i,j}:\mathcal{K}_{r_{i}}\hookrightarrow\mathcal{K}_{r_{j}}$ induces a
map $\iota^{i,j}_{*}:H_{*}(\mathcal{K}_{r_{i}},\mathbb{F})\to
H_{*}(\mathcal{K}_{r_{j}},\mathbb{F})$ on homology, where $\mathbb{F}$ is a
field. For the remainder of this paper, we compute homology over the field
$\mathbb{F}=\mathbb{Z}/2\mathbb{Z}$. The _$q$ th-persistent homology_ (PH) is
the pair

$\Big{(}\\{H_{q}(\mathcal{K}_{r_{i}},\mathbb{F})\\}_{1\leq i\leq
k}\,,\\{\iota_{*}^{i,j}\\}_{1\leq i\leq j\leq k}\Big{)}\,.$

A homology class is _born_ at $r_{i}$ if it is not in the image of
$\iota_{*}^{i,i-1}$. The homology class _dies_ at $r_{j}>r_{i}$ if $j$ is the
minimum index such that $\iota_{*}^{i,j}$ maps the homology class to zero.
(Such a $j$ may not exist; in that case, the homology class never dies.) The
Fundamental Theorem of Persistent Homology yields compatible choices of bases
for the vector spaces $H_{q}(\mathcal{K}_{r_{i}},\mathbb{F})$. The generators
in our definition of a persistence diagram, below, are the basis elements in
the decomposition given by the Fundamental Theorem of Persistent Homology.

Persistent homology is often visualized as a _persistence diagram_ (PD). The
$q$th persistence diagram $PD_{q}(f)$ is a multiset of points in the extended
plane $\overline{\mathbb{R}}^{2}$ that summarizes the $q$th persistent
homology. It contains the points on the diagonal with infinite multiplicity
(for technical reasons) and one point for every generator. If a generator is
born at $b$ and dies at $d$, then the coordinates of the corresponding point
in the PD are $(b,d)$, and if the generator is born at $b$ and never dies,
then the coordinates of the point are $(b,\infty)$. Note that off-diagonal
points in a PD may also appear with multiplicity.

One of the standard methods for computing PH is the algorithm introduced in
[6], which we review here. Let $f:\mathcal{K}\to\mathbb{R}$ be a filtration
function, and let $\sigma_{1},\ldots,\sigma_{N}$ be the simplices of
$\mathcal{K}$, indexed such that $i<j$ if $\sigma_{i}$ is a proper face of
$\sigma_{j}$. The algorithm requires a choice of _compatible simplex indexing_
$\text{idx}:\mathcal{K}\to\\{1,\ldots,N\\}$, where $N$ is the number of
simplices in $\mathcal{K}$.

######  Definition 2.1.

A _simplex indexing_ is a bijection
$\text{idx}:\mathcal{K}\to\\{1,\ldots,N\\}$.

######  Definition 2.2.

A _compatible simplex indexing_ is a simplex indexing
$\text{idx}:\mathcal{K}\to\\{1,\ldots,N\\}$ such that
$\text{idx}(\sigma_{i})<\text{idx}(\sigma_{j})$ if
$f(\sigma_{i})<f(\sigma_{j})$ or $\sigma_{i}$ is a proper face of
$\sigma_{j}$.

A compatible simplex indexing exists because monotonicity ensures that
$f(\sigma_{i})\leq f(\sigma_{j})$ if $\sigma_{i}$ is a face of $\sigma_{j}$.
Because there may not be a unique compatible simplex indexing, we define the
_simplex indexing induced by $f$_ as follows.

######  Definition 2.3.

The _simplex indexing induced by $f$_ is the unique compatible simplex
indexing $\text{idx}_{f}:\mathcal{K}\to\\{1,\ldots,N\\}$ such that
$\text{idx}_{f}(\sigma_{i})<\text{idx}_{f}(\sigma_{j})$ if either
$f(\sigma_{i})<f(\sigma_{j})$ or if $f(\sigma_{i})=f(\sigma_{j})$ and $i<j$.

The function $\text{idx}_{f}$ is a compatible simplex indexing because the
sequence $\sigma_{1},\ldots,\sigma_{N}$ of simplices is ordered such that
$i<j$ if $\sigma_{i}$ is a proper face of $\sigma_{j}$.

Let $D$ be the boundary matrix compatible with the simplex indexing
$\text{idx}_{f}$. That is, let $D$ be the matrix whose $(i,j)$th entry is

$D_{ij}=\begin{cases}1\,,&\text{idx}_{f}^{-1}(i)\text{ is an
$(m-1)$-dimensional face of the $m$-dimensional simplex
}\text{idx}_{f}^{-1}(j)\\\ 0\,,&\text{otherwise.}\end{cases}$

We decompose the boundary matrix $D$ into a matrix product $D=RU$ such that
$U$ is upper triangular and $R$ is a binary matrix that is “reduced”. A binary
matrix $R$ is _reduced_ if $\text{low}_{R}(j)\neq\text{low}_{R}(j^{\prime})$
whenever $j\neq j^{\prime}$ are the indices of nonzero columns in $R$. The
quantity $\text{low}_{R}(j)$, which is called the _pairing function_ , is the
row index of the last $1$ in column $j$ if column $j$ is nonzero and undefined
if column $j$ is the zero vector. An RU decomposition can be computed in
$\mathcal{O}(N^{3})$ time [6, 5].

Cohen-Steiner et al. [3] showed that the pairing function $\text{low}_{R}(j)$
depends only on the boundary matrix $D$ and not on the particular reduced
binary matrix $R$ in the decomposition $D=RU$. A pair
$(\text{idx}_{f}^{-1}(i),\text{idx}_{f}^{-1}(j))$ of simplices with
$i=\text{low}_{R}(j)$ represents a persistent homology generator. The _birth
simplex_ $\text{idx}_{f}^{-1}(i)$ creates the homology class and the _death
simplex_ $\text{idx}_{f}^{-1}(j)$ destroys the homology class. The two
simplices in a pair have consecutive dimensions; that is, if
dim$(\text{idx}^{-1}(i))=q$, then dim$(\text{idx}^{-1}(j))=q+1$. If
dim$(\text{idx}_{f}^{-1}(i))=q$ and dim$(\text{idx}_{f}^{-1}(j))=q+1$, then a
point with coordinates $(f(\text{idx}_{f}^{-1}(i)),f(\text{idx}_{f}^{-1}(j)))$
is added to the $q$th persistence diagram. We refer to
$f(\text{idx}_{f}^{-1}(i))$ as its _birth_ and to $f(\text{idx}_{f}^{-1}(j))$
as its _death_. Some simplices are not paired. If $i\neq\text{low}_{R}(j)$ for
all $j$, then the simplex $\text{idx}_{f}^{-1}(i)$ is a birth simplex for a
homology class that never dies. Its birth is $f(\text{idx}_{f}^{-1}(i))$ and
its death is $\infty$. If $\dim(\text{idx}_{f}^{-1}(i))=q$, then a point with
coordinates $(f(\text{idx}_{f}^{-1}(i)),\infty)$ is added to the $q$th
persistence diagram.

###  2.2. Vineyards

Let $\mathcal{K}$ be a simplicial complex. A _$1$ -parameter filtration
function_ on $\mathcal{K}$ is a function $f:\mathcal{K}\times I\to\mathbb{R}$,
where $I=[t_{0},t_{1}]$ is an interval in $\mathbb{R}$, such that $f(\cdot,t)$
is a filtration function on $\mathcal{K}$ for all $t\in I$. For each $t\in I$,
the $r$-sublevel sets $\mathcal{K}_{r}^{t}=\\{\sigma\in\mathcal{K}\mid
f(\sigma,t)\leq r\\}$ are a filtration of $\mathcal{K}$. The set
$\\{\\{\mathcal{K}_{r}^{t}\\}_{r\in\mathbb{R}}\\}_{t\in I}$ is a set of
filtrations parameterized by $t\in I$. For each $t\in I$, one can compute the
persistence diagram $PD(f(\cdot,t))$. The associated _vineyard_ is the
1-parameter set $\\{PD(f(\cdot,t))\\}_{t\in I}$ of persistence diagrams. We
visualize the vineyard in $\overline{\mathbb{R}}^{2}\times I$ as a continuous
stack of PDs (see Figure 1). The points in the PDs trace out curves with time;
these curves are the _vines_.

An algorithm for computing vineyards is given by [3], and we review it here.
Let $f:\mathcal{K}\times I\to\mathbb{R}$ be a $1$-parameter filtration
function, and let $\sigma_{1},\ldots,\sigma_{N}$ be the simplices of
$\mathcal{K}$, indexed such that $i<j$ if $\sigma_{i}$ is a proper face of
$\sigma_{j}$. As in Section 2.1, we fix a compatible simplex indexing
$\text{idx}_{f}:\mathcal{K}\times I\to\\{1,\ldots,N\\}$ such that
$\text{idx}_{f}(\sigma_{i},t)<\text{idx}_{f}(\sigma_{j},t)$ if we have either
$f(\sigma_{i},t)<f(\sigma_{j},t)$ or we have $f(\sigma_{i},t)=f(\sigma_{j},t)$
and $i<j$. The simplex indexing can only change at times $t_{*}$ at which
$f(\sigma,t_{*})=f(\tau,t_{*})$ for some $\sigma,\tau\in\mathcal{K}$. At
$t_{*}$, the indices of $\sigma$ and $\tau$ may change. (If $\sigma$, $\tau$
are the unique pair of simplices whose indices change at $t_{*}$, then
$\sigma$ and $\tau$ are transposed in the simplex indexing and $\sigma$ and
$\tau$ have consecutive indices in the indexing. Otherwise, there is a
sequence of such transpositions.) Let $D(t)$ denote the boundary matrix
compatible with the indexing at time $t$ and let $\text{low}_{R}(\cdot,t)$
denote the corresponding pairing function. One computes $D(t_{0})$ at the
initial time $t_{0}$ and an RU decomposition $D(t_{0})=R(t_{0})U(t_{0})$. The
initial pairing function $\text{low}_{R}(\cdot,t_{0})$ is read off from the
initial RU decomposition. Then we sweep through the intersections $t_{*}$
(from left to right) at which the simplex indexing changes. At each $t_{*}$,
we update the simplex indexing, RU decomposition, and pairing function. This
updating procedure yields the birth and death simplices for the filtration
function $f(\cdot,t_{*})$, which one can use to obtain the new persistence
diagram.

The procedure above is an efficient way of computing the diagrams
PD$(f(\cdot,t))$ for all $t\in[t_{0},t_{1}]$. At worst, updating $R(t_{*})$
requires adding one column to another and adding one row to another—similarly
for $U(t_{*})$. The addition of columns and rows is an $\mathcal{O}(N)$
operation, although in experiments, the authors of [3] found that updating
$R(t_{*})$ and $U(t_{*})$ can be done in approximately constant time if one
uses the sparse matrix representations that are given in [3]. If there is a
single transposition of simplices at $t_{*}$, then at most two (birth, death)
simplex pairs are updated, and these updates occur in $\mathcal{O}(1)$ time.

A special type of vineyard is a _piecewise-linear vineyard_. If we are only
given $f(\sigma,t_{i})$ at discrete time steps $t_{i}$, then for all $i$ we
extend $f(\sigma,t)$ to $t\in[t_{i},t_{i+1}]$ by linear interpolation. In this
case, one can compute the transposition times $t_{*}$ by using the
Bentley–Ottman planesweep algorithm [4]. This is because computing when (if)
two simplices $\sigma,\tau$ get transposed in $[t_{i},t_{i+1}]$ is equivalent
to finding the intersection (if it exists) between the line segments

$\displaystyle
y=\frac{f(\sigma,t_{i+1})-f(\sigma,t_{i})}{t_{i+1}-t_{i}}(t-t_{i})+f(\sigma,t_{i})\,,$
$\displaystyle
y=\frac{f(\tau,t_{i+1})-f(\tau,t_{i})}{t_{i+1}-t_{i}}(t-t_{i})+f(\tau,t_{i})\,.$

###  2.3. PD bundles

PD bundles were introduced in [7] as a generalization of vineyards in which a
set of filtrations is parameterized by an arbitrary topological space $B$. A
vineyard is the special case in which $B$ is an interval in $\mathbb{R}$.

######  Definition 2.4.

A _fibered filtration function_ is a set
$\\{f_{p}:\mathcal{K}^{p}\to\mathbb{R}\\}_{p\in B}$ of filtration functions
parameterized by a topological space $B$ (the _base space_). When
$\mathcal{K}^{p}\equiv\mathcal{K}$ for all $p\in B$, we define
$f(\sigma,p):=f_{p}(\sigma)$ for all simplices $\sigma\in\mathcal{K}$ and
$p\in B$.

######  Definition 2.5.

Let $\\{f_{p}:\mathcal{K}^{p}\to\mathbb{R}\\}_{p\in B}$ be a fibered
filtration function. The topological space $B$ is the _base space_. The space
$E:=\\{(p,z)\mid z\in PD_{q}(f_{p})\,,p\in B\\}$ is the $q$th _total space_.
We give $E$ the subspace topology inherited from the inclusion
$E\hookrightarrow B\times\overline{\mathbb{R}}^{2}$. The associated _$q$ th PD
bundle_ is the triple $(E,B,\pi)$, where $\pi$ is the projection from $E$ to
$B$.

In [3], it was computationally easier to work with a piecewise-linear
vineyard, which is a vineyard for a fibered filtration function of the form
$f:\mathcal{K}\times[t_{0},t_{1}]\to\mathbb{R}$ such that $f(\sigma,\cdot)$ is
piecewise linear for all $\sigma\in\mathcal{K}$. (See the discussion at the
end of Section 2.2.) Below, we define an analog of piecewise-linear vineyards.

######  Definition 2.6 (Piecewise-linear PD bundles).

Let $\\{f_{p}:\mathcal{K}^{p}\to\mathbb{R}\\}_{p\in B}$ be a fibered
filtration function in which $\mathcal{K}^{p}\equiv\mathcal{K}$. As in
Definition 2.4, we define $f(\sigma,p):=f_{p}(\sigma)$ for all
$\sigma\in\mathcal{K}$ and $p\in B$. If $B$ is a simplicial complex and
$f(\sigma,\cdot)$ is linear on each simplex of $B$ for all simplices
$\sigma\in\mathcal{K}$, then $f$ is a _piecewise-linear fibered filtration
function._ The resulting PD bundle is a _piecewise-linear PD bundle_.

For the remainder of the paper, we only consider piecewise-linear PD bundles.

For example, in Section 1 we considered a point cloud $X(t,\mu)$ whose
coordinates depended on time $t\in\mathbb{R}$ and parameter
$\mu\in\mathbb{R}$. Given only the coordinates of the point cloud at a
discrete set $\\{t_{i}\\}$ of time steps and a discrete set $\\{\mu_{j}\\}$ of
parameter values, we obtained a filtration function $f_{(t_{i},\mu_{j})}$ for
every $(t_{i},\mu_{j})$. We extended this to a piecewise-linear fibered
filtration function on $B=[\min t_{i},\max
t_{i}]\times[\min\mu_{j},\max\mu_{j}]$ via linear interpolation of the
filtration values for each simplex $\sigma\in\mathcal{K}$.

More generally, suppose that we are given a fibered filtration function of the
form $f:\mathcal{K}\times\prod_{i=1}^{m}\mathcal{I}_{i}\to\mathbb{R}$, where
each $\mathcal{I}_{i}$ is a finite subset of $\mathbb{R}$, and we wish to
extend $f$ to a fibered filtration function whose base is
$B:=\prod_{i=1}^{m}[\min\mathcal{I}_{i},\max\mathcal{I}_{i}]$. First, we
construct a triangulation $B$ (i.e., an $m$-dimensional simplicial complex) of
$\prod_{i=1}^{m}[\min\mathcal{I}_{i},\max\mathcal{I}_{i}]$ whose set of
vertices is $\prod_{i=1}^{m}\mathcal{I}_{i}$. (See [10], for example, for a
method to triangulate a cubical complex.) We then extend $f$ to a piecewise-
linear fibered filtration function $f:\mathcal{K}\times B\to\mathbb{R}$ by
linearly interpolating $f(\sigma,\cdot)$ on each simplex $\Delta\in B$ for all
simplices $\sigma\in\mathcal{K}$.

In [7], it was shown that if $f:\mathcal{K}\times B\to\mathbb{R}$ is a
piecewise-linear fibered filtration function on an $n$-dimensional simplicial
complex $B$, then $B$ can be subdivided into $n$-dimensional polyhedra such
that within each polyhedron $P$, there is a “template” from which
$PD_{q}(f(\cdot,p))$ can be computed for any $p\in P$. The template is a list
of (birth, death) simplex pairs $(\sigma_{b},\sigma_{d})$.

The polyhedra are defined as follows. We define

(1) $I(\sigma,\tau):=\\{p\in B\mid f(\sigma,p)=f(\tau,p)\\}\,.$

For every $n$-simplex $\Delta$ in $B$, the intersection
$I(\sigma,\tau)\cap\Delta$ is $\emptyset$, $\Delta$, a vertex of $\Delta$, or
the intersection of an $(n-1)$-dimensional hyperplane with $\Delta$. The set

(2) $\bigcup_{\Delta\in
B}\partial\Delta\cup\Big{\\{}\Big{(}I(\sigma,\tau)\cap\Delta\Big{)}\,\Big{|}\,\emptyset\subset\Big{(}I(\sigma,\tau)\cap\Delta\Big{)}\subset\Delta\Big{\\}}$

determines the boundaries of a polyhedral decomposition of $B$, where $\Delta$
denotes an $n$-simplex of $B$ and $\partial\Delta$ denotes the boundary of
$\Delta$. Every polyhedron face is either a subset of $\partial\Delta$ or a
subset of $I(\sigma,\tau)\cap\Delta$ for some simplex $\Delta$ of $B$.

As in Sections 2.1 and 2.2, we define the _simplex indexing induced by $f$_ as
follows. Let $\sigma_{1},\ldots,\sigma_{N}$ be the simplices of $\mathcal{K}$,
indexed such that $i<j$ if $\sigma_{i}$ is a proper face of $\sigma_{j}$. We
define the simplex indexing function $\text{idx}_{f}:\mathcal{K}\times
B\to\\{1,\ldots,N\\}$ to be the unique function such that
$\text{idx}_{f}(\sigma_{i},p)<\text{idx}_{f}(\sigma_{j},p)$ if we have
$f(\sigma_{i},p)<f(\sigma_{j},p)$ or we have $f(\sigma_{i},p)=f(\sigma_{j},p)$
and $i<j$.

######  Proposition 2.7 ([7]).

If $f:\mathcal{K}\times B\to\mathbb{R}$ is a piecewise-linear fibered
filtration function, then the set in Equation 2 subdivides $B$ into polyhedra
$P$ on which the simplex indexing is constant (i.e.,
$\text{idx}_{f}(\sigma,\cdot)|_{P}$ is constant for all
$\sigma\in\mathcal{K}$).

Proposition 2.7 implies that within each polyhedron $P$, the set
$\\{(\sigma_{b},\sigma_{d})\\}$ of (birth, death) simplex pairs for
$f(\cdot,p)$ is constant with respect to $p$.

##  3\. Computing piecewise-linear PD bundles

The algorithm to compute piecewise-linear PD bundles has three main parts.

  1. (1)

Compute the polyhedra: First, we compute the polyhedra on which the simplex
indexing (and thus pairing function) is constant (see Proposition 2.7). For
every pair $P_{1}$, $P_{2}$ of adjacent polyhedra, the face $Q$ that they
share is a subset of a set of the form
$I(\sigma_{i_{1}},\sigma_{j_{1}})\cap\cdots\cap
I(\sigma_{i_{m}},\sigma_{j_{m}})$ for some $m$, where
$I(\sigma_{i_{1}},\sigma_{j_{1}})$ is defined by Equation 1. In the “generic
case,” defined below at the beginning of Section 3.1.1, we have $m=1$. We
compute and store a reference that is associated with $Q$ to the set
$\\{(\sigma_{i_{1}},\sigma_{j_{1}}),\ldots,(\sigma_{i_{m}},\sigma_{j_{m}})\\}_{k=1}^{m}$.

  2. (2)

Compute the pairing function: Within each polyhedron $P$, the set of (birth,
death) simplex pairs for $f(\cdot,p)$ is constant with respect to $p\in P$. We
compute the set of (birth, death) simplex pairs for each $P$. First, we choose
an initial point $p_{*}\in B$ at which we compute the simplex indexing at
$p_{*}$, the boundary matrix $D(p_{*})$, an RU decomposition
$D(p_{*})=R(p_{*})U(p_{*})$, and the pairing function at $p_{*}$. (This is the
pairing function for the entire polyhedron that contains $p_{*}$.) These
computations take $\mathcal{O}(N^{3})$ time. We then traverse the polyhedra,
starting with the polyhedron that contains $p_{*}$ and visiting each
polyhedron at least once. As we move from one polyhedron to the next via a
shared face $Q$, we use the set
$\\{(\sigma_{i_{1}},\sigma_{j_{1}}),\ldots,(\sigma_{i_{m}},\sigma_{j_{m}})\\}_{k=1}^{m}$
computed earlier for $Q$ to update the RU decomposition and pairing function
via the update rules that are used when computing vineyards (see [3]). For
each polyhedron, we store its pairing function (i.e., the pairs
$(\sigma_{b},\sigma_{d})$ of birth and death simplex pairs and also the
unpaired simplices $\sigma_{b}$, which are birth simplices for homology
classes that never die).

  3. (3)

Query the PD bundle: To see the $q$th persistence diagram
$\text{PD}_{q}(f(\cdot,p))$ associated with point $p\in B$, we first locate
the polyhedron $P$ that contains $p$. For each pair $(\sigma_{b},\sigma_{d})$
of simplices in the pairing function for $P$, the diagram
$\text{PD}_{q}(f(\cdot,p))$ has a point with coordinates
$(f(\sigma_{b},p),f(\sigma_{d},p))$ if $\dim(\sigma_{b})=q$. For every
$q$-dimensional simplex $\sigma_{b}$ that is unpaired in $P$, the diagram
$\text{PD}_{q}(f(\cdot,p))$ contains the point $(f(\sigma_{b},p),\infty)$.

In what follows, I elaborate on each step of the algorithm above. We focus on
the case in which $B$ is a triangulated surface.

###  3.1. Special case: $B$ is a triangulated surface

Let $\mathcal{K}$ be a simplicial complex with simplices
$\sigma_{1},\ldots,\sigma_{N}$, indexed such that $i<j$ if $\sigma_{i}$ is a
proper face of $\sigma_{j}$. Let $f:\mathcal{K}\times B\to\mathbb{R}$ be a
piecewise-linear fibered filtration function, and suppose that $B$ is a
triangulated surface. If $\Delta$ is a triangle in $B$, then
$I(\sigma,\tau)\cap\Delta$ is either $\emptyset$, $\Delta$, a vertex of
$\Delta$, or a line segment whose endpoints are on $\partial\Delta$. Figure 2
shows a few possible cases for $I(\sigma,\tau)$. The set in Equation 2 is a
set $L$ of line segments, and the polygonal subdivision induced by $L$ is a
_line arrangement_ 111Usually, a “line arrangement” refers to a planar
subdivision that is induced by a set of lines in the plane. However, here I am
using the term “line arrangement” slightly more generally. Every line segment
in $L$ lies in some triangle $\Delta$ of $B$, and the line segments subdivide
each triangle into polygons. $\mathcal{A}(L)$. For an example of a line
arrangement, see Figure 3.

(a)

(b)

(c)

(d)

(e)

(f)

Figure 2. A few possible cases for the set $I(\sigma,\tau)$, which is defined
in Equation 1. The black lines are the 1-skeleton of $B$. Figure 3. A line
arrangement $\mathcal{A}(L)$ that represents the subdivision of a triangulated
base space $B$ into polygons (see Proposition 2.7). Within each polygon, the
simplex indexing is constant. (This figure appeared originally in [7].)

To simplify the exposition, we will make two genericity assumptions for the
remainder of Section 3.1. The idea of the algorithm is not different in the
general case, but it requires some technical modifications, which we discuss
in Section A.1. The assumptions are as follows:

  1. (1)

For all distinct simplices $\sigma,\tau\in\mathcal{K}$ and all vertices $v\in
B$, we have that $f(\sigma,v)\neq f(\tau,v)$. This implies that for all
triangles $\Delta\in B$, the intersection $I(\sigma,\tau)\cap\Delta$ is either
$\emptyset$ or a line segment whose endpoints are not vertices of $\Delta$.
See the examples in Figures 2a and 2b.

  2. (2)

For all distinct simplices
$\sigma_{i_{1}},\sigma_{j_{1}},\sigma_{i_{2}},\sigma_{j_{2}}\in\mathcal{K}$
and every triangle $\Delta\in B$ such that
$I(\sigma_{i_{1}},\sigma_{j_{1}})\cap\Delta$ and
$I(\sigma_{i_{2}},\sigma_{j_{2}})\cap\Delta$ are nonempty, the line segments
$I(\sigma_{i_{1}},\sigma_{j_{1}})\cap\Delta$ and
$I(\sigma_{i_{2}},\sigma_{j_{2}})\cap\Delta$ do not share any endpoints.

Throughout Section 3.1, we use the following notation. Let $m$ denote the
number of triangles in $B$. The numbers of vertices and edges in $B$ are both
$\mathcal{O}(m)$. Let $N$ denote the number of simplices in $\mathcal{K}$. For
every triangle $\Delta$ of $B$, let $\kappa_{\Delta}$ denote the number of
line segments $\ell$ of the form $I(\sigma,\tau)\cap\Delta$ and let
$\kappa:=\sum_{\Delta\in B}\kappa_{\Delta}$. The worst case is
$\kappa_{\Delta}=\mathcal{O}(N^{2})$ and $\kappa=\mathcal{O}(mN^{2})$, but
these are very crude upper bounds. Let $\mu_{\Delta}$ denote the number of
vertices of $\mathcal{A}(L)$ in the interior of triangle $\Delta\in B$; the
quantity $\mu_{\Delta}$ is equal to the number of points of the form
$I(\sigma_{i_{1}},\sigma_{j_{1}})\cap
I(\sigma_{i_{2}},\sigma_{j_{2}})\cap\Delta$. Let $\mu$ denote the total number
of vertices in $\mathcal{A}(L)$, which is $\sum_{\Delta\in
B}\mu_{\Delta}+\mathcal{O}(m+\kappa)$. In the worst case,
$\mu_{\Delta}=\mathcal{O}(\kappa_{\Delta}^{2})=\mathcal{O}(N^{4})$ and
$\mu=\mathcal{O}(mN^{4})$, but these are again very crude upper bounds. The
numbers of edges and polygons in $\mathcal{A}(L)$ are
$\mathcal{O}(\mu)$.222Within each triangle $\Delta$ of $B$, the numbers of
edges and polygons in $\mathcal{A}(L)\cap\Delta$ are
$\mathcal{O}(\mu_{\Delta})$ by Euler’s formula, as noted in [4].

####  3.1.1. Computing the polygons

For a piecewise-linear vineyard, computing the intervals on which the simplex
indexing is constant can be reduced to finding the intersections between the
piecewise-linear functions $f(\sigma,t)$ and $f(\tau,t)$ for all pairs
$(\sigma,\tau)$ of simplices in $\mathcal{K}$. Likewise, for a piecewise-
linear PD bundle, computing the polygons on which the simplex indexing is
constant can be reduced to finding the intersections $I(\sigma,\tau)$ for all
pairs $(\sigma,\tau)$ of simplices.

We seek to compute the line arrangement $\mathcal{A}(L)$, where $L$ is the set
of line segments defined by Equation 2. (See Figure 3.) The polygons of
$\mathcal{A}(L)$ are the polygons on which the simplex indexing is constant.
We store $\mathcal{A}(L)$ using a doubly-connected edge list (DCEL) data
structure, which is a standard data structure for storing a polygonal
subdivision of the plane [4]. The DCEL data structure can be used without
modification to represent $\mathcal{A}(L)$, which is a polygonal subdivision
of a triangulated surface.333The primary reason to consider triangulated
surfaces, rather than any simplicial complex $B$ such that $\dim(B)\leq 2$, is
so that we can use a DCEL data structure to represent $\mathcal{A}(L)$.
Otherwise, what follows in Section 3.1 works just as well for any $2$D
simplicial complex. The space complexity of $\mathcal{A}(L)$ is
$\mathcal{O}(\mu)$. We compute $\mathcal{A}(L)$ using the following algorithm
(illustrated in Figure 4):

(a)

(b)

(c)

(d)

(e)

(f)

(g)

(h)

Figure 4. Computing the polygons. (A) The line arrangement $\mathcal{A}(L)$ is
initialized to represent the triangulated base space $B$, which in this case
consists of two triangles. (B) We find the vertices $v$ of $\mathcal{A}(L)$
that lie on the $1$-skeleton of $B$. (C)–(H) For each triangle $\Delta$ of
$B$, we incrementally add the line segments of the form
$I(\sigma,\tau)\cap\Delta$. The endpoints of a given line segment are a pair
$(v,w)$ of vertices in (B). In (D), an internal vertex (a vertex at the
intersection of two line segments) is created when the second line segment is
added. Three of the line segments intersect at the internal vertex.

  1. (1)

We initialize $\mathcal{A}(L)$ so that it represents the triangulation $B$.
(See Figure 4a.) In addition to the usual data that a DCEL stores, we
enumerate the triangles in $B$ and every half edge444An “edge” is a line
segment in $\mathcal{A}(L)$ that connects two vertices $u$ and $v$. One can
think of each edge as two _half edges_ , which are represented as the two
oriented edges $u\to v$ and $v\to u$. If an edge is adjacent to two faces
$F_{1}$ and $F_{2}$, then one of the half edges is associated with $F_{1}$ and
the other is associated with $F_{2}$. One chooses an orientation (clockwise or
counterclockwise) so that for every face $F$ in the DCEL, the half-edges
associated with $F$ have this orientation (clockwise or counterclockwise)
around the boundary of $F$. $e$ stores the index of the triangle in $B$ that
$e$ is on the boundary of.

  2. (2)

For every triangle $\Delta\in B$, we initialize an empty dictionary
$\mathcal{D}(\Delta)$.555A dictionary is a data structure for storing (key,
value) pairs. The keys will be pairs $(\sigma,\tau)$ of simplices for which
$I(\sigma,\tau)\cap\Delta$ is a line segment. The value of $(\sigma,\tau)$
will be a list of the endpoints $(v,w)$ of the line segment. We denote the
value of $(\sigma,\tau)$ by $\mathcal{D}(\Delta)[(\sigma,\tau)]$. These
dictionaries use $\mathcal{O}(\kappa)$ space.

  3. (3)

For each edge $e$ of $B$, we compute the vertices of $\mathcal{A}(L)$ that lie
on $e$. (See Figure 4b.) These are the vertices that equal $I(\sigma,\tau)\cap
e$ for some pair $(\sigma,\tau)$ of simplices. To do this, we consider the
restriction of $f$ to $e$; the restriction $f|_{e}$ is a $1$-parameter
filtration function (the input to a vineyard). For each
$\sigma\in\mathcal{K}$, the set $\\{(p,f(\sigma,p))\mid p\in e\\}$ is a line
segment $\ell_{\sigma,e}$ and $I(\sigma,\tau)\cap e$ is the point $p\in e$ at
which $\ell_{\sigma,e}$ and $\ell_{\tau,e}$ intersect. We use the
Bentley–Ottman planesweep algorithm [4] to compute these intersections, thus
obtaining the vertices $v$ of $\mathcal{A}(L)$ that lie on $e$. For a vertex
$v$ that equals $I(\sigma,\tau)\cap e$, we add $v$ to the list
$\mathcal{D}(\Delta)[(\sigma,\tau)]$ for each triangle $\Delta$ to which $e$
is adjacent. Completing step 3 takes $\mathcal{O}(N)$ space and
$\mathcal{O}((mN+\kappa)\log N)$ time in total for all edges in $B$, and can
be parallelized over the edges.666In some cases, in may be more efficient to
consider the restriction of $f$ to an Euler path $\gamma$ through the
$1$-skeleton of $B$, rather than the restriction of $f$ to each edge
separately. For example, if $B$ is of the form in Figure 5a, then an Euler
path is given in Figure 5b.

(a)

(b)

Figure 5. (A) A triangulated base space $B$. (B) An Euler path $\gamma$
through the $1$-skeleton of $B$, starting at the bottom-left vertical edge
(violet) and ending at the top-right vertical edge (red).

  4. (4)

For each triangle $\Delta\in B$ and for each key $(\sigma,\tau)$ in the
dictionary $\mathcal{D}(\Delta)$, there is an associated pair $(v,w)$ of
vertices that are the endpoints of the line segment
$I(\sigma,\tau)\cap\Delta$. We seek to add all of these line segments to the
DCEL that represents $\mathcal{A}(L)$. (See Figures 4C–H.) There are many
standard algorithms for doing this; one example is the incremental algorithm
(see, e.g., Chapter 8.3 of [4]), in which the line segments are incrementally
added one at a time. The worst-case run time of the incremental algorithm in
triangle $\Delta$ is $\mathcal{O}(\kappa_{\Delta}^{2})$, which yields a total
run time of $\mathcal{O}(\sum_{\Delta\in B}\kappa_{\Delta}^{2})$. For better
performance, this algorithm can be parallelized over the triangles $\Delta\in
B$.

In Figure 4, we illustrate the algorithm for computing the polygons.

As we add line segments to $\mathcal{A}(L)$, we keep track of the pairs
$(\sigma,\tau)$ of simplices that correspond to each edge of $\mathcal{A}(L)$.
If edge $e$ of $\mathcal{A}(L)$ is a subset of $I(\sigma,\tau)$, then $e$
stores a reference to the pair $(\sigma,\tau)$. We add the reference to
$(\sigma,\tau)$ at the time that edge $e$ is created in $\mathcal{A}(L)$. If
$P_{1}$, $P_{2}$ are adjacent polygons of $\mathcal{A}(L)$ that share edge
$e$, then the simplex indexings in $P_{1}$, $P_{2}$ are related via the
transposition of $\sigma$ and $\tau$.

####  3.1.2. Computing the pairing function

Let $G$ be the dual graph of the line arrangement $\mathcal{A}(L)$. The graph
$G$ has a vertex $v_{P}$ for every polygon $P$ of $\mathcal{A}(L)$ and an edge
between two vertices if the corresponding polygons are adjacent. We compute a
path $\Gamma$ that visits every vertex of $G$ at least once. For an example,
see Figure 6. One way to obtain such a path is via depth-first search, which
takes $\mathcal{O}(\mu)$ time because the number of nodes in $G$ (i.e.,
polygons of $\mathcal{A}(L)$) is $\mathcal{O}(\mu)$. This yields a path
$\Gamma$ whose length is $\mathcal{O}(\mu)$.

Figure 6. A path $\Gamma$ (in white) that visits every polygon in the line
arrangement $\mathcal{A}(L)$.

At the first vertex $v_{P}$ of $\Gamma$, we compute the simplex indexing in
polygon $P$, the RU decomposition of the boundary matrix in $P$, and the
(birth, death) simplex pairs in $P$. To store the RU decomposition, we use the
sparse matrix data structure from [3]. The polygon $P$ stores a reference to
its (birth, death) simplex pairs. To store the current simplex indexing, each
simplex stores a reference to its index in the current indexing (which we
initialize to the indexing in $P$).

We traverse the path $\Gamma$. As we walk from one polygon $P_{1}$ to the next
polygon $P_{2}$ by crossing an edge $e$ in $\mathcal{A}(L)$, we update the
simplex indexing, the RU decomposition, and the (birth, death) simplex pairs.
To update the simplex indexing, recall that edge $e$ stores a reference to the
simplex pair $(\sigma,\tau)$ such that $e\subseteq I(\sigma,\tau)$. This
implies that the simplex indexings in $P_{1}$ and $P_{2}$ are related via the
transposition of $\sigma$ and $\tau$ because we must have (without loss of
generality) $f(\sigma,p)>f(\tau,p)$ for $p\in P_{1}$ and
$f(\sigma,p)<f(\tau,p)$ for $p\in P_{2}$, with $f(\sigma,p)=f(\tau,p)$ on the
shared edge $e$. We update the simplex indexing by swapping the indices that
$\sigma$ and $\tau$ store. To update the RU decomposition and the (birth,
death) simplex pairs, we apply the update algorithm of [3]. This takes
$\mathcal{O}(N)$ time in the worse case, but often approximately constant time
in practice; see [3] and the earlier discussion in Section 2.2. In $P_{2}$, we
store the new (birth, death) simplex pairs. Storing the simplex pairs uses
$\mathcal{O}(N)$ space for each polygon of $\mathcal{A}(L)$, so we use
$\mathcal{O}(N\mu)$ space in total.

We can optimize the space requirements by recalling from [3] that most updates
of the simplex indexing do not change the (birth, death) simplex
pairs.777However, we note that we do not know in advance whether an update of
the simplex indexing will change the (birth, death) simplex pairs. If the
update from $P_{1}$ to $P_{2}$ does not change the simplex pairs, we can
delete the edge in $\mathcal{A}(L)$ that $P_{1}$ and $P_{2}$ share, thus
merging $P_{1}$ and $P_{2}$ into a single polygon and reducing the size of
$\mathcal{A}(L)$.

####  3.1.3. Querying the PD bundle

We consider the scenario in which a user seeks to query many points $p\in B$
in real time and see the $q$th persistence diagram $\text{PD}_{q}(f(\cdot,p))$
associated with each queried point $p$.

To compute the $q$th persistence diagram $\text{PD}_{q}(f(\cdot,p))$
associated with a given $p$, we first identify the polygon $P$ of
$\mathcal{A}(L)$ that contains $p$. We do this in two steps. The first step is
to identify the triangle $\Delta\in B$ that contains $p$. This takes
$\mathcal{O}(m)$ worst-case time because it takes constant time to test if a
given triangle contains $p$. In certain cases, one can identify the triangle
$\Delta$ more efficiently. For example, if $B$ is a triangulation of the form
in Figure 5a, then one can locate the triangle $\Delta$ in $\mathcal{O}(1)$
time by simply examining the coordinates of $p$. The second step is to locate
the polygon $P$ in $\Delta$ that contains $p$. This is a well-studied problem
in computational geometry; it is known as the _point-location problem_ [14].
When one is planning to perform many point-location queries on the same line
arrangement (i.e., if one is querying many points $p\in B$), the standard
strategy is to precompute a data structure so that the subsequent point-
location queries can be done efficiently. There are many strategies for doing
this (see, e.g., Chapter 38 in [14]). One method is the slab-and-persistence
method [13], in which one precomputes a “persistent search tree” for the line
arrangement.888If $B$ is a 2D triangulated subset of $\mathbb{R}^{2}$, then
one can also use the slab-and-persistence method to perform the first
step—locating the triangle $\Delta\in B$ that contains $p$. We construct a
persistent search tree for each triangle in $B$. Using separate persistent
search trees for the planar subdivisions in each triangle, the slab and
persistence method takes $\mathcal{O}(\sum_{\Delta\in
B}\mu_{\Delta}\log(\mu_{\Delta}))$ preprocessing time, $\mathcal{O}(\mu)$
space, and $\mathcal{O}(\max_{\Delta\in B}\log\mu_{\Delta})$ time per query.

We obtain $\text{PD}_{q}(f(\cdot,p))$ by evaluating $f(\cdot,p)$ on the
simplices in the pairing function for polygon $P$, which was precomputed in
the previous step (see Section 3.1.2). This takes $\mathcal{O}(N)$ time. For
every (birth, death) pair $(\sigma_{b},\sigma_{d})$ of simplices,
$\text{PD}_{q}(f(\cdot,p))$ has a point with coordinates
$(f(\sigma_{b},p),f(\sigma_{d},p))$ if $\dim(\sigma_{b})=q$. For every
unpaired $q$-dimensional simplex $\sigma_{b}$, the diagram
$\text{PD}_{q}(f(\cdot,p))$ has a point with coordinates
$(f(\sigma_{b},p),\infty)$.

###  3.2. Generalizing to higher-dimensional base spaces

The algorithm of Section 3.1 (as outlined at the beginning of Section 3) does
not require many modifications for higher-dimensional base spaces $B$. We
replace the subdivision of $B$ into polygons by a subdivision of $B$ into
$n$-dimensional polyhedra, where $n=\dim(B)$. In the third step (querying the
PD bundle), one uses a point-location algorithm for a hyperplane arrangement
(see, e.g., [11, 1]). Only the first step (computing the polyhedra) requires a
meaningful modification, which I describe below.

When $n=2$, the intersection of $I(\sigma,\tau)$ with a triangle $\Delta\in B$
is the intersection of a line with $\Delta$, which is a line segment
$L_{\sigma,\tau,\Delta}$. These line segments completely determine the
polygonal subdivision of $B$ because the line segments are the faces of the
polygons. In turn, each line segment $L_{\sigma,\tau,\Delta}$ is completely
determined by the intersection of $L_{\sigma,\tau,\Delta}$ with the
$1$-skeleton of $B$; this intersection is a pair
$(v_{\sigma,\tau,\Delta},w_{\sigma,\tau,\Delta})$ of points. In Section 3.1.1,
we computed the set
{$(v_{\sigma,\tau,\Delta},w_{\sigma,\tau,\Delta})\\}_{\sigma,\tau,\Delta}$ by
restricting the fibered filtration function $f$ to each edge of $B$ and
applying the Bentley–Ottman planesweep algorithm.

In general, the intersection of $I(\sigma,\tau)$ with an $n$-simplex
$\Delta\in B$ is the intersection of an $(n-1)$-dimensional hyperplane
$H_{\sigma,\tau,\Delta}$ with $\Delta$; the intersection is an
$(n-1)$-dimensional polyhedron $P_{\sigma,\tau,\Delta}$. The set
$\\{P_{\sigma,\tau,\Delta}\\}_{\sigma,\tau,\Delta}$ completely determines the
polyhedral subdivision of $B$ that is given by Proposition 2.7 because the
polyhedra $P_{\sigma,\tau,\Delta}$ are the $(n-1)$-dimensional faces of the
$n$-dimensional polyhedra in the subdivision. In turn, each polyhedron
$P_{\sigma,\tau,\Delta}$ is completely determined by its intersection with the
edges of $B$, as follows. The $m$-dimensional faces of
$P_{\sigma,\tau,\Delta}$ are the set

$\\{H_{\sigma,\tau,\Delta}\cap\Delta^{(m+1)}\mid\Delta^{(m+1)}\text{ is an
}(m+1)\text{-dimensional face of }\Delta\text{ and
}H_{\sigma,\tau,\Delta}\cap\Delta^{(m+1)}\neq\emptyset\\}\,.$

For every $(m+1)$-dimensional face $\Delta^{(m+1)}$ of $\Delta$ such that
$H_{\sigma,\tau,\Delta}\cap\Delta^{(m+1)}\neq\emptyset$, the intersection
$H_{\sigma,\tau,\Delta}\cap\Delta^{(m+1)}$ is the $m$-dimensional polyhedron
whose $(m-1)$-dimensional faces are the set

$\\{H_{\sigma,\tau,\Delta}\cap\Delta^{(m)}\mid\Delta^{(m)}\text{ is an
}m\text{-dimensional face of }\Delta^{(m+1)}\\}\,.$

By induction, the faces of $P_{\sigma,\tau,\Delta}$ are determined by

$\\{H_{\sigma,\tau,\Delta}\cap e\mid e\text{ is a }1\text{-dimensional face of
}\Delta\text{ (i.e., $e$ is an edge)}\\}\,,$

which are the vertices of $P_{\sigma,\tau,\Delta}$. Consequently, we can
compute each $P_{\sigma,\tau,\Delta}$ by determining its vertices. As in the
case in which $B$ is a triangulated surface, we do this by restricting the
fibered filtration function $f$ to each edge of $B$ and applying the
Bentley–Ottman planesweep algorithm.

##  4\. Conclusions and Discussion

I introduced an algorithm for efficiently computing persistence diagram (PD)
bundles when the fibered filtration function is piecewise linear. I gave full
details for the case in which the base space $B$ is a triangulated surface.
Additionally, in Section 3.2, I discussed how one can generalize the algorithm
to higher dimensions. I conclude with some questions and proposals for future
work:

  * •

What invariants can we use to summarize and analyze PD bundles in ways that do
not require exploratory data analysis?

The current algorithm asks a user to “query” the PD bundle at various points
in the base space. This is useful for qualitative analysis or if one has a
function whose input is a PD and one seeks to optimize that function over $B$.
However, other applications may require global invariants.

  * •

How do we generalize the algorithm to fibered filtration functions that are
not piecewise linear?

For piecewise-linear fibered filtration functions, we used the fact that the
base space $B$ can be subdivided into polyhedrons such that there is a single
PD “template” (a list of (birth, death) simplex pairs) for each polyhedron.
The template can then be used to obtain $\text{PD}_{q}(f_{p})$ at any point
$p$ in the polyhedron. For “generic” fibered filtration functions, it was
shown in [7] that the base space $B$ is stratified such that for each stratum,
there is a single PD template that can be used to obtain
$\text{PD}_{q}(f_{p})$ at any point $p$ in the stratum.

  * •

How do we generalize to the case where $\mathcal{K}^{p}$ is not constant with
respect to $p\in B$?

In this case, simplices are added and removed from the filtration as $p\in B$
varies, so the algorithm must be modified.

## Acknowledgements

I thank Michael Lesnick and Nina Otter for helpful discussions about computing
PD bundles.

## Appendix

###  A.1. Technical details of the algorithm

Let $\mathcal{K}$ be a simplicial complex with $N$ simplices, indexed
$\sigma_{1},\ldots,\sigma_{N}$ such that $i<j$ if $\sigma_{i}$ is a proper
face of $\sigma_{j}$. Let $B$ be a triangulated surface, and let
$f:\mathcal{K}\times B\to\mathbb{R}$ be a piecewise-linear fibered filtration
function. In Section 3.1, we made two generic assumptions to simplify the
exposition. If assumption (1) holds, then every nonempty
$I(\sigma_{i},\sigma_{j})\cap\Delta$ is a line segment $\ell$ that subdivides
triangle $\Delta$ into polygons $Q_{1}$, $Q_{2}$ such that

$(\text{idx}_{f}(\sigma_{i},p_{1})-\text{idx}_{f}(\sigma_{j},p_{1}))\times(\text{idx}_{f}(\sigma_{i},p_{2})-\text{idx}_{f}(\sigma_{j},p_{2}))<0$

for all $p_{1}\in Q_{1}$ and $p_{2}\in Q_{2}$ (i.e., $\sigma_{i}$ and
$\sigma_{j}$ have different relative orders in $Q_{1}$ and $Q_{2}$). We say
that $\sigma_{i}$ and $\sigma_{j}$ _swap along $\ell$_ because $\sigma_{i}$
and $\sigma_{j}$ have different relative orders on either side of $\ell$. If
assumption (1) does not hold, then for any triangle $\Delta$ in $B$ and pair
$(\sigma_{i},\sigma_{j})$ of simplices, it is possible that
$I(\sigma_{i},\sigma_{j})\cap\Delta$ equals either $\Delta$ or an edge of
$\Delta$. If $e$ is an edge of triangle $\Delta$ such that $e\subseteq
I(\sigma_{i},\sigma_{j})\cap\Delta$, then $\sigma_{i}$ and $\sigma_{j}$ _swap
along line segment $e$_ if

$(\text{idx}_{f}(\sigma_{i},p_{1})-\text{idx}_{f}(\sigma_{j},p_{1}))\times(\text{idx}_{f}(\sigma_{i},p_{2})-\text{idx}_{f}(\sigma_{j},p_{2}))<0$

for all $p_{1}\in\Delta_{1}$ and $p_{2}\in\Delta_{2}$, where $\Delta_{1}$ and
$\Delta_{2}$ are the triangles of $B$ that are adjacent to $e$ (i.e.,
$\sigma_{i}$ and $\sigma_{j}$ have different relative orders in $\Delta_{1}$
and $\Delta_{2}$). In Figure 7, we illustrate an example where assumption (1)
does not hold. We highlight the line segments that $\sigma_{i}$ and
$\sigma_{j}$ swap along.

If assumption (2) does not hold, then there may be a line segment $\ell$ in a
triangle $\Delta$ such that
$I(\sigma_{i_{1}},\sigma_{j_{1}})\cap\Delta=\ell=I(\sigma_{i_{2}},\sigma_{j_{2}})\cap\Delta$
for two distinct pairs $(\sigma_{i_{1}},\sigma_{j_{1}})$ and
$(\sigma_{i_{2}},\sigma_{j_{2}})$ of simplices in $\mathcal{K}$.

(a)

(b)

(c)

Figure 7. Examples of fibered filtration functions for which assumption (1) of
Section 3.1 does not hold. (A) The pair $(\sigma_{i},\sigma_{j})$ is a pair of
simplices such that $I(\sigma_{i},\sigma_{j})\cap\Delta=\Delta$ for every pink
triangle $\Delta$ and $I(\sigma_{i},\sigma_{j})\cap\Delta=\ell$ if
$\ell\subseteq\Delta$ is a pink line segment. Without loss of generality
$i<j$, so $\text{idx}_{f}(\sigma_{i},p)<\text{idx}_{f}(\sigma_{j},p)$ if
$f(\sigma_{i},p)=f(\sigma_{j},p)$. (B) Suppose that
$f(\sigma_{j},p)<f(\sigma_{i},p)$ on green triangles and
$f(\sigma_{i},p)<f(\sigma_{j},p)$ on yellow triangles. In blue, we draw the
line segments on which the pair $(\sigma_{i},\sigma_{j})$ swaps. (C) Suppose
that $f(\sigma_{i},p)<f(\sigma_{j},p)$ on green triangles and
$f(\sigma_{j},p)<f(\sigma_{i},p)$ on yellow triangles. In blue, we again draw
the line segments on which the pair $(\sigma_{i},\sigma_{j})$ swaps.

In Sections A.1.2 and A.1.3, I explain the modifications for the algorithm
when we do not make the assumptions of Section 3.1. It suffices to modify step
1 (see Section 3.1.1) and step 2 (see Section 3.1.2). Step 3 (Section 3.1.3)
remains the same.

####  A.1.1. Preliminaries

As in Section 3.1.1, we consider the restriction of $f$ to every edge $e$ of
$B$ to find the vertices of $\mathcal{A}(L)$ that lie on the $1$-skeleton of
$B$.

######  Definition A.1.

A vertex $v$ for a pair $(\sigma,\tau)$ of simplices is _detected along edge
$e$_ if, while traversing edge $e$ during the Bentley–Ottman algorithm, we
detect the point $v\in e$ as a point at which the relative order of $\sigma$
and $\tau$ changes.

A vertex $v$ for the pair $(\sigma,\tau)$ is detected along edge $e$ if and
only if $\sigma$ and $\tau$ have different relative orders at the endpoints of
$e$.

######  Definition A.2.

A line segment $(v,w)$ is _detected in triangle $\Delta$_ if there is a pair
$(\sigma,\tau)$ of simplices such that vertex $v$ is detected along an edge
$e_{1}$ of $\Delta$ for $(\sigma,\tau)$ and vertex $w$ is detected along an
edge $e_{2}$ of $\Delta$ for $(\sigma,\tau)$.

Lemma A.3 characterizes the conditions for which a pair of simplices swaps
along a line segment.

######  Lemma A.3.

Let $(v,w)$ be a line segment that is not on the boundary of $B$.

  1. (1)

If $v$ and $w$ are not the endpoints of an edge in $B$, let $\Delta$ be the
unique triangle that contains $(v,w)$. A pair $(\sigma_{i},\sigma_{j})$ of
simplices swaps along the line segment $(v,w)$ if and only if $(v,w)$ is
detected in triangle $\Delta$.

  2. (2)

If $v$ and $w$ are the endpoints of an edge in $B$, let $\Delta_{1}$,
$\Delta_{2}$ be the two triangles that are adjacent to that edge. A pair
$(\sigma_{i},\sigma_{j})$ of simplices swaps along the line segment $(v,w)$ if
and only if $(v,w)$ is detected in exactly one of $\Delta_{1}$, $\Delta_{2}$.

###### Proof.

Statement (1) was the situation in Section 3.1, so it remains only to prove
statement (2).

Suppose that $(v,w)$ are the endpoints of an edge $e$ in $B$ that is not on
the boundary of $B$. Let $\Delta_{1}$, $\Delta_{2}$ be the two triangles that
are adjacent to $\Delta$. As illustrated in Figure 8a, we denote the third
vertex of $\Delta_{1}$ by $u_{1}$, the third vertex of $\Delta_{2}$ by
$u_{2}$, the other two edges in $\Delta_{1}$ by $e_{2}$, $e_{3}$, and the
other two edges in $\Delta_{2}$ by $e_{4}$, $e_{5}$.

A pair $(\sigma_{i},\sigma_{j})$ swaps along $(v,w)$ only if $e\subseteq
I(\sigma_{i},\sigma_{j})\cap\Delta_{1}\cap\Delta_{2}$. If $e\subseteq
I(\sigma_{i},\sigma_{j})\cap\Delta$ for triangle $\Delta$, then either
$I(\sigma_{i},\sigma_{j})\cap\Delta=e$ or
$I(\sigma_{i},\sigma_{j})\cap\Delta=\Delta$. If both
$I(\sigma_{i},\sigma_{j})\cap\Delta_{1}=\Delta_{1}$ and
$I(\sigma_{i},\sigma_{j})\cap\Delta_{2}=\Delta_{2}$, then
$(\sigma_{i},\sigma_{j})$ does not swap along $(v,w)$ because $\sigma_{i}$ and
$\sigma_{j}$ have the same relative order in $\Delta_{1}$ and $\Delta_{2}$.
Therefore, the pair $(\sigma_{i},\sigma_{j})$ swaps along $(v,w)$ only if the
intersection of $I(\sigma_{i},\sigma_{j})$ with one triangle is $e$ and the
intersection with the other triangle is either $e$ or the entire triangle.
Without loss of generality, $I(\sigma_{i},\sigma_{j})\cap\Delta_{1}=e$ and
either $I(\sigma_{i},\sigma_{j})\cap\Delta_{2}=\Delta_{2}$ or
$I(\sigma_{i},\sigma_{j})\cap\Delta_{2}=e$.

Similarly, the line segment $(v,w)$ can be detected in triangle $\Delta_{k}$
only if $e=I(\sigma_{i},\sigma_{j})\cap\Delta_{k}$. If
$I(\sigma_{i},\sigma_{j})\cap\Delta_{1}=e$, then we must also have $e\subseteq
I(\sigma_{i},\sigma_{j})\cap\Delta_{2}$, so either
$I(\sigma_{i},\sigma_{j})\cap\Delta_{2}=\Delta_{2}$ or
$I(\sigma_{i},\sigma_{j})\cap\Delta_{2}=e$ (and vice versa if
$I(\sigma_{i},\sigma_{j})\cap\Delta_{2}=e$). Therefore, $(v,w)$ is detected in
$\Delta_{k}$ only if $I(\sigma_{i},\sigma_{j})=e$ and the intersection with
the other triangle is either $e$ or the whole triangle. Without loss of
generality, $\Delta_{k}=\Delta_{1}$.

In Figures 8b–8f, we illustrate the possible cases in which
$I(\sigma_{i},\sigma_{j})\cap\Delta_{1}=e$ and either
$I(\sigma_{i},\sigma_{j})\cap\Delta_{2}=\Delta_{2}$ or
$I(\sigma_{i},\sigma_{j})\cap\Delta_{2}=e$. We will show that statement (2)
holds in each of these cases. In all other cases, we have already shown that
neither $(\sigma_{i},\sigma_{j})$ swaps along the line segment $(v,w)$ nor is
$(v,w)$ detected in $\Delta_{1}$ or $\Delta_{2}$.

(a)

(b) Case 1.1

(c) Case 1.2

(d) Case 2.1

(e) Case 2.2

(f) Case 2.3

Figure 8. (A) The vertices, edges, and triangles that were defined in the
proof of Lemma A.3. (B–F) The cases in the proof of Lemma A.3. Pink regions
are regions on which $\sigma_{i}$ and $\sigma_{j}$ have equal filtration
values.

Without loss of generality, we assume that $i<j$, so
$\text{idx}_{f}(\sigma_{i},p)<\text{idx}_{f}(\sigma_{j},p)$ if $p\in
I(\sigma_{i},\sigma_{j})$.

Case 1: $I(\sigma_{i},\sigma_{j})\cap\Delta_{2}=\Delta_{2}$.

There are two subcases.

  1. (1)

Case 1.1:  (see Figure 8b) $f(\sigma_{i},p)<f(\sigma_{j},p)$ for all
$p\in\Delta_{1}\setminus e$.

In this subcase, $\text{idx}_{f}(\sigma_{i},p)<\text{idx}_{f}(\sigma_{j},p)$
for all $p\in\Delta_{1}\cup\Delta_{2}$. Therefore, the pair
$(\sigma_{i},\sigma_{j})$ does not swap along $(v,w)$. Neither $v$ nor $w$ is
detected along any edges of $\Delta_{1}$ or $\Delta_{2}$, so the line segment
$(v,w)$ is not detected in either $\Delta_{1}$ or $\Delta_{2}$.

  2. (2)

Case 1.2:  (see Figure 8c) $f(\sigma_{j},p)<f(\sigma_{i},p)$ for all
$p\in\Delta_{1}\setminus e$.

In this subcase, $\text{idx}_{f}(\sigma_{i},p)<\text{idx}_{f}(\sigma_{j},p)$
for all $p\in\Delta_{2}$ and
$\text{idx}_{f}(\sigma_{j},p)<\text{idx}_{f}(\sigma_{i},p)$ for all
$p\in\Delta_{1}\setminus e$. Therefore, the pair $(\sigma_{i},\sigma_{j})$
swaps along $(v,w)$. The vertex $v$ is detected along edge $e_{2}$ and the
vertex $w$ is detected along the edge $e_{3}$. Because $e_{2}$ and $e_{3}$ are
edges of $\Delta_{1}$, the line segment $(v,w)$ is detected in $\Delta_{1}$.
The vertices $v$ and $w$ are not detected along any edge of $\Delta_{2}$, so
$(v,w)$ is not detected in $\Delta_{2}$.

Case 2: $I(\sigma_{i},\sigma_{j})\cap\Delta_{2}=e$.

There are three subcases.

  1. (1)

Case 2.1: (see Figure 8d) $f(\sigma_{j},p)<f(\sigma_{i},p)$ for all
$p\in(\Delta_{1}\cup\Delta_{2})\setminus e$.

In this subcase, $\text{idx}_{f}(\sigma_{j},p)<\text{idx}_{f}(\sigma_{i},p)$
for all $p\in(\Delta_{1}\cup\Delta_{2})\setminus e.$ Therefore, the pair
$(\sigma_{i},\sigma_{j})$ does not swap along $(v,w)$. The vertex $w$ is
detected along edges $e_{3}$ and $e_{5}$. The vertex $v$ is detected along
edges $e_{2}$ and $e_{4}$. Therefore, $(v,w)$ is detected in both $\Delta_{1}$
and $\Delta_{2}$.

  2. (2)

Case 2.2: (see Figure 8e) Either we have $f(\sigma_{j},p)<f(\sigma_{i},p)$ for
all $p\in\Delta_{1}\setminus e$ and $f(\sigma_{i},p)<f(\sigma_{j},p)$ for all
$p\in\Delta_{2}\setminus e$ or we have $f(\sigma_{i},p)<f(\sigma_{j},p)$ for
all $p\in\Delta_{1}\setminus e$ and $f(\sigma_{j},p)<f(\sigma_{i},p)$ for all
$p\in\Delta_{1}\setminus e$. Without loss of generality, we assume the former.

In this subcase, $\text{idx}_{f}(\sigma_{i},p)<\text{idx}_{f}(\sigma_{j},p)$
for all $p\in\Delta_{2}$ and
$\text{idx}_{f}(\sigma_{j},p)<\text{idx}_{f}(\sigma_{i},p)$ for all
$p\in\Delta_{1}\setminus e$. Therefore, the pair $(\sigma_{i},\sigma_{j})$
swaps along $(v,w)$. The vertex $v$ is detected along $e_{2}$ and the vertex
$w$ is detected along $e_{3}$, so $(v,w)$ is detected in triangle
$\Delta_{1}$. Neither $v$ nor $w$ is detected along any edge of $\Delta_{2}$,
so $(v,w)$ is not detected in $\Delta_{2}$.

  3. (3)

Case 2.3: (see Figure 8f) $f(\sigma_{i},p)<f(\sigma_{j},p)$ for all
$p\in(\Delta_{1}\cup\Delta_{2})\setminus e$.

In this subcase, $\text{idx}_{f}(\sigma_{i},p)<\text{idx}_{f}(\sigma_{j},p)$
for all $p\in\Delta_{1}\cup\Delta_{2}$. Therefore, the pair
$(\sigma_{i},\sigma_{j})$ does not swap along $(v,w)$. Neither $v$ nor $w$ is
detected along any edge of $\Delta_{1}$ or $\Delta_{2}$, so $(v,w)$ is not
detected in either $\Delta_{1}$ or $\Delta_{2}$.

∎

We use Lemma A.4 to modify step 2 of the algorithm: computing the simplex
pairing function.

######  Lemma A.4.

Let $\text{idx}_{0},\text{idx}_{1}:\mathcal{K}\to\\{1,\ldots,N\\}$ denote two
different simplex indexings (not necessarily compatible with $f$; see
Definition 2.1), where $N$ is the number of simplices in $\mathcal{K}$. Let
$\\{(\sigma_{i_{k}},\sigma_{j_{k}})\\}_{k=1}^{m}$ be the set of pairs
$(\sigma_{i_{k}},\sigma_{j_{k}})$ such that

$(\text{idx}_{0}(\sigma_{i_{k}})-\text{idx}_{0}(\sigma_{j_{k}}))\times(\text{idx}_{1}(\sigma_{i_{k}})-\text{idx}_{1}(\sigma_{j_{k}}))<0\,.$

That is, $\sigma_{i_{k}}$ and $\sigma_{j_{k}}$ have different relative orders
in $\text{idx}_{0}$ and $\text{idx}_{1}$. Let $\zeta_{0}:=\text{idx}_{0}$, and
for $k\in\\{1,\ldots,m\\}$, let $\zeta_{k}:\mathcal{K}\to\\{1,\ldots,N\\}$ be
the simplex indexing obtained by transposing $(\sigma_{i_{k}},\sigma_{j_{k}})$
in the simplex indexing $\zeta_{k-1}$. If $\zeta_{k-1}(\sigma_{i_{k}})$ and
$\zeta_{k-1}(\sigma_{j_{k}})$ are consecutive integers for all $k$, then
$\zeta_{m}=\text{idx}_{1}$. Furthermore, the sequence
$\\{(\sigma_{i_{k}},\sigma_{j_{k}})\\}_{k=1}^{m}$ can be ordered so that this
conditions holds.

###### Proof.

First, we prove that there is at least one pair
$(\sigma_{i_{k}},\sigma_{j_{k}})$ such that $\text{idx}_{0}(\sigma_{i_{k}})$
and $\text{idx}_{0}(\sigma_{j_{k}})$ are consecutive integers. Let
$k_{*}=\operatorname*{arg\,min}_{k}\lVert\text{idx}_{0}(\sigma_{i_{k}})-\text{idx}_{0}(\sigma_{j_{k}})\rVert$.
To obtain a contradiction, suppose that
$s_{1}:=\text{idx}_{0}(\sigma_{i_{k_{*}}})$ and
$s_{2}:=\text{idx}_{0}(\sigma_{j_{k_{*}}})$ are not consecutive integers.
Without loss of generality, $s_{1}<s_{2}$. For
$r\in\\{1,\ldots,s_{2}-s_{1}-1\\}$, let
$\tau_{s_{1}+r}:=\text{idx}^{-1}_{0}(s_{1}+r)$. That is,
$\tau_{s_{1}+1},\ldots,\tau_{s_{2}-1}$ are the simplices between
$\sigma_{i_{k_{*}}}$ and $\sigma_{j_{k_{*}}}$. For all $r$, either
$(\sigma_{i_{k_{*}}},\tau_{s_{1}+r})\in\\{(\sigma_{i_{k}},\sigma_{j_{k}})\\}_{k=1}^{m}$
or
$(\sigma_{j_{k_{*}}},\tau_{s_{1}+r})\in\\{(\sigma_{i_{k}},\sigma_{j_{k}})\\}_{k=1}^{m}$.
That is, either the relative order of $\sigma_{i_{k_{*}}}$ and
$\tau_{s_{1}+r}$ changes or the relative order of $\sigma_{j_{k_{*}}}$ and
$\tau_{s_{1}+r}$ changes. By definition of $k_{*}$, none of the pairs in the
set $\\{(\tau_{s_{1}+r_{1}},\tau_{s_{1}+r_{2}})\\}_{r_{1},r_{2}}$ are in
$\\{(\sigma_{i_{k}},\sigma_{j_{k}})\\}_{k=1}^{m}$. That is, the relative order
of the simplices $\tau_{s_{1}+1},\ldots,\tau_{s_{2}-1}$ does not change.
Therefore, either
$(\sigma_{i_{k_{*}}},\tau_{s_{1}+1})\in\\{(\sigma_{i_{k}},\sigma_{j_{k}})\\}_{k=1}^{m}$
or
$(\sigma_{j_{k_{*}}},\tau_{s_{1}+r})\in\\{(\sigma_{i_{k}},\sigma_{j_{k}})\\}_{k=1}^{m}$
for all $r$. In either case, one of these is a transposition of simplices
whose indices in $\text{idx}_{0}$ are consecutive integers, which is a
contradiction.

Now we prove the lemma by induction on $m$. When $m=1$, we just showed that
$\text{idx}_{0}(\sigma_{i_{1}})$ and $\text{idx}_{0}(\sigma_{j_{1}})$ are
consecutive integers, so $\zeta_{1}=\text{idx}_{1}$. In the general case, we
can assume $\text{idx}_{0}(\sigma_{i_{1}})$ and
$\text{idx}_{0}(\sigma_{j_{1}})$ are consecutive integers without loss of
generality. The set $\\{(\sigma_{i_{k}},\sigma_{j_{k}})\\}_{k=2}^{m}$ is the
set of pairs $(\sigma_{i_{k}},\sigma_{j_{k}})$ such that

$(\zeta_{1}(\sigma_{i_{k}})-\zeta_{1}(\sigma_{j_{k}}))\times(\text{idx}_{1}(\sigma_{i_{k}})-\text{idx}_{1}(\sigma_{j_{k}}))<0\,.$

That is, $\sigma_{i_{k}}$ and $\sigma_{j_{k}}$ have different relative orders
in $\zeta_{1}$ and $\text{idx}_{1}$ for all $k\in\\{2,\ldots,m\\}$. By
induction, we can assume that the sequence
$\\{(\sigma_{i_{k}},\sigma_{j_{k}})\\}_{k=2}^{m}$ is ordered such that
$\zeta_{k-1}(\sigma_{i_{k}})$ and $\zeta_{k-1}(\sigma_{j_{k}})$ are
consecutive integers for $k\in\\{2,\ldots,m\\}$ and
$\zeta_{m}=\text{idx}_{1}$. ∎

####  A.1.2. Modifications to step 1: Computing the polygons

In Section 3.1.1, we maintained a dictionary $\mathcal{D}_{1}(\Delta)$ for
each triangle $\Delta\in B$. The keys were pairs $(\sigma,\tau)$ such that
$I(\sigma,\tau)\cap\Delta$ was a line segment in $\Delta$, and the value of a
key $(\sigma,\tau)$ was the list $[v,w]$ of vertices in $\mathcal{A}(L)$ that
were the endpoints of the line segment $I(\sigma,\tau)\cap\Delta$.

Now we maintain two additional dictionaries $\mathcal{D}_{2}(\Delta)$ and
$\mathcal{D}_{3}(\Delta)$ for each triangle $\Delta\in B$. We initialize these
dictionaries to be empty, and we update them as we traverse the edges of $B$.
At any time in this process, the keys of $\mathcal{D}_{2}(\Delta)$ are pairs
$(v,w)$ of vertices in $\mathcal{A}(L)$ such that

  1. (1)

the line segment $(v,w)$ has been detected in $\Delta$,

  2. (2)

the line segment $(v,w)$ is not an edge of $B$.

The value of $\mathcal{D}_{2}(\Delta)[(v,w)]$ is a list
$[(\sigma_{i_{1}},\sigma_{j_{1}}),\ldots,(\sigma_{i_{m}},\sigma_{j_{m}})]$ of
the simplex pairs that we have found so far such that $\sigma_{i_{k}}$ and
$\sigma_{j_{k}}$ swap along $(v,w)$. The keys of $\mathcal{D}_{3}(\Delta)$ are
vertices $v\in\mathcal{A}(L)$ such that

  1. (1)

vertex $v$ has been detected along an edge $e$ of triangle $\Delta$,

  2. (2)

there is a pair $(\sigma,\tau)$ of simplices such that $(\sigma,\tau)$ swaps
at $v$ and we have not yet found a vertex $w$ such that
$I(\sigma,\tau)\cap\Delta=(v,w)$.

We modify the algorithm of Section 3.1.1 as follows. Suppose that we detect a
vertex $v$ along edge $e$ in $\mathcal{A}(L)$ for the set
$\\{(\sigma_{i_{1}},\sigma_{j_{1}}),\ldots,(\sigma_{i_{m}},\sigma_{j_{m}})\\}$
of simplex pairs. We do the following:

  1. (1)

Update $\mathcal{D}_{1}$: For each triangle $\Delta\in B$ that is adjacent to
$e$, we append $v$ to the list of vertices for
$\mathcal{D}_{1}(\Delta)[(\sigma_{i_{k}},\sigma_{j_{k}})]$ for all $k$, as in
Section 3.1.1.

  2. (2)

Update $\mathcal{A}(L)$: If $v$ is not an endpoint of $e$, we split the edge
$e$ in $\mathcal{A}(L)$ and add an internal vertex in $e$, as in Section
3.1.1. If $v$ is an endpoint of $e$, we do not split the edge or create a new
vertex because $B$ already has a vertex at $v$.

  3. (3)

Update $\mathcal{D}_{2}$, $\mathcal{D}_{3}$, and the edge labels:  For each
triangle $\Delta$ that is adjacent to $e$ and each
$(\sigma_{i_{k}},\sigma_{j_{k}})$:

     * •

If $v$ is the only vertex in the list
$\mathcal{D}_{1}(\Delta)[(\sigma_{i_{k}},\sigma_{j_{k}})]$, then we have not
yet detected a line segment for $(\sigma_{i_{k}},\sigma_{j_{k}})$ of the form
$(v,w)$ for some vertex $w$. We do the following: If $v$ is not in
$\mathcal{D}_{3}(\Delta)$, add the key $v$ to $\mathcal{D}_{3}(\Delta)$ with
value $[(\sigma_{i_{k}},\sigma_{j_{k}})]$. Otherwise, append
$(\sigma_{i_{k}},\sigma_{j_{k}})$ to $\mathcal{D}_{3}(\Delta)[v]$.

     * •

Otherwise, there is another vertex
$w\in\mathcal{D}_{1}(\Delta)[(\sigma_{i_{k}},\sigma_{j_{k}})]$. This implies
that we have just detected a line segment $(v,w)$ in $\Delta$ for
$(\sigma_{i_{k}},\sigma_{j_{k}})$. We remove $(\sigma_{i_{k}},\sigma_{j_{k}})$
from $\mathcal{D}_{3}(\Delta)[w]$.

       * –

If $v$ and $w$ are not both vertices of $B$, then $(v,w)$ is not an edge of
$B$. We do the following: If $(v,w)$ is not in $\mathcal{D}_{2}(\Delta)$, then
add key $(v,w)$ to $\mathcal{D}_{2}(\Delta)$ with value
$[(\sigma_{i_{k}},\sigma_{j_{k}})]$. Otherwise, append
$(\sigma_{i_{k}},\sigma_{j_{k}})$ to $\mathcal{D}_{2}(\Delta)[(v,w)]$.

       * –

Otherwise, $v$ and $w$ are both vertices of triangle $\Delta$. This means that
we have detected a line segment $(v,w)$ in $\Delta$ in which $v$ and $w$ are
the endpoints of an edge $e^{\prime}$ in $B$. If $e^{\prime}$ is an edge on
the boundary of $B$, then we do nothing. Otherwise, let $\Delta_{2}$ be the
other triangle that is adjacent to $e^{\prime}$. By Lemma A.3, the pair
$(\sigma_{i_{k}},\sigma_{j_{k}})$ swaps along $(v,w)$ if and only if the line
segment is not detected in $\Delta_{2}$. If $e^{\prime}$ already stores a
reference to $(\sigma_{i_{k}},\sigma_{j_{k}})$, then we remove it because this
implies that $e^{\prime}$ was already detected in $\Delta_{2}$. Otherwise, we
add a reference to $(\sigma_{i_{k}},\sigma_{j_{k}})$ on $e^{\prime}$.

When the traversal of the $1$-skeleton is done, we add lines to
$\mathcal{A}(L)$. For every triangle $\Delta\in B$ and every key
$(v,w)\in\mathcal{D}_{2}(\Delta)$, we add a line segment with endpoints $v$,
$w$ to the DCEL that represents $\mathcal{A}(L)$. For every edge in the DCEL
that is a subset of the line segment $(v,w)$, we label the edge with a
reference to the list $\mathcal{D}_{2}(\Delta)[(v,w)]$, which is the list
$\\{(\sigma_{i_{1}},\sigma_{j_{1}}),\ldots,(\sigma_{i_{m}},\sigma_{j_{m}})\\}$
of simplex pairs that swap along the line segment.

####  A.1.3. Modifications to step 2: Computing the pairing function

We compute a path $\Gamma$ as in Section 3.1.2 and traverse $\Gamma$. At each
step, we walk from one polygon $P_{1}$ to the next polygon $P_{2}$ by crossing
an edge $e$ in $\mathcal{A}(L)$. The edge $e$ stores a list of simplex pairs
$(\sigma,\tau)$ such that $\sigma$ and $\tau$ have different relative orders
in the polygons $P_{1}$, $P_{2}$. We update the simplex indexing one
transposition at a time. Let
$\overline{\text{idx}}:\mathcal{K}\to\\{1,\ldots,N\\}$ denote the current
indexing, which we initialize to the simplex indexing
$\text{idx}_{f}(\cdot,P_{1})$ in $P_{1}$. While the list that $e$ stores is
nonempty, we do the following:

  1. (1)

Let $(\sigma,\tau)$ be the first element of the list.

  2. (2)

If $\overline{\text{idx}}(\sigma)$ and $\overline{\text{idx}}(\tau)$ are
consecutive integers, then we update $\overline{\text{idx}}$ by swapping the
order of $\sigma$ and $\tau$. As in Section 3.1.2, we also update the RU
decomposition, and the (birth, death) simplex pairs. We remove $(\sigma,\tau)$
from the list.

  3. (3)

Otherwise, we move $(\sigma,\tau)$ to the end of the list.

At the end of this algorithm, $\overline{\text{idx}}$ is the simplex indexing
$\text{idx}_{f}(\cdot,P_{2})$ in $P_{2}$ (by Lemma A.4), the RU decomposition
is an RU decomposition for $P_{2}$, and we have computed the (birth, death)
simplex pairs for $P_{2}$.

## References

  * [1] Bernard Chazelle and Joel Friedman.  Point location among hyperplanes and unidirectional ray-shooting.  Computational Geometry, 4(2):53–62, 1994.
  * [2] Yao-li Chuang, Maria R. D’Orsogna, Daniel Marthaler, Andrea L. Bertozzi, and Lincoln S. Chayes.  State transitions and the continuum limit for a 2D interacting, self-propelled particle system.  Physica D: Nonlinear Phenomena, 232(1):33–47, 2007.
  * [3] David Cohen-Steiner, Herbert Edelsbrunner, and Dmitriy Morozov.  Vines and vineyards by updating persistence in linear time.  In Proceedings of the Twenty-Second Annual Symposium on Computational Geometry, SCG ’06, pages 119–126, New York, NY, USA, 2006. Association for Computing Machinery.
  * [4] Mark de Berg, Otfried Cheong, Marc van Kreveld, and Mark Overmars.  Computational Geometry: Algorithms and Applications.  Springer, Berlin, Heidelberg, 3rd edition, 2008.
  * [5] Herbert Edelsbrunner and John Harer.  Computational Topology: An Introduction.  American Mathematical Society, Providence, RI, 2010.
  * [6] Herbert Edelsbrunner, David Letscher, and Afra Zomorodian.  Topological persistence and simplification.  Discrete & Computational Geometry, 28:511–533, 2002.
  * [7] Abigail Hickok.  Persistence diagram bundles: A multidimensional generalization of vineyards.  arXiv:2210.05124, 2022.
  * [8] Michael Lesnick and Matthew Wright.  Interactive visualization of 2-D persistence modules.  arXiv:1512.00180, 2015.
  * [9] Yanjie Li, Dingkang Wang, Giorgio A. Ascoli, Partha Mitra, and Yusu Wang.  Metrics for comparing neuronal tree shapes based on persistent homology.  PLoS ONE, 12(8):e0182184, 2017.
  * [10] Jesus De Loera, Joerg Rambau, and Francisco Santos.  Triangulations: Structures for Algorithms and Applications, chapter 6.3.2, pages 314–316.  Springer Science & Business Media, 2010.
  * [11] S. Meiser.  Point location in arrangements of hyperplanes.  Information and Computation, 106(2):286–303, 1993.
  * [12] Nina Otter, Mason A. Porter, Ulrike Tillmann, Peter Grindrod, and Heather A. Harrington.  A roadmap for the computation of persistent homology.  EPJ Data Science, 6:17, 2017.
  * [13] Neil Sarnak and Robert E. Tarjan.  Planar point location using persistent search trees.  Communications of the ACM, 29(7):669–679, 1986.
  * [14] Csaba D. Toth, Joseph O’Rourke, and Jacob E. Goodman.  Handbook of discrete and computational geometry.  CRC Press, Boca Raton, FL, 3rd edition, 2017.
  * [15] Katharine Turner, Sayan Mukherjee, and Doug M. Boyer.  Persistent homology transform for modeling shapes and surfaces.  Information and Inference: A Journal of the IMA, 3(4):310–344, 2014\.