Decomposing Multiparameter Persistence Modules

We present an algorithm inspired by the Dey-Xin algorithm [29], speeding up various steps and lifting the assumption that the module needs to be distinctly graded. Our methods are collected in the software package aida (Figure 1), which allows the decomposition of presentations arising from large data within seconds. We demonstrate the usefulness of our improvements in extensive experimental evaluation. Since the methods used in aida are of more general interest we implemented them in a stand-alone, header-only C++library Persistence-Algebra. Both this library and the aida code are publicly available.¹¹1https://github.com/JanJend/AIDA and https://github.com/JanJend/Persistence-Algebra

While the distinctly graded case allows for an incremental column-by-column approach to compute a decompositon, the general case requires to consider all columns of the same degree at once, say $k_i$ many at the $i$-th step. It does not suffice to consider all permutations of these columns [5]; instead arbitrary changes-of-base within the columns could be part of a correct solution. Our first approach requires a brute-force iteration over these change-of-basis matrices (which is only possible over a finite field) and is therefore not polynomial in $k_i$. Nevertheless, through a careful analysis, we limit the number of cases and provide a tractable solution if $k_i$ is small. Since we can then use the poly-time strategy of [29], this yields a decomposition algorithm which is fixed-parameter tractable with respect to $k_{\max }≔{\max }_i{k}_i$.

To further reduce the iteration – or even avoid it altogether – we analyse which kinds of matrix operations can actually help in the decomposition. This leads to another FPT-algorithm in a parameter ${\kappa }_i$ smaller than $k_i$ which promises to accelerate the computation for modules where $k_i$ is too large. A modification of this idea also yields an algorithm which, for interval-decomposable modules, runs in $𝒪(n^3)$ time, independent of ${\kappa }_i$.

Multi-parameter persistent homology is an active research area with many theoretic, computational, and applied results in recent years  [15]. In particular, efficient algorithms have been devised for the computation of a presentation of standard bifiltrations [4, 2, 17, 18, 25, 44], for minimizing presentations [29, 33, 47] and for efficient distance computation [10, 12, 28]. Our decomposition algorithm complements this pipeline.

In theory we can decompose any module over a finite dimensional algebra in polynomial time by decomposing its endomorphism algebra [23], but this algorithm is not practical. The algorithms in the C MeatAxe Package ([50], [57]), in the MAGMA system [14], and in the GAP system [34], [54], rely on the Holt-Rees extension [37] of the MeatAxe algorithm [53], both of which are Las Vegas algorithms and too slow for the large sparse matrices appearing in TDA. We remark that that the run-time complexity of these algorithms is sometimes not precisely cited in the TDA literature. For example the bounds in [15, Section 8.5] and [29, Introduction]) are a misinterpretation of the usage of the MeatAxe algorithm in this pipeline and use the lower bound given in [36] for the non-practical algorithm described in [55].

The only algorithm specifically designed for MPM is the aforementioned Dey-Xin algorithm [29] which has a time complexity of $𝒪(n^{2\omega +1})$, where is the matrix multiplication constant. It only works for uniquely graded modules, a property that is satisfied by some, but not all modules that appear in practice. Moreover, an algorithm to compute graphcodes [43] yields a partial decomposition

It is well known that decompositions are not stable; every persistence module is even arbitrarily close to an indecomposable one in interleaving distance [9], but geometric bifiltrations typically do not create such pathologies. Moreover, Bjerkevik [11, Thm 5.4] constructs a module which captures all decompositions of an entire $ϵ$-interleaving-neighbourhood. His construction creates many generators and relations of the same degree so its decomposition requires an efficient algorithm we provide with AIDA.

Complicated indecomposables can arise already in relatively simple geometric examples  [20, 19] and most modules in practice are not interval-decomposable (e.g., most examples from our benchmark set contain non-interval indecomposables; see also [3]). However, we expect that in practice most modules will mainly consist of small or structurally simple indecomposables (see [1]), and we designed our algorithm to work fast in such scenarios.

We review persistence modules and their presentations in Section 2. The algorithm by Dey and Xin [29], which is the basis of our work, is reviewed and modified in Sections  3 and 4. In Section 5, we describe how to decompose non-distinctly graded modules with an exhaustive approach. In Section 6, we present an algorithm that reduces the potentially restrictive iteration of the exhaustive algorithm. The detailed descriptions of our algorithms and proofs of correctness are deferred to the extended arXiv version [26].

A persistence module can be succinctly represented by a presentation. This is a set of generators and a set of relations among these generators. If these sets are finite, they form a presentation matrix, where the row indices correspond to the generators, the column indices to the relations and an entry at $(i,j)$ tells us which scalar coefficient in $𝕂$ the $i$-th generator appears with in the $j$-th relation.

Decompositions of modules correspond to decompositions of their presentations, but this correspondence is only surjective, see [26, Example 2.30] for details.

Admissible operations preserve the isomorphism type of the presented module.

In particular any decomposition of a module lifts to a change-of-basis between a presentation matrix and a block-decomposed presentation matrix.

If we can decompose $\alpha$-decomposed presentations, then by iterating over the relations we can decompose all finite presentations. If $N$ has only one column this means we have to zero out every sub-column $N_b$ of $N$ whenever possible to reach the block-decomposition with the most blocks. This is done by a subroutine which Dey and Xin called BlockReduce, a name we will keep.

The operations of type $(a)$, $(b)$, and $(c)$ have a pleasant feature: their effect on the matrix does not change if the order of their application changes. Dey-Xin call such sets of operations independent. In the setting of 10, when restricting the entries of $Q$ and $P$ to those belonging to independent operations, the system of equations becomes linear.

Dey and Xin observed that following this strategy, repeated application of the operations (a)–(c) actually suffices to decompose an $\alpha$-decomposed presentation with $k=1$, but assumed that the generators are also uniquely graded. We prove a more general statement, which therefore also lifts this assumption from the GP Algorithm.

At last assume again that $k=1$. After having zeroed out as much as possible of $N$, the blocks $b$ for which $N_b$ are still non-zero must all merge into a new block, since their non-zero entries now overlap in this last column. This finishes the loop of the GP Algorithm.

Consider an $\alpha$-decomposed presentation $\left[M_{ℬ}N\right]$ and let $k\ge 1$ now be arbitrary, where we recall that $k$ is the number of columns of $N$. At first, we are still interested in deleting the entire sub-matrices $N_b$ for every $b\in ℬ$ wherever this is possible. Consider now sets of admissible operations of the types $(a)$-$(c)$ in 12, encoded as graded matrices. Using $M_{b,c}≔M_{b_{\text{rows}}\times c_{\text{cols}}}$, their effect on $\left(M_{ℬ}N\right)$ can be described:

1. (a)

   A graded matrix $U\in {𝕂}^{R_b\times k\cdot \alpha }$ and $N_b\leftarrow N_b+M_{b}U$.

2. (b)

   A graded matrix $Q\in {𝕂}^{G_b\times G_c}$ and $N_b\leftarrow N_b+Q{N}_c,M_{b,c}\leftarrow M_{b,c}+Q{M}_c$.

3. (c)

   A graded matrix $P\in {𝕂}^{R_b\times R_c}$ and $M_{b,c}\leftarrow M_{b,c}+M_{b}P$.

Finding $Q,P$ such that after performing $(b)$ and $(c)$ the matrix $M_{ℬ}$ has not changed is equivalent to demanding that $M_{b,c}=0$, which is a linear system ($\ast$). Finding $U,Q$ such that performing $(a)$, $(b)$ zeroes out $N_b$ is another one ($\ast \ast$).

Input: $\left[M_{ℬ}N\right]\in K^{G\times (R\cup k\cdot \alpha )}$, $b\in ℬ$.

In [29], Dey and Xin write these equations in a large matrix which they call $S$-matrix.

Let $(M_{ℬ}N)$ be $\alpha$-decomposed with $M$ of size $m\times n$, then the size of every block is also bounded by $m\times n$. It follows that ($\ast$) contains at most $(m^2+n^2)$ variables and $mn$ equations and ($\ast \ast$) contains at most $m^2+nk$ variables and $mk$ equations. Let denote the exponent for matrix multiplication.

When processing a new set of columns $N$ of degreee ${\alpha }_i$, we first try to zero out each $N_b$ only with column-operations (i.e., ignoring row-operations first), which is often sufficient and, if not, will maybe reduce the number of non-zero entries. To this end we reduce the sub-matrix $M_b^{\le {\alpha }_i}$ of $M_b$ containing the columns of degree $\le {\alpha }_i$. Importantly, this matrix can be stored and sometimes even updated when the next set of columns is considered in the same way as in the standard persistence algorithm for one parameter.

The following change is subtle but crucial:

1. 1.

   Find a basis ${(Q_i,-P_i)}_{i\in I}$ of the solution-space of ($\ast$) for every $c$ and store it.

2. 2.

   Search for the solution of ($\ast \ast$) in the vector spaces spanned by ${(Q_i,-P_i)}_{i\in I}$ for each $c$.

This split makes no computational difference asymptotically, and its overhead is negligible. In subsequent algorithms, we will call BlockReduce recursively without changing $M_{ℬ}$ (thanks to Lemma 13). By solving and storing ($\ast$), we thus avoid recomputations during the recursion.

Let $M_c\in {𝕂}^{G_1\times R_1},M_b\in {𝕂}^{G_2\times R_2}$ present persistence modules $X_c,X_b$. Up to sign, the solutions of BlockReduce ($\ast$) form the $𝕂$-vector space

and there is a natural $𝕂$-linear surjection $\mathrm{Hom}(M_c,M_b)\twoheadrightarrow \mathrm{Hom}(X_c,X_b)$.

So when solving BlockReduce ($\ast$) for each pair $c\ne b\in ℬ$, we are actually computing, akin to elementary matrices, the generators of a subgroup of $\mathrm{Aut}\left(M_{ℬ}\right)$. Additionally, the property of a morphism $(Q,P)\in \mathrm{Hom}(M_c,M_b)$ to be part of a solution of BlockReduce ($\ast \ast$) depends only on the map $\tilde{Q}:X_c\to X_b$ which it induces.

Let $X,Y$ be persistence modules, then $\mathrm{Hom}(X,Y)$ can be computed faster if $Y$ is pointwise “small”. This is because we can restrict the target space of each generator of $X$. We will explore this thoroughly in future research and also point to the following known result for “smallest” modules, also due to Dey and Xin.

Furthermore, not all homomorphisms are needed to find a decomposition.

Thus, instead of finding a basis for the whole solution space of BlockReduce ($\ast$) we actually only need to compute a set of solutions which span ${\mathrm{Hom}}^{\alpha }(X_c,X_b)$. This is a simple reduction based on 16 and will decrease the size of the system ($\ast \ast$), but more importantly for Section 6 it will sometimes reduce this solution space to zero.

For $X,Y$ interval modules we can describe this reduction precisely. By [28, Proposition 16] the vector space $\mathrm{Hom}(X,Y)$ has a basis consisting of “valid” overlaps:

At most one of the valid overlaps can contain $\alpha$ and so ${dim}_K{\mathrm{Hom}}^{\alpha }(X,Y)\le 1$. Define a relation $X{\to }_{\alpha }Y:\iff {\mathrm{Hom}}^{\alpha }(X,Y)\ne 0$. This relation is a partial order on any set of isomorphism classes of interval modules ([26, Lemma 6.16]). Additionally, an interval module can have at most dimension $1$ at $\alpha$. Therefore, when we consider a set of new relations $N$ in the algorithm, if $M_b$ presents an interval module, we can reduce every sub-column $N_{b,i}$ to have a single non-zero entry. Together with the transitivity of $\alpha$-maps, this observation allows Gaussian (row-)elimination on $N$ avoiding BlockReduce ($\ast \ast$) ( 6, [26, Algorithm B.1]).

One immediate difference to the distinctly graded case is that a decomposition might introduce several new blocks. Instead of finding all of them at once, we will try to split the presentation into two summands and use recursion. We invoke 13 again to find invertible matrices, which leave $M_{ℬ}$ unchanged and such that

It follows that there is an invertible matrix $T\in {𝕂}^{k\times k}$ and a partition $[k]=\{1,\mathrm{\dots },j\}\cup \{j+1,\mathrm{\dots },k\}$ of the columns of $N$, which allow us to decompose the whole presentation by zeroeing out submatrices with BlockReduce: After multiplying $N$ with $T$ from the right, for every block $b\in ℬ$ we can delete either $N_{b,[j]}$ or $N_{b,{[j]}^{\mathrm{\complement }}}$. Here we see that we would have to solve ($\ast$) again, but since we stored its solution this is not necessary and ($\ast \ast$) is efficiently solved. Since we work over a finite field, this immediately implies a brute-force algorithm:

For each invertible matrix $T$ over ${𝔽}_q$ and $j\in [k]$ set $N_1={\left(NT\right)}_{\{1,\mathrm{\dots },j\}}$, $N_2={\left(NT\right)}_{\{j+1,\mathrm{\dots },k\}}$ and decompose each separately with BlockReduce.

This is computationally unfeasible even for smaller $k$: The maximal number of iterations of this loop for $q=2$ and $k\in \{2,\mathrm{\dots },6\}$ are $6,\mathrm{\hspace{0.17em}342},\mathrm{\hspace{0.17em}60\hspace{0.17em}822},\mathrm{\hspace{0.17em}40\hspace{0.17em}058\hspace{0.17em}262},\mathrm{\hspace{0.17em}100\hspace{0.17em}833\hspace{0.17em}607\hspace{0.17em}062}$ and more generally $|{\mathrm{𝐆𝐋}}_k({𝔽}_q)|\in q^{k^2+𝒪\left(k\right)}$. Fortunately, we can be much smarter about the matrices and partitions we use: Assume that $T$ is as above so that $NT$ decomposes into two parts, supported on $[j]$ and ${[j]}^{\mathrm{\complement }}$, after using BlockReduce ($\ast \ast$). The two parts also induce a partition of $T=\left[T_1{T}_2\right]$. If $T_1^{\prime },T_2^{\prime }$ are matrices with the same column-spaces as $T_1,T_2$ respectively then $N\left[T_1^{\prime }{T}_2^{\prime }\right]$ could be brought to the same decomposition with BlockReduce. That is, finding a decomposition with this strategy depends only on the decomposition of ${𝔽}_q^k$ as a vector space which $T=\left[T_1{T}_2\right]$ represents. Additionally, once we have deleted as many parts of $N{T}_1$ as possible, we can update the block structure and do the same again for $N{T}_2$, but now also being able to use column-operations from the columns of $N{T}_1$ to those of $N{T}_2$ without breaking linearity of the system. That is, we also do not need to consider a decomposition $\left[T_1{T}_2\right]$, if a decomposition has been considered already that is of the form $\left[T_1{T}_2^{\prime }+T_{1}S\right]$, for any matrix $S$ and any $T_2^{\prime }$ with the same column-space as $T_2$. But this is the case for every other invertible matrix of the form $\left[T_1{T}_2^{\prime }\right]$, so that we only need to iterate over all proper subspaces $\text{span}\left(T_1\right)$ of ${𝔽}_q^k$, i.e. the Grassmanian ${\mathrm{𝐆𝐫}}_l({𝔽}_q^k)$ for every . These matrices can be generated fast by leveraging the Schubert cell decomposition [56]:

For each set $P\in \left(\genfrac{}{}{0pt}{}{[k]}{l}\right)$ construct a $k\times l$ matrix $T_1$ in column-echelon form with pivots given by $P$. Then construct all subspaces in the corresponding cell by changing the entries that are irrelevant for the echelon form and not in the same row as a pivot. For each $T_1$ generated this way, choose any $T_2$ with opposite pivots and store the pairs $\left[T_1{T}_2\right]$ in a set ${\mathrm{𝐃𝐞𝐜}}_q(k,l)$.

Since exchanging $T_1$ and $T_2$ in ${\mathrm{𝐃𝐞𝐜}}_q(k,l)$ would not make a difference for the decomposition we can actually ignore the sets ${\mathrm{𝐃𝐞𝐜}}_q(k,l)$ for $l>k/2$. In the case where $k$ is even and $l=k/2$ we can avoid double counting even more. Consider the graph on ${\mathrm{𝐆𝐫}}_{k/2}({𝔽}_q^k)$ where $U\sim V$ if they form a decomposition (also called the $q$-Kneser graph). By symmetry we only need to find a minimal vertex cover of this graph and use these subspaces for the computation above. This can be done by considering in Generate ${\mathrm{𝐃𝐞𝐜}}_q(k,k/2)$ only those pivot sets $P\in \left(\genfrac{}{}{0pt}{}{[k]}{k/2}\right)$ which contain the first position. Let ${\mathrm{𝐃𝐞𝐜}}_q(k)≔\bigcup _{l\le k/2}{\mathrm{𝐃𝐞𝐜}}_q(k,l)$, then we can iterate over ${\mathrm{𝐃𝐞𝐜}}_q(k)$ instead of $[k]\times {\mathrm{𝐆𝐋}}_k\left({𝔽}_q\right)$ for the exhaustive decomposition. The first terms of $|{\mathrm{𝐃𝐞𝐜}}_2(k)|$ are $1,\mathrm{\hspace{0.17em}2},\mathrm{\hspace{0.17em}7},\mathrm{\hspace{0.17em}43},\mathrm{\hspace{0.17em}186},\mathrm{\hspace{0.17em}1965},\mathrm{\hspace{0.17em}14605},297181,\mathrm{\dots }$ and the asymptotic growth is in $q^{k^2/4+𝒪\left(k\right)}$.

At last, when calling the recursion on a sub-matrix we can pull back decompositions we have already tried to the corresponding sub-space to avoid double counting again.

For at most $n$ times we need to compute BlockReduce ($\ast$) and then solve BlockReduce ($\ast \ast$) at most $q^{k^2/4+𝒪(k)}$ many times. Using 14 gives the term above. $◀$

We have seen in 19, that if $ℬ$ is acyclic with respect to ${\mathrm{Hom}}^{\alpha }\ne 0$, it is possible to process $N$ in a way so that all admissible operations we need to consider only flow in the direction of the graph. Consequently there is no iteration over subspaces needed. In particular the algorithm would then also work for any field $𝕂$, not only a finite one. Let us formalise the structure of the matrix we need to decompose:

When $ℬ$ is acyclic, we have seen in 19 that $𝒞$ and $𝒟$ contain just a single block and there are no ($\alpha$-)maps from $M_{𝒞}$ to $M_{𝒟}$. We argued that in this case there are also no column operations from $N_{𝒟}$ to the left which need to be considered and so we can use BlockReduce ($\ast \ast$) (on the extended presentations) to decide if $C$ can be zeroed out.

The correctness of the algorithm follows from our main theorem by lifting isomorphisms:

[26, Theorem 7.25]. $◀$

Notice that the condition means that $X$ is almost invariant under the action of $\mathrm{Aut}(Z)$: No automorphism could alter $X_{\alpha }$ and therefore also not the cocycle $C:A{(-\alpha )}^{\kappa }\to X$.

In general there is no order in which we can traverse the blocks, such that maps to the current block come only from parts of the matrix which already form an indecomposable decomposition. The solution here is to contract all strongly connected components of $ℬ$ forming the condensation $𝒯$ of this digraph.

The $\alpha$-decomposed extension we need to decompose in each step may then have many blocks in $𝒞$ and $𝒟$ and thus cannot be reduced by a linear system anymore. We employ another algorithm ExtensionDecomp ([26, Algorithm B.4]) which again loops over subspaces, but of ${𝕂}^{\kappa }$, the column-space of $C$. The resulting algorithm is FPT in ${\kappa }_{\max }$, the largest occurrence of $\kappa$, and called the Automorphism-invariant Iterative Decomposition Algorithm.

Assume that the input module is interval-decomposable. Then in Section 4 we have seen that $ℬ$ is transitive and each $N_{b,i}$ can be assumed to be a single entry. Also, for every $\alpha$-decomposed extension we need to reduce, the matrices $\left[M_{𝒞}{N}_{𝒞}\right]$ and $\left[M_{𝒟}{N}_{𝒟}\right]$ must each present a set of isomorphic intervals. Then all row and column operations on $C$ are possible and we can use normal matrix reduction ([26, Algorithm B.1]).

We considered a large collection of data files [41] including chain complex data generated by various bifiltrations, including function-Rips, lower-star bifiltrations, and multi-covers – see [33]. Also included are chain complexes of function-alpha bifiltrations for various choices of point clouds and density functions as described in the paper of Alonso et al. [4]. We additionally generated test cases which are not based on bifiltration data by simple random processes that can create interval-decomposable or random presentation inspired by the Erdös-Renyi-model of random graphs ([26, Appendix A]). In total, the benchmark set contains around 1300 different minimal presentation files. Our benchmarks scripts, the version of aida that we have used, and the results of our experiments are available at [40].

Most instances in the above collection are distinctly graded, that is $k_{\max }=1$. This is perhaps surprising because, for instance, function-alpha bifiltrations are generally not distinctly graded simplex-wise. Yet apparently their homology classes have unique degrees. A few instances also have $k=2$ or $k=3$, but the only instances with large values of $k$ are the one coming from multi-cover bilfiltrations as described in [25].

We implemented our code such that some optimizations can be disabled to gauge their importance. Our first finding is that the column sweep optimization as described at the beginning of Section 3 has the most visible effect on practical performance. In Table 1, we display the running time for some representative instances. We observe that the improvement factor depends on the instance type and even varies a lot for instances of the same type (as the third and fourth row of the table indicate). However, the improvement was consistent over all considered examples.

We next consider the effect of the optimized computation of the system ($\ast$). In many instances, this technique does not speed up the computation or even slightly slows it down, due to the technical overhead of solving ($\ast$) first. However, the slow-down was never more than a factor of $2$ in the experiments. On the other hand, for some instances, a speed-up is achieved; most significantly, this is the case for the off-datasets, which are triangular meshes in 3D bi-filtered along two of the coordinate axis (bottom half of Table 1). This happens because this filtration tends to create the large but pointwise low-dimensional summands which we have optimized the computation for.

We demonstrate in Table 2 that aida is able to decompose presentations with hundreds of thousands of generators and relations within seconds in many practical cases. This is also true for non-distinctly graded cases: the bottom half of the table displays multi-cover instances where $k_{\max }$ is larger than $1$ ($k_{\max }=4$ on the left, $k_{\max }=6$ on the right of the table). We also display the progression of aida’s time and memory consumption when the input size (roughly) doubles. We can see that in some cases, aida scales almost linearly with the size of the input presentation (the upper-left and lower-right part of the table) whereas in other cases the time grows rather quadratically (the upper-right and lower-left part). We also observed instances with a super-quadratic growth in our experiments; nevertheless, the empirical behavior is much better than what the asymptotic bound predicts. This is not surprising for realistic sparse input matrices, in analogy to the persistence algorithm for one parameter, because matrix reduction is generally faster. Also the summands typically stay small, so that ($\ast$) can be solved fast.

We are not aware of any other implementation for the special case of multi-graded persistence modules, but there exist some libraries that include a decomposition algorithm for modules over finite dimensional algebras, for example the algorithm by Lux and Szöke [50] in the C-MeatAxe library. Since using this algorithm would require too many preprocessing steps, we instead used the decomposition algorithm for quiver representations in the QPA pacakge [54] for GAP [34] and implemented (also in the aida package) an algorithm to convert MPM to quiver representations. Our tests show that the QPA algorithm already needs several minutes for presentations with around $100$ generators and relations and times out for all presentations we have considered for testing aida.

We demonstrate how decomposing a module can speed up computational tasks in multi-parameter persistence: barcode templates are a common tool to capture all combinatorial barcodes that arise from a $2$-parameter module through restrictions to one parameter. Such a barcode template is computed via the arrangement of lines in the plane, arising from the joins of generators and relations of a module via point-line duality. It is evident that this construction is compatible with decomposition, meaning that if a module $M=M_1\oplus M_2$ decomposes, the arrangements of $M_1$ and $M_2$ can be computed separately and then be overlaid.

We implemented a small test program that computes the barcode template using the arrangement package of Cgal [59] once for the undecomposed modules and once after applying aida. Table 3 shows that the arrangement size as well as the computation time drastically reduces when using aida. We remark that the variance of the aida arrangement in size was very large, leading to the effect that the average arrangment size even decreases in our sample. In either case, we speculate that incorporating aida would speed-up the visualization of persistence modules in Rivet and is a necessary preprocessing step in the exact computation of matching distances in practice.