Exact computation of the matching distance on 2-parameter persistence modules

TL;DR

This paper presents a polynomial-time algorithm for exactly computing the matching distance between 2-parameter persistence modules, enhancing the tool's practicality for topological data analysis.

Contribution

It introduces a novel method to compute the matching distance exactly in polynomial time for 2-parameter modules, using a subdivision of the line space into regions.

Findings

Matching distance can be computed exactly in polynomial time for 2-parameter modules.

The matching distance is a rational number when input modules have rational bigrades.

The approach involves subdividing the space of affine lines into regions for analysis.

Abstract

The matching distance is a pseudometric on multi-parameter persistence modules, defined in terms of the weighted bottleneck distance on the restriction of the modules to affine lines. It is known that this distance is stable in a reasonable sense, and can be efficiently approximated, which makes it a promising tool for practical applications. In this work, we show that in the 2-parameter setting, the matching distance can be computed exactly in polynomial time. Our approach subdivides the space of affine lines into regions, via a line arrangement. In each region, the matching distance restricts to a simple analytic function, whose maximum is easily computed. As a byproduct, our analysis establishes that the matching distance is a rational number, if the bigrades of the input modules are rational.

Figures

Equations

$$X_1\subseteq X_2\subseteq \dots \subseteq X_n$$

$$H_p(X_1,\mathbb{K})\to H_p(X_2,\mathbb{K})\to \dots \to H_p(X_n,\mathbb{K})$$

$${\mathrm{Gen}}_x:=\{e_i\mid \mathrm{gr}({\mathrm{row}}_i)\le x\}$$

$${\mathrm{Rel}}_x:=\{{\mathrm{col}}_j\mid \mathrm{gr}({\mathrm{col}}_j)\le x\}$$

$$M_x^P:=\mathrm{span}({\mathrm{Gen}}_x)/\mathrm{span}({\mathrm{Rel}}_x)$$

$$A$$

$$B$$

$$d_B(D_1,D_2)=\min \{ϵ\mid \text{there exists a matching of cost ϵ between D1 and D2}\}.$$

$$M_p=\{\begin{array}{ll}\mathbb{K}& \text{ if p∈[0,1)×[0,1)}\\ 0& \text{ otherwise,}\end{array}{M}_{p\to q}=\{\begin{array}{ll}{\mathrm{Id}}_{\mathbb{K}}& \text{ if p,q∈[0,1)×[0,1)}\\ 0& \text{otherwise.}\end{array}.$$

$$\blockarray{ccc}&(1,0)(0,1)\\ \block{c[cc]}(0,0)11\\ .$$

$$w(\ell ):=\{\begin{array}{ll}\frac{1}{\sqrt{1+s^2}}& s\ge 1\\ \frac{1}{\sqrt{1+\frac{1}{s^2}}}& 0<s<1\end{array}$$

$$F^{M,N}(\ell ):=w(\ell )\cdot d_B(M^{\ell },N^{\ell }).$$

$$\{L_v\mid v\in \mathrm{gr}(P^M)\cup \mathrm{gr}(P^N)\}.$$

$$\mathrm{push}(p,\ell ):=\{\begin{array}{ll}(p_x,s\cdot p_x+t)& \text{if p lies below ℓ}\\ (\frac{p_y-t}{s},p_y)& \text{if p lies on or above ℓ}\end{array}$$

$${\delta }_{p,q}(s,t):=\Vert \mathrm{push}(p,\ell )-\mathrm{push}(q,\ell ){\Vert }_2=|{\lambda }^{-1}\circ \mathrm{push}(p,\ell )-{\lambda }^{-1}\circ \mathrm{push}(q,\ell )|$$

$${\delta }_{p,q}(s,t)=\Vert (\frac{p_y-t}{s},p_y)-(\frac{q_y-t}{s},q_y){\Vert }_2=\sqrt{{\left(\frac{p_y-q_y}{s}\right)}^2+(p_y-q_y{)}^2}=|p_y-q_y|\sqrt{1+\frac{1}{s^2}}.$$

$${\delta }_{p,q}(s,t)=\Vert (p_x,s{p}_x+t)-(q_x,s{q}_x+t){\Vert }_2=\sqrt{(p_x-q_x{)}^2+(s{p}_x-s{q}_x{)}^2}=|p_x-q_x|\sqrt{1+s^2}.$$

$${\delta }_{p,q}(s,t)$$

$$=\sqrt{\frac{1}{s^2}{\left(p_y-t-s{q}_x\right)}^2+(p_y-s{q}_x-t{)}^2}=|p_y-t-s{q}_x|\sqrt{1+\frac{1}{s^2}}.$$

$$\begin{gathered}\begin{array}[]{rcl}\delta_{p,q}(s,t)&=&0\qquad\qquad\qquad\ \ \ \ \,\textup{for $p,q\in\mathrm{gr}(P^{M})$ or $p,q\in\mathrm{gr}(P^{N})$},\\[5.0pt] c_{pq}\delta_{p,q}(s,t)&=&c_{p^{\prime}q^{\prime}}\delta_{p^{\prime},q^{\prime}}(s,t)\qquad\textup{for $p,q,p^{\prime},q^{\prime}\in\mathrm{gr}(P^{M})\sqcup\mathrm{gr}(P^{N})$},\end{array}\end{gathered}$$

$$c_{pq}:=\{\begin{array}{ll}\frac{1}{2}& \text{if p,q∈gr(PM) or p,q∈gr(PN)},\\ 1& \text{otherwise}.\end{array}$$

$$|p_y-q_y|\sqrt{1+\frac{1}{s^2}}=0$$

$$|p_x-q_x|\sqrt{1+s^2}=0.$$

$$|p_y-t-s{q}_x|\sqrt{1+\frac{1}{s^2}}=0$$

$$|p_y-q_y|\sqrt{1+\frac{1}{s^2}}=|p_y^{\prime }-q_y^{\prime }|\sqrt{1+\frac{1}{s^2}}$$

$$|p_y-q_y|\sqrt{1+\frac{1}{s^2}}=|p_x^{\prime }-q_x^{\prime }|\sqrt{1+s^2}.$$

$$|p_y-q_y|=s|p_x^{\prime }-q_x^{\prime }|$$

$$|p_y-t-s{q}_x|\sqrt{1+\frac{1}{s^2}}=|p_y^{\prime }-t-s{q}_x^{\prime }|\sqrt{1+\frac{1}{s^2}}$$

$$|p_y-q_y|\sqrt{1+\frac{1}{s^2}}=|p_y-t-s{q}_x|\sqrt{1+\frac{1}{s^2}}.$$

$$|p_x-q_x|\sqrt{1+s^2}=|p_y^{\prime }-t-s{q}_x^{\prime }|\sqrt{1+\frac{1}{s^2}}$$

$$s|p_x-q_x|=|p_y^{\prime }-t-s{q}_x^{\prime }|$$

Full text

\ccsdesc[100]Mathematics of computing Algebraic topology \ccsdesc[100]Mathematics of computing Mathematical optimization \EventEditorsGill Barequet and Yusu Wang \EventNoEds2 \EventLongTitle35th International Symposium on Computational Geometry (SoCG 2019) \EventShortTitleSoCG 2019 \EventAcronymSoCG \EventYear2019 \EventDateJune 18–21, 2019 \EventLocationPortland, United States \EventLogosocg-logo \SeriesVolume129 \ArticleNo46 \hideLIPIcs

Exact Computation of the Matching Distance on 2-Parameter Persistence Modules

Michael Kerber

  

Michael Lesnick

  

Steve Oudot

1 Introduction

Multi-parameter persistent homology is receiving growing attention, both from the theoretical and computational points of view. Its motivation lies in the possibility of extending the success of topological data analysis to settings where the structure of data is best captured by 2-parameter rather than 1-parameter constructions. The basic algebraic objects of study in multi-parameter persistence are certain commutative diagrams of vector spaces called persistence modules. In the 1-parameter setting, persistence modules decompose in an essentially unique way into simple summands called interval modules. The decomposition is specified by a discrete invariant called a persistence diagram. In contrast, the algebraic structure of a 2-parameter persistence module (henceforth, bipersistence module) can be far more complex. As a result, a good definition of persistence diagram is unavailable for bipersistence modules [5].

Nevertheless, it is still possible to define meaningful notions of distance between multi-parameter persistence modules. Distances on 1-parameter persistence modules play an essential role in both theory and applications. To extend such theory and applications to the multi-parameter setting, one needs to select a suitable distance on multi-parameter persistence modules. However, progress on finding well-behaved, efficiently computable distances on multi-parameter persistence modules has been slow, and this is been an impediment to progress in practical applications.

The most widely studied and applied distances in the 1-parameter setting are the bottleneck distance and the Wasserstein distance [11, 20]. Both can be efficiently computed via publicly available code [14]. In the multi-parameter setting, the distance that has received the most attention is the interleaving distance. On 1-parameter persistence modules, the interleaving and bottleneck distances are equal [16]. The interleaving distance is theoretically well-behaved; in particular, on modules over prime fields, it is the most discriminative distance among the ones satisfying a certain stability condition [16]. However, it was proven recently to be NP-hard to compute or even approximate to any constant factor less than three on bipersistence modules, even under the restriction that the modules decompose as direct sums of interval modules [4]. While the problem lies in $P$ when further restricting the focus to the sub-class of interval modules themselves [10], in practice one often encounters modules that are not interval modules (nor even direct sums thereof).

This motivates the search for a more computable surrogate for the interleaving distance. One natural candidate is the matching distance, introduced by Cerri et al. [6]. It is a lower bound for the interleaving distance; this is implicit in [6] and shown explicitly in [15]. In the 2-parameter setting, the matching distance is defined as follows: Given a pair of bipersistence modules, we call an affine line $\ell$ in parameter space with positive slope a slice. Restricting the modules to $\ell$ yields a pair of 1-parameter persistence modules, which we call slice modules. These slice modules have a well-defined bottleneck distance, which we multiply by a positive weight to ensure that the matching distance will be a lower bound for the interleaving distance. The matching distance is defined as the supremum of these weighted bottleneck distances over all slices. See Section 3 for the precise definition. The definition generalizes readily to $n$-parameter persistence modules, for any $n\geq 1$; when $n=1$, the matching distance is equal to the bottleneck distance.

As Cerri et al. have observed, to approximate the matching distance up to any (absolute) precision, it suffices to sample the space of slices sufficiently densely and to return the maximum weighted bottleneck distance encountered. For a constant number of scale parameters and approximation quality $\epsilon$, a polynomial number of slices are sufficient in terms of module size and $\frac{1}{\epsilon}$, yielding a polynomial time approximation algorithm. [3]. This approach has been recently applied to the virtual ligand screening problem in computational chemistry [13]. To the best of our knowledge, there is no other previous work in which the problem of computing the matching distance has been considered.

Our contribution

We give an algorithm that computes the exact matching distance between a pair of bipersistence modules in time polynomial with respect to the size of the input. We assume that each persistence module is specified by a presentation. Concretely, this means that the module is specified by a matrix encoding the generators and relations, with each row and each column labeled by a point in $\mathbb{R}^{2}$; see Section 2.

To explain our strategy for computing the matching distance, consider the function $F$ that assigns a slice to its weighted bottleneck distance. The matching distance is then simply the supremum of $F$, taken over all slices. $F$ has a rather complicated structure, since it depends on the longest edge of a perfect matching in a bipartite graph whose edges lengths depend on both the slice and the two modules given as input. When the slice changes, the matching realizing the bottleneck distance undergoes combinatorial changes, making the function $F$ difficult to treat analytically.

We show, however, that the space of slices can be divided into polynomially many regions so that the restriction of $F$ to each region takes a simple closed from. Perhaps surprisingly, if we parameterize the space of slices by the right half plane $\Omega\subset\mathbb{R}^{2}$, the boundary between these regions can be expressed by the union of polynomially many lines in $\Omega$, making each region convex and bounded by (possibly unbounded) line segments. (This is analogous to the observation of [17] that for a single persistence module, the locus of lines where the combinatorial structure underlying the slice module can change is described by a line arrangement.) Moreover, the restriction of $F$ to each cell attains its supremum at a boundary vertex of the cell, or as the limit of an unbounded line segment in $\Omega$; this follows from a straightforward case analysis. These observations together lead to a simple polynomial time algorithm to compute the matching distance.

The characterization of the matching distance underlying our algorithm also makes clear that if the row and column labels of the presentations of the input modules have rational coordinates, then the matching distance is rational as well. We are not aware of a simpler argument for this property.

Outline

We introduce the underlying topological concepts in Section 2, and introduce the matching distance in Section 3. We define the line arrangement subdividing the slice space in Section 4 and give the algorithm to maximize each cell of the arrangement in Section 5. We conclude in Section 6.

2 Persistence modules

Single-parameter modules

Let $\mathbb{K}$ be a fixed finite field throughout. A persistence module $M$ over $\mathbb{R}$ is an assignment of $\mathbb{K}$-vector spaces $M_{x}$ to real numbers $x$, and linear maps $M_{x\to y}:M_{x}\to M_{y}$ to a pair of real numbers $x\leq y$, such that $M_{x\to x}$ is the identity and $M_{y\to z}\circ M_{x\to y}=M_{x\to z}$. Equivalently, in categorical terms, a persistence module is a functor from $\mathbb{R}$ (considered as a poset category) to the category of $\mathbb{K}$-vector spaces.

A common way to arrive at a persistence module is to consider a nested sequence

[TABLE]

of simplicial complexes and to apply homology with respect to a fixed dimension and base field $\mathbb{K}$. This yields a sequence

[TABLE]

of vector spaces and linear maps. To obtain a persistence module over $\mathbb{R}$, we pick grades $s_{1}<s_{2}<\ldots<s_{n}$ and set $M_{x}:=0$ if $x<s_{1}$ and $M_{x}:=H_{p}(X_{i},\mathbb{K})$ with $i=\max\{j\mid s_{j}\leq x\}$ otherwise. For $y\geq x$ and $M_{y}=H_{p}(X_{j},\mathbb{K})$, we define $M_{x\to y}$ as the map $H_{p}(X_{i},\mathbb{K})\to H_{p}(X_{j},\mathbb{K})$ induced by the inclusion map $X_{i}\to X_{j}$.

Finite presentations

In this work, we restrict our attention to persistence modules that are finitely presented in the following sense. A finite presentation is an $k\times m$ matrix $P$ over $\mathbb{K}$, where each row and each column is labeled by a number in $\mathbb{R}$, called the grade, such that if $P_{ij}\neq 0$, then $\mathrm{gr}(\mathrm{row}_{i})\leq\mathrm{gr}(\mathrm{col}_{j})$; here $\mathrm{gr}(-)$ denotes the grade of a row or column. We refer to the multiset of grades of all rows and columns of $P$ simply as the *set of grades of *$P$, and denote this set as $\mathrm{gr}(P)$. A finite presentation gives rise to a persistence module, as we describe next. The rows of $P$ represent the generators of the module, while the columns of $P$ encode relations (or syzygies) on the generators. Concretely, let $e_{1},\ldots,e_{k}$ denote the standard basis of $\mathbb{K}^{k}$, and for $x\in\mathbb{R}$ define the subset

[TABLE]

Likewise, define

[TABLE]

Then, we define

[TABLE]

and $M^{P}_{x\to y}$ simply as the map induced by the inclusion map $\mathrm{span}(\mathrm{Gen}_{x})\to\mathrm{span}(\mathrm{Gen}_{y})$. It can be checked easily that this indeed defines a persistence module $M^{P}$. If a persistence module $N$ is isomorphic to $M^{P}$, we say that $P$ is a presentation of $N$. We call a persistence module $N$ finitely presented if there exists a finite presentation of $N$. For instance, persistence modules as above arising from a finite simplicial filtration are finitely presented. Also, the representation theorem of persistence [21, 7] states that the category of persistence modules over $\mathbb{K}$ is isomorphic to the category of graded $\mathcal{R}$-modules with an appropriately chosen ring $\mathcal{R}$. With that, a persistence module is finitely presented if and only if the corresponding $\mathcal{R}$-module is finitely presented (in the sense of a module).

We assume for concreteness that a finite presentation is encoded in terms of a sparse matrix representation [11]; the size of the presentation is understood to be the bitsize of this sparse matrix. In what follows, for finite a presentation with $k$ generators and $m$ relations, we will express complexity bounds in terms of $n:=k+m$, the number of generators and relations. Note that an algorithm polynomial in $n$ is also polynomial in the size of the presentation.

Persistence diagrams

A persistence diagram is a finite multi-set of points of the form $(b,d)\in\mathbb{R}\times(\mathbb{R}\cup\{\infty\})$ with $b<d$. A well-known structure theorem tells us that we can associate to any finitely presented persistence module $M$ a persistence diagram $D(M)$, and that this determines $M$ up to isomorphism [8].

Given a presentation of $M$, the persistence diagram can be computed by bringing the presentation matrix into echelon form. This process takes cubic time in $n$ using Gaussian elimination [11, 21], or $O(n^{\omega})$ time using fast matrix multiplication, where $\omega\leq 2.373$ [19]. Elimination-based approaches to computing persistent homology perform very well in practice, and are routinely used to study real data.

Lemma 2.1**.**

For $P$ a finite presentation of a persistence module $M$ and $(b,d)\in D(M)$,

$b$* is the grade of some row of $P$,* 2. 2.

if $d<\infty$, then $d$ is the grade of some column of $P$.

Proof 2.2**.**

This follows from the correctness of the basic matrix reduction algorithm for computing persistent homology, as described in [21].

Bottleneck distance

Consider two persistence diagrams $D_{1}$ and $D_{2}$ and a bijection $\sigma:D^{\prime}_{1}\to D^{\prime}_{2}$ for $D^{\prime}_{1}\subseteq D_{1}$ and $D^{\prime}_{2}\subseteq D_{2}$. For $\delta>0$, we define $\mathrm{cost}(\sigma):=\max(A,B)$, where

[TABLE]

and it is understood that $\infty-\infty=0$. We define the bottleneck distance $d_{B}$ by

[TABLE]

For persistence modules $M$ and $N$, we write $d_{B}(D(M),D(N))$ simply as $d_{B}(M,N)$.

Lemma 2.3**.**

Let $P^{M}$ and $P^{N}$ be finite presentations of persistence modules $M$ and $N$, respectively. $d_{B}(M,N)$ is equal to one of the following:

The difference of a grade of $P^{M}$ and a grade of $P^{N}$, 2. 2.

half the difference of two grades in $P^{M}$, 3. 3.

half the difference of two grades in $P^{N}$.

Proof 2.4**.**

This follows immediately from Lemma 2.1 and the definition of $d_{B}$.

Given two finite persistence diagrams $D$, $D^{\prime}$ with $m$ points each, we can compute $d_{B}(D,D^{\prime})$ in time to $O(m^{1.5}\log m)$ [12]; see [14] for details, including a report on practical efficiency. The number of points in a diagram is at most the number of generators of the corresponding module, and hence upper-bounded by $n$. Thus, the complexity of computing the bottleneck distance of two persistence modules is dominated by the computation of the persistence diagrams, and has worst-case complexity $O(n^{\omega})$.

Bipersistence modules

The definitions of persistence modules and presentations extend to higher dimensions without difficulty. In the 2-parameter setting, this goes as follows: Define a partial order $\leq$ on $\mathbb{R}^{2}$ by $p\leq q$ if $p_{x}\leq q_{x}$ and $p_{y}\leq q_{y}$. A bipersistence module is an assignment of $\mathbb{K}$-vector spaces $M_{p}$ to points $p\in\mathbb{R}^{2}$, and linear maps $M_{p\to q}:M_{p}\to M_{q}$ to pairs of points $p\leq q\in\mathbb{R}^{2}$, such that $M_{p\to p}$ is the identity and $M_{q\to r}\circ M_{p\to q}=M_{p\to r}$ whenever $p\leq q\leq r$. In topological data analysis, $2$-dimensional persistence modules typically arise by applying homology to a bifiltered simplicial complex

A finite presentation of a bipersistence module $\mathbb{R}^{2}$ is defined in the same way as for one-parameter persistence modules, except that the labels of each row/column are elements of $\mathbb{R}^{2}$, and the $\leq$ relation appearing in the definition means the partial order over $\mathbb{R}^{2}$. From now on, we will assume that all bipersistence modules considered are finitely presented.

Example**.**

Let $M$ be the bipersistence module given by

[TABLE]

Then a presentation of $M$ is given by

[TABLE]

In topological data analysis, we do not typically have immediate access to a presentation of a bipersistence module $M$, but rather to a chain complex of bipersistence modules for which $M$ is a homology module. However, it has recently been observed that from such a chain complex, a (minimal) presentation of $M$ can be computed in cubic time [18]. The algorithm for this is practical, and has been implemented in the software package RIVET [9].

3 The matching distance

Slices

We define a slice as a line $\ell:y=sx+t$ where $s$ and $t$ are real numbers with $s>0$. Let $\lambda:\mathbb{R}\to\ell$ be an order-preserving, isometric parameterization of the slice. Concretely, such a parameterization is given by $\lambda(x)=\frac{1}{\sqrt{1+s^{2}}}(x,sx+t)$. Given a bipersistence module $M$ and slice $\ell$, we define a (1-parameter) persistence module $M^{\ell}$ via $M^{\ell}_{x}:=M_{\lambda(x)}$, with its linear maps induced by $M$. We call $M^{\ell}$ a slice module.

Matching distance

For a slice $\ell:y=sx+t$, we define a weight

[TABLE]

$w(\ell)$ is maximized for slices with slope $1$, and gets smaller when the slope goes to [math] or to $\infty$.

Let $\Lambda$ denote the set of all slices. Writing $\Omega:=(0,\infty)\times\mathbb{R}$, the map $\alpha:\Omega\to\Lambda$ that sends $(s,t)$ to the line $y=sx+t$ is clearly a bijective parameterization of $\Lambda$. We equip $\Lambda$ with the topology induced by $\alpha$, using the subset topology $\Omega\subset\mathbb{R}^{2}$.

For two persistence modules $M$, $N$ over $\mathbb{R}^{2}$, we define a function $F^{M,N}:\Lambda\to[0,\infty)$ by

[TABLE]

and we define the matching distance between $M$ and $N$ as $d_{\mathrm{match}}(M,N):=\sup F^{M,N}$. As noted in the introduction, the weights $w(\ell)$ are chosen to ensure that $d_{\mathrm{match}}$ is a lower bound for the interleaving distance.

Lemma 3.1**.**

Given two bipersistence modules $M$, $N$, the map $F^{M,N}$ is continuous.

Proof 3.2**.**

$w$* is clearly continuous, so it suffices to show that the function $\ell\mapsto d_{B}(M^{\ell},N^{\ell})$ is continuous. Let $\mathcal{D}$ denote the metric space of all finite persistence diagrams, with metric the bottleneck distance. It follows from [15, Theorem 2] that the map sending a slice $\ell$ to the persistence diagram $D(M^{\ell})$ is continuous with respect to the topology on $\mathcal{D}$. Thus the map sending $\ell$ to the pair $(D(M^{\ell}),D(N^{\ell}))$ is also continuous. Moreover, the bottleneck distance is clearly continuous as a map $\mathcal{D}\times\mathcal{D}\to[0,\infty)$, thanks to the triangle inequality. Thus, $\ell\mapsto d_{B}(M^{\ell},N^{\ell})$ is a composition of continuous functions, hence continuous.*

4 The arrangement

In what follows, we fix two bipersistence modules $M$, $N$ (each with at most $n$ generators and relations) and write $F:=F^{M,N}$. Using the parameterization $\alpha$ from above, we also have a function $F\circ\alpha:\Omega\to[0,\infty)$; slightly abusing notation, we will also denote this map by $F$.

In this section, we construct a line arrangement in $\Omega$ in such a way that it is simple to compute $\sup F$ on each face. Recall that a line arrangement of $\Omega$ is the subdivision of $\Omega$ into vertices, edges, and faces induced by a finite set of distinct lines $L_{1},\ldots,L_{n}$. The vertices of the arrangement are the intersection points of (at least) two lines, the edges are maximal connected subsets of lines not containing any vertex, and the faces are the connected components of $\Omega\setminus\bigcup_{i=1}^{n}L_{i}$. Clearly, each vertex, edge, and face of the arrangement is a convex set. The boundary of each face consists of a finite number of edges and vertices.

A first line arrangement

For $v\in\mathbb{R}^{2}$, let $L_{v}$ denote the line $y=-v_{x}\,x+v_{y}$. Note that $L_{v}\cap\Omega$ is exactly the set of parameterizations of slices containing $v$. Now, fix presentations $P^{M}$ and $P^{N}$ of $M$ and $N$. Let $\mathcal{A}_{0}$ denote the arrangement in $\Omega$ induced by the set of lines

[TABLE]

In what follows, we will refine $\mathcal{A}_{0}$ by adding more lines into the arrangement. For this we first need to introduce some definitions.

Pushes

For a point $p=(p_{x},p_{y})\in\mathbb{R}^{2}$ and a slice $\ell:y=sx+t$, we define the push of $p$ onto $\ell$ as

[TABLE]

Geometrically, $\mathrm{push}(p,\ell)$ gives the intersection point of $\ell$ and a vertical upward ray emanating from $p$ in the first case, and the intersection of $\ell$ with a horizontal right ray emanating from $p$ in the second case. See Figure 1 for an illustration.

Remark 4.1**.**

A finite presentation of $M$ induces a finite presentation of $M^{\ell}$ with the same underlying matrix, and each row or column grade $p\in\mathbb{R}^{2}$ replaced with $\lambda^{-1}(\mathrm{push}(p,\ell))$. Clearly, this presentation can be obtained in time linear in the size of the presentation of $M$.

For $p,q\in\mathbb{R}^{2}$, define $\delta_{p,q}:\Omega\to[0,\infty)$ by

[TABLE]

with $\ell$ the slice defined by $s$ and $t$. Again, see Figure 1 for an illustration.

We now give piecewise analytic formulae for $\delta_{p,q}$, which depend on whether the slice $\ell$ is below or above $p$ and $q$.

(I) slice is below both $p$ and $q$:

[TABLE]

(II) slice is above both $p$ and $q$:

[TABLE]

(III) slice is between $p$ and $q$:

There are two subcases, which we will call (IIIa) and (IIIb): Assuming $p$ lies above the slice (IIIa), the formula is

[TABLE]

If $p$ lies below the slice (IIIb), the formula is the same, with the roles of $p$ and $q$ exchanged.

The push map is easily seen to be continuous with respect to the slice $\ell$, so these formulae also extend to boundaries of the cases, i.e., when the slice contains $p$ or $q$.

Lemma 4.2**.**

If $p,q\in\mathrm{gr}(P^{M})\cup\mathrm{gr}(P^{N})$, then in each face of $\mathcal{A}_{0}$, exactly one of the conditions (I), (II), (IIIa), (IIIb) holds everywhere. Hence, $\delta_{p,q}$ can be expressed on the entire face by one of the analytic formulae above.

Proof 4.3**.**

Assume for a contradiction that two points $x$, $x^{\prime}$ in a face $f$ of $\mathcal{A}_{0}$ are of two different types among (I), (II), (IIIa), (IIIb). Let $\ell=\alpha(x)$, $\ell^{\prime}=\alpha(x^{\prime})$ denote the corresponding slices. Then either $p$ or $q$ (or both) switch sides from $\ell$ to $\ell^{\prime}$. Assume that $p$ lies above $\ell$ and $p$ lies below $\ell^{\prime}$ – all other cases are analogous. Then, on the line segment from $x$ to $x^{\prime}$, there is some point $x^{\prime\prime}$ such that $p$ lies on the slice $\alpha(x^{\prime\prime})$. Hence, $x^{\prime\prime}$ lies on $L_{p}$ and is therefore not in $f$. This contradicts the convexity of the faces.

In view of Lemma 4.2, for $p,q\in\mathrm{gr}(P^{M})\cup\mathrm{gr}(P^{N})$ we may define the type of $\delta_{p,q}$ on a face of $\mathcal{A}_{0}$ to be the case (I), (II), (IIIa), or (IIIb) which holds on that face.

Refinement of the arrangement

Now we further subdivide the arrangement $\mathcal{A}_{0}$. For that, consider the set of equations of the form

[TABLE]

where

[TABLE]

Lemma 4.4**.**

The solution set of each of the above equations restricted to a face $f$ is either the empty set, the entire face, the intersection of $f$ with a line, or intersection of $f$ with the union of two lines.

Proof 4.5**.**

First we show the statement for equations of type $\delta_{p,q}(s,t)=0$. There are $3$ cases:

$\delta_{p,q}$** is of type (I):**

the equation becomes

[TABLE]

for which either all $(s,t)\in f$ are a solution (if $p_{y}=q_{y}$), or no $(s,t)$ is a solution.

$\delta_{p,q}$** is of type (II):**

the same argument holds for the equation

[TABLE]

$\delta_{p,q}$** is of type (III):**

Swapping $p$ and $q$ if necessary, we obtain the equation

[TABLE]

and the solution set is made of all $(s,t)\in f$ for which $p_{y}-t-sq_{x}=0$, which is the equation of a line.

For the remaining equations, we give the proof in the special case that $c_{pq}=c_{p^{\prime}q^{\prime}}$; the proof in the other cases is essentially the same. For equations of the form $\delta_{p,q}(s,t)=\delta_{p^{\prime},q^{\prime}}(s,t)$, there are six cases to check, depending on the type of the $\delta$-functions on the left and right sides of the equation:

Both $\delta_{p,q}(s,t)$ and $\delta_{p^{\prime},q^{\prime}}(s,t)$ are of type (I):

the equation is

[TABLE]

and the equation is satisfied if and only if $|p_{y}-q_{y}|=|p^{\prime}_{y}-q^{\prime}_{y}|$, independent of $s$ and $t$. Hence, the solution set is either $f$ or $\emptyset$.

Both $\delta_{p,q}(s,t)$ and $\delta_{p^{\prime},q^{\prime}}(s,t)$ are of type (II):

the same argument as in the previous case applies, so the solution set is either $f$ or $\emptyset$.

$\delta_{p,q}$** is of type (I) and $\delta_{p^{\prime},q^{\prime}}$ of type (II):**

we get the equation

[TABLE]

Since $1+\frac{1}{s^{2}}=\frac{1+s^{2}}{s^{2}}$, this simplifies to

[TABLE]

and the solution set is either all of $f$ (if both absolute values vanish), the empty set (if only $p^{\prime}_{x}-q^{\prime}_{x}=0$), or the intersection of $f$ with the vertical line $s=\frac{|p_{y}-q_{y}|}{|p^{\prime}_{x}-q^{\prime}_{x}|}$ (otherwise).

Both $\delta_{p,q}$ and $\delta_{p^{\prime},q^{\prime}}$ are of type (III):

Swapping $p,q$ or $p^{\prime},q^{\prime}$ if necessary, we get

[TABLE]

hence, $(s,t)$ is a solution if and only if $p_{y}-t-sq_{x}=p^{\prime}_{y}-t-sq^{\prime}_{x}$ or $p_{y}-t-sq_{x}=-(p^{\prime}_{y}-t-sq^{\prime}_{x})$. The first equations yields again either $f$, $\emptyset$, or a vertical line as solution set, the second equation always defines a line. Exchanging the roles of $p$ and $q$, or the roles of $p^{\prime}$ and $q^{\prime}$, or both, does not change the conclusion.

$\delta_{p,q}$** is of type (I) and $\delta_{p^{\prime},q^{\prime}}$ is of type (III):**

Swapping $p^{\prime}$ and $q^{\prime}$ if necessary, the formula is

[TABLE]

$(s,t)\in f$* is a solution if $p_{y}-t-sq_{x}=p_{y}-q_{y}$ or $p_{y}-t-sq_{x}=q_{y}-p_{y}$, which is a line equation in both cases.*

$\delta_{p,q}$** is of type (II) and $\delta_{p^{\prime},q^{\prime}}$ is of type (III):**

Swapping $p^{\prime}$ and $q^{\prime}$ if necessary, we get

[TABLE]

which simplifies to

[TABLE]

Here $(s,t)\in f$ is a solution if and only if $s(p_{x}-q_{x})=p^{\prime}_{y}-t-sq^{\prime}_{x}$ or $s(p_{x}-q_{x})=-(p^{\prime}_{y}-t-sq^{\prime}_{x})$. Again, we obtain a line in both cases.

Definition 4.6**.**

Let $\mathcal{A}$ denote the line arrangement in $\Omega$ formed by the lines in $\mathcal{A}_{0}$, all lines from the case analysis above, and the vertical line $s=1$.

Lemma 4.7**.**

The arrangement $\mathcal{A}$ consists of $O(n^{4})$ lines.

Proof 4.8**.**

The case analysis in the proof of Lemma 4.4 was performed relative to a choice of face $f$ in $\mathcal{A}_{0}$. However, for a fixed choice of an equation in the set of equations (1), the lines which arise in the case analysis depend only on the types of $\delta_{p,q}$ and $\delta_{p^{\prime},q^{\prime}}$ on $f$. There are at most $4\times 4=16$ possible ways of jointly choosing those types, and for a given choice, at most two lines are added to the arrangement. Hence, each of the $O(n^{4})$ equations in the set of equations (1) contributes at most a constant number of lines to $\mathcal{A}$. The result now follows easily.

Note that the arrangement $\mathcal{A}$ depends on the choice of presentations for $M$ and $N$. The next statement says that within a face of the arrangement, the bottleneck distance is realized by the difference of the pushes of two fixed grades of the presentations.

Theorem 4.9**.**

For any face $f$ of $\mathcal{A}$, there is some choice of $p,q\in\mathrm{gr}(P^{M})\cup\mathrm{gr}(P^{N})$ such that $d_{B}(M^{\ell},N^{\ell})=c_{pq}\delta_{p,q}(\ell)$ for all $\ell\in\alpha(f)$.

The formal proof of this result is deferred to Appendix A. It can be summarized as follows: By Remark 4.1 and Lemma 2.1, for each $\ell\in f$ there is a collection $T_{M}^{\ell}$ of pairs

[TABLE]

such that $D(M^{\ell})$ is obtained by pushing the elements of $T_{M}^{\ell}$ onto $\ell$. Call such $T_{M}^{\ell}$ a template for $D(M^{\ell})$. For $\ell$ and $\ell^{\prime}$ such that the grades of $P^{M}$ push onto $\ell$ and $\ell^{\prime}$ in the same order, the templates for $D(M^{\ell})$ and $D(M^{\ell^{\prime}})$ are the same. But in fact, the grades of $P^{M}$ push onto all lines $\ell\in f$ in the same order, because whenever the order changes, we need to cross one of the lines of the arrangement $\mathcal{A}$. Hence, there exists $T_{M}^{f}$ that is a template for $D(M^{\ell})$, for all $\ell\in f$. Similarly for $N$.

By Lemma 2.3, for any fixed $\ell^{\prime}\in f$ we have that $d_{B}(M^{\ell^{\prime}},N^{\ell^{\prime}})=c_{pq}\delta_{p,q}(\ell^{\prime})$, for $p$ and $q$ each some coordinate of a pair in $T_{M}^{f}\cup T_{N}^{f}$. The order of the values of the functions

[TABLE]

remains the same across $f$, because any change in the order will result again in crossing one of the lines of the arrangement $\mathcal{A}$. Thus, in fact $d_{B}(M^{\ell},N^{\ell})=c_{pq}\delta_{p,q}({\ell})$ for all $\ell\in f$.

5 Maximization

We define a region of $\mathcal{A}$ as the closure of a face of $\mathcal{A}$ within $\Omega$. We can compute the matching distance by determining $\sup F(s,t)$ separately for each region in $\mathcal{A}$. We will show now that in each region, $\sup F(s,t)$ is either realized at a boundary vertex, or as the limit of an unbounded boundary edge, which can be computed easily.

We fix the following notation: A region $R\subseteq\Omega$ is an inner region if it is bounded as a set in $\mathbb{R}^{2}$ and it has a positive distance to the vertical line $s=0$ in $\mathbb{R}^{2}$ (in other words, $R$ does not reach the boundary of $\Omega$). An inner region is a convex polygon. Regions that are not inner region are called outer regions. Outer regions have exactly two outer segments in their boundary, which are infinite or converge to a point on the vertical line $s=0$ (except for the empty arrangement, which only consists of one outer region). See Figure 2 for an illustration.

For a fixed region $R$ of $\mathcal{A}$, Theorem 4.9 ensures that there is a pair of grades $(p,q)$ whose $\delta$-function realizes $d_{B}$ within the interior of $R$ (a face of $\mathcal{A}$). By continuity, this implies that $\delta_{p,q}$ also realizes $d_{B}$ on the entire region.

Lemma 5.1**.**

The supremum of $F$ within $R$ is attained either at a vertex on the boundary of $R$, or as the limit of $F$ along an outer segment. In the latter case, the limit can be expressed in simple terms based on the equation of the line segment and one of the functions $\delta_{p,q}$.

Proof 5.2**.**

We distinguish $6$ cases based on the type of the $\delta$-function and on whether $s\leq 1$ or $s\geq 1$ (note that each cell belongs to one of these cases, because the line $s=1$ is in $\mathcal{A}$).

$\delta_{p,q}$** of type (I), $s\leq 1$**

In that case,

[TABLE]

a constant function. Clearly, the supremum is attained everywhere, in particular at the boundary vertices of $R$.

$\delta_{p,q}$** of type (I), $s\geq 1$**

We get

[TABLE]

Clearly, this function becomes larger when $s$ gets smaller. Moreover, because $s\geq 1$ within the cell, there is a leftmost vertex on the boundary, which minimizes $s$ and therefore attains the supremum within the cell.

$\delta_{p,q}$** of type (II), $s\leq 1$**

We obtain

[TABLE]

Similarly to the previous case, there exists a rightmost boundary vertex in the cell (because $s\leq 1$), which realizes the supremum.

$\delta_{p,q}$** of type (II), $s\geq 1$**

The function simplifies to

[TABLE]

a constant function, which attains its supremum at any boundary vertex.

$\delta_{p,q}$** of type (III), $s\leq 1$**

Assuming that $p$ lies above the slice, we get

[TABLE]

If $R$ is an inner region, we are maximizing the above function over a closed convex polygon, and the maximum is achieved at a boundary vertex, because $|p_{y}-t-sq_{x}|$ is the maximum of two linear functions in $s$ and $t$.

It remains to analyze the case that $R$ is an outer region. We argue first that $R$ is bounded in $t$-direction from above and below: Since $\delta_{p,q}$ is of type (III), with $p$ lying above $\ell$, $(s,t)$ must be below the (non-vertical) line $t=-sp_{x}+p_{y}$ in the dual space. Likewise, since $q$ is below $\ell$, $(s,t)$ must be above $t=-sq_{x}+q_{y}$. Moreover, we have $0<s\leq 1$. If the above lines intersect at a point $r$ with $s$-value in $(0,1)$, $R$ is contained in the triangle spanned by the two lines and the vertical line $s=0$. Otherwise, $R$ is contained in the trapezoid induced by these two lines and the vertical lines $s=0$ and $s=1$.

It follows that the two outer segments of $R$ converge to the vertical line $s=0$. Let $(0,t_{1})$ denote the limit of the lower outer segment and $(0,t_{2})$ the limit of the upper outer segment. Clearly $t_{1}\leq t_{2}$. Let $\bar{R}$ denote the union of $R$ with the vertical line segment from $(0,t_{1})$ to $(0,t_{2})$; note that $\bar{R}$ is the closure of $R$ considered as a subset of $\mathbb{R}^{2}$. Observe that $|p_{y}-t-sq_{x}|$ is continuous over $\mathbb{R}^{2}$; therefore $F$ can be continuously extended to $\bar{R}$. It follows that the supremum of $F$ over $\bar{R}$ is attained at a boundary vertex, since $\bar{R}$ is a convex closed polygon. There are two cases: either the maximum is attained at a vertex of $\mathcal{A}$, or at $(0,t_{1})$ or $(0,t_{2})$. As we can readily see, the function values at the latter two points are $|p_{y}-t_{1}|$ and $|p_{y}-t_{2}|$, respectively. The case where $p$ is below the slice and $q$ is above is analyzed in the same way, with the roles of $p$ and $q$ swapped.

$\delta_{p,q}$** of type (III), $s\geq 1$**

Assuming that $p$ lies above the slice, we get

[TABLE]

We first consider the case where $R$ is an inner region, and we show that the function is maximized at a boundary vertex of $R$. The function $|\frac{p_{y}}{s}-\frac{t}{s}-q_{x}|$ has no local maximum over $\mathbb{R}^{2}$ since it is the absolute value of a linear function in $t$ for any fixed $s$. Hence, the supremum over $R$ must be attained on the boundary. We have to exclude the case that the maximum lies in the interior of an edge. For vertical edges this is obvious, because for a constant $s$, the function simplifies to the absolute value of a linear function in $t$ which must be maximized at a boundary vertex. For a non-vertical line of the form $t=as+b$, plugging in this equation for $t$ yields a function of the form

[TABLE]

This is the absolute value of a monotone function in $s$ and hence has no local maximum. Again, this implies that it is maximized at a boundary vertex.

Consider now the case where $R$ is an outer region. As in the previous case, $R$ is upper and lower bounded by two non-vertical lines, because we assume type (III). Hence, the two outer segments of $R$ cannot be vertical; the lower outer segment has a slope $r_{1}$ and the upper outer segment has a slope $r_{2}$ with $r_{1}<r_{2}$. We argue next that the supremum of $\frac{p_{y}}{s}-\frac{t}{s}-q_{x}$ is either attained at a boundary vertex, or equal to $|r_{1}+q_{x}|$, or equal to $|r_{2}+q_{x}|$. Let $(s_{i},t_{i})$ denote a sequence of points in $R$ such that $F(s_{i},t_{i})$ converges to the supremum. If $(s_{i},t_{i})$ converges to a point in $R$ (or has at least a convergent subsequence), it follows (similarly to the case of an inner region) that the limit point is a boundary vertex. Otherwise, we can assume (by passing to a subsequence) that the sequence $s_{i}$ is unbounded. Moreover, the sequence $\frac{t_{i}}{s_{i}}$ is bounded by $[r_{1},r_{2}]$ and therefore has a convergent subsequence with limit $r^{\prime}$. Passing to this subsequence, we obtain that

[TABLE]

Hence, the supremum must be of the form $|r^{\prime}+q_{x}|$ for some $r^{\prime}\in[r_{1},r_{2}]$. On the other hand, this expression is clearly maximized for either $|r_{1}+q_{x}|$ or $|r_{2}+q_{x}|$, and there exist sequences attaining these values, for instance when choosing $(s_{i},t_{i})$ on either of the outer segments. The case where $p$ lies below the slice and $q$ lies above is treated similarly, with the roles of $p$ and $q$ exchanged.

The algorithm

We now give the algorithm to compute the matching distance:

• •

Compute the arrangement induced by $\mathcal{A}$ from Definition 4.6.

• •

For each vertex $(s,t)$ in the arrangement, compute $F(s,t)$. Let $m$ be the maximum among all the values.

• •

For each outer region $R$, pick a point $(s,t)$ in the interior. Compute the bottleneck distance and identify a pair $(p,q)$ of grades that realizes the bottleneck. Determine whether $p$ and $q$ are above or below the slice $(s,t)$. If the region is of type (IIIa) with respect to $p$ and $q$, then do the following:

• –

If $R$ is to the left of $s=1$, compute the intersections $(0,t_{1})$, $(0,t_{2})$ of the outer segments of $R$ with the vertical line $s=0$. Set $m\leftarrow\max\{m,\ |p_{y}-t_{1}|,\ |p_{y}-t_{2}|\}$.

• –

If $R$ is to the right of $s=1$, let $r_{1}$, $r_{2}$ denote the slopes of the outer segments of $R$. Set $m\leftarrow\max\{m,\ |r_{1}+q_{x}|,\ |r_{2}+q_{x}|\}$.

If the region is of type (IIIb), proceed analogously, with the roles of $p$ and $q$ exchanged.

• •

Return $m$.

By “computing the arrangement”, we mean to store the planar subdivision induced by the lines of the arrangement (e.g.[2, Ch.2]). In fact, it is possible to implement the algorithm without explicitly constructing the line arrangement $\mathcal{A}$ or even storing its entire set of vertices. This reduces the space complexity of the algorithm. See Appendix B for the details.

Theorem 5.3**.**

The above algorithm computes the matching distance in polynomial time.

Proof 5.4**.**

Correctness follows from Lemma 5.1: as we check all vertices of the arrangement, we cover the supremum of all inner regions. The outer regions are handled separately in the last steps of the algorithm.

Running time: recall from Lemma 4.7 that we have $O(n^{4})$ lines in the arrangement $\mathcal{A}$. Hence, the arrangement has $O(n^{8})$ vertices, $O(n^{4})$ outer regions, and can be computed in $O(n^{8}\log n)$ time using an extension of the Bentley-Ottman sweep-line algorithm [1]. For each vertex and each outer region, we have to compute two persistence diagrams, which can be done in $O(n^{3})$ time, and a bottleneck distance whose complexity can be neglected. The remaining computations are negligible. Hence, we arrive at a $O(n^{11})$ algorithm.

We remark that the algorithm can be realized entirely with rational arithmetic if all grades are rational numbers. Indeed, all lines in the arrangement have rational coefficients, and so do their intersection points. An intersection points corresponds to a slice along which we are required to compute the bottleneck distance. Recall from Section 3 that the definition of the slice $\ell$ module introduces a grade of $\lambda^{-1}(\mathrm{push}(p,\ell))$ where $\lambda^{-1}(p_{x},p_{y})=\sqrt{1+s^{2}}p_{x}$. Hence, these grades are not rational numbers. However, the bottleneck distance is multiplied with the weight $w(\ell)$ of the slice afterwards, and instead of doing so, one can as well scale all grades with $w(\ell)$ in advance. A simple calculation shows that this indeed turns the grades into rational values.

A simple analysis also reveals that if the input coordinates are rational and of bitsize $\leq b$, all intermediate computations in the algorithm can be performed with a bitsize of $\leq cb$, with $c$ a (small) constant. Hence, the algorithm is strongly polynomial.

6 Discussion

We have presented the first polynomial time algorithm to exactly compute the matching distance for 2-parameter persistence modules. It is natural to ask about practicality of our approach. The large exponent of $n^{11}$ seems discouraging at first, but we mention first that the worst-case running time of $O(n^{3})$ for persistent homology is usually not appearing for real instances; indeed an almost linear behavior can be expected. Still, the large number of $O(n^{4})$ lines in the arrangement constitutes a computational barrier in practice. There are several possibilities, however, to reduce this effect:

• •

Instead of computing the arrangement $\mathcal{A}$ globally, we could compute the intermediate arrangement $\mathcal{A}_{0}$ and refine each face of it separately, using only those lines that affect the $\delta$-functions within this face.

• •

As a follow-up to the previous point, it might be possible to compute a smaller arrangement per face adaptively. The idea is to start at some interior point in a face of $\mathcal{A}_{0}$, identifying a pair $(p,q)$ that realizes the bottleneck distance and then to determine the boundary of the region where $(p,q)$ realizes the bottleneck distance.

• •

As a preprocessing step, we can minimize the input presentations, yielding presentations for the same modules with the smallest possible number of generators and relations (hence minimizing $n$) [18].

We pose the question of whether an implementation realizing the above ideas is competitive to an approximative, sampling-based approach for computing the matching distance.

Our algorithm needs to treat the outer edges of the arrangement $\mathcal{A}$ separately since our analysis does not rule out the possibility that the supremum is realized at the boundary of $\Omega$. On the other hand, we are not aware of an example of two finite presentations whose matching distance is not realized by a particular slice in $\Omega$. A proof that the supremum in the definition of the matching distance is in fact a maximum would greatly simplify our algorithm, since it would boil down to computing the intersection points of all lines and searching for the maximal $F$-value among them.

Finally, we have restricted attention to the case of two-parameter persistence modules. It is natural to conjecture that our algorithm extends to more parameters by constructing a hyperplane arrangement. It would be worthwhile to check this conjecture in future work.

A Appendix: proof of Theorem 4.9

Let $G^{M}=\mathrm{gr}(P^{M})$ and let $G^{N}=\mathrm{gr}(P^{N})$. As an intermediate result, we prove that the combinatorial structure of each persistence diagram stays the same across the face $f$:

Lemma A.1**.**

For each face $f$ of $\mathcal{A}$, there exist multisets

[TABLE]

such that for any $\ell\in\alpha(f)$,

[TABLE]

where by convention $\mathrm{push}((\infty,\infty),\ell):=(\infty,\infty)$ and $\lambda^{-1}((\infty,\infty)):=\infty$. We call $\mathcal{T}_{M}^{f}$ and $\mathcal{T}_{N}^{f}$ diagram templates.

Proof A.2**.**

This is a variant of [17, Theorem 4.1]. We give a succinct, algorithmically flavored proof here. We prove the result for $M$; the proof for $N$ is the same.

First, note that the restriction of the partial order on $\mathbb{R}^{2}$ to any $\ell\in\alpha(f)$ is a total order. This total order pulls back under the map $\mathrm{push}(-,\ell):\mathrm{gr}(P^{M})\to\ell$ to a totally ordered partition of $\mathrm{gr}(P^{M})$, (i.e., elements of the partition are level sets). It can be checked that, because $\mathcal{A}$ refines $\mathcal{A}_{0}$ and also contains all lines of the form $\delta_{p.q}=0$ for $p,q\in\mathrm{gr}(P^{M})$ with $p$ and $q$ incomparable in the partial order on $\mathbb{R}^{2}$, this totally ordered partition is the same for all $\ell\in\alpha(f)$; see [17, Corollary 3.4]. Thus, we obtain a unique totally ordered partition of $\mathrm{gr}(P^{M})$ associated to all of $f$. Let us refine this to a fixed total order on $\mathrm{gr}(P^{M})$.

Permuting the row or columns of a presentation for $M$ yields another presentation for $M$, so we may assume without loss of generality that the order of rows and columns for $P^{M}$ is consistent with our total order on $\mathrm{gr}(P^{M})$. Applying the matrix reduction underlying the standard persistence algorithm [21] to the matrix underlying $P^{M}$ yields a matching $\sigma$ of row and column indices, where a non-zero column in the reduced matrix is matched to the row of its pivot. We may define

[TABLE]

By Remark 4.1 and the correctness of the standard algorithm for computing persistent homology [21], it follows that

[TABLE]

as desired.

Suppose we are given two persistence diagrams $D$ and $D^{\prime}$. For each $(x,y)\in D\cup D^{\prime}$, let $w_{0}(x,y):=\frac{y-x}{2}$, and for each $(x,y),(x^{\prime},y^{\prime})\in D\times D^{\prime}$, let $w_{1}(x,x^{\prime}):=\max(|x-x^{\prime}|)$, and $w_{2}(y,y^{\prime}):=\max(|y-y^{\prime}|)$. We call these (possibly infinite) numbers the weights of the pair $(D,D^{\prime})$. Denoting the set of pairs indexing these weights as $I(D,D^{\prime})$, the weights define a function $w^{D,D^{\prime}}:I^{D,D}\to[0,\infty)$. Note that $d_{B}(D,D^{\prime})=w^{D,D^{\prime}}(x,y)$ for some $(x,y)\in I^{D,D^{\prime}}$.

Clearly, a pair of bijections of persistence diagrams $\zeta_{1}:D\to E$ and $\zeta_{2}:D^{\prime}\to E^{\prime}$ induces a bijection

[TABLE]

The following is an easy consequence of the definition of $d_{B}$.

Lemma A.3**.**

If $d_{B}(D,D^{\prime})=w^{D,D^{\prime}}(x,y)$ and $\zeta_{1}$, $\zeta_{2}$ preserve the order of the weights, in the sense that

[TABLE]

then $d_{B}(E,E^{\prime})=w^{E,E^{\prime}}(\zeta(x,y))$.

Proof A.4** (Proof of Theorem 4.9).**

Let $\mathcal{T}_{M}^{f}$ and $\mathcal{T}_{N}^{f}$ be diagram templates for the cell $f$. For any $\ell\in\alpha(f)$, Lemma A.1 gives us distinguished bijections

[TABLE]

It is easily checked that for $(x,y)\in D(M^{\ell})$, and $(p,q)=\gamma^{-1}(x,y)$, we have

[TABLE]

Similarly, for

[TABLE]

we have

[TABLE]

Note that as we move the slice $\ell$ inside $f$, the order of the values taken by the functions

[TABLE]

cannot change. Indeed, the functions $\ell\to\delta_{p,q}(\ell)$ are continuous, and by Lemma 4.4 and the definition of the arrangement $\mathcal{A}$, the intersection of $f$ with the solution set of each equation $c_{pq}\delta_{p,q}=c_{p^{\prime}q^{\prime}}\delta_{p^{\prime},q^{\prime}}$ is either $f$ or $\emptyset$.

The diagram templates provide bijections $D(M^{\ell})\to D(M^{\ell^{\prime}})$ and $D(N^{\ell})\to D(N^{\ell^{\prime}})$ for all $\ell,\ell^{\prime}\in\alpha(f)$. It now follows from equations (2-4) that for each $\ell,\ell^{\prime}$, these preserve the order on weights. Thus, by Lemma A.3, if

[TABLE]

for some $\ell\in\alpha(f)$, then

[TABLE]

for all $\ell^{\prime}\in\alpha(f)$, where

[TABLE]

is the bijection induced by the diagram templates. The result now follows from equations (2-4).

B A more space-efficient algorithm

The algorithm described in Section 5 consists of two major steps: we first find the maximal value of $F$ over all intersection points of lines in the arrangement $\mathcal{A}$, and then check the convergence along outer segments for certain outer regions of $\mathcal{A}$. To decide whether an outer region needs to be considered, we require one interior point of that region, to decide the type of the region. The arrangement $\mathcal{A}$ (stored as a planar subdivision as in [2, Ch.2]) contains enough information to access all necessary data conveniently. However, its space complexity is $O(n^{8})$ because the arrangement is induced by $O(n^{4})$ lines. We show next how to implement the algorithm without storing $\mathcal{A}$ in memory, yielding a space complexity of $O(n^{4})$:

Reporting the intersection points of a set of lines in the plane is one of the oldest problems in computational geometry. The Bentley-Ottmann sweep-line algorithm [1] reports all intersection points in time proportional to the number of such points. Note that there are up to $O(n^{8})$ intersections, so storing all such points is space consuming. However, there is no need to do so – whenever the sweep-line algorithm reports a new point, we compute $F$ at this point and compare the value with the maximal $F$-value seen before, updating the maximum if the newly encountered value is larger. There is no need to store the intersection beyond this moment. In this way, we obtain the maximal $F$-value among all intersection points using only $O(n^{4})$ space.

For handling outer regions, note that we can extend the sweep-line algorithm such that it also returns the leftmost intersection point with positive $s$-coordinate, and the rightmost intersection point, among all pairs of lines in $\mathcal{A}$. Let $0<s_{min}\leq s_{max}$ denote the $s$-coordinates of these intersection points. Now, consider the vertical line $s=\frac{s_{min}}{2}$, which is partitioned into open intervals

[TABLE]

where $b_{1},\ldots,b_{m}$ are the intersection points of the non-vertical lines of $\mathcal{A}$ with the vertical line. Note that $m=O(n^{4})$. The intervals in the sequence can be computed in $O(n^{4}\log n)$ time, just by sorting the $b_{i}$’s. The intervals in the sequence (5) are in one-to-one correspondence with the outer regions of $\mathcal{A}$ whose boundary in $\mathbb{R}^{2}$ has nonempty intersection with the line $s=0$. More precisely, each interval $(b_{i},b_{i+1})$ is contained in a unique outer region $R$; for $i\in\{1,\ldots,m-1\}$, the point $(\frac{s_{min}}{2},\frac{b_{i}\,+\,b_{i+1}}{2})$ is an interior point of $R$; and the outer segments of $R$ contain the points $b_{i}$ and $b_{i+1}$. This information suffices to determine the type of $R$, and to compute the limit values of the $F$-function if necessary. The two infinite regions corresponding to $(-\infty,b_{1})$ and $(b_{m},\infty)$ can be ignored, because these regions must be of type (I) or type (II). The same construction provides interior points of outer regions that are unbounded in the $s$-direction, considering the vertical line at $s=2s_{max}$.

With this variant, we compute the matching distance with space complexity $O(n^{4})$. The time complexity remains $O(n^{11})$ as in the original algorithm, because we still iterate over $O(n^{8})$ vertices and compute the bottleneck distance in each position in $O(n^{3})$ time.