Chunk Reduction for Multi-Parameter Persistent Homology

TL;DR

This paper introduces a new algorithm for simplifying multi-parameter persistent homology computations by reducing input complex size while preserving homological features, improving efficiency and scalability.

Contribution

It extends the chunk algorithm to multi-parameter settings, producing the smallest quasi-isomorphic chain complex and enabling parallelization.

Findings

The algorithm produces minimal multi-filtered chain complexes.

Experimental results show improved performance over existing methods.

Parallelization scheme enhances scalability for large datasets.

Abstract

The extension of persistent homology to multi-parameter setups is an algorithmic challenge. Since most computation tasks scale badly with the size of the input complex, an important pre-processing step consists of simplifying the input while maintaining the homological information. We present an algorithm that drastically reduces the size of an input. Our approach is an extension of the chunk algorithm for persistent homology (Bauer et al., Topological Methods in Data Analysis and Visualization III, 2014). We show that our construction produces the smallest multi-filtered chain complex among all the complexes quasi-isomorphic to the input, improving on the guarantees of previous work in the context of discrete Morse theory. Our algorithm also offers an immediate parallelization scheme in shared memory. Already its sequential version compares favorably with existing simplification schemes, as we show by experimental evaluation.

Figures

Tables

Table 1. Table 1: Results and performances obtained by the chunk and the DMT-based algorithms. The columns from left to right indicate: the name of the dataset ( Dataset ), the number of cells before and after the simplification algorithms ( Size ), the time ( Time ), expressed in seconds, needed to perform the preparation step ( Prep. ) and the simplification step ( Simpl. ), the maximum peak of memory, expressed in GB, required by the two algorithms ( Memory Usage ).

\multirow3*Dataset	Size		Time (sec.)				Memory Usage (GB)
	\multirow2*Input	\multirow2*Output	Chunk		DMT		\multirow2*  Chunk	\multirow2*DMT
			Prep.	Simpl.	Prep.	Simpl.
Eros	2.9 M	202 K	1.7	0.8	2.7	15.8	0.36	0.46
Donna	3.0 M	217 K	1.8	0.8	2.8	16.9	0.38	0.48
Chinese Dragon	3.9 M	321 K	2.5	1.1	3.9	22.3	0.52	0.64
Circular Box	4.2 M	365 K	2.9	1.2	4.3	24.0	0.68	0.68
Ramesses	5.0 M	407 K	3.4	1.3	5.5	29.4	0.68	0.81
Pensatore	6.0 M	369 K	3.8	1.6	6.8	34.3	0.76	0.97
Raptor	6.0 M	260 K	4.4	1.7	5.4	32.3	0.73	0.93
Neptune	12.0 M	893 K	8.4	4.4	14.9	69.2	1.52	1.94
Cube 1	590 K	67 K	0.5	0.3	0.7	3.2	0.09	0.10
Cube 2	2.4 M	264 K	1.8	1.1	2.6	13.1	0.35	0.40
Cube 3	9.4 M	1.0 M	7.6	4.8	11.0	53.1	1.37	1.58
Cube 4	37.7 M	4.2 M	31.9	19.4	44.9	216.0	5.50	6.32

Equations

$$(p,k,\Delta )$$

$${\bar{\partial }}_l(c)=\{\begin{array}{ll}{\partial }_k(c)+\lambda ⟨c,c_2⟩{\partial }_k(c_1)& \text{ if }l=k,\\ {\partial }_{k+1}(c)-\lambda ⟨{\partial }_{k+1}(c),c_2⟩c_1& \text{ if }l=k+1,\\ {\partial }_l(c)& \text{ otherwise.}\end{array}$$

$$f_l(c)=\{\begin{array}{ll}c-\lambda ⟨c,c_2⟩c_1& \text{ if }l=k,\\ c& \text{ otherwise;}\end{array}$$

$$g_l(c)=\{\begin{array}{ll}c+\lambda ⟨c,c_2⟩c_1& \text{ if }l=k,\\ c& \text{ otherwise.}\end{array}$$

$${\bar{C}}_l=\{\begin{array}{ll}\{c\in C_k|⟨c,c_1⟩=0\}& \text{ if }l=k,\\ \{c\in C_{k+1}|⟨c,c_2⟩=0\}& \text{ if }l=k+1,\\ C_l& \text{ otherwise;}\end{array}$$

$${\bar{\partial }}_l(c)=\{\begin{array}{ll}{\partial }_{k+1}(c)-{\lambda }^{-1}⟨{\partial }_{k+1}(c),c_1⟩{\partial }_{k+1}(c_2)& \text{ if }l=k+1,\\ {\partial }_{k+2}(c)-⟨{\partial }_{k+2}(c),c_2⟩c_2& \text{ if }l=k+2,\\ {\partial }_l(c)& \text{ otherwise.}\end{array}$$

$$r_l(c)=\{\begin{array}{ll}c-{\lambda }^{-1}⟨c,c_1⟩{\partial }_{k+1}(c_2)& \text{ if }l=k,\\ c-⟨c,c_2⟩c_2& \text{ if }l=k+1,\\ c& \text{ otherwise;}\end{array}$$

$$s_l(c)=\{\begin{array}{ll}c-{\lambda }^{-1}⟨{\partial }_{k+1}(c),c_1⟩c_2& \text{ if }l=k+1,\\ c& \text{ otherwise;}\end{array}$$

$$D_{\ast }^{<p}:=\sum _{q<p}{D}_{\ast }^q,$$

$${\delta }_k^p(D):=\dim \ker {\eta }_{k-1}^p+\dim \mathrm{coker}{\eta }_k^p,$$

$${\gamma }_k^p(D):=\dim D_k^p-\dim D_k^{<p}.$$

$${\delta }_k^p(D)\le {\gamma }_k^p(D)$$

$${\delta }_k^p(C)={\delta }_k^p(D)$$

$${\bar{\partial }}_{k-1}{\bar{\partial }}_k(c)$$

$${\bar{\partial }}_k{\bar{\partial }}_{k+1}(c)$$

$${\bar{\partial }}_{k+1}{\bar{\partial }}_{k+2}(c)$$

$${\bar{\partial }}_k{f}_k(c)$$

$${\bar{\partial }}_{k+1}{f}_{k+1}(c)$$

$$g_{k-1}{\bar{\partial }}_k(c)$$

$$g_k{\bar{\partial }}_{k+1}(c)$$

$$g_k{f}_k(c)$$

$$f_k{g}_k(c)$$

$$⟨{\bar{\partial }}_{k+1}(c),c_1⟩$$

$$⟨{\bar{\partial }}_{k+2}(c),c_2⟩$$

$${\bar{\partial }}_k{\bar{\partial }}_{k+1}(c)$$

$${\bar{\partial }}_{k+1}{\bar{\partial }}_{k+2}(c)$$

$${\bar{\partial }}_{k+2}{\bar{\partial }}_{k+3}(c)$$

$${\bar{\partial }}_k{r}_k(c)$$

$${\bar{\partial }}_{k+1}{r}_{k+1}(c)$$

$${\bar{\partial }}_{k+2}{r}_{k+2}(c)$$

Full text

Chunk Reduction for Multi-Parameter Persistent Homology

Ulderico Fugacci

  

Michael Kerber

1 Introduction

In the last decades, topology-based tools are gaining a more and more relevant role in the analysis and in the extraction of the core information of unorganized, high-dimensional and potentially large datasets. Thanks to its capability of keeping track of the changes in the homological features of a dataset which evolves with respect to a parameter, persistent homology has represented a real game-changer in this field. Recently, an extension of persistent homology called multi-parameter persistent homology is drawing the attention of a growing number of researchers. In a nutshell, multi-parameter persistence generalizes the classic persistent homology by studying multivariate datasets which are filtered by two or more (independent) scale parameters. Multi-parameter persistent homology of a dataset cannot be captured by complete discrete invariant [5], but this has not prevented the researchers from defining several descriptors based on multi-parameter persistence [6, 13, 16]. Since these descriptors tend to have high algorithmic complexity, it is a natural pre-processing step to pass to a smaller but equivalent representation of the input complex and to invoke any demanding computation on this smaller representation.

Contribution.

Inspired by the chunk algorithm [3] for persistent homology, we present a reduction algorithm which for a filtered dataset returns a filtered chain complex having the same multi-parameter persistence but a drastically smaller size. Our approach is based on the observation that even in the absence of a global persistence diagram, local matrix reductions yield pairs of simplices which can be eliminated from the boundary matrix without affecting the homological information. Our approach proceeds in two steps, first identifying such local pairs of simplices, and in a second step manipulating the non-local columns of the boundary matrix such that all indices of locally paired simplices disappear. Both steps permit a parallel implementation with shared memory. We prove that our simplification scheme is optimal in the sense that any chain complex quasi-isomorphic to the input complex must contain at least as many generators as our output complex. Our algorithm yields similar time and space complexity bounds as its one-parameter counterpart. We implemented our algorithm, making use of various techniques that have proven effective for persistence computation, such as the twist reduction [8] and efficient data structure for the columns of a boundary matrix [4]. We experimentally show that our implementation is effective (see next paragraph). For the sake of clarity, our approach is described for the two-parameter case. No constraint prevents the generalization of the proposed method to an arbitrary number of parameters.

Comparison with related work.

Our work is motivated by a line of research on complex simplification using discrete Morse theory (DMT) [12, 1, 2, 18]. In these works, the idea is to build a discrete gradient locally and to return the resulting Morse complex on critical simplices as a simplification. In analogy to the one-parameter case, our chunk algorithm is an attempt to realize this simplification scheme using persistent homology instead of DMT. This gives more flexibility, as in DMT, the paired cells are constrained to be incident, while this restriction is not present for persistence pairs. Consequently, we are able to prove optimality of our output size in general, while DMT-based approaches only succeeded to give guarantees for special cases such as multi-filtrations of 3D regular grids and of abstract simplicial 2-complexes [17]. In practice, the theoretical benefit in terms of output size is rather marginal, as our experiments show; however, the timings show that our improved theoretical guarantee comes without performance penalty; on the contrary, our algorithm is always faster than the DMT-based approach presented in [18] on all tested examples. We remark, however, that the DMT-based algorithm returns a complex endowed with a Forman gradient as well as the corresponding discrete Morse complex, which can be of potential use for other application domains than computing persistent homology.

A related question is the computation of a minimal presentation of a persistence module induced by a simplicial or general chain complex. Roughly speaking, a presentation consists of a finite set of (graded) generators and relations, and the minimal presentation is one with the minimal possible number of generators and relations. Our algorithm does not yield a minimal presentation, however, since more computations are needed to identify generators of the chain complex as generators or as relations. An algorithm by Lesnick and Wright111https://www.ima.umn.edu/2017-2018/SW8.13-15.18/27428 [14] computes such a minimal presentation through matrix reduction in cubic time, carefully choosing a column order to reduce the number of reduction instances. Their algorithm is used in their software package RIVET [15]. Our contribution can serve as a pre-processing step for their algorithm, reducing the size of the matrix through efficient and possibly parallelized computation and invoking their global reduction step on a much smaller chain complex.

2 Background

Homology.

Fixed a base field $\mathbb{F}$ and a positive integer $d$, let us consider, for each $0\leq k\leq d$, a finite collection of elements denoted as $k$-generators (or, equivalently, generators of dimension $k$). A finitely generated chain complex $C_{*}=(C_{k},\partial_{k})$ over $\mathbb{F}$ is a collection of pairs $(C_{k},\partial_{k})$ where:

• •

$C_{k}$ is the $\mathbb{F}$-vector space spanned by the generators of dimension $k$,

• •

$\partial_{k}:C_{k}\rightarrow C_{k-1}$ is called boundary map and it satisfies the property that $\partial_{k-1}\partial_{k}=0$.

The elements of $C_{k}$ are called $k$-chains and, by definition, are $\mathbb{F}$-linear combinations of the generators of dimension $k$. The support of a $k$-chain is the set of $k$-generators whose coefficient in the chain is not zero. In the following, we will simply use the term chain complex in place of finitely generated chain complex. Moreover, we will assume that, given a chain complex, an explicit set of generators is provided.

Given a chain complex $C_{*}$, we denote as $Z_{k}:=\mathrm{ker}\,\partial_{k}$ the space of the $k$-cycles of $C_{*}$, and as $B_{k}:=\operatorname{Im}\,\partial_{k+1}$ the space of the $k$-boundaries of $C_{*}$. The $k^{th}$* homology space* of $C_{*}$ is defined as the vector space $H_{k}(C_{*}):=Z_{k}/B_{k}$. The rank $\beta_{k}$ of the $k^{th}$ homology space of a chain complex $C_{*}$ is called the $k^{th}$* Betti number* of $C_{*}$.

A chain map $f_{*}:C_{*}\to D_{*}$ is a collection of linear maps $f_{k}:C_{k}\to D_{k}$ which commutes with the boundary operators of $C_{*}$ and of $D_{*}$. A simple example of a chain map is the inclusion map, if $C_{*}$ is a subcomplex of $D_{*}$. In general, a chain map $f_{*}$ induces linear maps $H_{k}(C_{*})\to H_{k}(D_{*})$ for every $k$.

Chain complexes allow for capturing the combinatorial and the topological structure of a discretized topological space. Given a finite simplicial complex $K$, the chain complex $C_{*}$ associated to $K$ is defined by setting $C_{k}$ as the $\mathbb{F}$-vector space generated by the $k$-simplices of $K$ and the boundary $\partial_{k}(c)$ of a $k$-chain $c$ corresponding to a $k$-simplex $\sigma$ as the collection of the $(k-1)$-simplices lying on the geometrical boundary of $\sigma$. Consequently, homology of a finite simplicial complex $K$ is defined as the homology of the chain complex $C_{*}$ associated to $K$. Intuitively, homology spaces of $K$ reveal the presence of “holes” in the simplicial complex. The non-null elements of each homology space are cycles, which do not represent the boundary of any collection of simplices of $K$. Specifically, $\beta_{0}$ counts the number of connected components of $K$, $\beta_{1}$ its tunnels and holes, and $\beta_{2}$ the shells surrounding voids or cavities.

Multi-parameter persistent homology.

In the following, we will focus for simplicity on datasets filtered by two independent scale parameters. All definitions and results in this paper can be generalized to more parameters without problems. Let $p=(p_{x},p_{y}),q=(q_{x},q_{y})\in\mathbb{R}^{2}$, we will write throughout $p\leq q$ if $p_{x}\leq q_{x}$ and $p_{y}\leq q_{y}$. Given a chain complex $C_{*}=(C_{k},\partial_{k})$, let us assume to have an assignment which, for each generator $g\in C_{k}$, returns a value $v(g)\in\mathbb{R}^{2}$ such that, for any generator $g^{\prime}\in C_{k-1}$, if $\langle\partial_{k}g,g^{\prime}\rangle\neq 0$, then $v(g^{\prime})\leq v(g)$. The assignment $v$ can be extended to a function $v:C_{*}\rightarrow\mathbb{R}^{2}$ assigning to each chain $c$ the least common upper bound $v(c)$ of the values $v(g)$ of the generators in the support of $c$. In the following, $v(c)$ will be called the value of $c$.

Given a chain complex $C_{*}=(C_{k},\partial_{k})$ endowed with a value function $v:C_{*}\rightarrow\mathbb{R}^{2}$ and fixed a value $p\in\mathbb{R}^{2}$ we define $C_{*}^{p}=(C_{k}^{p},\partial_{k}^{p})$ as the chain complex for which:

• •

$C_{k}^{p}$ is the space of the $k$-chains of $C_{k}$ having value lower than or equal to $p$,

• •

$\partial_{k}^{p}$ is the restriction of $\partial_{k}$ to $C_{k}^{p}$.

By the definition of $v$, for any chain $c$ of dimension $k$, we have that $v(\partial_{k}(c))\leq v(c)$. So, $C_{*}^{p}$ is well-defined. For $p,q\in\mathbb{R}^{2}$ with $p\leq q$, $C_{*}^{p}$ is a chain subcomplex of $C_{*}^{q}$. For this reason, we denote the collection $C$ of the chain complexes $C_{*}^{p}$ with $p\in\mathbb{R}^{2}$ as bifiltered chain complex.

Given $p,q\in\mathbb{R}^{2}$ with $p\leq q$, the inclusion map from $C_{*}^{p}$ to $C_{*}^{q}$ induces a linear map between the corresponding homology spaces $H_{k}(C_{*}^{p})$ and $H_{k}(C_{*}^{p})$. The multi-parameter persistence $k^{th}$ module $H_{k}(C)$ of a bifiltered chain complex $C$ is the collection of the homology spaces $H_{k}(C_{*}^{p})$ with $p$ varying in $\mathbb{R}^{2}$ along with all the linear maps induced by the inclusion maps.

Let $C,D$ be two bifiltered chain complexes. We call $C$ and $D$ homology-equivalent if, for any fixed $k\in\mathbb{N}$, $H_{k}(C_{*}^{p})$ and $H_{k}(D_{*}^{p})$ are isomorphic via a map $\psi^{p}_{k}$ and, for any $p,q\in\mathbb{R}^{2}$ with $p\leq q$, the diagram {diagram} commutes where horizontal maps are induced by inclusion maps. Moreover, we call $C$ and $D$ quasi-isomorphic if the isomorphisms of the above diagram are induced by a collection of chain maps $f^{p}_{*}:C^{p}_{*}\rightarrow D^{p}_{*}$ satisfying, for any $p\leq q$ and any $k$, the commutative diagram {diagram} in which horizontal maps are the inclusion maps. By definition, quasi-isomorphic bifiltered chain complexes are homology-equivalent. The converse, as depicted in Figure 1, does not hold in general.

Two chain maps $\alpha_{*},\beta_{*}:C_{*}\rightarrow D_{*}$ are called chain-homotopic if there exists a collection of maps $\phi_{k}:C_{k}\rightarrow D_{k+1}$ such that, for any $k$, $\partial^{D}_{k+1}\phi_{k}+\phi_{k-1}\partial^{C}_{k}=\alpha_{k}-\beta_{k}$. Two bifiltered chain complexes $C$ and $D$ are called homotopy-equivalent if, there exist two collections of chain maps $f^{p}_{*}:C^{p}_{*}\rightarrow D^{p}_{*}$ and $g^{p}_{*}:D^{p}_{*}\rightarrow C^{p}_{*}$ such that, for any $p\in\mathbb{R}^{2}$, $f^{p}_{*}$ and $g^{p}_{*}$ are homotopy-inverse one with respect to the other (i.e., $g^{p}_{*}f^{p}_{*}$ and $f^{p}_{*}g^{p}_{*}$ are chain-homotopic to $\mathrm{id}_{C^{p}_{*}}$ and $\mathrm{id}_{D^{p}_{*}}$, respectively) and they satisfy, for any $p\leq q$, the following commutative diagram {diagram}

in which horizontal maps are the inclusion maps. Homotopy-equivalent bifiltered chain complexes are necessarily quasi-isomorphic.

Representation of bifiltered chain complexes.

We assume that the bifiltered chain complex in input is provided as a finite sequence of triples of the form

[TABLE]

where $p\in\mathbb{R}^{2}$, $k$ is a positive integer, and $\Delta$ is a finite list of pairs $(l,\lambda)$ with $\lambda\in\mathbb{F}$ and $l\in\mathbb{N}$. In order to retrieve the bifiltered chain complex $C$ represented by this list, each triple $(p,k,\Delta)$ has to be interpreted as a $k$-generator $c$ with value $v(c)=p$ and whose boundary is encoded through $\Delta$. Specifically, by defining as $c_{l}$ the generator stored in position $l$ of the sequence of triples, the pair $(l,\lambda)$ in $\Delta$ represents that $c_{l}$ appears in the boundary of $c$ with coefficient $\lambda$. In order to ensure that the above-described interpretation of the input sequence actually returns a valid bifiltered chain complex, we require that, if $c^{\prime}$ is a generator appearing in the boundary of $c$, then the value $v$ associated to $c^{\prime}$ is lower than or equal to the one taken by $c$.

Given a chain complex $C_{*}$ endowed with a value function $v$, let us consider an injective function $i$ providing the set of generators of $C_{*}$ with a total order which is consistent with the partial order induced by the value funtion $v$. More precisely, given two generators $c_{1},c_{2}$ of $C_{*}$, the index function $i$ has to satisfy that $v(c_{1})\leq v(c_{2})$ implies $i(c_{1})\leq i(c_{2})$.

We uniquely represent the bifiltered chain complex $C$ as follows. Each generator $c$ of $C_{*}$ is stored in a data structure that we call the column of the generator. A column stores:

• •

the index $i(c)$,

• •

the value $v(c)$,

• •

the dimension $k$ of $g$,

• •

the boundary of $c$.

The boundary is stored as a (possibly empty) container of entries of the form $(i(c^{\prime}),\lambda)$, where $c^{\prime}$ is a generator of dimension $k-1$ such that $\lambda:=\langle\partial_{k}(c),c^{\prime}\rangle\neq 0$. In such a case, we say that $c^{\prime}$ is in the boundary of $c$. The boundary of $c$ is then the linear combination induced by this container. The name “column” comes from the idea that the collection of all columns can be visualized as a matrix. If we decorate the columns of the matrix with extra information (index, value, and dimension), we obtain the data structure from above. We will also write $k$-column if the column represents a generator of dimension $k$. Due to the fact that a column is nothing but an encoding of a generator, in the following, with a small abuse of notation we will often use interchangeably the two terms. Given a $k$-column $c$, we define the local pivot of $c$ as the generator in the boundary of $c$ with maximal index such that $v(c^{\prime})=v(c)$. If no such generator exists, we simply say that $c$ has no local pivot.

Arguably, the most common case of input data is simplicial complexes. Let $K$ denote a finite simplicial complex. A bifiltration is a collection $(K^{p})_{p\in\mathbb{R}^{2}}$ of subcomplexes of $K$ such that $K^{p}\subseteq K^{q}$ whenever $p\leq q$. Fixing a simplex $\sigma\in K$, we say that $p\in\mathbb{R}^{2}$ is critical for $\sigma$, if $\sigma\in K^{p}$, but $\sigma\notin K^{p-(\epsilon,0)}\cup K^{p-(0,\epsilon)}$ for every $\epsilon>0$. The bifiltration is called $h$-critical with $h\geq 1$, if for every $\sigma\in K$, there are at most $h$ critical positions in $\mathbb{R}^{2}$. For instance, $1$-critical means that every simplex in $K$ enters the filtration at a unique minimal value. Such 1-critical bifiltrations can easily be described by a sequence of the form $(p,k,\Delta)$ as above by adding one line per simplex which defines its critical position, its dimensions, and its boundary. The case of $h$-critical bifiltrations can be handled with the following trick (see also [7]). Letting $p_{1},\ldots,p_{h}$ denote the critical positions of $\sigma$, sorted by $x$-coordinate, we introduce $h$ distinct $k$-generators $c_{1},\ldots,c_{h}$ of the form $(p_{i},k,\Delta)$, that is, $h$ copies of the same simplex. For two consecutive positions $p_{i}=(x_{i},y_{i})$ and $p_{i+1}=(x_{i+1},y_{i+1})$, we add an additional $(k+1)$-generator at $(x_{i+1},y_{i})$ (which is the smallest point $q$ such that $p_{i}\leq q$ and $p_{i+1}\leq q$), whose boundary is equal to $c_{i}-c_{i+1}$.

3 Algorithm

Given a bifiltered chain complex $C$ encoded as a collection of columns, we define the chunk algorithm that returns as output a bifiltered chain complex as a collection of columns that is homotopy-equivalent to $C$ and has fewer generators. The algorithm works in three phases:

• •

local reduction;

• •

compression;

• •

removal of local columns.

Phase I: Local reduction.

The goal of this phase is to label columns as global, local positive or local negative columns. Initially, all columns are unlabeled. We proceed in decreasing dimensions, from the top-dimensions of $C$ down to [math]. For dimension $k$, the algorithm traverses the $k$-columns in increasing order with respect to $i$ and performs the following operations on a $k$-column $c$. If $c$ is already labeled, do nothing. Otherwise, as long as $c$ has a local pivot and there is a $k$-column $c^{\prime}$ with $i(c^{\prime})<i(c)$ and the same local pivot as $c$, perform the column addition $c\leftarrow c+\lambda c^{\prime}$, where $\lambda$ is chosen such that the local pivot of $c$ disappers. If at the end of this loop, the column $c$ does not have a local pivot, we label the column as global and proceed. Otherwise, we label $c$ as local negative and its local pivot $c^{\prime}$ as local positive. We call $(c^{\prime},c)$ a local pair in this case. This ends the description of Phase I of the algorithm.

Note that within the local reduction, any column addition of the form $c\leftarrow c+\lambda c^{\prime}$ implies that $v(c)=v(c^{\prime})$. Hence, the local reduction operates independently on columns of the same value. We call blocks of columns with the same value chunks; hence the name of the algorithm. Operations on one chunk do not affect columns on any other chunk, hence the local reduction phase can be readily invoked in parallel on the chunks of the chain complex.

Finally, note that by proceeding in decreasing dimension, we avoid performing any column additions on local positive columns. That is reminiscent of the clearing optimization in the one-parameter version [3, 8].

Phase II: Compression.

In the second phase, the algorithm removes local (positive or negative) generators from the boundary of all global columns in the matrix:

For each global $k$-column $c$, while the boundary of the column contains a generator that is local positive or local negative, the algorithm picks the local $(k-1)$-generator $c^{\prime}$ with maximal index.

• •

if $c^{\prime}$ is negative, remove $c^{\prime}$ from the boundary of $c$;

• •

if $c^{\prime}$ is positive, denote $c^{\prime\prime}$ as the (unique) local negative $k$-column with $c^{\prime}$ as local pivot and perform the column addition $c\leftarrow c+\lambda c^{\prime\prime}$, where $\lambda$ is chosen such that $c^{\prime}$ disappears in the boundary of $c$.

This ends the description of the compression phase. On termination, all columns in the boundary of a global $k$-column are global $(k-1)$-columns.

The above process terminates for a column $c$ because the index of the maximal local generator in the boundary of $c$ is strictly decreasing in each step. That is clear for the case that $c^{\prime}$ is local negative. If $c^{\prime}$ is local positive, then $c^{\prime}$ is the generator in the boundary of $c^{\prime\prime}$ with the maximal index, so the column addition does not introduce in the boundary of $c$ any generators with a larger index.

Note that the compression of a global column does not affect the result on any other global column. Thus, the phase can be parallelized as well.

Phase III: Removal of local pairs.

In this step, the chain complex becomes smaller. The procedure is simple: traverse all columns, and remove all columns labeled as local (positive or negative). Return the remaining (global) columns as resulting chain complex. This finishes the description of the phase and the entire chunk algorithm.

4 Correctness

In this section, we prove that the presented algorithm returns a bifiltered chain complex that is homotopy-equivalent to the input. For that, we define two elementary operations on chain complexes:

Order-preserving column addition.

An operation of the form $c\leftarrow c+\lambda c^{\prime}$ is called order-preserving if $v(c^{\prime})\leq v(c)$. Note that such an operation maintains the property that any generator $c^{\prime\prime}$ in the boundary of $c$ satisfies $v(c^{\prime\prime})\leq v(c)$, by transitivity of $v$. We remark that order-preserving column additions are the generalization of left-to-right column additions in the one-parameter case.

Removal of a local pair.

Fix a local pair $(c_{1},c_{2})$, that means, $c_{1}$ is a $k$-column, $c_{2}$ is a $(k+1)$-column, $v(c_{1})=v(c_{2})$ and $c_{1}$ is the local pivot of $c_{2}$. We call removal of the local pair $(c_{1},c_{2})$ the operation $Del(c_{1},c_{2})$ which acts on the columns as follows.

• •

For every $(k+1)$-column $c$, replace its boundary $\partial_{k+1}(c)$ with $\partial_{k+1}(c)-\lambda^{-1}\mu\partial_{k+1}(c_{2})$, where $\lambda$ and $\mu$ are the coefficients of $c_{1}$ in $\partial_{k+1}(c_{2})$ and in $\partial_{k+1}(c)$, respectively. In particular, after this operation, $c_{1}$ disappeared from the boundary of any $(k+1)$-columns.

• •

For every $(k+2)$-column $c$, update its boundary by setting the coefficient of $c_{2}$ in $\partial_{k+2}(c)$ to [math]. Visualizing the chain complex as a matrix, this corresponds to removing the row corresponding to $c_{2}$.

• •

Delete the columns $c_{1}$ and $c_{2}$.

Note that all column additions performed in the first step are order-preserving because the pair $(c_{1},c_{2})$ is local, that is, $v(c_{1})=v(c_{2})$.

We will show at the end of this section that both elementary operations leave the homotopy type the bifiltered chain complex unchanged.

Theorem 4.1**.**

Let $\bar{C}$ denote the bifiltered chain complex computed by the chunk algorithm from the previous section on an input bifiltered chain complex $C$. Then, $\bar{C}$ and $C$ are homotopy-equivalent.

Proof 4.2**.**

The idea is to express the chunk algorithm by a sequence of order-preserving column additions and removals of local pairs. Because every column addition in Phase I is between columns of the same value, all column additions are order-preserving. Hence, after Phase I, the chain complex is equivalent to the input.

In Phase II, note that all column additions performed are order-preserving. Indeed, if $c^{\prime}$ is in the boundary of column $c$, then $v(c^{\prime})\leq v(c)$ holds. If $c^{\prime}$ is local positive, it triggers a column addition of the form $c\leftarrow c+\lambda c^{\prime\prime}$ with its local negative counterpart $c^{\prime\prime}$. Since $v(c^{\prime})=v(c^{\prime\prime})$, $v(c^{\prime\prime})\leq v(c)$ as well.

A further manipulation in Phase II is the removal of local negative columns from the boundary of global columns. These removals cannot be directly expressed in terms of the two elementary operations from above. Instead, we define a slight variation of our algorithm: in Phase II, when we encounter a local negative $c^{\prime}$, we do nothing. In other words, the compression only removes the local positive generators from the boundary $c$, and keeps local negative and global generators. In Phase III, instead of removing local columns, we perform a removal of a local pair $(c_{1},c_{2})$ whenever we encounter a local negative column $c_{2}$ with local pivot $c_{1}$. We call that algorithm the modified chunk algorithm. Note that this modified algorithm is a sequence of order-preserving column additions, followed by a sequence of local pair removals, and thus produces a chain complex that is equivalent to the input $C$.

We argue next that the chunk algorithm and the modified chunk algorithm yield the same output. Since both versions eventually remove all local columns, it suffices to show that they yield the same global columns. Fix an index of a global column, and let $c$ denote the column of that index returned by the original chunk algorithm. Let $c^{\ast}$ denote the column of the same index produced by the modified algorithm after the modified Phase II. The difference of $c^{\ast}$ and $c$ lies in the presence of local negative generators in the boundary of $c^{\ast}$ which have been removed in $c$. The modified Phase III affects $c^{\ast}$ in the following way: when a local pair $(c_{1},c_{2})$ is removed, the local negative $c_{2}$ is, if it is present, removed from the boundary of $c^{\ast}$. There is no column addition during the modified Phase III involving $c^{\ast}$ because all local positive columns have been eliminated. Hence, the effect of the modified Phase III on $c^{\ast}$ is that all local negative columns are removed from its boundary which turns $c^{\ast}$ to be equal to $c$ at the end of the algorithm. Hence, the output of both algorithms is the same, proving the theorem.

We proceed with the proofs that both elementary operations yield homotopy-equivalent bifiltered chain complexes.

Order-preserving column addition.

Given two $k$-columns $c_{1},c_{2}$ such that $v(c_{1})\leq v(c_{2})$ and a scalar value $\lambda\in\mathbb{F}$, the algorithm we propose allows for adding $\lambda$ copies of the boundary of column $c_{1}$ to the boundary of column $c_{2}$. In this section, we formalize that in terms of modifications of chain complexes and we prove that this operation does not affect multi-parameter persistent homology.

Given a chain complex $C_{*}=(C_{l},\partial_{l})$ endowed with a value function $v:C_{*}\rightarrow\mathbb{R}^{2}$, let us consider two generators $c_{1},c_{2}$ among the ones of the space of the $k$-chains $C_{k}$ such that $\langle c_{1},c_{2}\rangle=0$ and $v(c_{1})\leq v(c_{2})$. Chosen a scalar value $\lambda\in\mathbb{F}$, let us define $\bar{C}_{*}=(\bar{C}_{l},\bar{\partial}_{l})$ by setting:

• •

$\bar{C}_{l}=C_{l}$,

• •

for any $c\in\bar{C}_{l}$,

[TABLE]

As formally proven in Appendix A, $\bar{C}_{*}$ is a chain complex.

Let us define the maps $f_{*}:C_{*}\rightarrow\bar{C}_{*}$, $g_{*}:\bar{C}_{*}\rightarrow C_{*}$ as follows:

• •

for any $c\in C_{l}$,

[TABLE]

• •

for any $c\in\bar{C}_{l}$,

[TABLE]

The just defined maps enable to prove the following result (see Appendix A for detailed proofs).

Proposition 4.3**.**

$C_{*}$* and $\bar{C}_{*}$ are isomorphic via the chain map $f_{*}$ and its inverse $g_{*}$.*

Since the function $v$ is a valid value function for $\bar{C}_{*}$ inducing a bifiltered chain complex $\bar{C}$ and, for any $p\in\mathbb{R}^{2}$, the restrictions $f^{p}_{l}:C^{p}_{l}\rightarrow\bar{C}^{p}_{l}$ and $g^{p}_{l}:\bar{C}^{p}_{l}\rightarrow C^{p}_{l}$ of the maps $f_{l}$ and $g_{l}$, respectively, are well-defined (see Appendix A), we are ready to prove the equivalence between the two bifiltered chain complexes $C$ and $\bar{C}$.

Theorem 4.4**.**

$C$* and $\bar{C}$ are homotopy-equivalent.*

Proof 4.5**.**

The previous results enables us to claim that, for any $p,q\in\mathbb{R}^{2}$ with $p\leq q$, the following diagram (in which horizontal maps are the inclusion maps) {diagram} commutes and the maps $f^{p}_{l}$, $f^{q}_{l}$ are isomorphisms.

Removal of a local pair.

Given a local pair $(c_{1},c_{2})$ of columns, the algorithm we propose allows for deleting them from the boundary matrix. In this section, we formalize that in terms of modifications of chain complexes and we prove that this operation does not affect multi-parameter persistent homology.

Given a chain complex $C_{*}=(C_{l},\partial_{l})$ endowed with a value function $v:C_{*}\rightarrow\mathbb{R}^{2}$, let us consider two generators $c_{1},c_{2}$ among the ones of the space of the $k$-chains $C_{k}$ and of the space of the $(k+1)$-chains $C_{k+1}$, respectively, such that $\lambda:=\langle\partial_{k+1}(c_{2}),c_{1}\rangle\neq 0$ and $v(c_{1})=v(c_{2})$. Let us define $\bar{C}_{*}=(\bar{C}_{l},\bar{\partial}_{l})$ by setting:

• •

the space of the $l$-chains $\bar{C}_{l}$ as

[TABLE]

• •

for any $c\in\bar{C}_{l}$,

[TABLE]

As formally proven in Appendix A, $\bar{C}_{*}$ is a chain complex.

Let us define the maps $r_{*}:C_{*}\rightarrow\bar{C}_{*}$, $s_{*}:\bar{C}_{*}\rightarrow C_{*}$ as follows:

• •

for any $c\in C_{l}$,

[TABLE]

• •

for any $c\in\bar{C}_{l}$,

[TABLE]

The just defined maps enable to prove the following result (see Appendix A for detailed proofs).

Proposition 4.6**.**

The maps $r_{*}$ and $s_{*}$ are chain maps which are homotopy-inverse one with respect to the other.

This latter combined with the fact that the function $v$ is a valid value function for $\bar{C}_{*}$ inducing a bifiltered chain complex $\bar{C}$ and, for any $p\in\mathbb{R}^{2}$, the restrictions $r^{p}_{l}:C^{p}_{l}\rightarrow\bar{C}^{p}_{l}$, $s^{p}_{l}:\bar{C}^{p}_{l}\rightarrow C^{p}_{l}$ of the maps $r_{l}$ and $s_{l}$, as well as the restrictions of the maps ensuring that $r_{*}$ and $s_{*}$ are homotopy-inverse, are well-defined (see Appendix A), guarantees the homotopy-equivalence between the two bifiltered chain complexes $C$ and $\bar{C}$.

Theorem 4.7**.**

$C$* and $\bar{C}$ are homotopy-equivalent.*

5 Optimality and complexity

Optimality.

Let $D$ be a bifiltered chain complex. For $p\in\mathbb{R}^{2}$, we define

[TABLE]

where the sum of chain complexes, analogously to the sum of the vector spaces, is the chain complex spanned by the union of the generators of the summands. Moreover, let $\eta^{p}_{k}$ be the homology map in dimension $k$ induced by the inclusion of $D^{<p}_{*}$ into $D^{p}_{*}$. We denote the number of variations in the $k^{th}$ homology space occurred at value $p$ as

[TABLE]

and the number of $k$-generators added at value $p$ in $D$ as

[TABLE]

Theorem 5.1**.**

Let $\bar{C}$ be the bifiltered chain complex obtained by applying the chunk algorithm to the bifiltered chain complex $C$. We have that $\delta^{p}_{k}(C)=\gamma^{p}_{k}(\bar{C})$.

The full proof is given in Appendix B and can be summarized as follows. Every global column with value $p$ either destroys a homology class of $H(C^{<p})$, or it creates a new homology class in $H(C^{p})$, which is not destroyed by any other column of value $p$. Hence, each global column contributes a generator to $\mathrm{ker}\,\iota^{p}_{k-1}$ or to $\mathrm{coker}\,\iota^{p}_{k}$, where $\iota^{p}_{l}$ is the map between the $l^{th}$ homology spaces induced by the inclusion of $C^{<p}_{*}$ into $C^{p}_{*}$. Local columns do not contribute to either of these two spaces. The result follows from the fact that the number of global columns at value $p$ is precisely the number of generators added at $\bar{C}$.

The next statement shows that our construction is optimal in the sense that any bifiltered chain complex $D$ that is quasi-isomorphic to $C$ must have at least as many generators as $\bar{C}$.

Theorem 5.2**.**

Any bifiltered chain complex $D$ quasi-isomorphic to $C$ has to add at least $\delta^{p}_{k}(C)$ $k$-generators at value $p$. I.e., $\delta^{p}_{k}(C)\leq\gamma^{p}_{k}(D)$.

The full proof is also given in Appendix B. To summarize it, it is not too hard to see that, for any bifiltered chain complex $D$,

[TABLE]

holds. Moreover, the quasi-isomorphism of $C$ and $D$ implies that $\mathrm{dim}\,\mathrm{ker}\,\eta^{p}_{k-1}=\mathrm{dim}\,\mathrm{ker}\,\iota^{p}_{k-1}$ and $\mathrm{dim}\,\mathrm{coker}\,\eta^{p}_{k}=\mathrm{dim}\,\mathrm{coker}\,\iota^{p}_{k}$. So,

[TABLE]

which implies that claim. The equality of the dimensions is formally verified using the Mayer-Vietoris sequence and the $5$-lemma to establish an isomorphism from $H_{k}(C^{<p}_{*})$ to $H_{k}(D^{<p}_{*})$ that commutes with the isomorphism at value $p$.

Complexity.

In order to properly express the time and the space complexity of the proposed algorithm, let us introduce the following parameters. Given a bifiltered chain complex $C$, we denote as $n$ the number of generators of $C$, as $m$ the number of chunks (i.e., the number of different values assumed by $v$), as $\ell$ the maximal size of a chunk, and as $g$ the number of global columns. Moreover, we assume the maximal size of the support of the boundary of the generators of $C$ as a constant. The latter condition is always ensured for bifiltered simplicial complex of fixed dimension.

Theorem 5.3**.**

The chunk algorithm has time complexity $O(m\ell^{3}\log\ell+g\ell n\log\ell)$ and space complexity $O(n\ell+g^{2})$.

Proof 5.4**.**

Due to the similarity between the two algorithms, the analysis of complexity of the proposed algorithm is analogous to the first two steps of the one-parameter chunk algorithm [3]. The additional factor of $\log\ell$ comes from our choice of using priority queues as column type and could be removed by using list representations as in [3].222We remark, however, that this would result in a performance penalty in practice. See [4].

On the space complexity, during the Phase I, the generators in the boundary of any column can be at most $O(\ell)$. So, $O(n\ell)$ is a bound on the accumulated size of all columns after Phase I. During Phase II, the boundary of any global column can consist of up to $n$ generators, but reduces to $g$ generators at the end of the compression of the column because all local entries have been removed. Hence, the final chain complex has at most $g$ entries in each of its $g$ columns, and requires $O(g^{2})$ space. Hence, the total space complexity is $O(n\ell+n+g^{2})$, where the second summand is redundant.333We remark that this bound only holds for the sequential version of the algorithm. In a parallelized version, it can happen that several compressed columns achieve a size of $O(n)$ at the same time.

6 Implementation and experimental results

We briefly introduce the developed implementation of the chunk algorithm and we evaluate its performance. The current C++ implementation of the chunk algorithm consists of approximately 350 lines of code. It takes as input a finite bifiltered chain complex expressed accordingly to the representation described in Section 2 and encodes each column as a std::priority_queue. Even if the chunk algorithm permits a parallelization in shared memory, this first implementation does not exploit this potential. Our experiments have been performed on a configuration Intel Xeon CPU E5-1650 v3 at 3.50 GHz with 64 GB of RAM.

For testing the implemented chunk algorithm, we have compared its performances with the simplification process based on discrete Morse theory proposed in [18] (DMT-based algorithm) whose implementation is publicly available [11]. We remark once more that the DMT-based approach yields a somewhat richer output than solely a simplified chain complex with the same homotopy; however, we only compare the size of the resulting structure (in terms of the number of generators per chain group) in this experimental comparison.

In our experiments, we have considered both synthetic and real datasets represented as simplicial complexes. The latter ones are courtesy of the AIM@SHAPE data repository [9]. Most of the datasets are of dimension 2 or 3 and are embedded in a 3D environment. For our experiments, we have adopted as filtering functions the ones obtained by extending to all the simplices of the complex the $x$ and the $y$ coordinates assigned to the vertices of the datasets. Table 1 displays the achieved results.

The data sets are given in off file format, that means, as a list of triangles specified by boundary vertices. In order to apply our algorithm, we first have to convert the data into a boundary matrix representation. Hence, we enlist in Table 1 the preparation time to create the boundary matrix and the simplification time to perform our chunk algorithm. In contrast, the DMT-based approach avoids the initial construction of the boundary matrix but transforms the input into a different data structure before starting its simplification step. We also list the running times of these separate steps in Table 1. In both cases, we do not list the time for reading the input file into memory and writing the output structure on disk.

The column Size in Table 1 collects the sizes (in terms of the number of simplices/columns) of the bifiltered chain complexes in input and in output. The compression factor achieved by both algorithms varies between 9 and 23. As average, the reduced chain complex in output is approximately 13 times smaller than the original one. Despite of the theoretical advantage that our approach yields an optimal size in all situations, the size of the output returned by the two algorithms is nearly the same in all listed examples. A difference can be noticed just for datasets including shapes that can be considered as “pathological”. For instance, an example of such a dataset is the conification of a dunce hat.

The column Time shows the computation times of both algorithms. In all tested examples, the chunk algorithm takes a small fraction of time with respect to the DMT-based approach. We also see that the creation of the boundary matrix (preparation step) takes more time than the chunk algorithm (simplification step) by a factor of $2$ to $3$, but even these two steps combined are faster than the simplification step of the DMT-based approach. It is worth investigating whether the boundary matrix creation becomes a more severe bottleneck for other datasets (e.g., in higher dimension), where the DMT-based approach might have advantages by not creating the boundary matrix explicitly.

The column Memory Usage shows the memory consumption of both approaches. The chunk algorithm does not need to encode auxiliary structures like the discrete Morse gradient stored in the DMT-based approach resulting in a slightly less amount of required space. Depending on the dataset, the chunk approach requires between 0.09 and 5.50 GB of memory and it is, on average, 1.2 times more compact than the DMT-based one.

In summary, the chunk algorithm satisfies a stronger optimality condition than [18], but still is an order of magnitude faster and uses comparable memory.

7 Conclusion

We have presented a pre-processing procedure for improving the computation of multi-parameter persistent homology and we have provided theoretical and experimental evidence of its effectiveness. In the future, we want to further improve the proposed strategy as follows. First, we would like to develop a parallel implementation with shared memory of the chunk algorithm and compare its performances with the ones obtained by the current version. Similarly to [18], we also want to evaluate the impact of the chunk algorithm for the computation of the persistence module and of the persistence space.

Our optimality criterion is formulated for the class of chain complexes quasi-isomorphic to the input. It is possible to produce counterexamples showing that this statement cannot be generalized to the class of the chain complexes homology-equivalent to the input. Nevertheless, we want to further investigate the optimality properties satisfied by our algorithm and compare our algorithm with the minimal presentation algorithm for persistence modules from RIVET.

Appendix A Proofs from Section 4

Order-preserving column additions.

In accordance with the notations adopted in the paragraph of Section 4 devoted to the analysis of order-preserving column additions, we show here the detailed proofs of the mentioned results.

Proposition A.1**.**

$\bar{C}_{*}$* is a chain complex.*

Proof A.2**.**

*It is enough to prove that $\bar{\partial}_{l-1}\bar{\partial}_{l}=0$ for $l=k,k+1,k+2$.

If $l=k$, then, thanks to the fact that $C_{*}$ is a chain complex,*

[TABLE]

If $l=k+1$, then, since $C_{*}$ is a chain complex and $\langle c_{1},c_{2}\rangle=0$,

[TABLE]

If $l=k+2$, then, since $C_{*}$ is a chain complex,

[TABLE]

Lemma A.3**.**

$f_{*}$* and $g_{*}$ are chain maps.*

Proof A.4**.**

*For proving that $f_{*}$ is a chain map is enough to show that $\bar{\partial}_{l}f_{l}=f_{l-1}\partial_{l}$ for $l=k,k+1$.

If $l=k$, then*

[TABLE]

If $l=k+1$, then

[TABLE]

*For proving that $g_{*}$ is a chain map is enough to show that $g_{l-1}\bar{\partial}_{l}=\partial_{l}g_{l}$ for $l=k,k+1$.

If $l=k$, then*

[TABLE]

If $l=k+1$, then

[TABLE]

Proposition A.5**.**

$C_{*}$* and $\bar{C}_{*}$ are isomorphic.*

Proof A.6**.**

It is enough to show that $g_{k}f_{k}=\mathrm{id}_{C_{k}}$ and $f_{k}g_{k}=\mathrm{id}_{\bar{C}_{k}}$. We have that

[TABLE]

Analogously,

[TABLE]

Lemma A.7**.**

The function $v$ is a valid value function for $\bar{C}_{*}$ inducing a bifiltered chain complex $\bar{C}$.

Proof A.8**.**

We have to show, that, for any $c\in\bar{C}_{l}$, $v(\bar{\partial}_{l}(c))\leq v(c)$. For $l\neq k,k+1$, that is trivial. Let us prove the other cases, recalling that:

• •

$v(c_{1})\leq v(c_{2})$,

• •

for any $c\in C_{l}$, $v(\partial_{l}(c))\leq v(c)$.

If $l=k$, $\bar{\partial}_{k}(c)=\partial_{k}(c)+\lambda\langle c,c_{2}\rangle\partial_{k}(c_{1})$. If $\langle c,c_{2}\rangle=0$, then $v(\bar{\partial}_{k}(c))=v(\partial_{k}(c))\leq v(c)$. Otherwise,

• •

$v(\partial_{k}(c_{1}))\leq v(c_{1})\leq v(c_{2})\leq v(c)$,

• •

$v(\partial_{k}(c))\leq v(c)$.

*Then, $v(\bar{\partial}_{k}(c))\leq v(c)$.

If $l=k+1$, $\bar{\partial}_{k+1}(c)=\partial_{k+1}(c)-\lambda\langle\partial_{k+1}(c),c_{2}\rangle c_{1}$. If $\langle\partial_{k+1}(c),c_{2}\rangle=0$, then $v(\bar{\partial}_{k+1}(c))=v(\partial_{k+1}(c))\leq v(c)$. Otherwise, $v(c_{1})\leq v(c_{2})\leq v(\partial_{k+1}(c))\leq v(c)$. Then, $v(\bar{\partial}_{k+1}(c))\leq v(c)$.*

Lemma A.9**.**

For any $p\in\mathbb{R}^{2}$, the restrictions $f^{p}_{l}:C^{p}_{l}\rightarrow\bar{C}^{p}_{l}$ and $g^{p}_{l}:\bar{C}^{p}_{l}\rightarrow C^{p}_{l}$ of the maps $f_{l}$ and $g_{l}$, respectively, are well-defined.

Proof A.10**.**

*For $f^{p}_{l}$, we have to prove that, for any $c\in C_{l}$, $v(f_{l}(c))\leq v(c)$. Let us focus on the only non-trivial case occurring for $l=k$.

If $l=k$, $f_{k}(c)=c-\lambda\langle c,c_{2}\rangle c_{1}$. If $\langle c,c_{2}\rangle=0$, then $v(f_{k}(c))=v(c)$. Otherwise, $v(c_{1})\leq v(c_{2})\leq v(c)$. Then, $v(f_{k}(c))\leq v(c)$.

Except for a sign, the proof for $g^{p}_{l}$ is identical.*

Removal of local pairs.

In accordance with the notations adopted in the paragraph of Section 4 devoted to the analysis of removal of local pairs, we show here the detailed proofs of the mentioned results.

Proposition A.11**.**

$\bar{C}_{*}$* is a chain complex.*

Proof A.12**.**

*First of all, we have to prove that $\bar{\partial}_{l}$ is well-defined. In order to do that, we have to show that, for $l=k+1,k+2$, if $c\in\bar{C}_{l}$ then $\bar{\partial}_{l}(c)\in\bar{C}_{l-1}$.

If $l=k+1$,*

[TABLE]

*Then, $\bar{\partial}_{k+1}(c)\in\bar{C}_{k}$.

If $l=k+2$,*

[TABLE]

*Then, $\bar{\partial}_{k+2}(c)\in\bar{C}_{k+1}$.

Second, we have to prove that $\bar{\partial}_{l}\bar{\partial}_{l+1}=0$. It is enough to show that for $l=k,k+1,k+2$.

If $l=k$, then, thanks to the fact that $C_{*}$ is a chain complex,*

[TABLE]

If $l=k+1$, then, since $C_{*}$ is a chain complex,

[TABLE]

If $l=k+2$, then, since $C_{*}$ is a chain complex,

[TABLE]

Lemma A.13**.**

$r_{*}$* and $s_{*}$ are chain maps.*

Proof A.14**.**

*For proving that $r_{*}$ is a chain map is enough to show that $\bar{\partial}_{l}r_{l}=r_{l-1}\partial_{l}$ for $l=k,k+1,k+2$.

If $l=k$, then*

[TABLE]

If $l=k+1$, then

[TABLE]

If $l=k+2$, then

[TABLE]

*For proving that $s_{*}$ is a chain map is enough to show that $s_{l-1}\bar{\partial}_{l}=\partial_{l}s_{l}$ for $l=,k+1,k+2$.

If $l=k+1$, then*

[TABLE]

If $l=k+2$, then

[TABLE]

In order to prove the following result, let us define a collection of maps $\phi_{l}:C_{l}\rightarrow C_{l+1}$, for any $l\in\mathbb{N}$. Given $c\in C_{l}$, the map $\phi_{l}$ is defined as:

[TABLE]

Proposition A.15**.**

The maps $r_{*}$ and $s_{*}$ are chain maps which are homotopy-inverse one with respect to the other.

Proof A.16**.**

In order to prove the thesis, we show that, for any $l$:

• (1)

$s_{l}r_{l}$* is chain-homotopic to $\mathrm{id}_{C_{l}}$ via $\phi_{l}$,*

• (2)

$r_{l}s_{l}=\mathrm{id}_{\bar{C}_{l}}$.

Let us start with proving $(1)$. We have to show that

[TABLE]

*The only non-trivial cases occur for $l=k,k+1$.

If $l=k$, we have that*

[TABLE]

If $l=k+1$, we have that

[TABLE]

*Now, let us prove $(2)$. The only non-trivial cases occur for $l=k,k+1$.

If $l=k$, since $c\in\bar{C}_{k}$ implies $\langle c,c_{1}\rangle=0$, we have that*

[TABLE]

If $l=k+1$, since $c\in\bar{C}_{k+1}$ implies $\langle c,c_{2}\rangle=0$, we have that

[TABLE]

Lemma A.17**.**

The function $v$ is a valid value function for $\bar{C}_{*}$ inducing a bifiltered chain complex $\bar{C}$.

Proof A.18**.**

We have to show, that, for any $c\in\bar{C}_{l}$, $v(\bar{\partial}_{l}(c))\leq v(c)$. For $l\neq k+1,k+2$, that is trivial. Let us prove the other cases, recalling that:

• •

$v(c_{1})=v(c_{2})$,

• •

for any $c\in C_{l}$, $v(\partial_{l}(c))\leq v(c)$.

If $l=k+1$, $\bar{\partial}_{k+1}(c)=\partial_{k+1}(c)-\lambda^{-1}\langle\partial_{k+1}(c),c_{1}\rangle\partial_{k+1}(c_{2})$. If $\langle\partial_{k+1}(c),c_{1}\rangle=0$, then $v(\bar{\partial}_{k+1}(c))=v(\partial_{k+1}(c))\leq v(c)$. Otherwise,

$v(\partial_{k+1}(c_{2}))\leq v(c_{2})=v(c_{1})\leq v(\partial_{k+1}(c))\leq v(c)$*. Then, $v(\bar{\partial}_{k+1}(c))\leq v(c)$.

If $l=k+2$, $\bar{\partial}_{k+2}(c)=\partial_{k+2}(c)-\langle\partial_{k+2}(c),c_{2}\rangle c_{2}$. If $\langle\partial_{k+2}(c),c_{2}\rangle=0$, then $v(\bar{\partial}_{k+2}(c))=v(\partial_{k+2}(c))\leq v(c)$. Otherwise, $v(c_{2})\leq v(\partial_{k+2}(c))\leq v(c)$. Then, $v(\bar{\partial}_{k+2}(c))\leq v(c)$.*

Lemma A.19**.**

For any $p\in\mathbb{R}^{2}$, the restrictions $r^{p}_{l}:C^{p}_{l}\rightarrow\bar{C}^{p}_{l}$, $s^{p}_{l}:\bar{C}^{p}_{l}\rightarrow C^{p}_{l}$, and $\phi^{p}_{l}:C^{p}_{l}\rightarrow C^{p}_{l+1}$ of the maps $r_{l}$, $s_{l}$, and $\phi_{l}$, respectively, are well-defined.

Proof A.20**.**

*For $r^{p}_{l}$, we have to prove that, for any $c\in C_{l}$, $v(r_{l}(c))\leq v(c)$. Let us focus on the only non-trivial cases occurring for $l=k,k+1$.

If $l=k$, $r_{k}(c)=c-\lambda^{-1}\langle c,c_{1}\rangle\partial_{k+1}(c_{2})$. If $\langle c,c_{1}\rangle=0$, then $v(r_{k}(c))=v(c)$. Otherwise, $v(\partial_{k+1}(c_{2}))\leq v(c_{2})=v(c_{1})\leq v(c)$. Then, $v(r_{k}(c))\leq v(c)$.

If $l=k+1$, $r_{k+1}(c)=c-\langle c,c_{2}\rangle c_{2}$. If $\langle c,c_{2}\rangle=0$, then $v(r_{k+1}(c))=v(c)$. Otherwise, $v(c_{2})\leq v(c)$. Then, $v(r_{k+1}(c))\leq v(c)$.

For $s^{p}_{l}$, we have to prove that, for any $c\in\bar{C}_{l}$, $v(s_{l}(c))\leq v(c)$. Let us focus on the only non-trivial case occurring for $l=k+1$.

If $l=k+1$, $s_{k+1}(c)=c-\lambda^{-1}\langle\partial_{k+1}(c),c_{1}\rangle c_{2}$. If $\langle\partial_{k+1}(c),c_{1}\rangle=0$, then $v(s_{k+1}(c))=v(c)$. Otherwise, $v(c_{2})=v(c_{1})\leq v(\partial_{k+1}(c))\leq v(c)$. Then, $v(s_{k+1}(c))\leq v(c)$.

For $\phi^{p}_{l}$, we have to prove that, for any $c\in C_{l}$, $v(\phi_{l}(c))\leq v(c)$. Let us focus on the only non-trivial case occurring for $l=k$ and for $c\in C_{k}$ such that $\langle c,c_{1}\rangle\neq 0$.

If this is the case, then $v(c_{2})=v(c_{1})\leq v(c)$. So, $v(\phi_{k}(c))=v(\lambda^{-1}\langle c,c_{1}\rangle c_{2})=v(c_{2})\leq v(c)$.*

Appendix B Proofs from Section 5

In accordance with the notations adopted in Section 5, we show here the detailed proofs of the mentioned results.

Let $A_{*},B_{*}$ be two chain complexes such that $A_{*}$ is a subcomplex of $B_{*}$ and let $\eta_{k}:H_{k}(A_{*})\rightarrow H_{k}(B_{*})$ be the map in homology induced by this inclusion. Let us consider the boundary matrix of $B_{*}$ in which the generators of $A_{*}$ are placed on the left with respect to the generators of $B_{*}$ not belonging to $A_{*}$. By applying the standard reduction algorithm for computing persistent homology [10], we obtain a reduced matrix $R$. To avoid a clash in terminology with previous parts of the paper, we refer to the entry of lowest position in the matrix among the non-null ones in column $c^{\prime}$ as the lowest index of $c^{\prime}$. Let us denote the columns of $R$ corresponding to the generators of $B_{*}$ not belonging to $A_{*}$ as the columns of $B_{*}\setminus A_{*}$ and the others as the columns of $A_{*}$.

Lemma B.1**.**

We have that:

• •

the number of the $k$-columns of $B_{*}\setminus A_{*}$ having as lowest index a column of $A_{*}$ is equal to $\mathrm{dim}\,\mathrm{ker}\,\eta^{p}_{k-1}$,

• •

the number of the $k$-columns of $B_{*}\setminus A_{*}$ with null boundary for which there is no column in $R$ having them as lowest index is equal to $\mathrm{dim}\,\mathrm{coker}\,\eta^{p}_{k}$.

Proof B.2**.**

We can partition the $k$-columns of $B_{*}\setminus A_{*}$ in three distinct classes. The first class includes the $k$-columns of $B_{*}\setminus A_{*}$ having as lowest index a column of $A_{*}$. According to the standard reduction algorithm, these columns are in one-to-one correspondence with the non-null homology classes of $H_{k-1}(A_{*})$ that become trivial in $H_{k-1}(B_{*})$. So, their number coincides with the dimension of $\mathrm{ker}\,\eta^{p}_{k-1}$. The second class includes the $k$-columns of $B_{*}\setminus A_{*}$ with null boundary for which there is no column in $R$ having them as lowest index. According to the standard reduction algorithm, these columns are in one-to-one correspondence with the homology classes that are born in $H_{k}(B_{*})$. So, their number coincides with the dimension of $\mathrm{coker}\,\eta^{p}_{k}$. The third class of $k$-columns of $B_{*}\setminus A_{*}$ are the ones with lowest index in $B_{*}\setminus A_{*}$ or the columns with null boundary which are the lowest index of a column of $B_{*}\setminus A_{*}$.

Theorem B.3**.**

Let $\bar{C}$ be the bifiltered chain complex obtained by applying the proposed reduction algorithm to the bifiltered chain complex $C$. The number $\gamma^{p}_{k}(\bar{C})$ of $k$-generators added at value $p$ in $\bar{C}$ is equal to $\delta^{p}_{k}(C)$.

Proof B.4**.**

Let us apply Lemma B.1 to the case $A_{*}=C^{<p}_{*}$ and $B_{*}=C^{p}_{*}$ (considered as unfiltered chain complexes). Note that the local reduction phase at value $p$ in the chunk algorithm corresponds to a partial reduction of the matrix of $C^{p}_{*}$, because all column additions are left-to-right, and it is known that the lowest indices of the reduced matrix are invariant under left-to-right column additions. At the end of the local reduction phase, each column is labeled as local (positive or negative), or global. A local column is paired with another local column of value $p$, and this pairing will not change when completing the partial reduction to a completely reduced matrix for $C^{p}_{*}$. The reason is simply that all columns in $C^{<p}_{*}$ have a lowest index with smaller index than any column of value $p$, since the boundary matrix is upper triangular. By the above reasoning, the local columns will neither contribute to kernel nor cokernel of $\eta^{p}_{k-1}$.

On the other hand, consider a global column with value $p$. Its lowest index corresponds to a column in $C^{<p}_{*}$ after the local reduction phase. Since the lowest index can only decrease when further reducing the matrix, it follows that in the reduced matrix, the global column has still a lowest index in $C^{<p}_{*}$, in which case it contributes to $\mathrm{ker}\,\eta^{p}_{k-1}$, or its boundary becomes null. In the latter case, it cannot be paired with another column at value $p$ since then, this pairing would have been identified already at the local reduction phase. Hence, the global column contributes to $\mathrm{coker}\,\eta^{p}_{k}$. After all, each global column contributes to the kernel or the cokernel, and therefore, the number of global columns is equal to $\delta^{p}_{k}(C)$. Moreover, the number of global columns at value $p$ is exactly the number of generators added at value $p$ to the output complex $\bar{C}$.

Given an arbitrary bifiltered chain complex $D$ and $p=(p_{x},p_{y})\in\mathbb{R}^{2}$, let us consider the chain complex $D^{<p}_{*}$ defined as $\sum_{q<p}D^{q}_{*}$. Since we are considering just finitely generated chain complexes, there exists a sufficiently small $\epsilon$ for which, by defining $p_{1}=(p_{x}-\epsilon,p_{y})$ and $p_{2}=(p_{x},p_{y}-\epsilon)$, we have that, for any $p_{1}\leq q\leq p$, $D^{q}_{*}=D^{p_{1}}_{*}$, and, for any $p_{2}\leq q\leq p$, $D^{q}_{*}=D^{p_{2}}_{*}$. Then, $D^{<p}_{*}$ can be expressed as the finite sum $D^{p_{1}}_{*}+D^{p_{2}}_{*}$. Moreover, defining $p_{0}$ as the value $(p_{x}-\epsilon,p_{y}-\epsilon)$, the 1-criticality of the bifiltration $D$ ensures that $D^{p_{0}}_{*}=D^{p_{1}}_{*}\cap D^{p_{2}}_{*}$. For any arbitrary bifiltered chain complex $D$, let $\eta^{p}_{k}$ denote the homology map in dimension $k$ induced by the inclusion of $D^{<p}_{*}$ into $D^{p}_{*}$. Since each $k$-generator newly introduced in $D^{p}_{*}$ can act on the homology by just increasing the dimension of $\mathrm{ker}\,\eta^{p}_{k-1}$, or alternatively of $\mathrm{coker}\,\eta^{p}_{k}$, at most by 1, we have that

[TABLE]

Theorem B.5**.**

Let $D$ be a bifiltered chain complex quasi-isomorphic to $C$ (through the chain maps $f^{p}_{*}:C^{p}_{*}\rightarrow D^{p}_{*}$). Then, $\delta^{p}_{k}(C)\leq\gamma^{p}_{k}(D)$.

Proof B.6**.**

Let $f^{p}_{*}:C^{p}_{*}\rightarrow D^{p}_{*}$ be the collection of chain maps ensuring the quasi-isomorphism between $C$ and $D$. Thanks to their commutativity with the inclusion maps, we have that the chain map $f^{<p}_{*}:C^{<p}_{*}\rightarrow D^{<p}_{*}$ is well-defined. This ensures for any $p\in\mathbb{R}^{2}$, the construction of the following commutative diagram {diagram} in which each row is exact. Thanks to the quasi-isomorphism established by $f$ and to the Mayer-Vietoris sequence, we can derive from the previous one the following commutative diagram {diagram} in which each row is exact and the horizontal and the vertical maps are induced by the corresponding ones between chain complexes. Since all the other four vertical maps are isomorphisms, the $5$-lemma ensures that map $\psi^{<p}_{k}$ is also an isomorphism. So, we have the following commutative diagram {diagram} in which vertical maps are isomorphisms. Then, for any $k$, $\mathrm{dim}\,\mathrm{ker}\,\eta^{p}_{k}=\mathrm{dim}\,\mathrm{ker}\,\iota^{p}_{k}$ and $\mathrm{dim}\,\mathrm{coker}\,\eta^{p}_{k}=\mathrm{dim}\,\mathrm{coker}\,\iota^{p}_{k}$. So, by applying Equation (1), we have that

[TABLE]