Clear and Compress: Computing Persistent Homology in Chunks

We present a parallelizable algorithm for computing the persistent homology of a filtered chain complex. Our approach differs from the commonly used reduction algorithm by first computing persistence pairs within local chunks, then simplifying the unpaired columns, and finally applying standard reduction on the simplified matrix. The approach generalizes a technique by G\"unther et al., which uses discrete Morse Theory to compute persistence; we derive the same worst-case complexity bound in a more general context. The algorithm employs several practical optimization techniques which are of independent interest. Our sequential implementation of the algorithm is competitive with state-of-the-art methods, and we improve the performance through parallelized computation.


Introduction
Persistent homology has developed from a theoretical idea to an entire research area within the field of computational topology. One of its core features is its multi-scale approach to analyzing and quantifying topological features in data. Recent examples of application areas are shape classification [2], topological denoising [1], or developmental biology [7].
A second major feature of persistent homology is the existence of a simple yet efficient computation method: The standard reduction algorithm as described in [8,17] computes the persistence pairs by a simple sequence of column operations; Algorithm 1 gives a complete description in just 10 lines of pseudo-code. The worst-case complexity is cubic in the input size of the complex, but the practical behavior has been observed to be closer to linear on average. Various variants have been proposed in order to improve the theoretical bounds [14,3] or the practical behavior [5].
Our first contribution consists of two simple optimization techniques of the standard reduction algorithm, which we call clearing and compression. Both approaches exploit the special structure of a filtered chain complex in order to significantly reduce the number of operations on real-world instances. However, the two methods cannot be easily combined because they require the columns of the boundary matrix to be processed in different orders.
Our second contribution is a novel algorithm that incorporates both of the above optimization techniques, and is also suitable for parallelization. It proceeds in three steps: In the first step, the columns of the matrix are partitioned into consecutive chunks. Each chunk is reduced independently, applying the clearing optimization mentioned above. In this step, the algorithm finds at least all persistence pairs with (index) persistence less than the size of the smallest chunk; let g be the number of columns not paired within this first step. In the second step, the g unpaired columns are compressed using the method mentioned above. After compression, each column has at most g non-zero entries and the unpaired columns form a nested (g × g)-matrix. In the third and final step, this nested matrix is reduced, again applying the clearing optimization.
The chunk algorithm is closely related to two other methods for computing persistence. First of all, the spectral sequence algorithm [6,§VII.4] decomposes the matrix into blocks and proceeds in several phases, computing in phase r the persistence pairs lying r − 1 blocks apart. The first step of the chunk algorithm is equivalent to applying the first two phases of the spectral sequence approach. Furthermore, the three step chunk algorithm is inspired by the approach of Günther et al. [10], which combines persistence computation and discrete Morse theory for 3D image data. The first step of that algorithm consists in constructing a discrete gradient field consistent with the input function; such a gradient field can be interpreted as a set of persistence pairs that are incident in the complex and have persistence 0. We replace this method by local persistence computations, allowing us to find pairs with are not incident in the complex.
We analyze the chunk algorithm in terms of time complexity. Let n be the number of generators (simplices, cells) of the chain complex, m the number of chunks, the maximal size of a chunk, and g as above. We obtain a worst-case bound of where the three terms reflect the worst-case running times of the three steps. For the filtration of a cubical complex induced by a d-dimensional grayscale image (with d some fixed constant), this bound simplifies to if the chunks are given by the cells appearing simultaneously in the filtration. This bound improves on the previous general bound of O(g 2 n log n) from [4]. Moreover, it matches the bound in [10], but applies to arbitrary dimensions. Of course, the bound is still cubic if expressed only in terms of n because g ∈ O(n) in the worst case. We implemented a sequential and a parallelized version of the chunk algorithm; both are publicly available in our new PHAT library (http://phat.googlecode.com/). The sequential code already outperforms the standard reduction algorithm and is competitive to other variants with a good practical behavior. The parallelized version using 12 cores yields a speed-up factor between 3 and 11 (depending on the example) in our tests, making the implementation the fastest among the considered choices. This is the first result where the usefulness of parallelization is shown for the problem of persistence computation through practical experiments.

Background
This section summarizes the theoretical foundations of persistent homology as needed in this work. We limit our scope to simplicial homology over Z 2 just for the sake of simplicity in the description; our methods generalize to chain complexes over arbitrary fields.
Homology Homology is an algebraic tool for analyzing the connectivity of topological spaces. Let K be a simplicial complex of dimension d. In any dimension p, we call a pchain a formal sum of the p-simplices of K with Z 2 coefficients. The p-chains form a group called the pth chain group C p . The boundary of a p-simplex σ is the (p − 1)-chain formed by the sum of all faces of σ of codimension 1. This operation extends linearly to a boundary operator δ : C p → C p−1 . A p-chain γ is a p-cycle if δ (γ) = 0. The p-cycles form a subgroup of the p-chains, which we call the pth cycle group Z p . A p-chain γ is called a p-boundary if γ = δ (ξ ) for some (p + 1)-chain ξ . Again, the p-boundaries form a group B p , and since δ (δ (ξ )) = 0 for any chain ξ , p-boundaries are p-cycles, and so B p is a subgroup of Z p . The pth homology group H p is defined as the quotient group Z p /B p . The rank of H p is denoted by β p and is called the pth Betti number. In our case of Z 2 coefficients, the homology group is a vector space isomorphic to Z β p 2 , hence it is completely determined by the Betti number. Roughly speaking, the Betti numbers in dimension 0, 1, and 2 yield the number of connected components, tunnels, and voids of K, respectively.
Persistence Let {σ 1 , . . . , σ n } denote the simplices of K. We assume that for each i ≤ n, For every dimension p and every K i , we have a homology group H p (K i ); we usually consider all dimensions at once and write H(K i ) for the direct sum of the homology groups of K i in all dimensions. The inclusion K i → K i+1 induces a homomorphism g i+1 i : H(K i ) → H(K i+1 ) on the homology groups. These homomorphisms compose and we can define g j i : H(K i ) → H(K j ) for any i ≤ j. We say that a class α ∈ H(K ) is born at (index) i if α ∈ im g i but α / ∈ im g i−1 . A class α born at index i dies entering (index) j if g j i (α) ∈ im g j i−1 but g j−1 i (α) / ∈ im g j−1 i−1 . In this case, the index pair (i, j) is called a persistence pair, and the difference j − i is the (index) persistence of the pair. The transition from K i−1 to K i either causes the birth or the death of an homology class. We call the added simplex σ i positive if it causes a birth and negative if it causes a death. Note that homology classes of the full complex K do not die during the filtration. We call a simplex σ i that gives birth to such a class essential. All other simplices are called inessential.
Boundary matrix For a matrix M ∈ Z n×n 2 , we let M j denote its j-th column, M i its ith row, and M i j ∈ Z 2 its entry in row i and column j. For a non-zero column 0 = M j = (m 1 , . . . , m n ) ∈ Z n 2 , we set pivot(M j ) := max{i = 1, . . . , n | m i = 1} and call it the pivot index of that column.
The boundary matrix D ∈ (Z 2 ) n×n of a simplexwise filtration (K i ) i is a n × n matrix with D i j = 1 if and only if σ i is a face of σ j of codimension 1. In other words, the jth column of D encodes the boundary of σ j . D is an upper-triangular matrix because any face of σ j must precede σ j in the filtration. Since the jth row and column of D corresponds to the jth simplex σ j of the filtration, we can talk about positive columns, negative columns, and essential columns in a natural way, and similarly for rows.
The reduction algorithm A column operation of the form M j ← M j + M k is called leftto-right if k < j. We call a matrix M derived from M if M can be transformed into M by left-to-right operations. Note that in a derivation M of M, the jth column can be expressed as a linear combination of the columns 1, . . . , j of M, and this linear combination includes M j . We call a matrix R reduced if no two non-zero columns have the same pivot index. If R is derived from M, we call it a reduction of M. In this case, we define Although the reduction matrix R is not unique, the sets P R and E R are the same for any choice of reduction; therefore, we can define P M and E M to be equal to P R and E R for any reduction R of M. We call the set P the persistence pairs of M. When obvious from the context, we omit the subscripts and simply write P for the persistence pairs. For the boundary matrix D of K, the pairs (i, j) ∈ P are the persistence pairs of the filtration (K i ) 0≤i≤n , and the indices in E correspond to the essential simplices of the complex. Note that E is uniquely determined by P and n as the indices between 1 and n that do not appear in any pair of P.
The simplest way of reducing D is to process columns from left to right; for every column, other columns are added from the left until the pivot index is unique (Algorithm 1). A lookup table can be used to identify the next column to be added in constant time. A flag is used for every column denoting whether a persistence pair with the column index has already been found. After termination, the unpaired columns correspond to the essential columns. The running time is at most cubic in n, and this bound is actually tight for certain input filtrations, as demonstrated in [16].
Mark columns i and j as paired and add (i, j) to P 10: With this definition, a matrix M is a reduction of D if and only if every column is reduced.

Speed-ups
Algorithm 1 describes the simplest way of reducing the boundary matrix, but it performs more operations than actually necessary to compute the persistence pairs. We now present two simple techniques which both lead to a significant decrease in the number of required operations.
Clearing positive columns The key insight behind our first optimization is the following fact: if i appears as the pivot in a reduced column of M, the index i is positive and hence there exists a sequence of left-to-right operations on M i that turn it to zero. Instead of explicitly executing this sequence of operations, we define the clear operation by setting column M i to zero directly. Informally speaking, a clear is a shortcut to avoid some column operations in the reduction when the result is evident.
In order to apply this optimization, we change the traversal order in the reduction by first reducing the columns corresponding to simplices with dimension d (from left to right), then all columns with dimension d − 1, and so on. After having reduced all columns with dimension δ , we have found all positive inessential columns with dimension δ − 1 and clear them before continuing with δ − 1. This way all positive inessential columns of the complex are cleared without performing any column additions on them. See [4] for a more detailed description.
Compression Alternatively, we can try to save arithmetic operations by reducing the number of non-zero rows among the unpaired columns. A useful observation in this context is given next. Proof. The statement is clearly true if M j is reduced, because in this case (i, j) is a persistence pair. If M j is not reduced, this means that after applying some sequence of left-to-right column operations, some reduced column has i as pivot index. As a consequence, whenever a negative column with index j has been reduced, row j can be set to zero before further reducing. We can even do more: let i be the pivot index of the reduced column M j and assume that the submatrix of M with column indices {1, . . . , j} and row indices {i, . . . , n} is reduced, i.e., the pivot indices are unique in this submatrix. By adding column j to each unreduced column in the matrix that has a non-zero entry at row i, we can eliminate all non-zero entries in row i from the unreduced columns. Note that if k < j and M i k = 0, then pivot(M k ) ≥ i and thus, by assumption, M k must be a reduced negative column. Therefore, for each unreduced column M k , the operation M k ← M k + M j is a left-to-right addition and thus does not affect the pairs.

Reduction in chunks
The two optimization techniques from Section 3 both yield significant speed-ups, but they are not easily combinable, because clearing requires to process a simplex before its faces, whereas compression works in the opposite direction. In this section, we present an algorithm which combines both optimization techniques.
Let m ∈ N. Fix m + 1 numbers 0 = t 0 < t 1 < . . . < t m−1 < t m = n and define the ith chunk of D to be the columns of D with indices {t i−1 + 1, . . . ,t i }. We call a column D j local if it forms a persistence pair with another column in the same chunk or in one of the adjacent chunks. In this case, we also call the persistence pair local. Non-local columns (and pairs) are called global. If is a lower bound on the size of each chuck, then every global persistence pair has index persistence at least . We also call an index j local if the jth column of D is local, and the same for global. We denote the number of global columns in D by g. The high-level description of our new algorithm consists of three steps: 1. Partially reduce every chunk independently, applying the clearing optimization, so that all local columns are completely reduced.
2. Independently compress every global column such that all its non-zero entries are global.
3. Reduce the submatrix consisting only of the global rows and columns.
We give details about the three steps in the rest of this section. The first two steps can be performed in parallel, whereas the third step only needs to reduce a matrix of size g × g instead of n × n. In many situations, g is significantly smaller than n.
Local chunk reduction The first step of our algorithm computes the local pairs by performing two phases of the spectral sequence algorithm [6,§VII.4]. Concretely, we apply left-to-right operations as usual, but in the first phase we only add columns from the same chunk, and in the second phase we only add columns from both the same chunk and its left neighbor. After phase r, for each b ∈ {r, . . . , m} the submatrix with column indices {1, . . . ,t b } and row indices {t b−2 +1, . . . , n} is reduced. If the reduction of column j stops at a pivot index i > t b−r , row j cannot be reduced any further by adding any column, so we identify (i, j) as a local persistence pair. Conversely, any local pair (i, j) is detected by this method after two phases. We incorporate the clearing operation for efficiency, that is, we proceed in decreasing dimension and set detected local positive columns to zero; see Algorithm 2. After its execution, L[i] contains the index of the local negative column with pivot index i for any local positive column i, and the resulting matrix R is a derivation of D in which all local columns are reduced. for r = 1, 2 do Perform two phases of the spectral sequence algorithm 5: for b = r, . . . , m do Loop is parallelizable 6: for if j is not marked as paired then 8: 10: if R j = 0 then 11: i ← pivot(R j ) Mark i and j as paired and add (i, j) to P 16: return (R, L, P) if k is marked as active/inactive then return true/false 6: for each non-zero entry index i of R k do 7: if is unpaired then 8: mark k as active and return true 9: else if i is positive then 10: (i, j) is persistence pair 11: if j = k and MARK COLUMN(R, j) then 12: mark k as active and return true 13: mark k as inactive and return false Global column compression Let R be the matrix returned by Algorithm 2. Before computing the global persistence pairs, we first compress the global columns, using the ideas from Section 3; recall that negative rows can simply be set to zero, while entries in positive rows can be eliminated by an appropriate column addition. Note, however, that a full column addition might actually be unnecessary: for instance, if all non-zero row indices in the added column belong to negative columns (except for the pivot), the entry in the local positive row could just have been zeroed out in the same way as in Corollary 3. Speaking more generally, it is more efficient to avoid column additions that have no consequences for global indices, neither directly nor indirectly.
In the spirit of this observation, we call an index i inactive if either it is a local negative index or if (i, j) is a local pair and all indices of non-zero entries in column R j apart from i are inactive. Otherwise, the index is called active. By induction and Corollary 3, we can show: The compression proceeds in two steps: first, every non-zero entry of a global column is classified as active or inactive (using depth-first search; see Algorithm 3). Then, we iterate over the global columns, set all entries with inactive index to zero, and eliminate any non-zero entry with a local positive index by column addition with L[ ] (see Algorithm 4). After this process, we obtain a matrix R with the same persistence pairs as R, such that the global columns of R have non-zero entries only in the global rows.
Submatrix reduction After having compressed all global columns, these form a g × g matrix "nested" in R (recall that g is the number of global columns). To complete the computation of the persistence pairs, we simply perform standard reduction on the remaining matrix. For efficiency, we perform steps 2 and 3 alternatingly for all dimensions in decreasing order and apply the clearing optimization; this way, we avoid the compression of positive global columns. Algorithm 5 summarizes the whole method. for each non-zero entry index of R k in decreasing order do 4: if is paired then 5: if is inactive then for j = 1, . . . , n with dim σ j = δ do Loop is parallelizable 7: if column j is not paired then for j = 1, . . . , n with dim σ j = δ do step 3: reduce global columns 10:

12:
if R j = 0 then Mark i and j as paired and add (i, j) to P 17: return P

Analysis
Algorithm 5 permits a complexity analysis depending on the following parameters: n, the number of simplices; m, the number of chunks; , the maximal size of a chunk; and g, the number of global columns. We assume that for any simplex, the number of non-zero faces of codimension 1 is bounded by a constant (this is equivalent to assuming that the dimension of the complex is a constant).

General complexity We show that the complexity of Algorithm 5 is bounded by
The three summands correspond to the running times of the three steps 1 . Note that g ∈ O(n) in the worst case. For the complexity of Algorithm 2, we consider the complexity of reducing one chunk, which consists of up to columns. Within the local chunk reduction, every column is only added with columns of the same or the previous chunk, so there are only up to 2 column additions per column. Moreover, since the number of non-zero entries per column in D is assumed to be constant, there are only O( ) many entries that can possibly become non-zero during the local chunk reduction. It follows that the local chunk reduction can be considered as a reduction on a matrix with columns and O( ) rows. If we represent columns by linked lists (containing the non-zero indices in sorted order), one column operation can be done in O( ) primitive operations, which leads to a total complexity of O( 3 ) per chunk.
The computation of active columns in Algorithm 3 is done by depth-first search on a graph whose vertices are given by the columns and whose edges correspond to their non-zero entries. The number of edges is O(n ), so we obtain a running time of O(n ).
Next, we consider the cost of compressing a global column with index j. After the previous step, the column has at most O( ) non-zero entries. We transform the presentation of the column from a linked list into a bit vector of size n. In this representation, adding another column in list representation with v entries to column j takes time proportional to v. In the worst case, we need to add all columns with indices 1, . . . , j − 1 to j. Each such column has O( ) entries. At the end, we transform the bit vector back into a linked list representation. The total cost is O(n + ( j − 1) + n) = O(n ) per global column.
Finally, the complexity of the global reduction is O(g 3 ), as in the standard reduction.
Choosing chunks We discuss different choices of chunk size and their complexities. A generic choice for an arbitrary complex is to choose O( √ n) chunks of size O( √ n) each. With that, the complexity of (1) becomes O(n 2 + g 1 n √ n + g 3 1 ).
Alternatively, choosing O( n log n ) chunks of size O(log n), the complexity becomes O(n log 2 n + g 2 n log n + g 3 2 ).
We replaced g by g 1 and g 2 to express that the number of global columns is different in both variants. In general, choosing larger chunks is likely to produce less global columns, since every global persistence pair has index persistence at least (the size of the smallest chunk) .
Cubical complexes We consider an important special case of boundary matrices: consider a d-dimensional image with p hypercubes, where each vertex contains a grayscale value. We assume that the cubes are triangulated conformally in order to get simplicial input -the argument also works, however, for the case of cubical cells. We assign function values inductively, assigning to each simplex the maximal value of its faces. Assuming that all vertex values are distinct, the lower star of vertex v is the set of all simplices which have the same function value as v. Filtering the simplices in a way that respects the order of the function values, we get a lower star filtration of the image. Now choose the lower stars as the chunks in our reduction algorithm. Note that the lower star is a subset of the star of the corresponding vertex, which is of constant size (assuming that the dimension d is constant). Therefore, the complexity bound (1) reduces to O(n + gn + g 3 ) = O(gn + g 3 ).
Note that global columns with large index persistence might still have very small, or even zero, persistence with respect to the function values, for instance in the presence of a flat region in the image where many vertices have similar values.

Experiments
We implemented two versions of the algorithm presented in Section 4: a sequential and a parallel version (using OPENMP), in which the first two steps of the algorithm are performed simultaneously on each chunk and on each global column, respectively. In both cases, we use √ n as the chunk size. For a fair comparison, we also re-implemented the algorithms introduced in [4,8] in the same framework, that means, using the same data representations and low-level operations such as column additions. Our implementation is publicly available in our new PHAT library for computing persistence homology, available at http://phat.googlecode.com/. Additionally, we compare to the memory efficient algorithm [10] based on discrete Morse theory [9] and to the implementation of the persistent cohomology algorithm [5] found in the DIONYSUS library [15].
To find out how these algorithms behave in practice, we apply them to five representative data sets. The first three are 3D image data sets with a resolution of 128 3 . The first of these is given by a Fourier sum with random coefficients and is representative of smooth data. The second is uniform noise. The third is the sum of the first two and represents largescale structures with some small-scale noise. These data sets are illustrated in Figure 1 by an isosurface. a) b) c) Figure 1: A single isosurface of the representative data sets used to benchmark and compare our algorithm: a) a smooth data set, b) uniformly distributed noise, c) the sum of a) and b).
In addition to the lower star filtrations of these image data sets, we also consider an alpha shape filtration defined by 10000 samples of a torus embedded in R 3 , and the 4-skeleton of the Rips filtration given by 50 points randomly chosen from the Mumford data set [11].
As pointed out in [5], the pairs of persistent cohomology are the same as those of persistent homology. We therefore also applied all algorithms to the corresponding cochain filtration, given in matrix form by transposing the boundary matrix and reversing the order of the columns; this operation is denoted by (·) ⊥ . When reducing such a coboundary matrix with the clearing optimization, columns are processed in order of increasing instead of decreasing dimension. Table 1 contains the running times of the above algorithms applied to filtrations of these five data sets run on a PC with two Intel Xeon E5645 CPUs. We can observe a huge speed-up caused by the clearing optimization as already reported in [4]. We can see that our chunk algorithm performs slightly worse than the one of [4] when executed sequentially, but faster when parallelized. We also observe that standard, twist, and chunk generally behave worse when computing persistent cohomology, except for the Rips filtration of the Mumford data set, where the converse is true.

Dataset
n · 10 −6 std. [8] twist [4] cohom. [ Table 1: Running time comparison of various persistent homology algorithms applied to the data sets described in Section 6. The last three columns contain information of the algorithm presented in this paper: the fraction of global columns g/n, and the running times using one and twelve cores, respectively.

Conclusion and Outlook
We have presented an algorithm for persistent homology that includes two simple optimization techniques into the reduction algorithm. It can be fully parallelized, except for the reduction of compressed global columns, whose number is often small. Besides our asymptotic complexity bounds, which give a detailed dependence on the parameters of the algorithm, our experiments show that significant speed-ups can be achieved through parallelized persistence computation. Similar observations have been made recently by Lewis and Zomorodian [12] for the computation of (non-persistent) homology; see also [13]. We plan a more extensive discussion of the practical effects of our optimizations and parallelization in an extended version of this paper.