Jean-Daniel Boissonnat
email: jean-daniel.boissonnat@inria.fr

Clément Maria
email: clement.maria@inria.fr

Computing Persistent Homology with Various Coefficient Fields in a Single Pass

.

Introduction

Persistent homology [START_REF] Edelsbrunner | Computational Topology -an Introduction[END_REF][START_REF] Zomorodian | Computing persistent homology[END_REF] is an algebraic method for measuring the topological features of the sublevel sets of a function defined on a topological space. Its generality and stability [START_REF] Cohen-Steiner | Stability of persistence diagrams[END_REF] with regard to noise have made it a widely used tool for the study of data. At the algebraic level [START_REF] Zomorodian | Computing persistent homology[END_REF], it admits a decomposition -represented by mean of a persistence diagram -only when considered with field coefficients (by opposition to integer coefficients). The persistence diagram contains a rich information about the topology of the studied space and very efficient methods exist to compute it. However, the integral homology groups of a topological space are strictly more informative than the homology groups with field coefficients, in particular because they convey information about "torsion". Torsion can be pictured geometrically as a "twisting" of the shape and happens frequently as global topological feature in topological data analysis where, for example, Klein bottles appear naturally [START_REF] Carlsson | On the local behavior of spaces of natural images[END_REF][START_REF] Martin | Topology of cyclo-octane energy landscape[END_REF]. Algebraically, torsion is characterized by cyclic subgroups of the integral homology groups. When computed with field coefficients, these subgroups may either vanish or appear as "infinite", and consequently obfuscate the study of the topology of data. A simple solution is to compute persistent homology with different coefficient fields and track the differences in the persistence diagrams.

We build on this idea and describe an algorithm to compute persistent homology with various coefficient fields Z q1 , • • • , Z qr in a single pass of the matrix reduction algorithm, where Z q denotes the finite field Z/qZ for a prime q. To do so, we introduce a method we call modular reconstruction consisting in using the Chinese Remainder Isomorphism to encode an element of Z q1 ו • •×Z qr with an element of Z q1•••qr . We describe algorithms to perform elementary row/column operations in a matrix with Z q1•••qr coefficients, corresponding to simultaneous elementary row/column operations in matrices with coefficients in the fields Z q1 , • • • , Z qr . The method results in an algorithm with an output-sensitive complexity in the total number of distinct pairs in the echelon forms of the matrices with Z q1 , • • • , Z qr coefficients, plus an overhead due to arithmetic operations on big numbers in Z q1•••qr . The method is generic and applies to every algorithm for persistent homology. Finally, we describe how to infer the torsion coefficients of the integral homology using the Universal Coefficient Theorem for Homology.

We provide detailed experimental analysis of the algorithm and show, in particular, that on practical examples our method is substancially faster than the brute-force approach consisting in reducing separately r matrices with coefficients in Z q1 , • • • , Z qr . It is important to note that the method does not pretend to scale to very large r, as the arithmetic complexity of operations in Z q1•••qr becomes problematic. Experiments show, however, that for very large r (up to 100000) our approach is still substancially faster than brute-force.

Computing persistent homology with different coefficients has been mentioned in the literature [START_REF] Zomorodian | Computing persistent homology[END_REF] in order to verify if a persisting feature was due to an actual "hole" (or high-dimensional equivalent) or to torsion (and consequently existed only for a certain coefficient field). However, to the best of our knowledge, this is the first work formalizing the inference of torsion coefficients in the framework of persistent homology and describing an efficient algorithm to compute persistence with various coefficient fields.

Multi-Field Persistent Homology

For simplicity, we focus in the following on simplicial homology. However, the approach applies to any type of boundary matrix (defined below).

Background on Simplicial Homology and Persistence

: A simplicial com- plex K on a set of vertices V = {1, • • • , n} is a collection of simplices {σ}, σ ⊆ V , such that τ ⊆ σ ∈ K ⇒ τ ∈ K. The dimension d = |σ| -1 of σ is its number of elements minus 1. For a ring R, the group of d-chains, denoted C d (K, R), of K is the group of formal sums of d-simplices with R coefficients. The boundary operator is a linear operator ∂ d : C d (K, R) → C d-1 (K, R) such that ∂ d σ = ∂ d [v 0 , • • • , v d ] = d i=0 (-1) i [v 0 , • • • , v i , • • • , v d ],
where v i means v i is deleted from the list. It will be convenient to consider later the endomorphism

∂ * : d C d (K, R) → d C d (K,
H d (K, R) of K by the quotient H d (K, R) = Z d (K, R)/B d (K, R).
If R is the ring of integers Z, H d (K, Z) is an abelian group and, according to the fundamental theorem of finitely generated abelian groups [START_REF] Munkres | Elements of algebraic topology[END_REF], admits a primary decomposition:

H d (K, Z) ∼ = Z β d (Z) q prime Z q k 1 ⊕ • • • ⊕ Z q k t(d,q
) for uniquely defined integer β d (Z), called the d th integral Betti number, and integers t(d, q) ≥ 0 and k i > 0 for every prime number q. If t(d, q) > 0, the integers q k1 , • • • , q k t(d,q) are called torsion coefficients, and they admit q as unique prime divisor. Intuitively, in dimension 0, 1 and 2, the integral Betti numbers count the number of connected components, the number of holes and the number of voids respectively. The torsion coefficients represent non-orientable twisting of different order of the shape. If R is a field F, H d (K, F) is a vector-space and decomposes into F) , where β d (F) is the d th field Betti number. The field Betti numbers (β d (F)) d are entirely determined by the characteristic of F and the inegral homology (see Section 5); hence, the integral homology is more informative than homology in F.

H d (K, F) ∼ = F β d (
A filtration of a complex is a function f :

K → R satisfying f (τ ) ≤ f (σ) whenever τ ⊆ σ. The sequence [σ i ] i=1,••• ,m sorted according to increasing f val- ues induces the sequence of inclusions ∅ = K 0 K 1 • • • K m-1 K m = K, K i = K i-1 ∪ {σ i }, and the sequence of d-homology groups 0 = H d (K 0 , R) → H d (K 1 , R) → • • • → H d (K m-1 , R) → H d (K m , R) = H d (K, R) connected by homomorphisms.
When R is a field, the later sequence admits a decomposition in terms of intervals {(i, j)}, called an indexed persistence diagram, where a pair (i, j) is interpreted as a homology feature that is born at index i and dies at index j. Computing persistent homology consists in computing the interval decomposition and hence the persistence diagram. We refer to [START_REF] Munkres | Elements of algebraic topology[END_REF] for an introduction to homology and to [START_REF] Edelsbrunner | Computational Topology -an Introduction[END_REF] for an introduction to persistent homology.

We call the algorithmic problem of computing persistent homology with various coefficient fields multi-field persistent homology. Computing multi-field persistence allows us to infer a more informative description of the topology of a space, compared to persistence in a single field (see Section 5 for details). For a complex of size m, we know that the persistence diagram for any coefficient field contains at most m pairs. When computing multi-field persistent homology for r coefficient fields, denote by m the total number of distinct pairs in all persistence diagrams. In practice, the fields are Z q1 , • • • , Z qr for the r first prime numbers q 1 , • • • , q r , where q r is an upper bound on the prime divisors of the torsion coefficients of the integral homology of the space, which are usually small (see Section 5). The quantity m satisfies m ≤ r × m but in practice we observe that m ≈ m. We design in the following an algorithm for the multi-field persistence problem whose complexity depends mostly on m . It is however an interesting open problem to exhibit a "natural" example where m is must larger than m and/or the prime divisors of the torsion coefficients are big (for a fix m).

Algorithm for Multi-Field Persistent Homology

For clarity, we focus in this section on the persistent homology algorithm as presented in [START_REF] Edelsbrunner | Computational Topology -an Introduction[END_REF], which consists in a reduction to column echelon form (defined later) of a matrix. All other persistent homology algorithms are based on similar reductions, and our approach adapts directly to them. In the following, Z n denotes the ring (Z n , +, ×) for any integer n ≥ 1. and Z × n the subset of invertible elements for ×. If exists, we denote the inverse of x ∈ Z n by x -1 .

In computer algebra, working modulo small prime numbers is usually desireable in order to reduce coefficient growth. Our work goes the otherway around: we introduce tools to reduce a family of r matrices with coefficients in the fields Z q1 , • • • , Z qr respectively, by means of a single reduction of a matrix with coefficients in Z q1•••qr . We give theoretical and experimental evidence that, for reasonable values of r, our algorithm is significantly more efficient than the brute-force approach consisting in reducing the r matrices separately.

Persistent Homology Algorithm: For an m × m matrix M, denote by C j the j th column of M, 1 ≤ j ≤ m, and C j [k] its k th coefficient. Let low(j) denote the row index of the lowest non-zero coefficient of C j . If the column j is entirely zero, low(j) is undefined. We say that M is in reduced column echelon form if low(j) = low(j ) for every non-zero columns C j and C j with j = j .

Let K = [σ i ] i=1•••m be a filtered complex. Its boundary matrix M ∂ is the m× m matrix, with F coefficients, of the endomorphism ∂ * in the basis {σ 1 , • • • , σ m } of d C d (K, F).
The basis is ordered according to the filtration. It is a matrix with {-1, 0, 1} coefficients, where 0 and 1 are the identities for + and × in F respectively, and -1 is the inverse of 1 in F. The persistent homology algorithm consists in a left-to-right reduction to column echelon form of M ∂ : we denote by R the matrix we reduce, with columns C j , which is initially equal to M ∂ . The algorithm returns the (indexed) persistence diagram, which is the set of pairs {(low(j), j)} in the reduced column echelon form of the matrix.

Data:

Boundary matrix R ← M ∂ , persistence diagram P ← ∅ Output: Persistence diagram P = {(i, j)} 1 for j = 1, • • • , m do 2
while there exists j < j with low(j ) = low(j) do

3 k ← low(j); 4 C j ← C j -C j [k] × C j [k] -1 • C j ; 5 end 6 if C j = 0 then P ← P ∪ {(low(j), j)} 7 end
The reduced form of the matrix is not unique, but the pairs (i, j) such that i = low(j) in the column echelon form are [START_REF] Edelsbrunner | Computational Topology -an Introduction[END_REF]. The algorithm requires O(m 3 ) arithmetic operations in F.

Modular Reconstruction for Elementary Matrix Operations:

We present a particular case of the Chinese Remainder Theorem [START_REF] Von | Modern Computer Algebra[END_REF]: For a family {q 1 , • • • , q r } of r distinct prime numbers, there exists a ring isomorphism ψ :

Z q1 × • • • × Z qr → Z q1•••qr . The isomorphisms ψ and ψ -1 can be computed in O(r) arithmetic operations in Z q1•••qr .
Let [r] refer to the set {1, • • • , r}. For a family of r distinct prime numbers {q 1 , • • • , q r }, and a subset of indices S ⊆ [r], Q S refers to s∈S q s , and we write simply

Q = Q [r]
. We define the function ψ S : s∈S Z qs → Z Q S realizing the isomorphism of the Chinese Remainder Theorem for the subset {q s } s∈S of primes, and we write simply ψ for ψ [r] . For a family of elements u s ∈ Z qs , s ∈ S, we denote the corresponding |S|-uplet (u s ) s∈S ∈ s∈S Z qs .

Finally, we recall Bezout's lemma [START_REF] Von | Modern Computer Algebra[END_REF]: For two integers a and b, not both 0, there exist integers v and w such that va + wb = gcd(a, b), the greatest common divisor of a and b, with |v| < |b/ gcd(a, b)| and |w| < |a/ gcd(a, b)|. The Bezout's coefficients (v, w) can be computed with the extended Euclidean algorithm [START_REF] Von | Modern Computer Algebra[END_REF]. Elementary Column Operations: We are given a family of distinct prime numbers {q 1 , • • • , q r }, and their product

Q = q 1 • • • q r . Let M Q be a matrix with coefficients in the ring Z Q . Denoting ψ -1 : Z Q → Z q1 × • • • × Z qr
the isomorphism of the Chinese Remainder Theorem, and π s :

Z q1 × • • • × Z qr → Z qs the projection on the s th coordinate, we call projection of M Q onto Z qs , denoted M Q (Z qs ), the matrix M qs with Z qs coefficients, obtained by applying π s • ψ -1 to each coefficient of M Q . Conversely, given r (m × m)-matrices M q1 , • • • , M qr with coefficients in Z q1 , • • • , Z qr respectively,
there exists a unique matrix M Q with Z Q coefficients such that for every s the projection of M Q onto Z qs is M qs . This is simply a matrix version of the Chinese remainder theorem. Elementary column operations on M q with Z q coefficients are of three kinds:

(i) exchange Col k and Col (ii) multiply Col k by -1 ∈ Z q (iii) replace Col k by (Col k)+ x×(Col ), for x ∈ Z q .
For an elementary column operation ( * ) (i.e. an operation of type (i), (ii) or (iii) applied to columns k (and )), we denote by ( * ) • M q the result of applying ( * ) to M q . In this section, we introduce algorithms to run elementary column operations simultaneously on the matrices (M qs ) s=1,••• ,r by performing "partial column operations" on M Q . Specifically, for an elementary column operation ( * ) and a subset of indices S ⊆ [r], we call partial column operation on M Q the operation transforming M Q into M Q such that: for every s / ∈ S, the projection onto

Z qs satisfies M Q (Z qs ) = M Q (Z qs ) = M qs and for every s ∈ S, the projection onto Z qs satisfies M Q (Z qs ) = ( * ) • M qs . As the correspondence ψ : Z q1 × • • • × Z qr → Z Q is a ring homomor- phism, it satisfies the properties: ψ(u 1 , • • • , u r )+ψ(v 1 , • • • , v r )×ψ(w 1 , • • • , w r ) = ψ(u 1 +v 1 ×w 1 , • • • , u r +v r ×w r
) and we can compute addition and multiplication componentwise in Z q1 × • • • × Z qr using addition and multiplication in Z Q . In order to compute partial column operations, we first introduce the set of partial identities, which are coefficients that allow us to proceed to the partial column operations of type (i) and (ii). Secondly, as the rings Z qs are fields, we need to compute the multiplicative inverse of an element, that is used as multiplicative coefficient x in elementary column operation (iii). As Z Q is not a field, inversion is not possible, and we introduce the concept of partial inverse to overcome this difficulty. In the following, the term "arithmetic operation" refers to any op-eration {+, -, ×, gcd(•, •), • mod Q S , Extended Euclidean algorithm} on integer smaller than Q. Note they do not have constant time complexity for large Q. Partial Identity and Partial Inverse: Given a subset of indices S ⊆ [r], we define the partial identities w.r.t. S, denoted L S , by L S = ψ(δ 1,S , • • • , δ r,S ) where the symbol δ t,S ∈ Z qt is equal to 1 if t ∈ S and to 0 otherwise. For any S ⊆ [r], the partial identity L S can be constructed in O(r) arithmetic operations in Z Q by evaluating ψ on (δ 1,S , • • • , δ r,S ). However, it is important to notice that if We define now the partial inverse of an element in the ring Z Q :

S = [r], L [r] = ψ(1, • • • , 1) = 1,
Definition 1 (Partial Inverse) Given a set S ⊆ [r] of indices, the partial inverse of x = ψ(u 1 , • • • , u r ) with regard to S is the element x S ∈ Z Q : x S = ψ(u 1 S , • • • , u r S ), with u s S = u -1 s if s ∈ S and u s ∈ Z × qs 0 otherwise
Using elementary algebra (see [START_REF] Boissonnat | Computing persistent homology with various coefficient fields in a single pass[END_REF] for details) we prove:

Proposition 2 (Partial Inverse Construction) For x = ψ(u 1 , • • • , u r ) ∈ Z Q and S ⊆ [r],
(1) gcd(x, Q S ) = Q R for some R ⊆ S and for all s ∈ S, u s is invertible in Z qs iff s / ∈ R; we denote T = S \ R. (2) The Bezout's identity for x and Q T gives vx + wQ T = 1, where v satisfies

v mod Q T = ψ T ((u -1 s ) s∈T ) (3) x S = ψ T ((u -1 s ) s∈T ) × L T mod Q ∈ Z Q
, where L T is the partial identity w.r.t T .

We deduce directly an algorithm to compute the partial inverse of x w.r.t S if Q S is given: compute Q R = gcd(x, Q S ) and Q T = Q S /Q R , then v using the extended Euclidean algorithm and finally

x S = (v mod Q T ) × L T mod Q. Computing the partial identity L T requires O(r) arithmetic operations in Z Q , but is constant if T = [
r], which happens iff S = [r] and x is invertible in Z Q . Consequently, computing x S requires O(r) arithmetic operations in general, but only O(1) arithmetic operations in the later case.

Modular Reconstruction for Multi-Field Persistent Homology

Let K be a filtered complex with m simplices. Define M ∂ (Z qs ) to be the (m×m) boundary matrix of K with Z qs coefficients. Define M to be the (m × m) matrix with Z Q coefficients such that the projection of M onto Z qs is equal to M ∂ (Z qs ), for all s ∈ [r]. Note that the matrices M and M ∂ (Z qs ), for any s, are "identical" matrices in the sense that they contain 0, 1 and -1 coefficients at the same positions, where 0, 1 and -1 refer respectively to elements of Z Q and Z qs .

We reduce a matrix R which is initially equal to M. Denote by C j the j th column of R. Define low(j, Q S ) to be the index of the lowest element of C j such that C j [low(j, Q S )] mod Q S = 0. In particular, low(j, q s ) is equal to the index of the lowest non-zero element of column j in the projection R(Z qs ). After iteration j, we say that the columns C 1 , • • • , C j are reduced. We maintain, for every reduced column C j , the collection of "lowest indices" i as a set L(j) = {(i, Q S )} satisfying:

-For every (i, Q S ) ∈ L(j), i = low(j, Q S ) -For every (i, Q S ), (i , Q S ) ∈ L(j), either i = i and S = S , or i = i and

S ∩ S = ∅ -∪ (i,Q S )∈L(j) S = [r]
The algorithm returns the set of triplets P = {(i, j, Q S )} such that i = low(j) in the column echelon form of the matrix M ∂ (Z qs ) iff s ∈ S, or, equivalently, (i, Q S ) ∈ L(j) once C j has been reduced. This is a compact encoding of the persistence diagrams of the filtered complex in persistent homology with all coefficient fields. We call it multi-field persistence diagram.

Data: Matrix R = M Output: Multi-field persistence diagram P = {(i, j, Q S )} 1 for j = 1, • • • , m do 2 Q S ← Q [r] ; 3 while low(j, Q S ) is defined do 4 k ← low(j, Q S ); Q T ← Q S / gcd(C j [k], Q S ) ; 5
while there exists j < j with (i, Q T ) ∈ L(j ) satisfying

6 [i = low(j, Q S ) and gcd(Q T , Q T ) > 1] do 7 C j ← C j -C j [k] × C j [k] T • C j ; 8 Q T ← Q T / gcd(Q T , Q T ); 9 end 10 if Q T = 1 then P ← P ∪ {(k, j, Q T )}; Q S ← Q S /Q T ; 11 end 12 end
The {L(j)} j form an index table that we maintain implicitely. At iteration j of the for loop, we use Q S for the product of all prime numbers s∈S q s for which the column j in R(Z qs ) has not yet been reduced.

Correctness: First, note that all operations processed on R correspond to leftto-right elementary column operations in the matrices R(Z qs ) for all s ∈ [r]. By definition of the partial inverse, the column operation in line 7 can only reduce the value of low(j, Q S ). Moreover, one iteration of the while loop in line 3 either strictly reduces Q S by dividing it by Q T (in line 10) or set (C j [k] mod Q S ) to zero, hence reducing strictly low(j, Q S ). The later case happens when Q T is set to 1 in line 8. Consequently, the algorithm terminates.

We prove recursively, on the numbers of columns, that each of the matrix R(Z qs ) gets reduced to column echelon form. We fix an arbitrary field Z qs : suppose that the j -1 first columns of R(Z qs ) have been reduced at the end of iteration j -1 of the for loop in line 1. We prove that at the end of the j th iteration of the for loop in line 1, the j first columns of the matrix R(Z qs ) are reduced. Consider two cases. First suppose there is a triplet (i, j, Q T ) ∈ P for some i < j and Q T satisfying q s | Q T . This implies that the algorithm exits the while loop in line 5 with

q s | Q S (because by definition of Q T , in line 4, Q T | Q S ) and there is no j < j such that [low(j , Q T ) = low(j, Q S ) and gcd(Q T , Q T ) > 1].
This in particular implies that there is no j < j such that low(j , q s ) = low(j, q s ) and column j is reduced in R(Z qs ).

Secondly, suppose that there is no such pair (i, j, Q T ) in P, with q s dividing Q T . Consequently, during all the computation of the while loop in line 3, q s | Q S . When exiting this while loop, low(j, Q S ) is undefined, implying in particular that low(j, q s ) is undefined and column j of R(Z qs ) is zero, and hence reduced.

Output-Sensitive Complexity Analysis

Arithmetic Complexity Model for Large Integers: During the reduction algorithm we perform arithmetic operations on big integers, for which we describe a complexity model [START_REF] Von | Modern Computer Algebra[END_REF]. Suppose that on our architecture, a memory word is encoded on w bits (on modern architectures, w is usually 64). Computer chips contains Arithmetic Logic Units that allow arithmetic operations on a 1-memory word integer in O(1) machine cycles. Let the length of an integer z be defined by: λ(z) = log 2 z/w + 1, i.e. by the number of memory words necessary to encode z. We express the arithmetic complexity as a function of the length. For any positive integer z of length λ(z) = B, operations in Z z cost A + (z) = O(B) for addition, A × (z) = O(M (B)) for multiplication and A ÷ (z) = O(M (B) log B) for (extended) Euclidean algorithm, inversion and division, where M (n) is a monotonic upper bound on the number of word operations necessary to multiply two integers of length B. By a result of [START_REF] Fürer | Faster integer multiplication[END_REF], M (B) = O(B log B 2 O(log * B) ), where log * n is the iterated logarithm of n. In the following, we write A(z) for a bound on the complexity of arithmetic operations on integers smaller than z.

In the case of multi-field persistent homology, we are interested in the value of λ for an element in Z Q , Q = q 1 • • • q r , in the case where {q 1 , • • • , q r } are the first r prime numbers. We know [START_REF] Rosser | Approximate formulas for some functions of prime numbers[END_REF] that ln Q < 1.01624q r and q r < r ln(r ln r) for r ≥ 6. Consequently, λ(Q) < 1.46613r ln(r ln r)/w + 1. Note that λ(Q) r for r ln r e w , which is a reasonable assumption.

Complexity of the Modular Reconstruction Algorithm: Let K be a filtered complex of size m. The persistent homology algorithm described in Section 3, applied on K with coefficients in a field F, requires O(m 3 ) operations in F. For a field Z q these operations take constant time and the algorithm has complexity O(m 3 ). The output of the algorithm is the persistence diagram, which has size O(m) for any field.

For a set of prime numbers {q 1 , • • • , q r }, let m be the total number of distinct pairs in all persistence diagrams for the persistent homology of K with coefficient fields Z q1 , • • • , Z qr . We express the complexity of the modular reconstruction algorithm in terms of the size of its output (i.e. the multi-field persistence diagram of size m ), the number of fields r and the arithmetic complexity A(Q). First, note that, for a column j in the reduced form of R, the size of L(j ) is equal to the number of triplets of the multi-field persistence diagram with death index j ; denote this quantity by |L(j )|. Hence, when reducing column j > j , the column C j is involved in a column operation

C j ← C j + α • C j at most |L(j )| times. Consequently, reducing C j requires O( j <j |L(j )|) = O(m ) column op- erations.
There is a total number of O(m × m ) column operations to reduce the matrix, each of them being computed in time O(m × A(Q)).

Computing the partial inverse of an element x ∈ Z Q takes time O(r×A(Q)) in the general case, and only

O(A(Q)) if x is invertible in Z Q . The partial inverse of an element x = C j [k] is computed only if there is a pair (k, Q T ) ∈ L(j). This element is not invertible in Z Q iff |L(j)| > 1.
There are consequently O(|m -m|) non-invertible elements x that are at index low(j, Q T ) in some column j, for some Q T . If we store the partial inverses when we compute them, the total complexity for computing all partial inverses in the modular reconstruction algorithm is O((m + r × (m -m) × A(Q)). We conclude that the total cost of the modular reconstruction algorithm for multi-field persistent homol-

ogy is O( r × (m -m) + m 2 m × A(Q)) = O( r × (m -m) + (m ) 3 × A(Q)),
while the brute-force algorithm, consisting in computing persistence separately for every field Z q1 , • • • , Z qr requires O(r × m 3 ) operations.

Note that asymptotically in r, one arithmetic operation in Z Q [r] becomes more costly than r distinct arithmetic operations in Z q1 , • • • , Z qr , in which case the modular reconstruction approach developed in this article becomes worse than brute-force (even when m and m are close). This however happens for extremely big values of r (see Section 5) and has no incidence on practical cases.

Remark: On all datasets considered in our experiments, we have found no example where m was significantly bigger than m. However, it is unclear whether many "short-lived torsion" might appear in general. We prove in a long version of the paper [START_REF] Boissonnat | Computing persistent homology with various coefficient fields in a single pass[END_REF] that this is not an issue, by giving a finer complexity analysis of the algorithm in terms of index persistence |j -i| of the pairs (i, j) in the persistence diagram.

Experiments

In this section, we report on the performance of the modular reconstruction algorithm for multi-field persistent homology. Our implementation is in C++, and we use the GMP library for storing large integers. All timings are measured on a 64 bits Linux machine with 3.00 GHz processor and 32 GB RAM. All timings are averaged over 10 independent runs. We compute the persistent homology of Rips complexes [START_REF] Edelsbrunner | Computational Topology -an Introduction[END_REF] built on a variety of both real and synthetic datasets. We use the compressed annotation matrix implementation of persistence [START_REF] Boissonnat | The compressed annotation matrix: An efficient data structure for computing persistent cohomology[END_REF] efficiency and stability. Bud is a set of points sampled from the surface of the Stanford Buddha in R 3 . Bro is a set of 5 × 5 high-contrast patches derived from natural images, interpreted as vectors in R 25 , from the Brown database (with parameter k = 300 and cut 30%) [START_REF] Carlsson | On the local behavior of spaces of natural images[END_REF]. Cy8 is a set of points in R 24 , sampled from the space of conformations of the cyclo-octane molecule [START_REF] Martin | Topology of cyclo-octane energy landscape[END_REF], which is the union of two intersecting surfaces. Kl is a set of points sampled from the surface of the figure eight Klein Bottle embedded in R 5 . Finally S3 is a set of points distributed on the unit 3-sphere in R 4 . Datasets are listed in Figure 1 with the size of points sets |P|, the ambient dimension D and intrinsinc dimension d of the sample points (if known), the parameter r for the Rips complex and the size of the complex |K|. The values T r for r ∈ {1, 50, 100, 200} refers to the running time of the modular reconstruction algorithm for the r first prime numbers, and R r refers to the ratio between the brute-force approach and the modular reconstruction algorithm.

Interpretation of the Results: Surprisingly, we have observed that, on all experiments, the number of differences between persistence diagrams with various coefficient fields was extremely small. As a consequence, m -m can be considered as a very small constant in our experiments (≤ 10). We have also observed that these differences appeared for small prime numbers q s . Figure 1 presents the timings of the modular reconstruction approach for a variety of simplicial complexes ranging between 114 and 193 million simplices. We note that from r = 1 to r = 200 prime numbers, the time for computing multi-field persistence using the modular reconstruction approach only increases by 23 to 41%, when the brute-force approach requires about 200 times more time. This difference appears in the speedup expressed by the ratio R r . For r = 1, the modular reconstruction approach is about twice slower than the standard persistent homology algorithm in one field, because modular reconstruction is a more complex procedure and deals, in our implementation, with GMP integers that are slower than the classic int used in the standard persistent homology algorithm. However, this difference fades away as soon as r > 1 and the modular reconstruction is significantly more efficient than brute-force: it is, in particular, between 50.3 and 102.2 times faster for r = 200.

Figures 2 and3 present the evolution of the running time of the modular reconstruction approach and the brute-force approach for an increasing number of fields r (using the first r prime numbers). Persistence is computed for a Rips complex built on a set of 10000 points sampling a Klein bottle, which contains torsion in its integral homology, resulting in a simplicial complex of 6.14 million simplices. We analyze the in terms of the complexity analysis of Section 4. The quantity m is fixed and m is fixed for r ≥ 2. The complexity of the brute-force algorithm is O(r × m 3 ) and we indeed observe a linear behavior when r increases. The complexity of the modular reconstruction approach is O( r × (m -m) + m 3 A(Q [r] )). The part r × (m -m) of the complexity is negligeable because m -m is extremely small. For medium values of r (≤ 150), like in Figure 2, the arithmetic complexity O(A(Q [r] )) increases very slowly because λ(Q [r] ) = log 2 Q [r] /w + 1 increases slowly. We consequently observe a very slow increasing of the time complexity compare to the one of brute-force.

Figure 3 describes the asymptotic behavior of the modular approach, where the arithmetic operations become costly. We observe that the timings for the modular reconstruction approach follow a convex curve. The convexity comes from the growth of λ(Q [r] ), which is asymptotically Θ(r log r)) [START_REF] Rosser | Approximate formulas for some functions of prime numbers[END_REF]. However, the increasing of the slope is very slow: all along this experiment, we have been unable to reach a value of r for which the modular approach is worse than the bruteforce approach. For readability, the timings for the brute-force approach are implicetly represented through their ratio with the modular approach: all along the experiment, for 10000 ≤ r ≤ 100000, the modular approach is between 55 and 90 times faster. Based on a linear interpolation of the timings for the bruteforce approach, and a polynomial interpolation of the modular reconstruction timings, we expect the modular reconstruction to become worse than bruteforce for a number of primes r bigger than 4.9 million. In the case of multi-field persistent homology however, there is no need to take r bigger than 200, because r is related to torsion coefficients (see Section 5), which are small in practice.

Back to Topology: Inference of Torsion For a topological space X, the Universal Coefficient Theorem for Homology [START_REF] Munkres | Elements of algebraic topology[END_REF] establishes the relationship between the homology groups H d (X, Z) with Z coefficients and the homology

  R) extended by linearity to the external sum of chain groups. Denote by Z d (K, R) and B d-1 (K, R) the kernel and the image of ∂ d respectively. Observing ∂ d • ∂ d+1 = 0, we define the d th homology group

  because ψ is a ring isomorphism, and L r is computed in time O(1). Knowing the partial identities, we can implement the partial column operations (i) and (ii) for a set of indices S. Partial column operation (i) is implemented by replacing column k by (Col k × L [r]\S + Col × L S ) and column by (Col × L [r]\S + Col k × L S ). Partial column operation (ii) is implemented by multiplying column k by L [r] -2 × L S .

Fig. 1 .

 1 Fig. 1. Timings of the modular reconstruction algorithm vs brute-force.

Fig. 2 .Fig. 3 .

 23 Fig. 2. Timings for the modular reconstruction algorithm and brute force.

  for its

	Data	|P| D d	r	|K|	T1 R1	T50 R50 T100 R100 T200 R200
	Bud 49,990 3 2 0.09 127 • 10 6 96.3 0.51 110.3 22.2 115.9 42.3 130.7 75.0
	Bro 15,000 25 ? 0.04 142 • 10 6 123.8 0.41 143.5 17.8 150.2 34.0 174.5 58.5
	Cy8 6,040 24 2 0.8 193 • 10 6 121.2 0.63 134.6 28.2 139.2 54.6 148.8 102.2
	Kl	90,000 5 2 0.25 114 • 10 6 78.6 0.52 89.3 23.0 93.0 44.1 105.2 78.0
	S3	50,000 4 3 0.65 134 • 10 6 125.9 0.40 145.7 17.2 152.6 32.8 177.6 50.3

This research has been partially supported by the European Research Council under Advanced Grant 339025 GUDHI (Algorithmic Foundations of Geometric Understanding in Higher Dimensions).

groups H d (X, Z q ) with coefficients in the field Z q (of characteristic q), for q prime. We use the following corollary: Corollary 3 (Universal Coefficient Theorem [START_REF] Munkres | Elements of algebraic topology[END_REF]) For β d (Z) and β d (Z q ) the Betti numbers of H d (X, Z) and H d (X, Z q ) respectively, and t(j, q) the number of Z q k i summands in the primary decomposition of H j (X, Z), we have:

are the first r prime numbers and q r is a strict upper bound on the prime divisors of the torsion coefficients of X. Consequently, according to Corollary 3, β d (Z qr ) = β d (Z) for all dimensions d. Moreover, we know [START_REF] Munkres | Elements of algebraic topology[END_REF] that there is no torsion in 0-homology (i.e. t(0, q) = 0 for all primes q). Given the Betti numbers of X in all fields Z qs , 1 ≤ s ≤ r, we deduce from Corollary 3 the recurrence formula t(d, q s ) = β d (Z qs ) -β d (Z qr ) -t(d -1, q s ), from which we compute the value of t(d, q) for every dimension d and prime q. For any dimension d, we consequently infer the integral Betti numbers and the number t(d, q) of Z q k i summands in the primary decomposition of H d (X, Z). We note however that the powers k i from the decomposition remain unknown.

We describe in a long version of this article [START_REF] Boissonnat | Computing persistent homology with various coefficient fields in a single pass[END_REF] a representation of multifield persistence diagrams with torsion coefficients. We also describe an efficient algorithm to compute distances between them.