An algorithm for calculating top-dimensional bounding chains

Introduction

Topological spaces are by and large characterized by the cycles in them (i.e., closed paths and their higher dimensional analogues) and the ways in which they can or cannot be deformed into each other. This idea had been recognized by Poincaré from the beginning of the study of topology. Consequently much of the study of topological spaces has been dedicated to understanding cycles, and these are the key features studied by the topological data analysis community (Edelsbrunner & Harer, 2010).

One key part of topological data analysis methods is to distinguish between different cycles; more precisely, to characterize different cycles according to their homology class. This can be done efficiently using cohomology (Pokorny, Hawasly & Ramamoorthy, 2016; Edelsbrunner, Letscher & Zomorodian, 2002). However, such methods only distinguish between non-homologous cycles, and do not quantify the difference between cycles. A possible way to quantify this difference is to solve the problem of finding the chain whose boundary is the union of the two cycles in question, as was proposed in Chambers & Vejdemo-Johansson (2015) by solving the underlying linear system, where the authors also take the question of optimality into account (with regards to the size of the resulting chain). This line of research ellaborates on Dey, Hirani & Krishnamoorthy (2011) where the authors considered the problem of finding an optimal chain in the same homology class as a given chain. More recently, in Rodríguez et al. (2017) the authors have taken a similar approach to ours using a combinatorial method to compute bounding chains of 1-cycles on 3-dimensional simplicial complexes.

In this paper, we explore the geometric properties of simplicial n-manifolds to provide an algorithm that is able to calculate a chain whose boundary is some prescribed (n − 1)-dimensional cycle, and we show that the proposed algorithm has a complexity which is linear in the number of (n − 1)-faces of the complex. This is substantially faster than calculating a solution to the linear system as considered in Chambers & Vejdemo-Johansson (2015), for which the complexity would be, at least, quadratic on the number of n-dimensional faces (Gloub & Van Loan, 1996).

Methodology

For any simplicial complex S, and any pair of simplices σ, τ ∈ S such that, τ ◃ σ, we define the index of τ with respect to σ as 〈τ, ∂σ〉 = 〈e_τ, de_σ〉 (Forman, 1998). Note that the index corresponds to the orientation induced by the boundary of σ on τ and can be computed in O(d) time by the following algorithm:

 
___________________________________________________________________________ 
Index(σ,τ): 
     param: σ - k-simplex represented as a sorted list of indices of points. 
     param: τ - (k − 1)-face of σ represented as a sorted list of indices. 
1: for each   i ← 0...dimτ : 
2:      if   τi ⁄= σi: 
3:         orientation ← (−1)i 
4:         break loop 
5: return orientation 
___________________________________________________________________________    

By inspecting the main loop we can see that Index (σ, τ) returns (−1)ⁱ where i is the index of the first element at which σ and τ differ. We assume τ is a top dimensional face of σ, so if σ = [s_0:d], then by definition τ = [s_0:(i−1), s_(i+1):d] for some i, and so the coefficient of τ in the boundary of σ is (−1)ⁱ as per the definition of the boundary operator. This is also the index of the first element at which τ and σ differ, since they are represented as sorted lists of indices.

The following is an intuitive result, that will be the basis for the Coefficient-Flow algorithm that we will present in the sequel.

Proof. If we expand the equation ∂c = p, we get p_τ = ∑_τ◃ω〈e_τ, d(c_ωe_ω)〉, recall that by definition d(c_ωe_ω) = ∑_ν◃ω〈ν, ∂ω〉c_ωe_ν; and so we get p_τ = ∑_τ◃ω〈τ, ∂ω〉c_ωe_τ.

Now since S is a manifold-like simplicial complex and τ = ∂σ∩∂σ′, then σ, σ′ are the only cofaces of τ, and hence we have: ${\displaystyle p_{\tau }=〈\tau ,\partial \sigma 〉c_{\sigma }+〈\tau ,\partial {\sigma }^{\prime }〉c_{{\sigma }^{\prime }}}$

which can be reorganized to $c_{\sigma }=\frac{p_{\tau }-〈\tau ,\partial {\sigma }^{\prime }〉c_{{\sigma }^{\prime }}}{〈\tau ,\partial \sigma 〉}$. Finally, since the index 〈τ, ∂σ〉 is either 1 or −1, we can rewrite this equation as: ${\displaystyle c_{\sigma }=〈\tau ,\partial \sigma 〉\left(p_{\tau }-〈\tau ,\partial {\sigma }^{\prime }〉c_{{\sigma }^{\prime }}\right)}$

Next, we present an algorithm to calculate a bounding chain for a (n − 1)-cycle in an n-manifold-like simplicial complex. The algorithm proceeds by checking every top-dimensional face σ, and calculating the value of the chain on adjacent top-dimensional faces, using Proposition 1.

In order to prove that the Coefficient-Flow algorithm solves problem (1), will use the fact that we can see the algorithm as a traversal of the dual graph.

 
________________________________________________________________________________________________________________________________________________ 
 
Coefficient-Flow(p,σ0,v0): 
     param: p — an n − 1 boundary of the simplicial complex S 
    param: σ0 — an n–simplex from where the calculation will start 
     param: v0 — the value to assign the bounding chain at sigma 
     return: c — a bounding chain of p satisfying cσ0  = v0 (if it exists) 
6:   initialize c and mark every n- and (n − 1)-cell of S as not seen. 
7:   initialize an empty queue Q 
8:   let τ0 ∈ ∂σ0 
9:   enqueue (σ0,τ0,v0) into Q. 
10:   while   Q is non-empty: 
11:       (σ,τ,v) ← pop first element from Q 
12:       if   σ has been marked seen: 
13:            if   v ⁄= cσ: the problem has no solution 
14:       else: 
15:            if   τ has been marked seen: skip 
16:            cσ ← v 
17:            mark τ and σ as seen 
18:            for each   τ′ ∈ ∂σ: 
19:                if   σ is the only coface of τ′: 
20:                    mark τ′ as seen 
21:                    if   pτ′ ⁄= 〈∂σ,τ′〉v : the problem has no solution 
22:                else: 
23:                    if   τ′ has not been marked as seen: 
24:                        σ′ ← other coface of τ′ 
25:                        v′ ←〈∂σ′,τ〉(pτ −〈∂σ,τ〉v) 
26:                        enqueue (σ′,τ′,v′) into Q 
27:   return c 
__________________________________________________________________________________________________________________________________    

Proof. We start by proving the bound on the number of executions of the main loop. This is guaranteed by the fact that the loop is executed while the queue is non-empty, and a triple (σ, τ, v) can only be inserted in line (21) if τ has not been marked as seen. Furthermore, since τ has at most two cofaces, say, σ, σ′, we can only enqueue τ if we are analyzing σ or σ′ and so for each τ, at most two elements are placed in the queue, and hence the main loop gets executed at most 2|S⁽ⁿ⁻¹⁾| times.

To prove correctness of the algorithm, we have to prove that it outputs an error if and only if the problem has no solution, and otherwise it outputs a bounding chain of p with c_σ0 = v.

First, we note that if a face σ is marked as seen, the value of c_σ can never be reassigned. This is because the program branches on whether or not σ has been marked as seen, and c_σ can only be assigned on line (11) which is bypassed if σ has been previously marked as seen. From this fact we conclude that c_σ0 = v₀ as it is assigned on line (11) and σ₀ is marked as seen in the first iteration of the main loop.

Second, note that there is an edge between two n-faces in the dual graph if and only if they share an (n − 1)-face. This implies that as we execute the algorithm, analyze a new n-face σ and successfully add the other cofaces of elements of the boundary of σ, we add the vertices neighboring σ in the dual graph. Since the dual graph is connected all of the nodes in the graph are eventually added, and hence all of the n-faces are analyzed.

Third, we note that for any pair (τ, σ) with dimσ = n and τ ◃ σ, either σ is the only coface of τ, or τ has another coface, σ′. In the first case, if p_τ ≠ c_σ〈τ, ∂σ〉 an error is detected on line (16). In the second case, assuming that the triple (σ, τ, v) is enqueued before (σ′, τ, v′) we have v′ = 〈∂σ′, τ〉(p_τ − 〈∂, σ〉v) as is assigned in line (20) then

${\displaystyle {\left(dc\right)}_{\tau }=〈\partial \sigma ,\tau 〉v+〈\partial {\sigma }^{\prime },\tau 〉v^{\prime }=〈\partial \sigma ,\tau 〉v+〈\partial {\sigma }^{\prime },\tau 〉\left(〈\partial {\sigma }^{\prime },\tau 〉\left(p_{\tau }-〈\partial ,\sigma 〉v\right)\right)=〈\partial \sigma ,\tau 〉v+p_{\tau }-〈\partial \sigma ,\tau 〉v=p_{\tau }.}$

Finally, since upon the successful return of the algorithm, this equation must be satisfied by every pair τ◃σ, it must be the case that dc = p. If this is not the case, then there will be an error in line (08) and the algorithm will abort.□

Note that the connectivity condition can be removed if we instead require a value for one cell in each connected component of the graph G(S), and throw an error in case there is an (n − 1)-simplex τ with no cofaces, such that p_τ ≠ 0. Furthermore, the algorithm can be easily parallelized using a thread pool that iteratively processes elements from the queue.

Finally, in the case where it is known that S has an (n − 1)-face τ with a single coface, we do not need to specify σ or v in Coefficient-Flow, and instead use the fact that we know the relationship between the coefficient p_τ and that of its coface in a bounding chain c of p, i.e., p_τ = 〈∂σ, τ〉c_σ. This proves Corollary 1.

 
____________________________________________________________________________________________________________________________________________________ 
 
Bounding-Chain(p): 
     param: p — an n − 1 boundary of the simplicial complex S 
    return: c — a bounding chain of p 
28:   for each   τ ∈ Sn−1: 
29:       if   τ has a single coface: 
30:            σ ← the coface of τ 
31:            v ←〈∂σ,τ〉pτ 
32:            break loop 
33:   c ← Coefficient-Flow(p,σ,v) 
34:   return c 
_______________________________________________________________________________________________________________________________________    

Example runs and tests

In Fig. 1 we provide an example of the output of the Coefficient-Flow algorithm for the mesh of the Stanford bunny (The Stanford 3D scanning repository, 1994) and the eulaema meriana bee model from the Smithsonian 3D model library (Smithsonian, 2013).

For comparison, we performed the same experiments using the Eigen linear algebra library (Eigen, 2017) to solve the underlying linear system,² and summarized the results in Table 1. This allowed us to see that even though both approaches remain feasible with relatively large meshes, solving the linear system consistently underperforms using the Coefficient-Flow algorithm.

Table 1:

Timing for computation of bounding chains using Coefficient-Flow, and using Eigen in several meshes.

“Bunny” and “Bee (Sample)/Bee (Full)” refer to the meshes in Figs. 1A and 1B, respectively. The mesh “Bee (Full)” is the one obtained from Smithsonian (2013), whereas the one labeled “Bee (Sample)” is a sub-sampled version of it.

Mesh	Faces	Edges	Vertices	Eigen	Coefficient-Flow
Bunny	69 663	104 496	34 834	2.073 (s)	0.48911 (s)
Bee (Sample)	499 932	749 898	249 926	116.553 (s)	3.04668 (s)
Bee (Full)	999 864	1 499 796	499 892	595.023 (s)	7.15014 (s)

DOI: 10.7717/peerjcs.153/table-1

Even though Coefficient-Flow is expected to outperform a linear system solver (an exact solution to a linear system has Ω(n²) time complexity), we wanted to test it against an approximate sparse system solver. Such solvers (e.g., conjugate gradient descent (Gloub & Van Loan, 1996)) rely on iterative matrix products, which in the case of boundary matrices of dimension d can be performed in O((d + 1)n) where n is the number of d-dimensional simplices, placing the complexity of the method in Ω((d + 1)n). In order to observe the difference in complexity class we performed randomized tests on both implementations. In this scenario we made a mesh on a unit square from a random sample. By varying the number of points sampled from the square, we effectively varied the resolution of the mesh. Finally, at each resolution level we snapped a cycle onto the mesh, and computed its bounding chain using both Coefficient-Flow and by solving the sparse linear system as before. We plotted the timings in Fig. 2 from where we can experimentally observe the difference in the complexity class between our algorithm and the solution to the sparse linear system.

Furthermore by analyzing the Log-Log plot in Fig. 2B, we can confirm our complexity estimates by analyzing the slope of the lines where the samples are distributed, i.e., solving the sparse linear system is approximately O(n^1.7) complexity,³ whereas Coefficient-Flow displays linear complexity.

Conclusion and future work

While the problem of finding a bounding chain for a given cycle in a simplicial complex remains a challenging one for large complexes, we showed that this problem can be solved efficiently for codimension-1 boundaries. We implemented and tested our algorithm and have provided complexity bounds for its run-time.

However, this leaves open the question of finding bounding chains for boundaries of higher codimension, for which solving a large sparse linear system is still, to the best of our knowledge, the only feasible approach, save for codimension 2 cycles in dimension 3 (Rodríguez et al., 2017). In the future we would like to generalize our algorithm to be able to work with cobounding cochains (i.e., in cohomology), as well as considering the optimality question (i.e., finding the minimal bounding chain w.r.t. some cost function).

Supplemental Information