Persistent (Co)Homology in Matrix Multiplication Time

For completeness, we recount the standard persistence algorithm, see [15], [27], and [16] for a more complete description. Let $𝐃$ denote the boundary operator, which we represent as an $n\times n$-matrix over some fixed (finite) field $𝔽$. The rows (top-down) and columns (left-right) are ordered by appearance in the filtration. Hence, the $i$-th row and column represents the simplex $\sigma =K_i-K_{i-1}$ and $𝐃[\cdot ,i]$ represents the boundary of $\sigma$ respectively. The goal is to compute the $\mathrm{𝐑𝐕}$-decomposition [8]. That is, find matrices $𝐑$ and $𝐕$ such that

where $𝐑$ is a matrix in reduced column echelon form and $𝐕$ is a full-rank upper-triangular matrix; we let $𝐔={𝐕}^{-1}$. We introduce the helper function

which is defined whenever $𝐑[\cdot ,j]\ne \mathrm{𝟎}$. When defined, it returns the largest row index $i$ where the column is non-zero; we refer to these lowest non-zero elements as pivots. Throughout this paper, $\mathrm{𝐥𝐨𝐰}$ always applies to $𝐑$. The persistence diagram is then given by:

The standard algorithm is a variant of Gaussian elimination. Algorithm 1, which we call the lazy reduction, reduces each column by considering pivots in the columns to the left. In each step of the loop, it removes the lowest non-zero entry until a new pivot is found or the entire column is zeroed out. A straightforward analysis gives a running time bound of $O(n^3)$. In [22], a filtration where this algorithm takes cubic time was given, showing that the bound is tight. The updates of matrix $𝐕$ follow those of matrix $𝐑$. The updates in matrix $𝐔$ undo the updates in $𝐕$ to maintain $𝐔\cdot 𝐕=𝐈$. Because the columns of $𝐑$ are processed from left to right, during the update of matrix $𝐔$, the row $𝐔[j,\cdot ]$ has a single (diagonal) element. Therefore, the update in $𝐔$ is equivalent to just setting the entry $𝐔[j^{\prime },j]$ to $\alpha$.

We introduce another variant of persistence computation used in applications. This one “looks ahead” and eliminates as many elements from the matrix as it can; so we call it the exhaustive reduction. In Algorithm 2, by the time a column is considered, it is already reduced. It is then applied to all columns to the right, zeroing out the entire row. In this way, when a column is processed, all previous pivots have already been applied. Because column $𝐑[\cdot ,j^{\prime }]$ may have already been used to reduce another column $𝐑[\cdot ,j^{\prime \prime }]$, when an update is applied to it, the row $𝐔[j^{\prime },\cdot ]$ need not consist of a single diagonal element. Hence, the row update in matrix $𝐔$ is not equivalent to just setting $𝐔[j,j^{\prime }]$ to $\alpha$, like in the lazy version of the algorithm.

Given an $𝐑=\mathrm{𝐃𝐕}$ decomposition, we distinguish between two types of simplices in the filtration: positive simplices create new homology classes; negative simplices destroy them. From the matrix decomposition, one can recover a set of cycle representatives of persistent homology. There are two sources of this information, although with different content.

1. 1.

   $𝐕$ matrix: The columns of positive simplices in $𝐕$ correspond to zero columns in the reduced matrix $𝐑$. These columns in $𝐕$ are cycles, by definition.

2. 2.

   $𝐑$ matrix: The columns of negative simplices in $𝐑$ are non-zero. They store linear combinations of boundaries, which are cycles (at the time of their birth) that eventually die in the filtration.

The cycles one recovers from matrix $𝐑$ differ between the lazy and exhaustive reductions. Exhaustive reduction was called total reduction by Cohen-Steiner et al. [9]. They show that these cycles form a lexicographically optimal basis, which can be used to triangulate point cloud data. Nigmetov and Morozov [24] use the lazy reduction – specifically, the columns and rows from matrices $𝐕$ and $𝐔$ – to recover “critical sets,” subsets of simplicial complexes affected by optimization with the loss formulated in terms of a persistence diagram.

We present a fast version of exhaustive reduction in Algorithm 3. The general structure follows Algorithm 2. However, rather than apply the pivots when we find them, we apply them in batches via the Schur complement through a standard recursion. We first augment the boundary matrix with an identity matrix, which is of size at least $n$, for reasons we shall describe below. To simplify notation, we assume $n$ is a power of $2$. If it is not, we can set the size of the identity matrix to $m+n$ with such that $m+2n=2^x$. Therefore, our initial matrices are:

Throughout the algorithm, we implicitly perform two types of permutations by taking appropriate subsets of rows and columns. The first are column permutations. To ensure that the appropriate submatrices are invertible, we require that after processing $k$ columns of $𝐑$, there are $k$ pivots, or, equivalently, that the matrix corresponding to the processed columns is full (column) rank. Because cycles reduce to zero columns, the if-statement in Algorithm 3 permutes a zero column from $𝐑$ with a non-zero column from the identity matrix ${𝐈}_{m+n}$, and records the indices of such columns in list $Z$, to eventually undo the permutation.

We must also consider row permutations. To apply Theorem 7, we must find the inverse ${𝐑}^{-1}[L,B]$ in Algorithm 3, which is used in the Schur complement. We construct a sequence $L=⟨\mathrm{𝐥𝐨𝐰}(\mathrm{\ell }):\mathrm{\ell }\in B⟩$ of the pivot rows in columns $B$. Denoting by $\overline{L}$ the sequence of rows outside of $L$, we denote a permutation that places rows $L$ below $\overline{L}$ as ${𝐏}_L$:

See Figure 2(b-c). The permutation is explicitly given by,

where $B(i)$ is the $i$-th entry of the indices in $B$.

We delay the proof until after the description of the algorithm. We now follow a recursion: we split matrix $𝐑$ into the first half of the columns denoted by $B$ and the second half denoted by $C$. We proceed to call the algorithm on the first half of the matrix. Once we reach a leaf in the recursion tree (shown in Figure 3), the matrix consists of a single column. The key invariant is that whenever we reach a leaf, the corresponding column is reduced with respect to all the columns which come prior to it in the filtration. For the first column, this is tautologically true and we delay the proof for the general case until later.

If we are not in the base case (a single column), we must apply the pivots in the columns indexed by $B$ (left half of the submatrix) to the the columns indexed by $C$ (right half of the submatrix). As described in Equation 3, we construct the permutation ${𝐏}_L$ such that

By construction of ${𝐏}_L$, $𝐑[L,B]$ is lower-triangular with non-zero diagonal entries and so is full-rank. We can then apply the Schur complement to zero out all the corresponding rows in the columns $C$, $𝐑[L,C]\leftarrow 0$ and update the rest of the rows:

To update matrix $𝐕$, we perform the same operations $𝚲=𝐑{[L,B]}^{-1}\cdot 𝐑[L,C]$ on $𝐕$:

Once the left half of the matrix has been applied to the right half, we recurse on the second half, i.e., $C$. Here again, because $𝐑[L,B]$ is lower-triangular full-rank and $𝐑[L,C]=0$, the first column of $𝐑[\cdot ,C]$ is reduced. We now prove correctness.

The fast version of Algorithm 4 is given in Algorithm 5. As above, it reduces the matrix row by row from the bottom up. This variant does not require matrix inversion and requires less padding. The trade-off is that it requires a more complex sequence of updates. Assuming $𝐃$ is an $n\times n$ matrix, we set such that $m+n=2^x$. We initialize our matrices as:

To simplify notation, from this point on, we assume $n$ is a power of 2 (as ), so the resulting matrix is . This ensures that $\mathrm{𝐥𝐟𝐭}(i)$ is defined for all $i$ in $𝐑$, i.e., the identity matrix ensures there are no zero rows.

It will be convenient to again consider permutations, but here they will be column permutations. We build up the permutation incrementally. In the $i$-th step, we update the permutation with the column transposition

which we use a shorthand for the two mappings $\mathrm{𝐥𝐟𝐭}(n-i+1)\mapsto i$ and $i\mapsto \mathrm{𝐥𝐟𝐭}(n-i+1)$. As $\mathrm{𝐥𝐟𝐭}$ is always defined, two columns are permuted in each step. The permutation after $i$ steps is then given by

Note that we apply the permutations in order, i.e., for $i=1,2,\mathrm{\dots }$ While the columns corresponding to $\mathrm{𝐥𝐟𝐭}(n-i+1)$ are fixed at $i$ for all permutations, the original column at $i$ may be permuted several times. In the algorithm, we omit the subscript; $𝐏$ refers to the accumulated permutation. As in the case of the column algorithm, we perform the permutations implicitly. We observe that in the permuted order, we arrange pivots on the left along the anti-diagonal, see Figure 4.

We keep track of the operations directly in the matrix $𝚲$ which is initialized as the identity matrix. Before delving into the details of the algorithm, assume that the all rows below the $i$-th row have been reduced. The updates in the $i$-th row are given by

We digress here to explain why this works. As mentioned above, if we consider the column-permuted matrix $𝐑\cdot 𝐏$, it has the pivots in the bottom-left along the anti-diagonal (Figure 4). The $i$-th column corresponds to the column given by $\mathrm{𝐥𝐟𝐭}(n-i+1)$ (in the original filtration order). Hence, it is the earliest column with a non-zero entry in the $i$-th row from the bottom such that there is no pivot below. For any non-zero entry to the right of $i$, it must occur to the right of $i$ in the original filtration order (before permuting), or it would have been returned by $\mathrm{𝐥𝐟𝐭}(n-i+1)$ rather than the current $i$-th column in the permuted order.

Finally, as we permute the columns we must also permute the rows of the $𝚲$ matrix so that the recorded column operations match. This represents the base case of the recursion. As $𝚲$ is initialized as the identity matrix, the case $n-i+1=j$ is taken care of implicitly.

Returning to the algorithm, we proceed by recursing as in the column algorithm, but on the rows rather than on the columns. We divide the rows of the sub-matrix into the top and bottom halves. We recurse first on the bottom half. Because the pivots are determined by rows below the current one, we can directly identify the pivot with $\mathrm{𝐥𝐟𝐭}$.

The base case is a single row. By assumption, rows below the current row have been updated. At row $n-i+1$, we first find the left-most non-zero entry and perform a column permutation setting the column $\mathrm{𝐥𝐟𝐭}(n-i+1)$ to the $i$-th column, making $𝐑[n-i,i]$ a pivot as required. By the reasoning above, we may zero out the row, recording the operations in $𝚲$ and we can update $(𝐑\cdot 𝐏)[n-i,\{i+1,\mathrm{\dots },n\}]\leftarrow 0$.

Though the recursion proceeds as in Algorithm 3, splitting the matrix into the indices in the “first half” and “second half,” a complication is that the columns are processed left to right (same as the indexing), while the rows are processed bottom-up while being indexed top-down. Hence, we introduce:

by which we mean for any $i\in C$, $i\mapsto n-i+1$. This indexes the rows from the bottom of the matrix as required by the algorithm. Returning to the recursion, once we processed the columns for $B$, we first apply the updates to the rows $C_r$,

shown in orange in Figure 5. Note that though $B_r$ is shown in Figure 5, it is not used explicitly in the algorithm as all the operations in those rows have already been performed. On the other hand, $B$ in $(𝐏\cdot 𝚲\cdot 𝐏)[B,\cdot ]$ and in $(𝐑\cdot 𝐏)[C_r,B]$ do not need to be reversed, as they both represent column operations which are in the standard order. After this, the first row in $C_r$ is now up-to-date with all operations from rows below it. This allows us to correctly identify $\mathrm{𝐥𝐟𝐭}$ for this row and perform the column transposition accordingly. Hence, when we recurse on $C$, the base case can be carried out.

Returning from $C$, we complete the recursion by performing the update operations on the remaining parts of the columns in $C$, shown in green in Figure 5. Note, the columns are indexed by $C$ not $C_r$. Let ${\overline{A}}_r=\{1,\mathrm{\dots },\min (C_r)-1\}.$ The column update is given as

Now all columns prior to and including $C$ and all rows below and including $C_r$ are up to date. We show the update sequence in Figure 5 (omitting the permutations). Note that in the column update, the rows in ${\overline{A}}_r$ are correctly indexed, so no reversal is necessary.