Fast Free Resolutions of Bifiltered Chain Complexes

Multi-parameter persistence is a branch of topological data analysis where a data set (e.g., a point cloud) is filtered with respect to two or more parameters and the topological evolution of the data when changing the parameters is analyzed. In this context, the first step of a computational pipeline for two parameters typically consists of the computation of a bifiltration of simplicial complexes, that is, a family of simplicial complexes indexed by ${\mathbb{N}}^2$ that grows when increasing the parameters. We refer to the parameter set ${\mathbb{N}}^2$ as the grades of the bifiltration. A bifiltration can be equivalently described by determining the support of each simplex, that is, the set of grades at which the simplex is part of the complex. A bifiltration is called free or $1$-critical if the support of every simplex is a principal upset, that is, the upward closure of a single element in the parameter space. More generally, a bifiltration is $k$-critical if the support of a simplex is the upward closure of at most $k$ elements; see Figure 1 for a visualization.

A bifiltration gives rise to a bifiltered chain complex $C_{\bullet }$

with the $k$-simplices at a given grade forming the basis of the $k$-chains at that grade. Free and $k$-critical filtered chain complexes are defined analogously.

Bifiltered chain complexes arise naturally in various constructions for topological data analysis, and several constructions produce bifiltrations that are not free. The most prominent example is given by the degree-Rips bifiltration [19, 23, 24, 7]. Other examples are the approximate multi-cover bifiltration [1] and Delaunay core bifiltration [6]. On the other hand, free bifiltrations are most suitable for algorithmic and computational purposes. There are fast algorithms for minimizing a bifiltered chain complex without changing its homology [14, 15, 22] as well as for computing the homology of a chain complex in terms of a minimal presentation [20, 18, 15, 4], but both approaches require free chain complexes as input.

In homological algebra, a standard technique is to replace a general chain complex by a free one, connected to the original chain complex by a chain map that induces an isomorphism in homology (a quasi-isomorphism). The free chain complex together with the chain map is called a free resolution of the original complex. The problem studied in this paper is how to efficiently compute a free resolution of a non-free bifiltered chain complex.

At the homology level, one way to address the case of a non-free bifiltration was proposed by Chacholski, Scolamiero and Vaccarino [8]: given a segment $C:C_{m+1}\to C_m\to C_{m-1}$ of a chain complex, the authors describe an algorithm to compute a free chain complex $F:F_{m+1}\to F_m\to F_{m-1}$ such that $H_m(C)$ and $H_m(F)$ are isomorphic, that is, a free implicit representation. This algorithm suffices if one is interested in a presentation of $H_m(C)$ for further processing. Methods to compute a projective implicit representation from families of simplicial complexes and general simplicial maps have been developed in [11, 12], where, as in our approach, graph theoretic methods are used to speed up computations.

There are good reasons to work on the level of entire chain complexes instead. First of all, there are potential computational advantages, especially if one is interested in multiple homology dimensions (see Section 6). Moreover, the chain complex structure can encode subtle information on the data that is lost at the homology level: for example, two chain complexes may have isomorphic homology in every dimension without being quasi-isomorphic.

Computing free resolutions of general chain complexes is a standard task in computational algebra, available in the computer algebra software Macaulay2 for a much wider range of chain complexes. In the context of applied topology, we are mostly interested in very large bifiltered chain complexes with millions of generators but a simple combinatorial structure (i.e., simplicial boundary maps). The goal of this paper is to develop specialized and highly optimized algorithms for this type of input data, which the general purpose algorithms implemented in existing computer algebra systems are not tailored for.

Our main contribution is to propose two algorithms to compute a free resolution of a bifiltered chain complex. Both algorithms rely on the same simple idea of expanding any $k$-critical simplex in dimension $p$ into a sequence of $k$ free copies, with consecutive copies related by $(p+1)$-dimensional elements at the join of their grades. In algebraic terms, this corresponds to a free resolution of the upset module associated to the simplex. In order to construct a valid total complex, the algorithms introduce further maps to establish the chain complex property while maintaining quasi-isomorphism to the original filtered complex. Finding these maps is computationally inexpensive and takes place at a purely combinatorial level.

The algorithms differ mainly in the choice of free resolution of the upset modules. In the first algorithm, the path algorithm, the free resolution corresponds to the chain complex of a filtered path graph (as a simple special case of a cellular resolution [5, 21]). We give an example of a family of simplicial bifiltrations so that the resulting free chain complex has a dense boundary matrix (for any choice of basis). This example shows tightness of the worst-case runtime $O(n^2)$, where $n$ is the number of input generators (description size).

The second algorithm, the log-path algorithm, extends the path algorithm by adding additional relations to the free resolution of a simplex, such that any pair of copies of a simplex is connected via a sequence of relations of logarithmic length. This ensures sparsity of the boundary matrices in the output chain complex, but also requires adding further higher relations (syzygies). Again, further maps are required to establish the chain complex property of the resulting total complex, allowing the algorithm to maintain sparsity and obtain a resolution with worst-case run time in $O(n{\log }^{2}n)$.

Our findings lead to an interesting observation: as shown by our worst case example, a minimal resolution may require dense boundary matrices, while a non-minimal resolution may actually admit a sparse matrix representation, with asymptotically fewer non-zero entries ($O(n{\log }^{2}n)$ instead of $O(n^2)$, where $n$ is the description size of the complex). This observation suggests that minimizing chain complexes does not necessarily speed-up subsequent algorithmic tasks, at least in certain worst-case examples.

We provide implementations of the path and log-path algorithms, in addition to the Chacholski–Scolamiero–Vaccarino algorithm for computing free implicit representations of homology. Systematic tests on various $k$-critical bifiltrations show that the overhead of the log-path algorithm over the path algorithm does not exceed a factor of 3 in run time and a factor of 2 in the number of non-zero entries for our examples, while showing the expected improvement on the mentioned worst-case examples. Subsequently minimizing the free chain complex yields a further significant reduction of the size. Furthermore, we consider the task of computing minimal presentations of homology in all degrees, comparing the approach of first computing a “global” free resolution with the approach of computing free implicit representations. Our results show a clear computational advantage for the global approach. Remarkably, for some instances, computing all minimal presentations using a free resolution is faster than computing a single minimal presentation using the Chacholski–Scolamiero–Vaccarino algorithm.

A simplicial bifiltration $𝒦$ is an abstract simplicial complex together with a collection of subcomplexes ${({𝒦}_s)}_{s\in {\mathbb{N}}^2}$ such that ${𝒦}_s\subseteq {𝒦}_t$ whenever $s\le t$ (which means that $s_1\le t_1$ and $s_2\le t_2$). As shown in Figure 1, each simplex $\sigma$ enters $𝒦$ along a staircase bounding the support of $\sigma$, denoted $\mathrm{supp}(\sigma )$, which is an upset (an upward closed subset of ${\mathbb{N}}^2)$. There is a unique minimal set of grades $𝒢(\sigma ):=\{x_1,\mathrm{\dots },x_m\}\subseteq {\mathbb{N}}^2$ such that $\mathrm{supp}(\sigma )=\{s\in {\mathbb{N}}^2\mid \exists x_i\le s\}$. If $m\le k$ for each $\sigma \in 𝒦$, then the bifiltration is $k$-critical.

Simplicial bifiltrations give rise to bifiltered chain complexes. We first describe the elementary building blocks. Each $p$-simplex $\sigma$ has a support $\mathrm{supp}(\sigma )$ with minimal generating set $𝒢(\sigma )=\{x_1,\mathrm{\dots },x_m\}$, determining an upset module $U_{\{x_1,\mathrm{\dots },x_m\}}$ (which we also denote by $U_{\sigma }$ for brevity) given by

as illustrated in Figure 2. This is a special case of a bipersistence module $M$, which is a family ${(M_s)}_{s\in {\mathbb{N}}^2}$ of vector spaces over a field $\mathrm{𝕜}$ together with a family of homomorphisms ${(M_{s,t}:M_s\to M_t)}_{s\le t\in {\mathbb{N}}^2}$. Although our results extend to arbitrary fields straightforwardly, we will stick to the case $\mathrm{𝕜}={\mathbb{Z}}_2$ in order to simplify the exposition. A morphism of bipersistence modules $\phi :M\to N$ is a natural transformation, that is, a family of linear maps ${(f_s)}_{s\in {\mathbb{N}}^2}$ that commute with the structure maps of $M$ and $N$: $f_t\circ M_{s,t}=N_{s,t}\circ f_s$ for $s\le t$.

A direct sum of upset modules, each generated by a single grade, is called free. A basis of a free bipersistence module $F$ is a set of elements $\{b^{x_1},\mathrm{\dots },b^{x_n}\}$ where each $b^{x_i}\in F_{x_i}$, such that every $v\in F_y$, for any grade $y$, can be written in unique way as a linear combination of images of the basis elements under the structure maps (which are inclusions into the vector space ${\bigcup }_z{F}_z$). Note that the grades appearing in this linear combination with a non-zero coefficient must be less or equal to $y$. A morphism $f$ between free modules is then determined by the images of the basis elements for the domain, and can be encoded by a matrix $[f]$.

The upset modules associated to the simplices assemble to

which is a bifiltered chain complex with $C_i:={⨁}_{\sigma \in {𝒦}^{(i)}}{U}_{\sigma }$ and boundary maps ${\partial }_i:C_i\to C_{i-1}$ inherited from the simplicial complex $𝒦$, satisfying ${\partial }_i\circ {\partial }_{i+1}=0$. The notion of a $k$-critical bifiltered chain complex is defined analogously to the setting of a simplicial bifiltration.

Even though a bifiltered chain complex has a simple combinatorial structure, a free chain complex (a chain complex $F_{\bullet }$ with each $F_i$ free) is preferred in computational (and algebraic) settings. This motivates the goal of finding a free resolution, that is, a free chain complex together with a chain map $F_{\bullet }\to C_{\bullet }$ that induces an isomorphism $H_i(F_{\bullet })\cong H_i(C_{\bullet })$ on homology for each $i\ge 0$. Such a chain map is called a quasi-isomorphism.

The list of all simplices of a bifiltration $𝒦$, each with a lexicographically ordered list of the minimal generating set of its support, and a list of its facets gives a full combinatorial description of $𝒦$ and will serve as the input to our algorithms. This data representation generalizes to regular cell complexes, and in our case to bifiltered chain complexes where $\mathrm{𝕜}={\mathbb{Z}}_2$. For the simplicial case, the number of facets is constant for each simplex dimension, ensuring sparsity. Therefore the description size of the input is given, up to a constant factor, by the total number of generators in this list. Our output is a free resolution $F_{\bullet }$ of our input chain complex $C_{\bullet }$, again represented in the same data format.

To compute the homology $H_m(C_{\bullet })$ of a $k$-critical chain complex in a chosen dimension $m$, Chacholski, Scolamiero, and Vaccarino [8] provide a procedure to construct a free chain complex segment $X\stackrel{𝑓}{\leftarrow }Y\stackrel{𝑔}{\leftarrow }Z$ from the input segment

such that $H_m(C_{\bullet })\cong \ker f/\mathrm{im}g$. The pair of graded matrices $([f],[g])$ is called a free implicit representation (FI-rep) of $H_m(C_{\bullet })$ [20], and it serves as input for the computation of minimal presentations of $H_m(C_{\bullet })$ [20, 18].

We consider an input complex $C_{\bullet }$ induced by a simplicial bifiltration $𝒦$. The first step is to extend each $U_{\sigma }$ in $C_{m-1}$, induced by an $(m-1)$-simplex $\sigma$, to a free upset module $U_{\{0\}}$. The resulting free module $D_{m-1}≔{⨁}_{\sigma \in {𝒦}^{(m-1)}}{U}_{\{0\}}$ , whose basis we denote by ${\{b_{\sigma }^0\}}_{\sigma \in {𝒦}^{(m-1)}}$, contains $C_{m-1}$ as a submodule, since the support of $U_{\sigma }$ is contained in the support of $U_{\{0\}}$. Thus, postcomposing ${\partial }_m$ with this submodule inclusion does not affect the kernel: we have $\ker {\partial }_m=\ker \iota \circ {\partial }_m$.

As a second step, we cover each $U_{\sigma }$ in $C_{m+1}$, induced by an $(m+1)$-simplex $\sigma$ by a free bipersistence module. This means that we replace each upset module $U_{\sigma }$ by the free module $G_{\sigma }≔{⨁}_{x_i\in 𝒢(\sigma )}{U}_{\{x_i\}}$, connected to the upset module by a surjection $G_{\sigma }\twoheadrightarrow U_{\sigma }$ since the support of $G_{\sigma }$ equals the support of $U_{\sigma }$. We call the basis elements $g_{\sigma }^{x_i}$ the generators of $U_{\sigma }$. The module $G_{m+1}$ arises from $C_{m+1}$ by covering $U_{\sigma }$ by $G_{\sigma }$ for each $(m+1)$-simplex via a canonical surjection . Hence, precomposing ${\partial }_{m+1}$ with ${\alpha }_{m+1}$ does not affect the image. The input (1) can now be replaced by (2) which has isomorphic homology.

As the third and final step it remains to replace $C_m$. We first construct $G_m$ analogously to $G_{m+1}$, by substituting each $U_{\sigma }$ in $C_m$ by $G_{\sigma }$. Precomposing $\iota \circ {\partial }_m$ with the surjection yields a map $\gamma ≔\iota \circ {\partial }_m\circ {\alpha }_m:G_m\to D_{m-1}$, as depicted in (3). For $\sigma \in {𝒦}^{(m)}$ with $\partial \sigma ={\tau }_0+\mathrm{\cdots }+{\tau }_m$, the resulting map $\gamma$ sends a generator $g_{\sigma }^x$ in $G_m$ to $b_{{\tau }_0}^0+\mathrm{\cdots }+b_{{\tau }_m}^0$ in $D_{m-1}$.

Replacing $C_m$ by $G_m$ makes it necessary to also replace the map ${\partial }_{m+1}\circ {\alpha }_{m+1}:G_{m+1}\to C_m$ by a map $f_{m+1}^0:G_{m+1}\to G_m$, representing ${\partial }_{m+1}$ on the generators and also making the diagram commute. This map $f_{m+1}^0$ sends each generator $g_{\sigma }^x$ of an $(m+1)$-simplex $\sigma$ with $\partial \sigma ={\tau }_0+\mathrm{\cdots }+{\tau }_m$ to $f_{m+1}^0(g_{\sigma }^x)=g_{{\tau }_0}^{y_0}+\mathrm{\cdots }+g_{{\tau }_m}^{y_m}$, where each $g_{{\tau }_i}^{y_i}$ is a generator of ${\tau }_i$ chosen such that its grade satisfies $y_i\le x$; such a $y_i$ always exists because the faces ${\tau }_i$ of a simplex $\sigma$ are present in the bifiltration whenever $\sigma$ itself is present. The map $f_{m+1}^0$ thus makes the square in (3) commute and is therefore called a lift of ${\partial }_{m+1}$. See Figure 3 for an illustration.

Note that the kernel of $\gamma$ may not be isomorphic to the kernel of ${\partial }_m$; the surjection ${\alpha }_m$ maps each generator $g_{\sigma }^x$ of an $m$-simplex in $G_m$ to the same potential cycle in $C_m$, thus increasing the dimension of the kernel. This can be resolved by relating these generators appropriately. Assume that the generators $g_{\sigma }^{x_1},\mathrm{\dots },g_{\sigma }^{x_k}$ of $U_{\sigma }$ are ordered w.r.t. the first coordinate of their grades. To represent $U_{\sigma }$ correctly, two consecutive generators $g_{\sigma }^{x_i}$ and $g_{\sigma }^{x_{i+1}}$ have to be identified at the join $y_i=x_i\vee x_{i+1}$ of their grades by a relation $r_{\sigma }^{y_i}$. With the free bipersistence module $R_{\sigma }≔{⨁}_{y_i}{U}_{\{y_i\}}$ with basis ${\{r_{\sigma }^{y_i}\}}_i$, we can present $U_{\sigma }$ by generators and relations via

where $p_{\sigma }^1(r_{\sigma }^{y_i})=g_{\sigma }^{x_i}+g_{\sigma }^{x_{i+1}}$. Since $p_{\sigma }^1$ is injective, ${\alpha }_{\sigma }$ is surjective, and $\ker {\alpha }_{\sigma }=\mathrm{im}{p}_{\sigma }^1$, the sequence in (4) is a short exact sequence. This means that $(G_{\sigma }\stackrel{p_{\sigma }^1}{\leftarrow }{R}_{\sigma },{\alpha }_{\sigma })$ already determines a free resolution of $U_{\sigma }$, as defined in Section 2, where $U_{\sigma }$ is considered as a chain complex concentrated in degree $0$.

The free resolution in Diagram 4 has a simple combinatorial structure, given by a path graph ${𝒫}_{\sigma }$ with $k$ vertices corresponding to the generators $g_{\sigma }^{x_i}$ and $k-1$ edges $\{g_{\sigma }^{x_i},g_{\sigma }^{x_{i+1}}\}$ corresponding to the relations $r_{\sigma }^{y_i}$. Note that vertices and edges implicitly carry the grades of the generators and relations. The morphism $p_{\sigma }^1$ is represented by the (graded) incidence matrix of the graph. We call this free resolution the path resolution of $U_{\sigma }$, see Figure 4. Note that a path resolution is a special case of a cellular resolution [5, 21].

We now use the free resolution to replace the map ${\partial }_{m+1}\circ {\alpha }_{m+1}$ in (2), assembling the following output complex from (3)

By construction, the homology of this complex segment in degree $m$ is isomorphic to that of (2), and hence to that of the input complex.

We now turn to the problem of finding a free resolution, that is, a chain complex of free modules that is quasi-isomorphic to the input $C_{\bullet }$, thus preserving homology in all dimensions. For this problem, extending the codomain of ${\partial }_i:C_i\to C_{i-1}$ to a free module generated in degree $0$, as done in the Chacholski–Scolamiero–Vaccarino algorithm, is not feasible anymore, since this changes $H_{i-1}(C_{\bullet })$. Instead, we now carry out the step of substituting $C_i$ by its free resolution $G_i\stackrel{p_i^1}{\leftarrow }{R}_i\leftarrow 0\leftarrow \mathrm{\cdots }$ in every dimension $i$ and replacing ${\partial }_{i+1}$ by the (zeroth) lifts $f_{i+1}^0:G_{i+1}\to G_i$, yielding maps as illustrated in (6).

However, the sequence of maps, arising from this construction, does not yield a chain complex. Given a relation $r_{\sigma }^{y_j}$ in $R_i$, we obtain

This term (also the term $f_{i-1}^0\circ f_i^0$) may not vanish, since the generators of $\sigma$ may be sent to different generators of its boundary simplices ${\tau }_j$ by $f_i^0$, i.e., $g_{{\tau }_j}^{z_j}\ne g_{{\tau }_j}^{z_j^{\prime }}$. However, it vanishes modulo relations, in the sense that its representative vanishes in $G_{i-1}/\mathrm{im}{p}_{i-1}^1$, see Lemma 1.

This can be seen combinatorially. The free resolution of $U_{\tau }$ corresponds to the path graph ${𝒫}_{\tau }$, in the sense that the restriction of $p_{i-1}^1$ to a map $R_{\tau }\to G_{\tau }$ corresponds to the boundary map of the path graph. The relations $r_{\tau }^{w_j},\mathrm{\dots },r_{\tau }^{w_{\mathrm{\ell }-1}}$ then define the unique path connecting the vertices $g_{\tau }^{x_j}$ and $g_{\tau }^{x_{\mathrm{\ell }}}$ in ${𝒫}_{\tau }$. $◀$

The map $p_i^1$ sends the edge $r_{\sigma }^{y_j}$ of ${𝒫}_{\sigma }$ to the vertices $g_{\sigma }^{x_j}$ and $g_{\sigma }^{x_{j+1}}$ which in turn are mapped to $g_{{\tau }_j}^{z_j}$ and $g_{{\tau }_j}^{z_j^{\prime }}$ on ${𝒫}_{{\tau }_j}$ by $f_i^0$. We can therefore modify the boundary map in (6) by additionally sending $r_{\sigma }^{y_j}$ to the path $r_{{\tau }_j}^{w_1}+\mathrm{\cdots }+r_{{\tau }_j}^{w_{\mathrm{\ell }}}$ guaranteed by Lemma 1, see Figure 5. This defines a correction term $f_i^1:R_i\to R_{i-1}$ such that $p_{i-1}^1\circ f_i^1=f_i^0\circ p_i^1$. The maps $f_i^1$ are called the first lifts of ${\partial }_i$.

Even after adding the correction term $f_i^1$ to (6), the resulting maps may still fail to be boundary maps, since the composition $f_{i-1}^0\circ f_i^0(g_{\sigma }^x)$ need not vanish for each generator $g_{\sigma }^x$ in $G_i$, as seen in Figure 6. Again, $f_{i-1}^0\circ f_i^0(g_{\sigma }^x)$ is a sum consisting of pairs $g_{\tau }^y,g_{\tau }^{y^{\prime }}$ of possibly different generators of the same simplex $\tau$, which can be connected by a path of relations as above. A second correction term $h_i^0:G_i\to R_{i-2}$ is defined analogously to the lifts $f_i^1$. By construction, $h_i^0$ satisfies $p_{i-2}^1\circ h_i^0=f_{i-1}^0\circ f_i^0$, making the following diagram commute:

Adding $h_i^0$ and $f_i^1$ to the boundary maps in (6) yields the desired free resolution of $C_{\bullet }$. The next theorem, proved in A.1 of the full version [3], shows correctness of the path algorithm.

For the description of the algorithm and its asymptotic time complexity, we assume that the input $C_{\bullet }$ is induced by a $k$-critical bifiltration of a simplicial complex with fixed dimension $d$, given as a list of simplices $\sigma$, each with its dimension, grades $𝒢(\sigma )$ and boundary $\partial \sigma$. The output is described by a list of generators and relations (each carrying a grade) as well as their images under the maps $p_i^1,f_i^0,f_i^1$ and $h_i^0$ (which together form the boundary maps of the free resolution). The output is initialized as an empty list.

1. 1.

   Compute the bases of $G_i\oplus R_{i-1}$ and the matrices $[p_i^1]$. For each $i$-simplex $\sigma$ and $x_j\in 𝒢(\sigma )$, add a generator $g_{\sigma }^{x_j}$ to $G_i$. For any two consecutive $x_j,x_{j+1}\in 𝒢(\sigma )$, add a relation $r_{\sigma }^{y_j}$ with $y_j=x_j\vee x_{j+1}$ to $R_i$. Moreover, define $p_i^1(r_{\sigma }^{y_j})≔g_{\sigma }^{x_j}+g_{\sigma }^{x_{j+1}}$.

2. 2.

   Compute the matrices $[f_i^0]$. For each generator $g_{\sigma }^x$ in $G_i$, find generators $g_{{\tau }_0}^{x_0},\mathrm{\dots },g_{{\tau }_i}^{x_i}$ for ${\tau }_0+\mathrm{\cdots }+{\tau }_i=\partial \sigma$ such that $x_j\le x$ and define $f_i^0(g_{\sigma }^x)≔g_{{\tau }_0}^{x_0}+\mathrm{\cdots }+g_{{\tau }_i}^{x_i}$.

3. 3.

   Compute $[f_i^1]$ and $[h_i^0]$. For $[f_i^1]$, perform a matrix multiplication $[f_i^0][p_i^1]$. The resulting columns are indexed by the relations $r_{\sigma }^y$ in $R_i$ and consist of one pair of generators $g_{\tau }^{z_j},g_{\tau }^{z_{\mathrm{\ell }}}$ for each facet $\tau$ of $\sigma$. Such a pair is connected by a path $c_{\tau }:=r_{\tau }^{w_j}+\mathrm{\cdots }+r_{\tau }^{w_{\mathrm{\ell }}}$ in ${𝒫}_{\tau }$, and we set $f_i^1(r_{\sigma }^y):={\sum }_{\tau }{c}_{\tau }$; see Figure 7. For $[h_i^0]$, proceed analogously to multiply $[f_{i-1}^0][f_i^0]$. The resulting columns are indexed by the generators $g_{\sigma }^x$ in $G_i$ and consist of one pair of generators $g_{\tau }^{z_j},g_{\tau }^{z_{\mathrm{\ell }}}$ in $G_{i-2}$ for each codimension $2$ face $\tau$ of $\sigma$. Such a pair is connected by a path $c_{\tau }:=r_{\tau }^{w_j}+\mathrm{\cdots }+r_{\tau }^{w_{\mathrm{\ell }}}$ in ${𝒫}_{\tau }$, and we set $h_i^0(g_{\sigma }^x):={\sum }_{\tau }{c}_{\tau }$.

The running time of the path algorithm is $O(nk)$ for a $k$-critical bifiltration of description size $n$; see A.2. in the full version [3]. This runtime is worst-case optimal, since the description size of a computed free resolution can be $\mathrm{\Omega }(nk)$ in the worst-case. The reason is that the paths in the computation of $f_i^1$ and $h_i^0$ can be of size up to $k$, and many simplices may require such long paths in the worst case. Figure 8 gives a construction for which the matrix $[h_2^0]$ is of size $\mathrm{\Omega }(nk)$; see Figure 15 in [3] for a similar construction for $[f_1^1]$.

A free chain complex is minimal if there is no quasi-isomorphic free chain complex with a smaller number of generators. In general, free resolutions computed by the path algorithm are not minimal. However, a slight modification of the construction in Figure 8 yields the following result (see A.3. in the full version [3] for details).

In the preceding section, we leveraged the fact that the upset modules induced by simplices of a bifiltration admit path-shaped resolutions, allowing us to compute the correction terms $f_i^1$ and $h_i^0$ in a simple way. In the worst case, many generators and relations are mapped to long paths, making the matrices $[f_i^1]$ and $[h_i^0]$ dense ($\mathrm{\Omega }(nk)$ size) and leading to a quadratic running time of $\mathrm{\Theta }(nk)$ for the path algorithm.

Such long paths can be avoided by introducing shortcuts in the path resolution ${𝒫}_{\sigma }$ (4), as illustrated in Figure 9. Recall that the grades $𝒢(\sigma )=\{x_1,\mathrm{\cdots },x_{n_{\sigma }}\}$ generating $\mathrm{supp}(\sigma )$ are totally ordered (by their first coordinate). Now any two vertices $g_{\sigma }^{x_j}$ and $g_{\sigma }^{x_{\mathrm{\ell }}}$ with are connected by an additional shortcut edge if there are numbers $s\ge 0,t>0$ with $j=s{2}^t$ and $\mathrm{\ell }=(s+1)2^t$. Extending ${𝒫}_{\sigma }$ by this edge corresponds to adding a relation $r_{\sigma }^z$ to $R_{\sigma }$, where $p_{\sigma }(r_{\sigma }^z)=g_{\sigma }^{x_j}+g_{\sigma }^{x_{\mathrm{\ell }}}$ and $z=x_j\vee x_{\mathrm{\ell }}$. In total, only $O(n_{\sigma })$ many edges are added.

This construction ensures that any two vertices $g_{\sigma }^{x_j}$ and $g_{\sigma }^{x_{\mathrm{\ell }}}$ with are connected by a monotone path of length logarithmic in $n_{\sigma }$. A shortest monotone path can be constructed in a greedy way. We start in $g_{\sigma }^{x_j}$ and take the longest possible edge in each step that does not overshoot $g_{\sigma }^{x_{\mathrm{\ell }}}$. We conclude the following extension of Lemma 1 (proved in [3, B.1]).

Lemma 4 appears to resolve the size issue: the same algorithm as in the previous section can be used, except that in the construction of $h_i^0$ and $f_i^1$ use Lemma 4 instead of Lemma 1. Then, the number of nonzero boundary coefficients for every generator and relation has only size $O(\log n_{\sigma })$ instead of $O(n_{\sigma })$, which would lead to an output boundary matrix with $O(n\log k)$ nonzero entries.

Now recall that a key reason for why the path construction in (7) yields a chain complex is the uniqueness of paths between vertices. This leads to the vanishing of the terms $h_i^0\circ f_{i+1}^0+f_{i-1}^1\circ h_{i+1}^0$ and $h_i^0\circ p_i^1+f_{i-1}^1\circ f_i^1$ when composing boundary maps in (7). However, after adding shortcut edges the paths between vertices are no longer unique. As illustrated in Figure 11, for a relation $r_{\sigma }^y$, the terms $h_i^0\circ p_i^1(r_{\sigma }^y)$ and $f_{i-1}^1\circ f_i^1(r_{\sigma }^y)$, when restricted to a component, may correspond to different paths connecting the same endpoints and may therefore enclose a cycle.

Any such cycle constitutes an obstruction to the vanishing of $h_i^0\circ p_i^1+f_{i-1}^1\circ f_i^1$. To eliminate these obstructions, we fill the cycles with faces, or, in algebraic terms, introduce syzygies, thus extending the modified relations to a free resolution, as depicted in Figure 10:

The syzygies $s_{\sigma }^z$ for $S_{\sigma }$ correspond to the triangles with vertex grades $(x_{(s-1)2^t},x_{s{2}^t},x_{(s+1)2^t})$ for any odd $s>0$; see Figure 9. The grade $z$ is given by the grade of the longest edge, that is, $x_{(s-1)2^t}\vee x_{(s+1)2^t}$. The map $p_{\sigma }^2$ is defined by sending $s_{\sigma }^z$ to the sum of the three edges of the triangle. We call (8) the log-path resolution of $U_{\sigma }$, noting that it is a bifiltered two-dimensional simplicial chain complex ${ℒ}_{\sigma }$. Its construction yields the following property:

To fill a cycle in ${ℒ}_{\sigma }$, we decompose it into simple cycles and use the observation that a simple cycle is filled with the triangles corresponding to the non-extremal vertices $x_{s{2}^t}$ of the cycle, see Figure 12. See B.2 of the full version [3] for details.

We now describe the algorithm, again assuming that the input is a $k$-critical simplicial bifiltration of size $n$ and constant dimension $d$ with chain complex $C_{\bullet }$.

Step 1: Compute the log-path resolution ${ℒ}_i:G_i\stackrel{p_i^1}{\leftarrow }{R}_i\stackrel{p_i^2}{\leftarrow }{S}_i$ of each module $C_i$. This requires iterating through the generators of each input simplex and adding a total of $O(n)$ relations and syzygies, which takes $O(n)$ time in total.

Step 2: Compute $f_i^0:G_i\to G_{i-1}$ in $O(n\log k)$ time as in Section 4.

Step 3: Compute the maps $f_i^1:R_i\to R_{i-1}$ and $h_i^0:G_i\to R_{i-2}$ as in the path algorithm, except that paths between $g_{\sigma }^{x_j},g_{\sigma }^{x_{\mathrm{\ell }}}$ coming from the same simplex $\sigma$ are computed via Lemma 4. We can find every path in $O(\log k)$ time and compute matrices $[f_i^1]$ and $[h_i^0]$ with at most $O(\log k)$ entries per column in a total running time of $O(n\log k)$ for this step.

Step 4. All maps computed in the previous steps yield the following extension of (7).

A construction as in Section 4 leads to boundary morphisms whose composition is non-zero. Again, the boundary defect is repaired by introducing correction maps $f_i^2$, $h_i^1$, and $H_i^0$:

Commutativity of Diagram (9) yields the following observation

The log-path resolution ${ℒ}_i$ of $C_i$ is the direct sum of the simplicial chain complexes ${ℒ}_{\sigma }$, taken over each simplex $\sigma$ in $C_i$. The map $p_i^1$ is given by the $1$-dimensional boundary map of this complex, and hence its kernel is the collection of cycles in the $1$-skeletons of ${ℒ}_{\sigma }$. More specifically, each triangle $s_{\sigma }^z$ in $S_{i+1}$ is sent by $p_{i+1}^2$ to its boundary, a sum of three edges.

By Lemma 6, applying $f_{i+1}^1$ to this boundary yields another cycle in $R_i$, which by construction of ${ℒ}_i$ decomposes as ${\sum }_{\tau }{c}_{\tau }$ for $\tau$ running over the facets of $\sigma$. For each such cycle $c_{\tau }$, by Lemma 5 there is a unique chain $b_{\tau }$ in $S_i$ such that $p_i^2(b_{\tau })=c_{\tau }$, as illustrated in Figure 13. The map $f_{i+1}^2$ is then defined via $f_{i+1}^2(s_{\sigma }^z):={\sum }_{\tau }{b}_{\tau }$, with $\tau$ running over the facets of $\sigma$. Thus $f_{i+1}^2$ repairs the boundary defect $f_{i+1}^1\circ p_{i+1}^2$ as $f_{i+1}^1\circ p_{i+1}^2=p_i^2\circ f_{i+1}^2$.

Similar arguments lead to correction terms $h_i^1$, and $H_i^0$ satisfying

Their construction is analogous to the above, and we omit the details.

Note that in each of the three cases $f_{i+1}^2$, $h_i^1$, and $H_i^0$, the construction involves matrix multiplication (evaluation of the left hand sides in (12)), but also cycle filling (invoking Lemma 5). The former is efficient, as paths of logarithmic length lead to sparse matrices $[h_i^0]$ and $[f_i^1]$. More precisely, the column sparsity of $[f_i^1]$ and $[h_i^0]$ is $O(\log k)$, while in $[p_i^2]$ every column has exactly three non-zero entries. All sparse matrix products in Step $4$ can thus be computed in $O(n{\log }^{2}k)$ and the worst case column sparsity of the results is $O({\log }^{2}k)$. The resulting matrices contain cycles of length $O({\log }^{2}k)$, which can be filled by triangles in $O({\log }^{2}k)$ time. Overall, Step $4$ takes $O(n{\log }^{2}k)$ time.

We can now assemble the output chain complex as

For a proof, see chapters B.3 and B.4 of the full version [3]. We summarize our observations on the running time in the following theorem.

Applying this result to the modified wheel example ([3, A.3]) yields a non-minimal free resolution with smaller description size than any minimal one. Thus, perhaps surprisingly, minimal free resolutions are not asymptotically optimal in terms of description size.

We ran our experiments on a total of 342 test instances from different sources. First, we generated degree-Rips bifiltrations for (noisy) point samples drawn from a torus embedded in ${ℝ}^5$, from the “swiss roll” embedded in ${ℝ}^3$, and from the $7$-dimensional unit cube, using the python package tadasets. In all cases, we generated the complex up to $3$-simplices.

Furthermore, we used a simple method for generating bifiltrations of a simplicial complex equipped with two non-negative real-valued functions $f_0$ and $f_1$ on its simplices, described in Figure 14. We employed this construction on meshes from the Aim@Shape repository²²2http://visionair.ge.imati.cnr.it/, using the squared mean curvature and the distance to the barycenter as the two functions. Furthermore, we bifiltered the Delaunay triangulations from the same point samples as above (torus, cube, swiss roll) using the minimum enclosing radius of a simplex and the average distance to the $\lfloor \log n\rfloor$-nearest neighbors as filter functions. In all cases, we generated $k$-critical instances with $k=2,4,8$. We computed the filtration values with CGAL [10, 26, 13].

Finally, we generated wheel bifiltrations as in Figure 8 for different sample sizes.

Table 1 shows some of the results obtained from both algorithms on these datasets; the results for the other instances are similar.

In almost all cases, the log-path algorithm is slower than the path algorithm by a factor of up to 3 and produces an output complex of around twice the size. Also, it uses around twice as much memory (not shown in the table). Moreover, minimizing the chain complex (rightmost column of Table 1) further reduces the complex size significantly, generally at a small cost – sometimes it even saves time because the time used for minimization is less than the cost of producing the larger output file.

The only exception is the wheel bifiltration, where our experiments show the expected asymptotic worst-case behavior, with the log-path algorithm outperforming the path algorithm. In this case, minimizing the chain complex produced by the log-path algorithm will necessarily introduce a dense matrix and destroy the advantage of the log-path algorithm – for instance, its running time in the instance in the last row of Table 1 is 33 seconds, with 14.4 seconds for minimization itself and 18 seconds to write the 2.3 GB output file.

For a multi-critical filtered chain complex, we investigate the task of computing minimal presentation matrices for homology. We provide two variants: In the first variant, we compute a free resolution of the input chain complex with the path algorithm (which is usually faster than the log-path algorithm, as seen in the previous experiment) and subsequently minimize using multi-chunk. We then compute minimal presentations for each homology degree separately using the mpfree library [18, 15].

The second variant operates by directly structuring the input complex into chain complex segments per degree. It then uses the Chacholski–Scolamiero–Vaccarino algorithm to compute a free implicit representation for each segment, before finally employing mpfree once more to generate a minimal presentation.

The results of this experiment are presented in Table 2. For the degree-Rips instances, the first variant provides a modest improvement over the second. The reason is that the vast majority of simplices are in the top dimension $3$, so that computing the presentation for $H_2$ is the bottleneck in the computation. This step, however, does not differ significantly in both approaches: most of the time is spent to determine which $3$-simplices are killing $2$-cycles, by reducing the boundary matrix for $2$- and $3$-simplices.

For the bifunction instances, the speed-up of the first variant using free resolutions is much more pronounced. The reason is that these instances have a more balanced distribution of simplices over different dimensions: while for the second variant, the algorithm for $H_k$ still has to find the bounding $(k+1)$-simplices for every $k$, the multi-chunk algorithm makes use of the clearing optimization [9] and hence avoids the reduction of large parts of the boundary matrices. Remarkably, this technique is so effective that computing all minimal presentations via free resolutions is faster than computing a single minimal presentation via the approach using the Chacholski–Scolamiero-Vaccarino algorithm, even though restricting to a single dimension $k$ allows this approach to disregard all chains in dimensions other than $k+1$, $k$, or $k-1$.

While the Macaulay2 software includes a general freeResolution routine [25], it is unsuitable for our purposes. Its use requires converting our bigraded input into a graded chain complex over ${\mathbb{Z}}_2[x,y]$, a conversion that we found to be slow itself. More importantly, the core computation in Macaulay2 was orders of magnitude slower than our method, even on much smaller instances. We infer that the software was not designed for the large inputs typical of TDA. The conversion script is available on request.