An Algorithmic Way to Generate Simplexes for Topological Data Analysis

Abstract

In this article we present a new algorithm for creating simplicial Vietoris-Rips complexes that is easily parallelizable using computation models like MapReduce and Apache Spark. The algorithm does not involve any computation in homology spaces.

Keywords

• data analysis
• topology
• simplicial complex

Introduction

The article of Silva and Grist [ 2 ] showed that the algebraic topology was a very practical tool. In this paper its authors analyse the coverage of an area with a network of sensors. This network is interpreted as a simplicial complex. Its homology type is then computed and exploited. If this complex is homotopy equivalent with point, the coverage is complete. Therefore, the authors translated a technical problem into a mathematical problem.

In the article [ 3 ] Grist searches for topology structures in big data sets. Again, this results in building simplicial complexes and analysing them.

The abovementioned works inspire to ask the question how to build simplicial complexes efficiently and how to analyse them. In this article w propose a novel algorithm to build the simplicial complex. This algorithm has the unique property that it can be implemented within massive parallel computational models like MapReduce and Apache Spark. have no insight into how the data actually look like and how it is embedded, but it would be enough to know their “shape” to say something about how it is are arranged in the original space.

In recent papers of Carlsson and collaborators [ 1 ] presented ideas to explore some properties of high-dimensional data sets. They use algebraic topology tools to find certain characteristics of the data or its simpler representation. This approach allows to discover objects and features that are immersed in highdimensional space.

We start presentation of basic facts from the notion of a simplex. Colloquially speaking, simplex is a generalization of a segment, a triangle and a tetrahedron to any dimension. More formally, suppose the k + 1 points p0; : : : ; pk 2 Rk are affinely independent, which means p1 p0; : : : ; pk p0 are linearly independent. Then, k-simplex determined by them is the set of points

( := 0p0 + + kpk j i 0; 0 i

k ) k; X i = 1 ; i=0 (1) i.e. k-simplex is a k-dimensional polytope which is the convex hull of its k + 1 vertices. We often write = [p0; : : : ; pk] for short. The points p0, : : :, pk are called vertices of the k-simplex. Each simplex is uniquely determined by its vertices. Any subset of vertices constitutes simplex with dimension smaller that the whole one, called a face of simplex.

We say two vertices v and w are neighbors if there is 1-simplex (edge) connecting them.

Definition 1. A simplicial complex K is a set of simplexes that satisfies the

following conditions:

1. Any face of a simplex from K is also in K.

2. The intersection of any two simplexes 1, 2 2 K is either ; or a face of

both 1 and 2.

Having data set P , we can use it to build a variety of complexes. We shall now formulate two most frequent definitions.

Definition 2. For a discrete set of points P

define the Cˇ ech complex of radius as Cˇ(P; ) by Rd and a parameter > 0, we Cˇ(P; ) := ( [p1; p2; : : : ; pk] j fp1; p2; : : : ; pkg

P; \ B(pi; =2) 6= ? i ) ; (2) where B(p; ) is the closed ball of radius centered at p.

The nerve theorem states that Cˇ ech complex Cˇ (P; ) has homotopy type of the union of closed balls with radius =2 around data P This means that the Cˇ ech complex gives faithful representation of thickened data (they are common ”shape”).

Vietoris-Rips complex is closely related to the Cˇ ech complex.

Definition 3. For a given > 0, the Vietoris-Rips complex (later called Rips

complex) R(P; ) is the largest simplicial complex that shares the same 1-skeleton (graph) as the Cˇech complex Cˇ(P; ). More explicitly,

R(P; ) := f := [p1; p2; : : : ; pk] j fp1; p2; : : : ; pkg P; diam( ) g ; (3) where diam( ) is the diameter of simplex .

Lemma 1 (Vietoris-Rips). For every finite set of points P

Rd and r Cˇ(P; )

R(P; )

Cˇ(P; 2 ): Thus, if the Cˇ ech complex at the same time for and 2 approximates the data in a good way, then Vietoris-Rips complex do it as well. This estimate can be further improved [ 2 ].

Remark 1. Cˇ ech complex is difficult to calculate, but it is quite small. However, Vietoris-Rips complex is easy to calculate, but is usually very big. Both Cˇ ech 0, (4) and Vietoris-Rips complexes can produce simplicial complex of dimension greater than the considered space. Therefore, there are considered many alternatives for them: Delaunay Complex, Alpha Complex, (Lazy) Witness Complex, Mapper Complex, Sparsified Rips complex, Graph Induced Complex (GIC).

For small > 0 we have disjoint set of balls, whereas for big the covering by balls is simply connected (without holes), so we are totally losing the knowledge about the data structure. Since (in general) we do not know which radius to take, we consider them all and obtain in this way a filtration of simplicial complexes, i.e.

Cˇ (P; 1)

Cˇ (P; 2) : : :

Cˇ (P; n); (5) for 0 < 1 2 : : : n. With the change of radius also changes the topology of the data: some holes are created, and some disappear. Roughly speaking, those holes that last longest represent the shape of data. It is usually represented in the form of the so-called barcode [ 3 ].

Let us recall the following definition from topology.

Definition 4. Two topological spaces X and Y are called homotopy equivalent if

there exist continuous maps f : X ! Y and g : Y ! X, such that the composition f g is homotopic (i.e. it can be deformed in a continuous way) to the identity idY on Y , and g f is homotopic to idX . 3

Related Work

Data anaylsis is a process of modeling data in order to discover useful information. Unfortunately, mining high-dimensional data is very challenging. Therefore, a number of methods was created so as to make it easier for researchers. One of the youngest but rapidly growing field now is topological data analysis (TDA).

TDA is an approach to analyse data sets using techniques from topology. This method relies on calculation of some properties of the space/data, which maintain (are invariant) under certain transformations. Although most of topological invariants are difficult to calculate, there is an algorithm which can calculate homology groups and barcodes [ 4, 5 ]. These groups are often used in applications [ 3 ].

Topological methods due to their nature can be used to cope with many problems where traditional methods fail. It has broad application in coverage problem in sensor network [ 2, 3 ], detecting breast cancer [ 6 ], computer vision [ 7 ], detection of periodicity in time series (behaviour classification) [ 8 ] etc.

Another issue addressed in this subject is the speed of algorithms. In paper [ 9 ] Dłotko et al. presented distributed algorithm for homology computation over a sensor network. Their technique involve reduction and coreduction of simplicial complexes.

Regardless of the algebraic approach, there are created algorithms that do not use calculations on homology groups to obtain some approximation of the shape of data [ 10 ]. Algorithm proposed in the last paper can have problems with more complicated structures, like the sphere, but it is shown that it does well with coverage in sensor networks.

Notice that apart from using Cˇ ech and Rips complexes there are other directions in computational topology, but we do not discuss them here [ 2 ]. 4

A Method to Build Complex using MapReduce

In this Section we present an efficient method to build Vietoris-Rips complexes. We make the following assumptions: – R is a Euclidean space. – A fixed number of dimensions n is given. – M Rn is the set of input data points. – We assume that each input data point has a unique identifier id (x). Those identifiers are items of an linearly ordered set. – For any i 2 f1; 2; : : : ; ng we define i : Rn ! R to be the projection to the i-th dimension. – In Rn we assume the Euclidean distance. For any two x; y 2 Rn let jx; yj denote the Euclidean distance of x and y. – Let us choose > 0. We are going to build the Vietoris-Rips complex R(M; )

The first step to build a Vietoris-Rips complex is to find all pair of points x; y 2 M such that jx; yj < . If the set M is large, computing distances between all pairs of points is too time consuming. If the number of dimensions n is equal to 1, we can assign all points to segments of the form hk ; (k + 1) ). Then, instead of computing distances for all pairs, we consider only those pairs of points that fall into the same segment or two neighbouring segments. As the result, each point is paired only with points from three segments. This can significantly accelerate the computation. Of course, in a one-dimensional space, a simpler method that encompasses sorting and sequential scanning will be faster. However, this method cannot be extended to more dimensions.

For a two-dimensional space we can apply a method similar to the former one described above. Instead of segments, we have then squares with sides of length . Eventually, each point will be paired for distance computation only with points from nine other squares. The three-dimensional case is similar, but we use cubes with side of length and each point is paired with points from 27 cubes.

Unfortunately, enormous growth of the number of dimensions prohibits direct usage of this method. In case of a 100-dimentional space, the “segments” will be 100-dimensional hypercubes with -side. Each point will be paired for distance computation with points from 3100 hypercubes. Obviously the na¨ive method will amount to be better for sure in this case.

However, we can apply the abovementioned three-dimensional method in multi-dimensional spaces indirectly. If we choose three dimensions i; j; k, we can consider projecting the whole space to them and further search for close pairs like in a three-dimensional space. Let us consider the projection = i j Rn ! R3. Obviously for each pair x; y 2 Rn the following inequality holds: k : j (x); (y)j3 jx; yj; 3 By j : : : j3 we denote the Euclidean distance in R .

We will assign all points to -sided cubes according to their projections onto the chosen dimensions. Then, it is enough to pair each point for distance computation only with points residing in the same cube or in 26 neighbouring cubes. If two points are not assigned to the same cube or neighbouring cubes, the distance of their projections onto the chosen dimensions is larger than . Therefore, so is their distance in Rn and their pair does not have to be considered.

The efficiency of this method notably depends on the choice of the three projected dimension. Chosen dimensions may reduce the number of distance computations required. Our first idea was to compare the projection onto each single dimension with the uniform distribution using the algorithm based on Kullback– Leibler divergence [ 11 ]. However, this approach required two scans of input data. Moreover, it is hard to distinguish more dense and more sparse distributions. The problem lies rather in sparseness (with respect to ) that uniformity. As the result, we invented the following quality measure of dimensions.

For a given dimension i and an integer t we define: (i; t) = jfx 2 M; t i(x) < (t + 1) gj The number (i; t) tells how many points will fall into the segment identified by t if we project to the dimension i. Note that if a point x satisfies the condition in the definition above, then t = b i(x)= c. bac is the largest integer not greater than a.

For a dimension i 2 f1; 2; : : : ; ng and an input data point x 2 M we define ti(x) to be the integer:

ti(x) = b i(x)= c:

Observe that the computation of all non-zero values of is possible using only a single scan of input data. Moreover, it is easily implementable in the MapReduce computation model. The mapper emits hi; ti(x)i for each input data point x 2 M and for each dimension i. The reducer simply counts the number of occurrences of pairs hi; ti.

Non-dispersion measure of a dimension i is the number

i = X t2Z

(i; t)2

Now, for each dimension we compute its non-dispersion measure by means of the MapReduce procedure presented above. Then we chose three dimensions with smallest non-dispersions.

Assume i; j; k to be the three dimensions with smallest non-dispersion. The algorithm to find all pairs of points closer that using these dimensions can be easily parallelized using the MapReduce computation model.

The mapper reads all input data points and for each point x emits 27 keyvalue pairs:

hht1; t2; t3i; xi where t1; t2; t3 satisfy the following conditions: t1 2 f ti(x) t2 2 f tj (x) t3 2 f tk(x) 1; ti(x); ti(x) + 1 g 1; tj (x); tj (x) + 1 g 1; tk(x); tk(x) + 1 g

A single reducer collects all pairs that came with a given key ht1; t2; t3i. If a pair x; y satisfies the conditions id (x) < id (y), t1 = ti(x), t2 = tj (x), t3 = tk(x), then the distance of points x; y is computed. The three dimension equalities assure that x lies in the very -sided cube ht1; t2; t3i. If jx; yj < , the pair hx; yi is added to the result.

The collection of all pairs hx; yi that satisfy jx; yj < is the most timeconsuming phase of the construction of the Vietoris-Rips complex. Those pairs are single-dimensional simplexes. The second phase computes triplets hx; y; zi such that the three points are pairwise closer than and id (x) < id (y) < id (z). Following phases consist in identifying longer and longer lists of points that satisfy analogous conditions. It can be done as in the Apriori algorithm to find frequent subsets, because if all points in a set are closer than that the same must hold for all its subsets. If is not too big, these phases are not timeconsuming. However, it is much simpler to build higher-dimensional simplexes using the following lemma.

Lemma 2 (Raising the dimension of simplexes). Let R be a Vietoris-Rips complex. A simplex a := hx1; :::; xn; xn+1; xn+2i belongs to R if and only if all the three following simplexes are also in R: b := hx1; :::; xn; xn+1i c := hx1; :::; xn; xn+2i d := hxn+1; xn+2i Proof. The “only if” part is straightforward. If a 2 R, then all its subsimplexes obviously belong to R. Thus, so do b; c; d.

The proof of the “if” part is based on the following property of Vietoris-Rips complexes. If a simplex a belongs to R, then all its one-dimensional subsimplexes also belong to R. Assume a one-dimensional simplex e := hxj ; xki a, where j; k 2 f1; : : : ; n + 2g. We will show that e 2 R.

We have to consider three cases. Firstly, if j = n + 1 ^ k = n + 2, then e = d 2 R. Secondly, if j < n + 1 ^ k = n + 2, then e c 2 R, thus e 2 R. Thirdly, if k < n + 2, then j < n + 1, thus e b 2 R and also e 2 R.

Since the choice of j and k is arbitrary, the “if” part has been proven. So has been the whole lemma. tu

As the result of the first phase of the algorithm, we have a collection pairs that contains all pairs hxi; xj i such that jxi; xjj <

^ id j < id k We recall that id i is the identifier of the point xi.

In the second phase of the algorithm, we build collections of simplexes of subsequent dimensions n that belong to the Vietoris-Rips complex R. The algorithm stops when for a given n we get the empty list of simplexes.

For a given n, we search for all b, c and d that satisfy the constraints of Lemma 2. The symbols introduced in Lemma 2 will prove to be handy again. Therefore, we perform the following steps: 1. We convert the collection of simplexes hx1; x2; : : : ; xn+1i of the dimension n to the collection of pairs hhx1; x2; : : : ; xni; xn+1i. 2. We compute a natural self-join of this set of pairs using the first coordinate. As the result we get the set of all triplets: hhx1; x2; : : : ; xni; xj; xki (a as in Lemma 2) such that there are pairs hhx1; x2; : : : ; xni; xji (b) and hhx1; x2; : : : ; xni; xki (c) in the previous step. 3. We leave (select) only those triples that satisfy id j < id k. 4. We equi-join the remaining triplets with the set of all pairs using the last two coordinates of the triplets. We do it to assure that the last two coordinates of triplets (d) form a one-dimensional simplex belonging to the complex. 5. We flatten remaining triplets to get complexes of the dimension n + 1 of the form hx1; x2; : : : ; xn; xj; xki.

Note, that this procedure contains only operations that are expressible in relational algebra (or calculus). Therefore, the whole algorithm is formulated in a highly abstract language that is known to be easy to optimise and execute in parallel. We exploited this possibility and implemented it in the Apache Spark for efficient execution on large computing infrastructures of commodity computers. An example Python code that implements the whole algorithm is shown below. komplex = [] komplex.append(points.map(lambda p: p[0])) komplex.append(pairs) # initial base = pairs pairsAsKey = pairs.map(lambda p: (p, 1)) noAnylizedObjects = pairs.count() #iterated while noAnylizedObjects > 0: base = base.map(lambda a:(a[:-1], a[-1])) base = base.join(base) \ .filter(lambda a: a[ 1 ][0] < a[ 1 ][ 1 ]) \ .map(lambda a: (a[ 1 ] , a[0])) \ .join(pairsAsKey) \ .map(lambda a:rigthTupleSum(a[ 1 ][0], a[0])) noAnylizedObjects = base.count() if noAnylizedObjects > 0:

komplex.append(base) 5

Conclusions

In this article we have presented an efficient method to build Vietoris-Rips complexes. Moreover, this method is easily parallelizable using the MapReduce computation model or Apache Spark. Further research will focus on algorithms that reduce constructed complexes. It will facilitate finding homology types of complexes. As a result it will also allow assessing of the correctness of the selection of the value.

References

1. Carlsson , G. , Ishkhanov , T. , De Silva , V. , Zomorodian , A. : On the local behavior of spaces of natural images . International journal of computer vision 76 ( 2008 ) 1 - 12

2. De Silva , V. , Ghrist , R.: Coverage in sensor networks via persistent homology . Algebraic & Geometric Topology 7 ( 2007 ) 339 - 358

3. Ghrist , R.: Barcodes: the persistent topology of data . Bulletin of the American Mathematical Society 45 ( 2008 ) 61 - 75

4. Edelsbrunner , H. , Letscher , D. , Zomorodian , A. : Topological persistence and simplification . Discrete and Computational Geometry 28 ( 2002 ) 511 - 533

5. Zomorodian , A. , Carlsson , G.: Computing persistent homology . Discrete & Computational Geometry 33 ( 2005 ) 249 - 274

6. Nicolau , M. , Levine , A.J. , Carlsson, G.: Topology based data analysis identifies a subgroup of breast cancers with a unique mutational profile and excellent survival . Proceedings of the National Academy of Sciences 108 ( 2011 ) 7265 - 7270

7. Freedman , D. , Chen , C. : Algebraic topology for computer vision . Computer Vision ( 2009 ) 239 - 268

8. Vejdemo-Johansson , M. , Pokorny , F.T. , Skraba , P. , Kragic , D. : Cohomological learning of periodic motion . Applicable Algebra in Engineering, Communication and Computing 26 ( 2015 ) 5 - 26

9. Dłotko , P. , Ghrist , R. , Juda , M. , Mrozek , M. : Distributed computation of coverage in sensor networks by homological methods . Applicable Algebra in Engineering, Communication and Computing 23 ( 2012 ) 29 - 58

10. Ćwiszewski , A. , Wiśniewski , P. : Coverage verification algorithm for sensor networks . In: Computer Applications for Bio-technology, Multimedia , and Ubiquitous City. Springer ( 2012 ) 397 - 405

11. Amari , S.I. , Cichocki , A. , Yang , H.H. , et al.: A new learning algorithm for blind signal separation . Advances in neural information processing systems ( 1996 ) 757 - 763