Topological Methods in Machine Learning: A Tutorial for Practitioners
Baris Coskunuzer Cüneyt Gürcan Akçora







coskunuz@utdallas.edu
0000-0001-7462-8819

  University of Texas at Dallas
  Richardson
  TX
  USA


cuney.akcora@ucf.edu
0000-0002-2882-6950

  University of Central Florida
  Orlando
  FL
  USA





§ ABSTRACT 
Topological Machine Learning (TML) is an emerging field that leverages techniques from algebraic topology to analyze complex data structures in ways that traditional machine learning methods may not capture. This tutorial provides a comprehensive introduction to two key TML techniques, persistent homology and the Mapper algorithm, with an emphasis on practical applications. Persistent homology captures multi-scale topological features such as clusters, loops, and voids, while the Mapper algorithm creates an interpretable graph summarizing high-dimensional data. To enhance accessibility, we adopt a data-centric approach, enabling readers to gain hands-on experience applying these techniques to relevant tasks. We provide step-by-step explanations, implementations, hands-on examples, and case studies to demonstrate how these tools can be applied to real-world problems. The goal is to equip researchers and practitioners with the knowledge and resources to incorporate TML into their work, revealing insights often hidden from conventional machine learning methods. The tutorial code is available at <https://github.com/cakcora/TopologyForML>.


 


<ccs2012>
   <concept>
       <concept_id>10002950.10003648.10003700</concept_id>
       <concept_desc>Mathematics of computing Algebraic topology</concept_desc>
       <concept_significance>500</concept_significance>
       </concept>
   <concept>
       <concept_id>10010147.10010257.10010293.10010294</concept_id>
       <concept_desc>Computing methodologies Learning paradigms</concept_desc>
       <concept_significance>500</concept_significance>
       </concept>
   <concept>
       <concept_id>10010147.10010257.10010282</concept_id>
       <concept_desc>Computing methodologies Machine learning algorithms</concept_desc>
       <concept_significance>500</concept_significance>
       </concept>
   <concept>
       <concept_id>10010147.10010257.10010293.10010319</concept_id>
       <concept_desc>Computing methodologies Machine learning applications</concept_desc>
       <concept_significance>500</concept_significance>
       </concept>
   <concept>
       <concept_id>10010405.10010497</concept_id>
       <concept_desc>Applied computing Life and medical sciences</concept_desc>
       <concept_significance>300</concept_significance>
       </concept>
   <concept>
       <concept_id>10010405.10010432</concept_id>
       <concept_desc>Applied computing Physical sciences and engineering</concept_desc>
       <concept_significance>300</concept_significance>
       </concept>
 </ccs2012>


[500]Mathematics of computing Algebraic topology
[500]Computing methodologies Learning paradigms
[500]Computing methodologies Machine learning algorithms
[500]Computing methodologies Machine learning applications
[300]Applied computing Life and medical sciences
[300]Applied computing Physical sciences and engineering





 


20 September 2024
[revised]12 March 20xx
[accepted]5 June 20xx

 
Topological Methods in Machine Learning: A Tutorial for Practitioners
    Cüneyt Gürcan Akçora
    Received XXXX; Accepted YYYY
=====================================================================







§ INTRODUCTION
 



As the complexity of datasets has grown in recent years, topological methods have emerged as powerful complements to state-of-the-art ML methods. The advent of ML has revolutionized the way we analyze and interpret complex data, yet there remain challenges in capturing the intrinsic topological structures inherent in such data. Traditional ML techniques, while powerful, often fall short in identifying and leveraging these structures, leading to the potential loss of valuable insights. Topological Machine Learning (TML) bridges this gap by integrating concepts from algebraic topology into ML workflows, enabling the discovery of patterns and features that are otherwise elusive. Despite this utility, much of the existing literature on topological methods in ML is highly technical, making it challenging for newcomers to grasp the direct connections to practical applications. 

In this tutorial, we introduce the fundamental concepts of topological methods to the machine learning (ML) community and a broader audience interested in integrating these novel approaches into their research. No prior knowledge of topology or ML is required. Our primary aim is to address this pressing need by providing a practical guide for non-experts looking to employ topological techniques in various ML contexts. To maintain accessibility, we will simplify the exposition and offer references to more detailed technical resources for those interested in further exploration.

In this paper, we teach two cornerstone techniques of TML: Persistent Homology and Mapper algorithm, and their effective utilization in ML. Persistent homology offers a robust, multi-scale analysis of topological features, allowing researchers to detect and quantify structures such as clusters, loops, and voids across different scales within the data. This capability is particularly beneficial for understanding the intricate relationships and hierarchies that may exist in complex datasets. From an ML perspective, PH offers a great feature extraction method for complex datasets, which was impossible with other methods. In this part, we will focus on this aspect of PH, giving hands-on instructions with illustrations on deriving effective topological vectors from complex data. On the other hand, the Mapper algorithm complements this by providing a visual and interpretable summary of high-dimensional data. By constructing a summary graph that mirrors the underlying topology of data, the Mapper algorithm facilitates the exploration and interpretation of data in an intuitive and informative way. This technique is instrumental in uncovering data's geometric and topological essence, making it accessible for practical analysis. 


Throughout the paper, we provide comprehensive explanations and step-by-step implementations of these techniques, supported by case studies spanning diverse applications such as cancer diagnosis, shape recognition, genotyping, and drug discovery. We aim to equip researchers and practitioners with the necessary knowledge and tools to integrate TML techniques into their studies, thereby unlocking new avenues for discovery and innovation. By demonstrating the practical utility of persistent homology and the Mapper algorithm, we highlight their potential to reveal insights that traditional methods may overlook, ultimately advancing the field of Machine Learning. 

We note that this paper is not intended to be a survey of recent advances in topological data analysis but rather a tutorial aimed at introducing the fundamentals of the topic to ML practitioners. For comprehensive surveys in topological data analysis, refer to <cit.>. For in-depth discussions of these topics, see the excellent textbooks on TDA and computational topology <cit.>.



 §.§ Roadmap for the Tutorial

  

We recommend reading the entire tutorial for a complete understanding, but readers may skip sections not pertinent to their needs. Here, we give a quick overview of the paper's structure.


In <Ref>, we provide the essential topological background needed to follow the concepts discussed throughout the rest of the paper. We aim to introduce these topological concepts, help build an intuition for TDA approaches, and demonstrate how they can be adapted and applied to various needs. For those unfamiliar with topology, we strongly encourage reading our introductory crash course in <Ref>. When discussing homology, we start with a non-technical, brief overview (<Ref>), which should suffice for following the rest of the paper. For readers interested in more technical details of homology computation, we offer a more in-depth explanation in <Ref>.

 

In <Ref>, we introduce Persistent Homology (PH) in three key sections: constructing filtrations (<Ref>), deriving persistence diagrams (PD) (<Ref>), and applying PDs to ML tasks, including vectorization (<Ref>) and neural networks (<Ref>). The filtration process is tailored to different data formats, as methods differ significantly based on the type of data—whether it’s point clouds (<Ref>), images (<Ref>), or networks (<Ref>). If you focus on a specific data type, you can skip sections on other formats. In each subsection, we also discuss the hyperparameter selection process. Lastly, <Ref> covers available software for PH. 

In <Ref>, we give a brief introduction to a niche subfield, Multiparameter Persistence (MP), an effective extension of PH. Again, we outline how to tailor the method to specific data formats for multifiltrations, i.e., point clouds (<Ref>), images (<Ref>), and networks (<Ref>). Next, we outline the state-of-the-art methods on how to integrate MP information in ML pipelines (<Ref>).

In <Ref>, we begin with a friendly introduction to the Mapper method (<Ref>), an effective TML tool for unsupervised learning. We then cover hyperparameter selection for Mapper in <Ref>, which is key for applications. Although the original Mapper algorithm is intended for point clouds, <Ref> explores recent advancements that extend its application to images and networks. Lastly, we provide an overview of available software libraries for Mapper in <Ref>.

In <Ref>, we summarize five real-life applications of these methods from published works, namely shape recognition for point clouds (<Ref>), anomaly forecasting for transaction networks (<Ref>), cancer detection from histopathological images (<Ref>), computer-aided drug discovery (<Ref>), and cancer genotyping from RNA sequencing (<Ref>).

Finally, in <Ref>, we outline potential future directions to advance TML methods, aiming to improve their practical use in ML and discuss strategies for broadening the application of topological methods to new and emerging fields.


 



§ BACKGROUND
 

In this section, we provide some background that will be used later. We first introduce several mathematical concepts that will be used in the second part, where we describe homology, which is essential for introducing the methods in subsequent sections.




 §.§ A Crash Course on Topology
 

Topology, a core discipline in mathematics, studies the properties of shapes and spaces that remain unchanged under continuous transformations such as stretching, crumpling, and bending, but not tearing or gluing. In machine learning, topology provides a powerful framework for examining complex data structures flexibly and intuitively. To provide clarity, we begin by defining several key terms essential for following the paper.



  §.§.§ Topological Space
 From a mathematical perspective, topology refers to the structure of a set. For instance, consider the set of real numbers between 0 and 1, i.e., = {x∈| 0 ≤ x ≤ 1}. Many readers might immediately consider the closed interval [0,1], whose shape resembles a stick. However,  is merely a set with no defined structure yet. This means we do not know which points in  are close to or far from others, as there is no concept of neighborhoods.

For example, if we define a distance (metric) on  such as d_1(x, y) = 1 for x ≠ y and d_1(x, y) = 0 for x = y, the "shape" of the set  would be entirely different, resembling an infinite number of points dispersed in a very high-dimensional space. Conversely, if we use a metric like d_2(x, y) = |x - y| on , we retrieve our familiar closed interval [0,1], a stick of length 1. Notice that due to the different metrics used, the neighborhood structures and shapes of (, d_1) and (, d_2) are completely different. This brings us to the first principle in Topology: the topology of a set is defined as the complete neighborhood information on the set.  

The principle is not without exceptions: a significant research area called point-set topology studying topological spaces that do not necessarily have a metric <cit.>. However, these topologies are beyond our scope. Thus, throughout the paper, a topological space refers to a set  (or dataset) equipped with a distance d(·,·), forming what is known as a metric space (,d). In other words, a set qualifies as a topological space if we can unambiguously identify the neighbors of its points. 



Although we might not explicitly reference distance, most datasets inherently possess some form of a metric for detailed analysis. For instance, point clouds use the metric of the space they are embedded in, providing neighboring information. In graphs, adjacent nodes are naturally considered neighbors. Similarly, in images, neighboring pixels are deemed neighbors. This paper treats these examples as topological spaces with their natural topologies. However, specific datasets, such as RNA-sequencing data, lack an inherent metric, requiring users to define a metric to determine which data points are considered close and distant, depending on the context.



  
*Simplicial Complexes. One special family of topological spaces used throughout our paper will be simplicial complexes. A simplicial complex is a mathematical structure used in topology and combinatorics to study the properties of shapes and spaces in a discrete manner. It consists of vertices, edges, and higher-dimensional simplices such as triangles, tetrahedra, and their higher-dimensional counterparts, which are glued together in a specific way to form a coherent whole. 



Each simplex is a generalization of the concept of a triangle to higher dimensions. By this, we mean that the properties and structure of a triangle (such as having vertices, edges, and faces) are extended into higher dimensions, even though the shapes look different as we go up in dimensions. In particular, a k-simplex is defined as the convex hull of affinely independent k+1 points in ^k. For example, 2-simplex is a triangle (with its inside filled), e.g., the convex hull of (smallest convex set containing) three points {(0,0), (1,0),(0,1)} in ^2. Similarly, 3-simplex is a tetrahedron (with inside filled), e.g., the convex hull of {(0,0,0), (1,0,0),(0,1,0),(0,0,1)} in ^3. A union of simplices is called a simplicial complex if any two simplices in the complex either do not intersect or intersect in a complete subsimplex (See <Ref>).


r3in


    
    < g r a p h i c s >

        Simplicial Complexes. Among the complexes, only b and f fail to be simplicial complexes, as their simplices do not intersect at complete subsimplices. All others are valid simplicial complexes. 
    





  
*Dimension in Topology One of the most confusing concepts for non-experts in topology is the notion of dimension. To study topological spaces more effectively, we focus on a specific family of spaces that meet certain regularity conditions, known as manifolds. A k-dimensional manifold is a topological space that locally resembles a ball in ^k, i.e., a small neighborhood of any point looks like (homeomorphic) a ball in ^k. 
Importantly, we do not concern ourselves with the ambient (i.e., surrounding) space in which our manifold resides; we focus solely on its intrinsic properties, e.g., the material the manifold is made of. For example, although a circle is often visualized in ^2 as a two-dimensional object, from a topological perspective, it is actually one-dimensional. This is because, at any point on the circle, its local neighborhood resembles a line segment, i.e., it is made of 1-dimensional material. Consider an ant walking on a circle; it has two options: moving forward or backward. While we might need two coordinates to describe the ant's position on the circle mathematically, from the ant's perspective, it only experiences a continuous path where it can move in either direction. To the ant, the circle feels like an endless line. Similarly, a sphere is a 2-dimensional manifold (surface) even though it is often visualized in ^3. This is because a small neighborhood of any point on the sphere resembles a piece of in ^2. Notice that with this definition, the dimension of the ambient space becomes irrelevant; only the topological properties define the dimension of a manifold, i.e., a loop (circle) in ^2  or ^3 is called 1-dimensional.  Similarly, a torus (a hollow donut) and a genus-2 surface (a surface with two holes) are examples of 2-dimensional manifolds. If this discussion leads you to wonder about manifolds with more than two dimensions, we have some discouraging news: they do exist, but visualizing them is not intuitive.





  
*Boundary. In topology, understanding the boundary of a manifold is crucial to grasping its structure. A k-dimensional manifold M is a space where each point has a neighborhood that resembles an open subset of ^k. The boundary ∂ M of M comprises points where these neighborhoods resemble an open subset of the half-space ^k_+ = { x ∈^k | x_k ≥ 0 }. This implies that the local structure is similar to part of ^k_+ near each boundary point.

To illustrate, consider a sphere, a 2-dimensional surface without any boundaries. An ant walking on the sphere can indefinitely move in any direction without encountering an edge. In contrast, the unit disk  in ℝ^2 is a 2-dimensional surface with a boundary, specifically the unit circle .  An ant residing in  could not continue its journey once it reaches the boundary , or the "border" of . This is denoted as ∂ =. Similarly, the unit ball  in ^3   (the solid ball) is a 3-dimensional manifold with a boundary, and its boundary is the unit sphere , denoted by ∂ =. In this context, we say that the sphere  bounds the ball .

If you're puzzled by the statement that  a sphere has no boundary, whereas the solid ball does, consider it from the perspective of "its inhabitants." An ant living on the sphere's surface can move freely in any direction without ever encountering a limit (no border), which is why we say the sphere has no boundary. However, a mouse living inside the ball will eventually hit the boundary (border)—the sphere's surface—beyond which it cannot go. That limiting surface is what we call the solid ball's boundary. 


A key principle in topology is that a boundary's boundary is always empty. Mathematically, this is expressed as ∂ (∂) = ∅. Furthermore, the boundary of a k-dimensional manifold is generally a (k-1)-dimensional manifold with no boundary. For example, the boundary of a disk  (a 2-dimensional manifold with a boundary) is a circle  (a 1-dimensional manifold without a boundary). This observation will soon be important when discussing homology. 




  §.§.§ Topological Equivalence
 Topology primarily focuses on the global properties of shapes rather than their local characteristics. For example, topology is concerned with whether a space is connected or has holes without regard to the size of the object or the holes. 

Topological equivalence can be intuitively understood as follows: two shapes,  and  are topologically equivalent if one can be continuously deformed into the other. Imagine  and  are made of Play-Doh. If you can reshape  into  without tearing, gluing, or collapsing any part of the shape, then they are considered to be continuously deformable into each other.  In this context, "collapsing" refers to reducing a part of the shape to a lower dimension, such as squishing a surface or line down to a point or flattening a 3-dimensional object into a 2-dimensional plane.

In mathematical terms, this concept corresponds to a homeomorphism (the prefix "homeo-" comes from the Greek word "homoios," which means "similar" or "like"). In particular, such a deformation from  to  represents a bijective map, which means it creates a one-to-one correspondence between elements of  and , ensuring that each element in  is paired with a unique element in  and vice versa. The map φ:→ keeps track of how each point x ∈ becomes a point y ∈ after the deformation, i.e., φ() = y. The condition of no tearing ensures that this map is continuous as nearby points in  must map to nearby points in . Gluing is the opposite of tearing, so we also require that the inverse map φ^-1 is continuous. In summary, φ:→ is called a homeomorphism if φ is a continuous bijection with a continuous inverse.

We previously noted that size does not affect the topological structure. For example, let =, the set of all real numbers, and = (-1, 1), the open interval containing all real numbers between -1 and 1, excluding the endpoints. These two spaces are topologically equivalent via the homomorphism φ:→(-1,1) with  φ(x) = x/1 + |x|. However,  has an infinite diameter, while (-1,1) has a diameter of only 2. Hence, this is a good example of how the space size does not matter in topology. On the other hand, consider = (0,1) ∪ (1,2) and = (0,2). While both spaces are similar as sets, they are not topologically equivalent because  consists of two separate connected components, whereas  is a single connected piece. A continuous deformation (homeomorphism) must preserve the number of components.

 


  
*Homotopy. Another important concept in topology is homotopy, which can be seen as a flexible deformation of one space into another. Unlike homeomorphism, homotopy allows for collapsing but not tearing or gluing. We call two spaces homotopic to each other if such a deformation exists from one to another. For instance, a disk and a point are homotopic because the entire disk can be collapsed into a point by pushing inwards from all directions towards the center. Similarly, a punctured disk (𝐃^2-{(0,0)}) and a circle are homotopic, as a punctured disk can be continuously deformed towards its boundary, starting from the puncture point (0,0). Homotopy provides a powerful tool for classifying topological spaces and understanding their fundamental properties, as most topological invariants are homotopy invariants, meaning they yield the same output for two homotopic spaces. For example, the connectivity and the count of holes or cavities do not change under homotopy.





  §.§.§ Topological Invariant
  If  two spaces  and  are topologically equivalent, showing that there is a homeomorphism φ: → would be sufficient. However, if they are topologically different, one must demonstrate that no such map can exist, which is highly challenging. Mathematicians use invariants to show such inequivalence. An invariant refers to a property or feature of an object that remains unchanged under certain transformations. A topological invariant is an invariant that is preserved under homeomorphism. A topological invariant can be a number (number of components, Euler characteristics, Betti numbers), a mathematical object (homology groups, fundamental group, cohomology ring), or a mathematical property (compactness, connectedness). 



r3in


    [Sphere]
    < g r a p h i c s >

    [Cube]
    < g r a p h i c s >

    [Torus]
    < g r a p h i c s >

    The sphere and the cube are topologically equivalent, whereas the torus is different from both.
    
    



In the example above, = (0,1) ∪ (1,2) versus = (0,2), we used the number of components as a topological invariant to show they are not equivalent. However, more subtle examples require more sophisticated invariants. For instance, a cube and a sphere are topologically equivalent because one can deform a cube into a sphere without tearing or gluing. In contrast, comparing a sphere and a torus (the surface of a donut as shown in  <ref>) is more complex. They are not topologically equivalent because the torus has "holes." By holes, we mean the loops on the surface, which cannot be shrunk down to a point without leaving the surface. From this perspective, a sphere has no hole, however, on a torus, some loops (meridian and longitude) can't be shrunk down to a point in the torus, representing the "holes" in torus (see <Ref>). Despite this difference, proving that a sphere and a torus are not homeomorphic is challenging. Therefore, subtle topological invariants that detect these holes are crucial for distinguishing such spaces. 



  
*Euler Characteristics One important and well-known example of such topological invariants is the Euler characteristic. If  is a simplicial complex, we define the Euler characteristic as the alternating sum of the count of k-simplices in , i.e., χ() = ∑_k (-1)^k n_k, where n_k is the count of k-simplices in . It turns out the Euler characteristic is a topological invariant. i.e., if  and  are homeomorphic, then χ() = χ(). The Euler characteristic is, in fact, a homotopy invariant as well (see <cit.>). For example, a sphere  is topologically equivalent to the surface of a hollow tetrahedron , which has 4 vertices (n_0), 6 edges (n_1), and 4 triangles (n_2). Therefore, we have χ() = χ() = n_0 - n_1 + n_2 = 4 - 6 + 4 = 2. Similarly, if one computes the Euler characteristic of a torus  through a simplicial complex, we find that χ() = 0. Since χ() ≠χ(), it follows that a sphere is not homeomorphic to a torus.

 



  §.§.§ Geometry vs. Topology
 Before proceeding, let's clarify a common confusion between the terms topology and geometry. Both are branches of mathematics, but they focus on different aspects of metric spaces.

Geometry focuses on the local study of shapes, sizes, and properties of space, as well as how these spaces embed (fit into) higher-dimensional spaces.  For example, a 2-dimensional sphere  can be embedded into ^3 (or ^10) in various shapes and sizes.
Geometry focuses on precise measurements (e.g., angles, distances, curvature) and the relationships between points, lines, surfaces, and solids. It explores how shapes bend and twist, providing tools to understand complex surfaces and spaces through concepts such as curvature and geodesics (i.e., shortest paths).

In contrast, topology examines the global properties of space that remain unchanged under continuous deformations. It focuses on qualitative aspects like connectedness and continuity rather than precise measurements. For instance, topology considers a cube and a sphere as topologically equivalent because these shapes can be transformed into one another through continuous deformation despite their geometric differences—one being flat and the other curved. For instance, no matter how we embed a 2-sphere into ℝ^10, its topological properties remain unchanged. However, the geometry varies significantly: a round sphere with a radius of 1 differs greatly from one with a radius of 100. Moreover, an irregular sphere would be geometrically distinct from both. 

In data science, geometry usually refers to the local shape characteristics of a dataset, such as distances, curvature, and angles, whereas topology pertains to global characteristics, such as connectedness and the number of holes/cavities. For example, dimension reduction methods like PCA <cit.> and UMAP <cit.> are considered geometric methods as they heavily depend on the distances and how the dataset sits in the high dimensional space. In contrast, methods that count the number of components or holes in the space are called topological methods.






 §.§ Homology
 

To distinguish topological spaces, the most common method is to use topological invariants such as Euler characteristics (<Ref>), the fundamental group <cit.>, or homology. Among these, homology is the most versatile and robust invariant that applies to a wide range of spaces such as surfaces (e.g., spheres, tori), simplicial complexes (e.g., triangulated shapes), and manifolds (e.g., higher-dimensional analogs of curves and surfaces). 


There are various ways to compute homology (cellular <cit.>, simplicial, Morse <cit.>), where the outputs are the same, but the computation methods applied are different. To utilize computational tools more effectively,  it's more efficient to use discrete representations of topological spaces, like simplicial complexes. Simplicial homology is particularly suited for TDA because it deals with simplicial complexes, which are sets of vertices, edges, triangles, and their higher-dimensional counterparts. Considering these aspects, TDA mostly employs simplicial homology to capture the topological patterns in data, 
although Morse or cellular homology are used in specific applications within TDA, such as cellular in image classification <cit.>.



 Here, we provide a brief overview (TLDR) of homology, followed by a formal yet accessible introduction. For a detailed, friendly introduction to simplicial homology, refer to <cit.>, or for an in-depth study, see <cit.>.

 



  §.§.§ TLDR
 
Homology is a fundamental invariant in topology that captures information about a space's structure by examining its holes/cavities of various dimensions. The focus on holes/cavities may surprise the reader, but holes are preferred because they are fundamental features that significantly influence the structure and properties of a space, as we outline below. To simplify the exposition, we use the concept of a k-hole in a topological space . Although this term slightly abuses notation, it refers to a k-dimensional submanifold  in  that cannot be continuously deformed into a point within the space. In reality, this k-hole corresponds to a (k+1)-dimensional "cavity" Ω in . Since this cavity represents a "missing" region within the space, we describe it by using its boundary ∂Ω =, which is non-contractible in . The k^th homology group _k() captures these k-holes, or k-dimensional manifolds in , that do not bound any (k+1)-dimensional region in the space.

 

  * k-holes are topological invariants, meaning they remain unchanged under continuous deformations, such as stretching or bending, that do not involve tearing or gluing. And homology can detect them.

  * 0-holes represent the connected components of a space. The dimension (or rank) of _0() corresponds to the number of these components in .

  * 1-holes correspond to loops in the space that cannot be contracted to a point without leaving the space. The dimension of _1() represents the number of such loops (1-holes) in the space.

  * 2-holes correspond to cavities within the space, which can be considered hollow regions enclosed by surfaces (e.g., the interior of a sphere or torus). The dimension of _2() represents the number of such cavities (2-holes) in the space.
 



r2.8in


    
    < g r a p h i c s >

        Toy examples for homology. We present the ranks of the homology groups of various topological spaces, i.e., for _0()=^k we write only _0=k for simplicity, representing the count of the k-dimensional holes in .      





To make the concept of homology groups more accessible, we will offer a simplified and visual explanation. 
We denote homology groups by _k(), where  represents the space under consideration, and k indicates the dimension being analyzed. Independent k-holes are generators for the homology group _k(). Therefore, the rank of _k() indicates the number of k-dimensional holes present in the space . These numbers are also called Betti numbers {β_k()}. e.g., if _2()=^3, then we say rank(_2())=3 , or β_2()=3, meaning  has three 2-holes.  Below, we will give some examples for dimensions 0,1 and 2 (See <Ref>).




 
    
∙ _0(): 0-dimensional homology represents the connected components in . If  has three components, then we say _0()=^3 (or _2^3 if used _2-coefficients), and rank(_0())=β_0()=3.
    

    

    
∙ _1(): 1-dimensional homology computes the non-contractible loops (holes) in .
    For example, a sphere  has no nontrivial loops, as any loop in the sphere bounds a disk in the sphere, i.e., β_1()=0.
    A torus  (hollow donut), on the other hand, has meridian and longitude circles (see Figure <ref>), both being non-contractible loops, making the total β_1()=2. On the other hand, a solid donut  would have only the longitude circle as noncontractible, β_1()=1.


    ∙ _2(): 2-dimensional homology group captures two-dimensional cavities (voids) in the space . To count these cavities, we utilize surfaces in the space which does not bound a 3-dimensional domain in the space. For example, in the unit ball  in ^3 (solid ball), any closed surface bounds a 3-dimensional domain within , and there are no cavities within it. Hence, we say rank(_2())=β_2()=0. Similarly, if one removes two disjoint smaller balls from , say ', we would have two different cavities, which can be represented by different spheres in , enclosing these balls. Then, we say β_2(')=2. In <Ref>, the sphere, the torus, or the genus-2 surface have one 2-dimesional cavities represented by themselves.
    
  

We conclude this section with an interesting fact. While we defined the Euler characteristic for simplicial complexes as the alternating sum of the number of k-simplices, there is an alternative way using Betti numbers. In particular, if  is a topological space, the Euler characteristic is given by χ()=∑_k(-1)^kβ_k().
 



  §.§.§ Computation of Homology
 
We now turn to a formal introduction. Homology is a mathematical operation whose inputs are topological spaces and outputs are groups. In particular, for a given space , _k() represents the k^th homology group of , summarizing the k-dimensional non-collapsible submanifolds in , each representing different (k+1)-dimensional "cavity" of . We will call these "k-holes" by abusing notation. 

The concept of homology groups stems from the idea that a (k+1)-dimensional hole or cavity in  is detected by the presence of its k-dimensional boundary, which cannot be continuously contracted within . While this method might seem indirect—using boundaries to infer the existence of cavities—the difficulty arises from the need to identify what is missing in . Here, a fundamental principle of topology comes into play: the boundary of a boundary is always empty. If Ω is a (k+1)-dimensional domain in  with boundary = ∂Ω, then the boundary of  is empty. i.e., ∂(∂Ω) = ∂=∅ Hence, to identify cavities, we first find all k-dimensional submanifolds with no boundary in  (k-cycles). Then, by eliminating those that bound a domain in  (k-boundaries), we are left with the true cavities. This process forms the core idea behind homology computation.


To keep the exposition focused, we will describe only simplicial homology with ℤ_2 coefficients, the most common version used in TDA. For other homology settings, refer to <cit.>. In particular, simplicial homology involves representing a given space 𝒳 as a simplicial complex (a collection of simplices) and performing computations on these simplices. This approach allows us to discretize the problem by focusing on the k-dimensional "building blocks" of the topological space, such as vertices for 0-dimension and edges for 1-dimension. By considering k-submanifolds  (or k-subcomplexes) as unions of k-simplices, we can identify the special ones that correspond to true cavities in 𝒳 by using computational tools. To formalize this concept, we now introduce the relevant mathematical notions.




  
*i. Representation of k-simplices In a simplicial complex , we describe k-simplices by listing their vertices. For example, a 1-simplex (an edge) e with endpoints v_2 and v_4 is denoted as e = [v_2, v_4]. Similarly, a 2-simplex (a triangle) τ with vertices v_1, v_5, and v_7 is written as τ = [v_1, v_5, v_7]. Since we are using ℤ_2 coefficients, the order of the vertices does not matter. However, in other versions, such as with ℤ-coefficients, the order of the vertices would be significant. In general, a k-simplex Δ is represented by its k+1 vertices, i.e., Δ = [v_i_0, v_i_1, …, v_i_k]. We will call a union of k-simplices a k-subcomplex of . 



  
*ii. k-chains _k() To describe all k-dimensional subcomplexes (which correspond to all k-submanifolds, with or without boundary) within a simplicial complex , we define a group _k(), known as the group of k-chains. Recall that any union of k-simplices forms a k-subcomplex. For instance, if the simplicial complex  consists of three edges {e_1, e_2, e_3}, then all possible 1-subcomplexes are {e_1, e_2, e_3, e_1 ∪ e_2, e_1 ∪ e_3, e_2 ∪ e_3, e_1 ∪ e_2 ∪ e_3}. We represent this collection using group elements, where each 1-simplex acts as a generator. In particular, the union e_1 ∪ e_3 is represented by σ_1 = (1,0,1), while e_1 ∪ e_2 is represented by σ_2 = (1,1,0). Hence, we obtain the group _1() = _2^3, where 0 element corresponds to the empty set. The group operation is addition, and since 1+1=0 in _2, we have σ_1 + σ_2 = (1,0,1) + (1,1,0) = (0,1,1), corresponding to the union e_2 ∪ e_3. Similarly, if  contains m k-simplices {Δ_1, Δ_2, …, Δ_m}, then _k() = _2^m, where a group element like (1,0,1,…,1) represents the k-subcomplex Δ_1 ∪Δ_3 ∪Δ_m within .





  
*iii. Boundary operator ∂_k To identify k-holes in a simplicial complex, we must first identify k-subcomplexes that have no boundary. This is accomplished by defining a boundary operator ∂_k: 𝒞_k() →𝒞_k-1(), which maps k-chains to their (k-1)-dimensional boundaries. In particular, ∂_k maps each k-simplex to its boundary. For example, if we have a 1-simplex (edge) e = [v_0, v_1], then ∂_1[e] = v_0 + v_1 ∈_0(), representing the boundary of the edge as the union of its end vertices {v_0}∪{v_1}. Similarly, for a 2-simplex (triangle) τ = [v_0, v_1, v_2], the boundary is given by ∂_2τ = [v_0, v_1] + [v_1, v_2] + [v_2, v_0] ∈_1(), representing the boundary of the triangle as the union [v_0, v_1] ∪ [v_1, v_2] ∪ [v_2, v_0]. 
To define ∂_k for general k-chains, we sum the boundaries of each k-simplex within the chain. For instance, if σ = Δ_1 + Δ_3 + Δ_7 is a k-chain, then ∂_kσ = ∂_kΔ_1 + ∂_kΔ_3 + ∂_kΔ_7. 
A k-chain σ is said to have no boundary if ∂_kσ = 0. In other words, a k-subcomplex with no boundary in  must map to 0 under the boundary operator.

The boundary operator ∂_k: 𝒞_k() →𝒞_k-1() is a linear operator and can be represented as a matrix. For example, if 𝒞_k() = ℤ_2^n and 𝒞_k-1() = ℤ_2^m, then ∂_k can be written as an m × n matrix 𝐀, with each column 𝐀_i corresponds to ∂_k Δ_i, where Δ_i is a k-simplex in . In <Ref>, we give an explicit example of a matrix representation of a boundary operator the simplicial complex in <Ref>.

r0.2


    
    < g r a p h i c s >

        Toy example for homology.       





  
*iv. k-cycles _k() We define a special subgroup _k(), named k-cycles, within _k() for k-subcomplexes that have no boundary. As previously mentioned, a k-subcomplex with no boundary in  must map to zero under the boundary operator. Hence, we define the subgroup _k() =  ∂_k, which includes all k-chains that map to zero.  Recall that 1-chains with no boundary correspond to loops, while 2-chains with no boundary correspond to closed surfaces, like a sphere or torus.

For example, consider a square-shaped simplicial complex  (<Ref>) with four vertices {v_0,v_1,v_2,v_3} and four edges {e_1,e_2,e_3,e_4}  where e_1=[v_0,v_1],e_2=[v_1,v_2],e_3=[v_2,v_3] and e_4=[v_3,v_0]. In this case, _1()=_2^4 and _0()=_2^4. Suppose σ_1=e_1+e_2; then ∂_1σ_1= (v_0+v_1)+(v_1+v_2)=v_0+v_2, meaning that σ_1 represents a 1-subcomplex with a boundary. However, if σ_2=e_1+e_2+e_3+e_4, then ∂_1σ_2= (v_0+v_1)+(v_1+v_2)+(v_2+v_3)+(v_3+v_0)=2(v_0+v_1+v_2+v_3)=0, indicating that σ_2 represents a 1-subcomplex with no boundary. Notice that σ_2 corresponds to a complete loop in . See <Ref> for more details.

r1.3in


∂_1 = 
[             ∂ e_1  ∂ e_2  ∂ e_3  ∂ e_4; [ 1 0 0 1; 1 1 0 0; 0 1 1 0; 0 0 1 1 ] ][    ; v_0; v_1; v_2; v_3 ]
∂_1:_1()→_0() is represented as a 4× 4 binary matrix. Columns represent the edges in _1(), and rows correspond to the vertices in _0(). For example,  ∂ e_3= v_2+v_3 can be read from the third column 





We can also describe the boundary map ∂_1: _1() →_0() as a 4 × 4 matrix, as shown in <Ref>. It is clear that the only element in _1() that maps to (0,0,0,0) ∈_0() is (1,1,1,1), which corresponds to the loop σ_2. This means  has only one 1-cycle, which is σ_2. Then, _k() is the subgroup of _k() generated by the single element (1,1,1,1). If you are unfamiliar with group theory, you can think of _1() as a four-dimensional vector space, where _1() is a one-dimensional subspace generated by the vector (1,1,1,1).




  
*v. k-boundaries _k() While we identified all k-subcomplexes with no boundary, none correspond to true cavities. We must eliminate the ones that bounds a domain in . To find them, we use again the boundary operator. As _k+1() represents all (k+1)-subcomplexes in , then the image of ∂_k+1, i.e., _k() = ∂_k+1𝒞_k+1() ⊂𝒞_k(), represents the ones which bounds a (k+1)-domain in . We call an element in _k() a k-boundary. Recall that ∂_k(∂_k+1σ)=0 from earlier discussion. This means for any k-boundary φ=∂_k+1σ in _k(), ∂_kφ=0. Therefore, we have _k()⊂_k(). For example, if  is a simplicial complex formed by only one 2-simplex τ with vertices {v_0,v_1,v_2}, then _1() would be a subgroup in _1(), generated by ∂_2 τ=σ=(1,1,1) while _1() would be the same group, i.e., _1()=_1(). This means there is only one loop σ in  but it does not represent a hole in  as it is filled by τ, i.e., σ=∂τ.



  
* vi. Homology group _k() Now we are ready to define homology. Notice that with the boundary operator, we obtained the following sequence of groups and maps. We can consider these as vector spaces and linear maps. For each k, k-chains _k() represent all k-subcomplexes in , k-cycles _k()⊂_k() represent k-subcomplexes with no boundary, and finally, _k()⊂_k() represent the k-subcomplexes which bounds a (k+1)-domain in .


    ⋯[r, "∂_k+2"]    _k+1() [r, "∂_k+1"]    _k() [r, "∂_k"]    _k-1() [r, "∂_k-1"]    ⋯[r, "∂_1"]    _0()[r]     0



 
Now, to identify k-holes (true cavities), we consider all k-cycles in  (_k()), which potentially represent a cavity of , and, among them, eliminate all k-boundaries in  (_k()), as they represent the "fake" cavities. Hence, we define k^th homology group as the quotient group 
    _k() = _k()/_k() = ker(∂_k)/image(∂_k+1)

In terms of the sequence above, this quotient _k() = _k()/_k()  effectively counts the k-dimensional cycles that correspond to actual holes or cavities in , not those that are merely boundaries of higher-dimensional regions. From a computational perspective, with this formulation, we only need to compute the kernels and images of a sequence of linear maps {∂_k} (binary matrices) to compute homology. 


r0.2


    
    < g r a p h i c s >

        Toy example <ref> for homology.     
        


To clarify these notions, we give two toy examples for explicit computation of homology.

  Consider the example in <Ref> where  is a square-shaped simplicial complex with four vertices {v_0, v_1, v_2, v_3} and four edges. For k ≥ 2, _k() = {0} since there are no k-simplices in . Both _1() and _0() are isomorphic to _2^4, as discussed earlier (k-cycles above). The group _1() has only one generator, the sum of all edges. Since _2() = {0}, we have _1() = {0} because _1() = ∂_2 _2(). Therefore, _1() = _1()/_1() = _2/{0} = _2, meaning that _1() has rank 1, corresponding to the entire square loop. The boundary map ∂_1 is given in <Ref>.


For _0(), we need to determine _0() and _0(). Since ∂_0 sends everything to 0, _0() = _0() = _2^4. The group _0() = ∂_1 _1() is generated by the boundaries of the edges: (v_0 + v_1), (v_1 + v_2), (v_2 + v_3), (v_3 + v_0). However, these boundaries are not linearly independent because v_3 + v_0 is the sum of the other three boundaries. Thus, _0() ≃_2^3. Consequently, _0() = _0()/_0() = _2^4/_2^3 = _2. Recall that the rank of _0() represents the number of connected components and the fact that rank(_0()) = 1 confirms that there is only one connected component.  



r3.2in


∂_2 = 
[                                    ∂τ_1 ∂τ_2  ∂τ_3  ∂τ_4; [ 1 0 1 0; 1 1 0 0; 0 1 0 1; 0 0 1 1; 1 0 0 1; 0 1 1 0 ] ][    ; e_1; e_2; e_3; e_4; e_5; e_6 ]  
∂_1 = 
[                  ∂ e_1  ∂ e_2  ∂ e_3  ∂ e_4  ∂ e_5  ∂ e_6; [ 1 0 0 1 1 0; 1 1 0 0 0 1; 0 1 1 0 1 0; 0 0 1 1 0 1;   ] ][    ; v_0; v_1; v_2; v_3 ]
Boundary maps for Ex. <ref>. The boundary map ∂_2:_2()→_1() is represented as a 6× 4 binary matrix (left). The top of each column corresponds to the 2-simplex in _2() whose image is represented by that column, while each row's corresponding edge is given next to the column. For example, the boundary of 2-simplex τ_2, ∂τ_2= e_2+e_3+e_6 by reading the second column in ∂_2 matrix. The boundary map ∂_1:_1()→_0() is similar (right). 



Second example is in one dimension higher.

  Consider the hollow tetrahedron  of <Ref>, composed of four triangular faces: τ_1 = [v_0, v_1, v_2], τ_2 = [v_1, v_2, v_3], τ_3 = [v_0, v_1, v_3], and τ_4 = [v_0, v_2, v_3]. The complex  contains four 2-simplices (triangles), six 1-simplices (edges), and four 0-simplices (vertices). Let e_1=[v_0,v_1],e_2=[v_1,v_2], e_3=[v_2,v_3],e_4=[v_3,v_0], e_5=[v_0,v_2] and e_6=[v_1,v_3]. Therefore, we have _2() = _2^4, _1() = _2^6, and _0() = _2^4. 



A straightforward computation shows that the kernel of ∂_2, denoted _2(), has a single generator: the sum of all the triangles, τ_1 + τ_2 + τ_3 + τ_4. Since there is no 3-simplex in , we find that _2() = {0}, leading to _2() = _2. This result aligns with the fact that rank(_2()) = 1, indicating the presence of one cavity in , consistent with  being a hollow tetrahedron. We 
Furthermore, calculating _1() = _2^3 and _1() = _2 shows that these cancel out, yielding _1() = {0}. In other words, any loop in  is filled, so there is no 1-hole in .
Similarly, following the same reasoning as in the example above, we obtain _0() = _2 as expected.  


















































 


 

























§ PERSISTENT HOMOLOGY
 
In this section, we introduce Persistent Homology (PH), a foundational technique that played a pivotal role in the emergence of TDA, as developed by Carlsson, Edelsbrunner, Zomorodian, and others in the early 2000s <cit.>. PH captures the underlying shape patterns within complex data sets by studying the evolution of topological features across multiple scales.


 

r3in


    
    < g r a p h i c s >

    PH Pipeline.  From data acquisition and filtration complex construction to generating persistence diagrams using software libraries. The final step highlights methods for integrating persistence diagrams into downstream ML tasks. 
    


PH first constructs a nested sequence of simplicial complexes, known as filtration, and tracks the birth and death of features, such as connected components, loops, and voids, in this sequence. The resulting multi-scale representation highlights significant features while filtering out noise, making PH a valuable tool in various fields, including medical imaging <cit.>, biomedicine <cit.>, time series analysis <cit.>, material science <cit.>, geography <cit.>, shape analysis <cit.> and finance <cit.>.

In the following, we explain PH in an accessible way, focusing on its key aspects relevant to ML applications. The main idea behind PH is to capture the hidden shape patterns in the data by using algebraic topology tools. PH achieves this by keeping track of the evolution of the topological features (k-holes, components, loops, and cavities) created in the data while looking at it in different resolutions.

In simple terms, PH can be summarized as a three-step procedure as follows.


    
  * Filtration: Generate a  nested sequence of simplicial complexes derived from the data.
    
  * Persistence Diagrams: Record the evolution of topological features across this sequence.
    
  * ML Integration: Transform the persistence diagrams into vectors for efficient use in ML models.



 

We provide details of these steps in the following sections. Although the second and third steps are conceptually similar across most settings, the first step, constructing filtrations, varies significantly depending on the data type, i.e., point clouds, images, and networks. See <Ref> for a visual summary of PH pipeline.


 

 



 §.§ Constructing Filtrations
 

In <Ref>, we introduced simplicial complexes which enable us to study the topological spaces via computational tools. To study the given data in different resolutions, PH generates a nested sequence of simplicial complexes _1⊂_2⊂…⊂_n induced from the data. Such a sequence is called a filtration. This step can be considered the most crucial for the effectiveness of PH in ML applications. The primary reason is that the filtration process involves examining the data at different resolutions by adjusting a "scale parameter". The choice of this "scale parameter" can greatly influence the performance of the method.

For each data type, there are well-established methods to construct filtrations that have proven highly effective in their respective contexts. These methods vary significantly depending on the data type. While point cloud <cit.> and image <cit.> settings use relatively common approaches in their settings, the construction of filtrations in graph settings is notably more versatile <cit.>.

r3in


    [ϵ=0.4]
        
    < g r a p h i c s >

        
    
    [ϵ=0.7]
        
    < g r a p h i c s >

        
    
    [ϵ=2]
        
    < g r a p h i c s >

        
    
     For 8-shaped point cloud , {_()} filtration steps for thresholds =0.4,=0.7 an =2.      
    





  §.§.§ Filtrations for Point Clouds
 
As the process is generally similar for various metric spaces, for simplicity, we will describe it for a point cloud  in ℝ^N using the Euclidean metric. Let = {x_1, x_2, …, x_m} be a point cloud in ℝ^N. We will define a nested sequence of simplicial complexes _1 ⊂_2 ⊂…⊂_n induced by .  The central idea here is to build a series of simplicial complexes that progressively capture the topological picture of the point cloud  in different resolutions as we move from _1 to _n. The "nested" part signifies that each simplicial complex is contained within the next one, like a series of Russian Matryoshka dolls.

 


Before moving on to simplicial complexes, we will define a simpler nested sequence for . Let _r(x)={y∈^N| d(x,y)≤ r} be the closed r-ball around x. Then, let r-neighborhood of , _r() = ⋃_i=1^m _r(x_i) be the union of r-balls around the points in . By declaring r as our scale parameter, we first fix a monotone sequence of threshold values 0 = r_1 < r_2 < … < r_n, where r_n = max_i,j{d(x_i, x_j)}, the diameter of . These values intuitively represent the resolution at which we observe the point cloud . In particular, a smaller value of r indicates a closer, more detailed examination of , while a larger value of r means observing the point cloud with a broader view from a greater distance, making it difficult to distinguish between points that are close together in . This naturally gives a nested sequence of topological spaces _r_1() ⊂_r_2() ⊂…⊂_r_n() (See <Ref>). While these neighborhoods  {_r()} give a natural sequence for the data, to effectively leverage computational tools, we need to induce a sequence by simplicial complexes. There are two common ways to achieve this while preserving the underlying topological information.






 
r3in


    [Cech Complex]
        
    < g r a p h i c s >

            
       [Rips Complex]
        
    < g r a p h i c s >

            
        [Cech Complex]
        
    < g r a p h i c s >

            
    Comparison of Rips and Čech Complexes. In panels <ref> and <ref>, the Čech and Rips complexes differ: the Čech complex does not form a 2-simplex among v_1, v_2, and v_3 due to insufficient ball overlap, while the Rips complex does. A larger radius, as in panel <ref>, is needed for the Čech complex to include this simplex.
    
    





  
*Method 1 - Rips complexes: For ⊂^N, and for r>0, the Rips complex (aka Vietoris-Rips complex) is the abstract simplicial complex _r() where a k-simplex σ=[x_i_0,x_i_1,…,x_i_k]∈_r()
if and only if d(x_i_m,x_i_n)< r for any 0≤ m,n≤ k. In other words, for r>0, if for k+1-points are pairwise r-close to each other, they form a k-simplex in _r(). 




  
*Method 2 - Čech complexes: Similarly,  for ⊂^N, and for r>0, the Čech complex is the abstract simplicial complex _r() where a k-simplex σ=[x_i_0,x_i_1,…,x_i_k]∈_r() if and only if ⋂_j=0^k B_r(x_i_j)≠∅. Here, the condition to build k-simplex is different as we ask for a nontrivial intersection of the r-balls of the k+1-points. 




 




The main relationship between our original filtration {_r_i()} and the simplicial complex filtrations arises from the Nerve Lemma. This lemma states that the Čech complex 𝒞̌_r() is homotopy equivalent to _r() for any r ≥ 0 <cit.>. Since homotopy equivalent spaces have the same homology, we have _k(𝒞̌_r()) ≃_k(_r()) for any k ≥ 0. Therefore, from a PH perspective, the filtrations {_r_i()}_1^n and {𝒞̌_r_i()}_1^n are essentially the same. Furthermore, for ⊂^N with the Euclidean metric, Rips and Čech complexes are closely related and produce similar topological information as  𝒞̌_r() ⊂_2r() ⊂𝒞̌_√(2)r() <cit.>.

While both complexes provide similar filtrations, Rips complexes are more commonly used in practice. This preference is due to the fact that Rips complexes only require the pairwise distances {d(x_i, x_j)} among points, which can be easily obtained at the beginning of the process. Hence, constructing Rips complexes is straightforward once these distances are known. In contrast, constructing Čech complexes requires checking whether collections of r-balls have nontrivial intersections. <Ref> depicts their differences in a toy scenario.

In particular, Rips complexes are the most common filtration type used for point clouds because of their computational practicality, and because of the Nerve Lemma, they capture very similar information produced by the simple neighborhood filtration {N_r()} described earlier.

Although we discuss the construction for point clouds in ^N, the same process can be applied to any metric space (any space with a distance function). Furthermore, in some cases, the point cloud  is given in an abstract setting, where only pairwise distances {d(x_i, x_j)} among the points are provided. The Rips complex filtration can be effectively applied to such point clouds. In <Ref>, we detail the real-life application of this approach for shape recognition. 




  
*Other complex types. One of the primary challenges in applying PH to point clouds is the computational cost. To mitigate this, witness complexes offer a valuable alternative for efficiently analyzing the topological features of point clouds, especially in high-dimensional spaces. Unlike the Rips complexes, which can be computationally expensive, witness complexes reduce complexity by utilizing a representative subset of the data points, known as witnesses, to construct the simplicial complex <cit.>. This method strikes a balance between computational efficiency and topological accuracy, making it particularly well-suited for large datasets where constructing full complexes would be impractical. By concentrating on representative data points (witnesses), witness complexes facilitate a more manageable and scalable computation of persistent homology, enabling the detection of the underlying shape and features of the point cloud across various scales. In addition to Rips and Čech complexes, other types of complexes, such as alpha complexes and Delaunay complexes, are particularly useful in lower-dimensional spaces. However, to maintain the focus of this paper, we refer the reader to <cit.> for an in-depth discussion of these complexes.
 











  §.§.§ Filtrations for Images
 

For images, constructing filtrations differs significantly due to the unique structure of image data. To keep the explanation simple, we will focus on 2D images, though the concepts can be extended to 3D and other types of images. Filtration for images typically involves nested sequences of binary images, known as cubical complexes <cit.>. Starting with a given color (or grayscale) image  of dimensions r × s, we first select a specific color channel (e.g., red, blue, green, or grayscale). The color values γ_ij∈ [0, 255] of individual pixels Δ_ij⊂ are used, where Δ_ij represents a closed box (square including its boundary) in the i^th row and j^th column of the image .



Next, we choose the number  of thresholds "n" to span the color interval [0, 255], i.e., 0 = t_1 < t_2 < … < t_n = 255. This determines the length of our filtration sequence, with a typical range being between 50 and 100. Using these thresholds, we define a nested sequence of binary images (cubical complexes) _1 ⊂_2 ⊂…⊂_n such that _m = {Δ_ij⊂|γ_ij≤ t_m} (see <Ref>).

In particular, this involves starting with a blank r × s image and progressively activating (coloring black) pixels as their grayscale values reach each specified threshold t_m. This process is known as sublevel filtration, applied to  relative to the chosen color channel (in this case, grayscale). Alternatively, pixels can be activated in descending order, referred to as superlevel filtration. In this context, let _m = {Δ_ij⊂|γ_ij≥ s_m}, where 255 = s_1 > s_2 > … > s_n = 0, and _1 ⊂_2 ⊂…⊂_n is called superlevel filtration. The persistence homology process involving such cubical complexes has a special name, called cubical persistence. In <Ref>, we detail a successful application of this method in cancer detection from histopathological images.

While these sublevel and superlevel filtrations are common choices for color or grayscale images, there are other filtration types used for binary images, e.g., erosion, dilation, and signed distance <cit.>. These filtrations are specific to binary images and effective in capturing the size and other properties of topological features. In <Ref>, we give an example of erosion filtration.


















  §.§.§ Filtrations for Graphs
 

Graphs are a widely used application domain for PH because they generate highly effective features that are not easily obtained through other methods. While filtration methods for point clouds and images are relatively standard, graph filtrations offer a variety of choices. The way you construct these filtrations can significantly impact the performance of ML models. Unlike other data formats, there are various methods to construct filtrations from graphs, where the details can be found in <cit.>. We categorize these methods into two groups:




  
*i. Filtrations through node/edge functions. This type of filtration is computationally efficient in most cases and is commonly used in applications. The main concept involves defining a filtration function on nodes or edges and using the order dictated by these functions to create a sequence of subgraphs and corresponding simplicial complexes. 

r3in



    < g r a p h i c s >

Graph Filtration. For 𝒢=𝒢_3 in both examples, the top figure illustrates a superlevel filtration using the node degree function with thresholds 3>2>1, where nodes of degree 3 are activated first, followed by those of lower degrees.  Similarly, the bottom figure illustrates a sublevel filtration based on edge weights with thresholds 1.5< 1.8< 2.1.  




Given a graph =(,), let f:→ be a node filtration function. Common examples of such functions include degree, centrality, betweenness, closeness, and heat kernel signatures (HKS). Additionally, functions may be derived from the domain of the graph, such as atomic number for molecular graphs or account balance for transaction networks. These functions establish a hierarchy among the nodes, ordering them from less important to more important within the graph within the context defined by the filtration function.


To define the resolution of our filtration, we choose a set of thresholds that cover the range of f, denoted as ℐ={ϵ_i}_1^n, where ϵ_1=min_v ∈ f(v) < ϵ_2 < … < ϵ_n = max_v ∈ f(v). For each threshold ϵ_i ∈ℐ, let _i = {v ∈| f(v) ≤ϵ_i}. Define _i as the subgraph of  induced by _i, i.e., _i = (_i, _i) where _i = {e_ij∈| v_i, v_j ∈_i}. In other words, at each threshold, we activate the nodes whose values reach that threshold, along with the edges in the graph between these activated nodes.


This process results in a nested sequence of subgraphs _1 ⊂_2 ⊂…⊂_n = (See <Ref>). For each _i, we define an abstract simplicial complex _i for 1 ≤ i ≤ n, creating a filtration of simplicial complexes _1 ⊆_2 ⊆…⊆_n. Typically, clique complexes are used, where the complex  is obtained by assigning a k-simplex to each (k+1)-complete subgraph ((k+1)-clique) in  <cit.>. The term clique refers to a group (clique) of k+1 vertices forming a k-simplex.  This is called sublevel filtration for the filtration function f. Similarly, one can reverse this process (using a decreasing order for activation of the nodes) to obtain a superlevel filtration, where _i = {v ∈| f(v) ≥ϵ_i'} where {_i'} is a monotone decreasing sequence. 

While we have explained the process for node filtration functions, edge filtration functions are also common. This approach is frequently employed in weighted graphs or graphs with key edge functions for downstream tasks. For a graph , let g: → be an edge filtration function. In the case of weighted graphs, the weight function w: → with w(e_ij) = ω_ij serves as an example of such a function. Other common examples include Forman Ricci and Ollivier Ricci curvatures, which assign the (weighted or unweighted) edge a curvature value by interpreting the geometry of the edge's neighborhood <cit.>. Furthermore, edge filtration functions can be derived from the graph's domain, such as transaction amounts in blockchain networks <cit.> or density in traffic networks <cit.>. These functions establish a hierarchy among the edges, similar to node filtration functions.

  

Again, by choosing the threshold set ℐ={_i}_i=1^n with _1=min_e_ij∈g(e_ij)<_2<…<_n=max_e_jk∈g(e_jk), we define the filtration as follows. For _i∈ℐ, let ℰ_i={e_jk∈ℰ| g(e_jk)≤_i}. 
Let 𝒢_i be a subgraph of 𝒢 induced by ℰ_i, i.e., 𝒢_i is the smallest subgraph of 𝒢 containing the edge set ℰ_i (Fig. <ref>). In other words, for each threshold, we activate the edges whose value reaches that threshold, along with the nodes attached to them. Again, this induces a nested sequence of subgraphs 𝒢_1⊂…⊂𝒢_n=𝒢. Then, one can follow the same steps by using clique complexes as before. 

In <Ref>, we give a toy example for superlevel (node) and sublevel (weight) filtrations. Here, we list the most common node and edge filtration functions used in applications:


    
  * Node Filtration Functions: Degree, betweenness, centrality, eccentricity, heat kernel signature (HKS), Fiedler values (spectral), node functions from the domain of the data.
    
  * Edge Filtration Functions: edge weights (weighted graphs), Ollivier Ricci, Forman Ricci curvature, 







 




  
*ii. Graph distance-based filtrations. Next to the sublevel filtrations above, a completely distinct filtration method for graphs uses the distances between nodes. Essentially, it treats the graph as a point cloud, where pairwise distances between nodes are defined by graph distance, and applies Rips filtration (<Ref>) to this point cloud. While this method is more computationally intensive than filtrations via node or edge functions, it can be highly effective for small graphs as it captures distance information and finer topological details.


To apply this method, we need to define graph distances between nodes. In an unweighted graph = (, ), a common approach is to set the distance between adjacent nodes to 1 (e_ij = 1) and define the distance between v_i and v_k as the length of the shortest path τ_ik in  with endpoints v_i and v_k, i.e. d(v_i, v_k) = minτ_ik. Thus, each pairwise distance is an integer representing the number of hops needed to travel from v_i to v_j in . The largest such distance is the diameter of , denoted as n.

After obtaining the pairwise distances {d(v_i, v_j)} between the nodes, we treat the graph as a point cloud = {v_1, v_2, …, v_m} with these distances. Since all distances are integers, we use integer thresholds {r_i = i} for 0 ≤ i ≤ n = diam(). Now, we are ready to define our filtration. For the filtration index, we use superscripts instead of subscripts, where the reason will become clear later. 




To make the exposition simpler, let ^k be the Rips complex corresponding to r=k, and let ^k be the graph corresponding to the 1-skeleton of ^k. Here 1-skeleton means the graph itself with its nodes and edges, without considering any higher-dimensional elements like faces or solids that might be part of a more complex simplicial structure. For the distance parameter r = 0, we have ^0 = since there are no edges in the Rips complex. It is easy to see that at the distance threshold of 1, 𝒢^1 = 𝒢 because all edges are automatically included in the complex when r = 1. Recall that ^1 is a Rips complex, hence any complete j-subgraph in ^1 generates a (j-1)-simplex (j-clique) in the Rips complex ^1. In particular, ^k is nothing but the clique (or flag) complex of ^k. Furthermore, ^k is the graph with node set , and edge set ^k = {e_ij| d(v_i, v_j) ≤ k}, i.e., ^k=(,^k). Therefore, in the filtration ^0 ⊂^1 ⊂…⊂^n, all simplicial complexes have the same node set, but at each step, we add new edges and cliques. Finally, for ^n, the 1-skeleton is a complete graph with m nodes, and hence the corresponding simplicial complex ^n would be an (m-1)-simplex. 

r3in

    < g r a p h i c s >

Graph Powers. A graph G=G^1 and its graph powers.
	Red edges are added in G^2, and green ones are added in G^3.
	Note 𝒢^3 is the complete graph on 7 vertices since D=diam(G)=3. Hence, 𝒢^3 would be a complete graph 𝒦_7, and all higher powers are the same, i.e. 𝒢^n=𝒢^3 for n≥ 3.  

The graphs {^k} are called graph powers (See <Ref>), and this filtration is called the power filtration (or Rips filtration), hence the superscripts. Observe that the power filtration calculates the distances between every pair of vertices in the graph. Therefore, even vertices that are not direct neighbors or appear distant in the graph can still form a simplex in the later stages of the filtration. Note that you don't need to complete the entire filtration. In most cases, the critical insights information lies in the first few steps (e.g., up to =5 or 10) for the power filtration. Given the high computational cost, it's both practical and reasonable to stop early.
 

For a weighted graph = (, , ), edge lengths can be defined using weights {ω_ij}, where e_ij = f(ω_ij) instead of the uniform e_ij = 1 used in unweighted graphs. The function f: →^+ assigns distances based on weights, where a smaller f(ω_ij) indicates that nodes v_i and v_j are closely related, and a larger distance indicates they are less related. This approach is particularly useful in financial networks where large edge weights (i.e., amounts) indicate stronger financial connections or dependencies. In <Ref>, we elaborate on this method in the context of a real-world application: Crypto-token anomaly forecasting.

Once edge lengths are defined, pairwise distances between any node pair can be calculated as the length of the shortest path as before, i.e., d(v_i, v_j) = minϵ_ij, where ϵ_ij is any  such path in the graph with endpoints v_i and v_j. After obtaining pairwise distances, the process is the same as for unweighted graphs.


r2in


    
    < g r a p h i c s >

    Coral Reduction. With CoralTDA, 2-core of  has the same persistence diagrams _k() for k≥ 1.
    
    
 




  
*Graph Reduction and PH Although PH generates highly effective topological embeddings for graph representation learning, scalability remains a challenge. To address this issue, Akcora et al. <cit.> introduced two key methods to significantly reduce the computational costs of the PH process. The first method, CoralTDA, leverages the observation that nodes with low graph-core do not contribute to persistence diagrams in higher dimensions. Specifically, they demonstrate that the (k+1)-core of a graph, which is a subgraph where each vertex has at least k neighbors, is sufficient to compute the _k() of the original graph. This improvement can be implemented as few as three lines of Python code (see the repository).  The second algorithm, PrunIT, introduces an efficient method to reduce the size of graphs without altering their persistence diagrams. In this approach, a vertex u is said to dominate another vertex v if the 1-neighborhood of u contains 1-neighborhood of v, i.e., (u)⊃(v). The authors show that removing (pruning) a dominated vertex from the graph does not affect the persistence diagrams at any level as long as the dominated vertex enters the filtration after the dominating vertex.




  
*Which filtration to use for graphs? We first note that filtrations using node/edge functions and filtrations based on graph distances are entirely different methods, producing distinct outputs. Distance-based filtrations are computationally demanding, but for datasets with small graphs (such as molecular graph datasets), power filtrations can be highly effective. Conversely, for larger graphs, the computational costs necessitate the use of filtration functions. In these cases, the choice of filtration function could be critical for performance <cit.>. Typically, relevant node or edge functions specific to the domain of the data are the best choices. If such functions are unavailable, heat kernel signatures (HKS) for node functions and Ollivier Ricci for edge functions are excellent alternatives. Choosing a filtration function can be seen as an outdated method, as letting ML algorithms select or construct the best filtration function for optimal performance is more reasonable. However, in many settings, model interpretability or greater control over the process is needed. In these instances, selecting a relevant filtration function is more suitable for the model.

If the goal is to achieve optimal performance with PH machinery, learnable filtration functions are a promising alternative. There are significant works, along with available code, that can be adapted for various tasks in graph representation learning <cit.>.




  §.§.§ Choosing Thresholds
  One of the key steps in effectively applying PH to ML problems is selecting the thresholds {_i}_i=1^m for constructing the filtration. The first decision involves determining the number of thresholds, i.e., m. This number can be viewed as the resolution of the filtration: a larger m implies higher resolution, while a smaller m indicates lower resolution. Then, a natural question arises: is there any disadvantage to choosing a large m? The answer is "Yes". Firstly, the computational cost increases with m. Secondly, the key information captured in the data may become diluted in higher dimensions, depending on the vectorization step. 

For graph filtrations, depending on the filtration function, selecting m between 10 and 20 generally yields good results. Once the number of thresholds is set, they can be chosen equally spaced if the filtration function f: → (or g: →) is appropriate. Another effective method is to examine the distribution of the value set {f(v_i)} for all vertices in the graphs across the dataset (or a random subsample), and then select the thresholds as the corresponding quantiles in this distribution. Similarly, for distance-based filtrations, analyzing the distribution of pairwise distances between vertices in a random subsample of the dataset can provide valuable insights for selecting appropriate thresholds.



For sublevel filtrations in the image setting, it’s important to choose thresholds that cover the color range [0,255]. Typically, using 50-100 thresholds yields good results. However, if the task only concerns a specific color interval [a,b], it’s advisable to concentrate most or all thresholds within that range. These thresholds don't need to be evenly spaced. In the case of binary image filtrations (e.g., erosion, height, signed distance), thresholds are usually evenly distributed up to the maximum value.


r2.8in


    [ϵ=0.4]
        
    < g r a p h i c s >

        
    
        [PB up to ϵ=0.4]
        
    < g r a p h i c s >

        
    
    
    [ϵ=0.7]
        
    < g r a p h i c s >

        
    
        [PB up to ϵ=0.7]
        
    < g r a p h i c s >

        
    
    
    [ϵ=2]
        
    < g r a p h i c s >

        
    
    [Persistence barcode]
        
    < g r a p h i c s >

        
    
    Evolution of Persistence Barcode (PB). In PBs shown on the right, red bars represent PB_0(), and blue bars represent PB_1(). The 8-shaped point cloud  with 35 points,  initially has 35 components, corresponding to 35 red bars. As ϵ increases, these components merge. At ϵ = 0.4, only the top 11 red bars remain (b), indicating the number of components in 𝒩_0.4() (a). By ϵ = 0.7, the space becomes fully connected (c), and all red bars in PB_0() terminate, except for the top ∞-bar (d). Regarding loops (1-holes), a small loop appears around ϵ = 0.3 (a-b) and persists in 𝒩_0.7() (c), while a larger loop emerges around ϵ = 0.6 (d). In PB_1() (f), two blue bars, [0.3, 0.9) and [0.6, 1.7), correspond to the small and large loops in the 8-shape.    






For point clouds, the number of thresholds directly determines the filtration's resolution. While thresholds are generally evenly spaced, in the presence of outliers, they can be more sparsely distributed after a certain value. Additionally, since computational cost is often a concern with point clouds, it's important to find a balance between the number of thresholds and the associated computational expense.

 



 §.§ Persistence Diagrams
 

 
After constructing the filtration _1 ⊂_2 ⊂…⊂_n for a data type , PH systematically tracks the evolution of topological features (k-holes) in the filtration sequence and records this information in a persistence diagram, which we define next. The nontrivial elements in the homology groups _k(_i) for 1≤ i≤ n represent the k-dimensional topological features (or k-holes) appearing in the filtration sequence. Furthermore, the inclusion map ι: _i ↪_i+1 allows us to determine if a k-hole σ in _i persists in _i+1 through the induced map ι_*: _k(_i) →_k(_i+1).



For each k-hole σ, PH records its first appearance in the filtration sequence, denoted _i_0, and its first disappearance in a later complex, _j_0. We define b_σ = ϵ_i_0 as the birth time of σ and d_σ = ϵ_j_0 as the death time of σ, where {ϵ_i}_1^n is the threshold set used for the filtration. The difference d_σ - b_σ is called the lifespan of σ. For example, if a k-hole τ first appears in _3 and disappears in _7, we mark the birth time as b_τ = ϵ_3 and the death time as d_τ = ϵ_7. The lifespan of τ is then ϵ_7 - ϵ_3.

 
 



 


For each nontrivial σ∈_k(_i) for 1 ≤ i ≤ n, we represent σ with a 2-tuple (b_σ, d_σ) to denote its birth and death times in the filtration. The collection of all such 2-tuples is called the persistence diagram (PD) as depicted in Figure <ref>. Note that a topological feature with the same birth and death time (0 lifespan) is considered a trivial topological feature, and they are represented with diagonal elements in the persistence diagrams. Therefore, the diagonal Δ={x=y}⊂^2 is always included in any persistence diagram. 
Then, k^th persistence diagram is defined as

    _k()={(b_σ, d_σ) |σ∈_k(_i)  for  b_σ≤ i < d_σ}∪Δ


Since trivial elements (t,t) in PDs can appear multiple times, we take the diagonal Δ with infinite multiplicity. Essentially, the persistence diagram is a subset of ^2 (_k() ⊂^2), where each point in _k() is either a pair of threshold values (ϵ_i, ϵ_j) for some i < j, or belongs to the diagonal (trivial element). Here, infinite multiplicity is a technical assumption, which will be important when discussing Wasserstein distance between PDs (<Ref>).

r3in


    [ϵ=2]
        
    < g r a p h i c s >

            
    [Persistence barcode]
        
    < g r a p h i c s >

        
       Persistence Barcode (left) and Persistence Diagram (right) for the 8-shaped point cloud  (<Ref>). Both representations convey the same information, while the PD plots each bar's birth time as its x-coordinate and death time as its y-coordinate. For example, the 35 red bars in PB_0() that begin at 0 are represented as 35 red points along the x=0 line in _0(), where the y-coordinates correspond to the bars' death times. Similarly, the blue points at (0.3, 0.9) and (0.6, 1.7) represent the blue bars [0.3, 0.9) and [0.6, 1.7), respectively.  
              

       

        


There is an equivalent concept called the persistence barcode (see Figure <ref>), which uses bars (half-open intervals) {[b_σ, d_σ)} instead of 2-tuples {(b_σ, d_σ)}. The bar notation [b_σ, d_σ) represents the entire lifespan of a topological feature σ as such an interval. In practice, persistence diagrams are more commonly used due to their practicality in applications.
 

Note that the choice between sublevel and superlevel filtration is irrelevant for the point cloud setting. In the image setting, this choice is not essential, as they essentially convey the same information due to Alexander duality <cit.>; when analyzing the topological features of images, the choice of focusing on either the binary image or its complement is not crucial because the topological information of one is directly related to the topological information of the other through this duality. The Alexander duality holds only when the ambient space is quite simple, such as when the space has a structure like an m × n rectangle, like images. However, sublevel and superlevel filtrations can yield significantly different results in the graph setting. Carriere et al. <cit.> introduced extended persistence diagrams to address this issue.  
This approach enables the birth and death pairs of the persistence diagram to be positioned below the diagonal (b > d). By doing so, it captures topological features identified through both superlevel and sublevel filtrations.

 
 



 



  §.§.§ Interpretation of Persistence Diagrams

A persistence diagram _k(𝒳) records the k-dimensional topological features of the data 𝒳 as a collection of points in ℝ^2. For example, _0(𝒳) records 0-holes (components), and _1(𝒳) records 1-holes (loops) appearing in the filtration sequence {𝒦_i}. Each point q_j = (x_j, y_j) represents a k-dimensional hole σ_j. For instance, consider a threshold set ℐ = {ϵ_i = i/10}_i=0^100 with 100 thresholds spanning the interval [0, 10]. Suppose we have two points, q_1 = (0.3, 9.7) and q_2 = (4.2, 4.6), in _2(𝒳). These points represent 2-holes (cavities) σ_1 and σ_2 that appear in the filtration 𝒦_0 ⊂𝒦_1 ⊂…⊂𝒦_100.

For q_1 = (0.3, 9.7), the cavity σ_1 first appears in 𝒦_3 and persists until 𝒦_97, indicating a long lifespan. Conversely, for q_2 = (4.2, 4.6), the cavity σ_2 first appears in 𝒦_42 and persists only until 𝒦_46, indicating a short lifespan. This suggests that σ_1 represents an important topological feature of the data 𝒳, while σ_2 is likely just topological noise.

In general, features with long lifespans (where y - x is large) are located far from the diagonal and are considered significant (or "big") features. Features with short lifespans (where y - x is small) are close to the diagonal and are considered insignificant (or "small") features. This is why many vectorization techniques try to involve the lifespan information in their computation to give different emphasis to small and big features in the output.



  §.§.§ Wasserstein Distance
 

After obtaining the persistence diagrams (PDs) for two datasets, 𝐗^+ and 𝐗^-, we can assess their topological similarities using these diagrams. One effective approach is to employ a metric in the space of PDs. If the distance between the two persistence diagrams is small, we can conclude that the two datasets exhibit similar topological characteristics; otherwise, they differ.
The most commonly used metric for this purpose is the Wasserstein distance (also known as the matching or earth mover distance) of the Optimal Transport Theory (see an ICML tutorial in <cit.>), which is defined as follows:

Let (^+) and (^-) be the persistence diagrams for the datasets ^+ and ^-, respectively (we omit the dimensions in PDs for simplicity). Denote (^+) = {q_j^+}∪Δ and (^-) = {q_l^-}∪Δ, where Δ represents the diagonal (indicating trivial cycles) with infinite multiplicity, and q_j^± = (b^±_j, d_j^±) ∈(^±) represents the birth and death times of a topological feature σ_j in ^±. Let ϕ: (^+) →(^-) represent a bijection (matching). The presence of the diagonal Δ on both sides ensures the existence of these bijections even if the cardinalities |{q_j^+}| and |{q_l^-}| differ. This means that for optimal matching, some points in {q_j^+} and {q_l^-} are matched to diagonal elements whenever necessary. Then, the x^th Wasserstein distance _p is defined as

    _p((^+), (^-)) = min_ϕ( ∑_j q_j^+ - ϕ(q_j^+)_∞^p )^1/p,    p ∈ℤ^+.


Here, (x_1,y_1)-(x_2,y_2)_∞=sup{|x_1-y_1|,|x_2-y_2|} is called supremum norm (or l^∞-norm). When p=∞, _∞ is called the bottleneck distance. i.e.,
_∞((^+), (^-)) = max_j {q_j^+ - ϕ(q_j^+)_∞} (See <Ref>).



r2.5in
    
    < g r a p h i c s >

    Wasserstein distances _1 and _∞ (bottleneck) between _1(8) and _1(9) for the point clouds of "8" and "9" on the left.
    
    


In applications, the most common choices for the Wasserstein distance are p=1, 2, and ∞. The bottleneck distance, unlike the p-Wasserstein distance, is insensitive to the number of points in _k(^±) and instead focuses on the distance between the farthest points in the optimal matching. For instance, ignoring the diagonal, if _k(^+) has 100 points and _k(^-) has only three points, the bottleneck distance is primarily determined by finding the three closest points in _k(^+) to those in _k(^-) and taking the maximum of these distances. Conversely, when using p=1 or p=2, the quantity of points in _k(^±) becomes significant as all points contribute to the distances, even if only slightly. Therefore, if there are only a few points in {_k(_j)} and the focus is on the distances of the most critical features, the bottleneck distance is preferable. However, if there are many points in {_k(_j)} and the primary topological patterns arise from the quantity and location of smaller features, then choosing p=1 (or p=2) for the Wasserstein distance would be more suitable. 
 

There are several ways to utilize the Wasserstein distance in applications. One common approach is to bypass the vectorization step, which requires some choices, and directly compare persistence diagrams. For instance, to determine if two datasets have similar shapes or structures, one can compute the persistence diagrams for both datasets and then calculate the Wasserstein distance between these diagrams. A smaller distance indicates more similar shapes. Another application lies in unsupervised learning, where datasets are clustered based on their topological similarity using persistence diagrams, utilizing the Wasserstein distance to measure distances between them.




 §.§ Integrating PDs to ML tasks 1: Vectorizations
 

To effectively use persistence diagrams (PDs) in ML applications, it is crucial to convert them into numerical or vector formats, allowing for the seamless integration of TDA outputs into standard ML workflows. Although PDs capture the birth and death of topological features within data, their variable size and structure make them challenging to handle directly. For example, in a binary graph classification task with 1000 graphs (400 in Class A and 600 in Class B), PDs might visually highlight differences between the classes. However, applying traditional statistical tools to PDs, which are subsets of ^2, is problematic; for instance, it is not possible to compute averages or confidence intervals for each class. Moreover, directly inputting PDs into ML models is impractical because their variable point count conflicts with the fixed-size input required by most ML algorithms.

Vectorization addresses this by making PDs compatible with traditional ML and statistical techniques. Below, we outline the most common vectorization methods used in practice. For a comprehensive review, refer to <cit.>.


 

In the following, for a fixed dataset 𝒳, a threshold set {ϵ_i}_i=1^n, and an induced filtration {𝒦_i}_i=1^n, we will introduce several vectorization methods. It is important to note that all these vectorization methods are applicable to any type of data since they simply transform a given PD into vectors. However, the effectiveness and common usage of certain vectorization methods can vary depending on the data type and the density or sparsity of the PDs. Although we present various vectorization methods and their configurations here, those who prefer to avoid the intricacies of vectorization selection and hyperparameter tuning can opt for automated approaches, as detailed in <Ref>.
 



  
*Function vs. Vector Format. In the following, we introduce several vectorization methods. For each method, we will describe how to convert a given PD into a vector or a function depending on the setting. Depending on the context, one of these formats might be preferable (visualization, ML input, etc.), however, both formats can be easily converted to each other. For example, if we have a function f:[ϵ_1,ϵ_n]→ℝ defined over an interval, we can convert f into an N-dimensional vector by sampling the values at f(t_j) where t_j = ϵ_1 + j/N(ϵ_n - ϵ_1) and j ranges from 1 to N. For each point t_j, we calculate f(t_j). The vector 𝐯_f is then  𝐯_f=[f(t_1), f(t_2), …, f(t_N)]. Conversely, for a given N-dimensional vector 𝐯=[v_1  v_2  … v_N], we can convert it to a function via a step function or linear spline interpolation, e.g., g_𝐯:[ϵ_1,ϵ_n]→ℝ such that g_𝐯(t_j)=v_j for t_j=ϵ_1+j/N(ϵ_n-ϵ_1). Then, for any point s∈[t_j,t_j+1], we define the linear extension g(s)=g(t_j)+g(t_j+1)-g(t_j)/t_j+1-t_j(s-t_j). Hence, any of the following vectorizations can be taken as a function or vector depending on the need.




  
*Topological Dimensions and Final Topological Vector. In each vectorization method, we convert a persistence diagram _k() into a vector (or array) 𝐯_k where k represents a topological dimension of k-holes. In particular, each persistence diagram produces a different vector. The user needs to decide which topological dimensions to be used in the method. After getting an N-dimensional vector for each dimension, a common method is to concatenate them to obtain the final topological vector. Hence, if we use m different topological dimensions, we have (m· N)-dimensional final vector 𝐯()=𝐯_0‖𝐯_1‖…‖𝐯_m-1. In most applications, the common dimensions used are k=0 and k=1, and hence the final topological vector is 2N-dimensional, i.e., 𝐯()=𝐯_0‖𝐯_1.




  
*Betti vectors One of the most straightforward and interpretable vectorization methods in TDA is the Betti function. The k^th Betti number of a topological space  essentially counts the total number of k-dimensional holes in . More formally, it is defined as β_k() = rank(H_k()), which is the rank of the k^th homology group of . For example, β_0() represents the number of connected components in , while β_1() indicates the number of 1-dimensional holes.



r0.4



    < g r a p h i c s >

Betti Function. The step function represents the first Betti number (β_1) over the figure-eight-shaped point cloud of Figure <ref> as a function of the threshold ϵ. The function shows the number of 1-dimensional holes (loops) in the point cloud. Initially, there are no loops (β_1 = 0), then a single loop appears at ϵ = 0.35 and persists until ϵ = 0.94, during which another loop appears at ϵ = 0.55 and vanishes at ϵ = 0.94. The final loop disappears at ϵ = 1.70, returning β_1 to 0. 
  
 


For given filtration {_i}_1^n, we define k^th  Betti vector β⃗_k()= [ β_k(_1)  β_k(_2) …β_k(_n)]. Therefore, for each topological dimension k, we obtain a Betti vector of the size of a number of thresholds. For example, k=0,1 with n=50 results in 50-dimensional vector for each dimension. In ML applications, typically, these vectors are concatenated to create a 100-dimensional topological vector as output.





  Consider the point cloud  in <Ref>. We define a filtration using five threshold values ϵ = [0, 0.25, 0.75, 1.5, 1.75], where _i = __i(). The Betti-0 vector β⃗_0() = [35, 20, 1, 1, 1] is obtained, with β_0() = 35, β_0(_0.25()) = 20, and β_0(_ϵ()) = 1 for ϵ≥ 0.7. This indicates that _0.25() consists of 20 connected components. The Betti-0 vector can be easily determined from the persistence barcode PB_0() (<Ref>) by counting the number of red bars intersected by a vertical line x=_i, representing the number of persistent topological features at the threshold _i. Thus, with 5 thresholds, we obtain a 5-dimensional vector, β⃗_0().




Similarly, we can determine the corresponding Betti-1 vector by examining the blue bars in the persistence barcode (<Ref>). We observe no blue bar at x = 0, 0.25, and 1.75. Furthermore, there are two blue bars at x = 0.75 and one blue bar at x = 1.5. Therefore, the Betti-1 vector is β⃗_1() = [0, 0, 2, 1, 0]. This indicates that there are no 1-dimensional holes (1-holes) in = _0(), _0.25(), or _1.75(), while _0.75() contains two 1-holes and _1.5() contains one. It is also common to use Betti vectors as Betti functions via step functions (See <Ref>).
Similarly, using 21 equally spaced thresholds = {0, 0.1, 0.2, …, 2} would yield 21-dimensional vectors β⃗_⃗0⃗ and β⃗_⃗1⃗.


We note that Betti vectors do not require the computation of persistence diagrams. Therefore, there are computationally more effective ways to produce Betti vectors <cit.>. Another favorable aspect of Betti vectors is their ease of interpretation. Simply put, β_k(ϵ_i) is equal to the number of k-dimensional topological features in 𝒦_i.



 r2in



    < g r a p h i c s >

Persistence Landscape (PL). The plot shows the PL for _1() for the 8-shaped point cloud filtration in <Ref>. In the plot, the blue function corresponds to λ^1, and the orange function to λ^2. 






  
*Persistence Landscapes Persistence Landscapes are one of the first vectorization methods in TDA, introduced by P. Bubenik, directly utilizing the lifespan information <cit.>. In particular, in this vectorization, the points away from the diagonal (large features) are easily distinguished and promoted. For a given persistence diagram _k()={(b_i,d_i)}, we first define generating functions Λ_i for each (b_i,d_i)∈_k(), i.e., Λ_i:[b_i,d_i]→ is a piecewise linear function obtained by two line segments starting from (b_i,0) and (d_i,0) connecting to the same point (b_i+d_i/2,b_i-d_i/2). Then, we define several piecewise-linear functions {λ^m} in the following way. For each t∈[ϵ_1,ϵ_n], we check all generating functions {Λ_i(t)}, and we mark m^th largest value. In particular, m^th Persistence Landscape function λ^m():[ϵ_1,ϵ_n]→ for t∈ [ϵ_1,ϵ_n] is defined as λ^m()(t)=m^th max_i{Λ_i(t)} (See <Ref>). Note that the first and second persistence landscapes are the most commonly used vectors, and they are used with concatenation in the applications. 
 


  r2.5in
 

    < g r a p h i c s >

Silhouette Functions. The plot shows the silhouette functions with tuning parameters p=0.5,1,2 _1() for the 8-shaped point cloud filtration given in <Ref>. 
 






  
*Silhouettes While persistence landscapes are among the first vectorizations to effectively utilize lifespan information, the need to consider m different maxima makes it difficult to use in ML applications.
In <cit.>, Chazal et al. proposed a practical modification called the Silhouette for persistence landscapes. This modification introduces a tuning parameter p to better utilize the lifespans of topological features. For a persistence diagram _k()={(b_i,d_i)}_i=1^N, let Λ_i be the generating function for (b_i,d_i) as defined in the persistence landscapes. The Silhouette function ψ is defined as: Ψ_k()=∑_i=1^N w_iΛ_i(t)∑_i=1^m w_i,  t∈[ϵ_1,ϵ_n],
where the weight w_i is usually chosen as the lifespan (d_i-b_i)^p. Thus, the p-silhouette function ψ^p_k():[ϵ_1,ϵ_n]→ is defined as: 

 Ψ^p_k()=∑_i=1^N (d_i-b_i)^pΛ_i(t)∑_i=1^N (d_i-b_i)^p




The tuning parameter p is crucial as it adjusts the silhouette function's emphasis on topological features with varying lifespans. When p<1, shorter lifespans (d_i-b_i) are given more weight, emphasizing smaller features. Conversely, shorter lifespans are down-weighted when p>1, highlighting larger features. Common choices for p are 1/2, 1, and 2 (See <Ref>). If the persistence diagram has a few points and the goal is to emphasize significant features, p=2 would be a good choice. If there are many points and the key information comes from smaller features, p=1/2 is considered more suitable. Similarly, Silhouette functions ψ^p_k() can be converted to vectors for each dimension k, and concatenation of these vectors can be used in the application as the topological vector of the dataset.



 r2.5in
 


    < g r a p h i c s >

Persistence Curves. The plot shows the persistence curves for _1() for the 8-shaped point cloud filtration in <Ref> by using two different summary statistics: sum and mean. The dashed line represents the persistence curve using the sum of lifespans, which emphasizes the cumulative contribution of all loops. The solid line represents the persistence curve using the mean of lifespans, which highlights the average significance of the loops. This visualization illustrates how different summary statistics influence the interpretation of topological features in the data.

 




  
*Persistence Curves The Persistence Curve (PC) framework, introduced by Chung et al. <cit.>, offers a structured approach to the vectorization process. The core idea revolves around utilizing a generating function ψ, a summary statistic T, and a filtration parameter ϵ. From these elements, a PC is defined as: 

PC(, ψ, T )(t) = T ([ψ(; b, d, ϵ) | (b, d) ∈_ϵ ]), ϵ∈. 

Here, _ϵ represents the points within the persistence diagram  that fall inside a region Δ_ϵ varying with the filtration value ϵ.  
The function ψ takes as input a point from the persistence diagram ,  and the filtration parameter, ϵ and outputs a real number.  This function can be chosen to prioritize certain points in the PD or to emphasize specific features. Finally, the statistic T  acts upon a multiset of values, aggregating them into a single real number.  Common examples include sum, mean, or max.  

 
PCs provide a general and unifying framework for vectorization methods.
By selecting different combinations of ψ and T, one can generate various functional summaries of the PD, each potentially highlighting different aspects of the data. Many established vectorizations of PDs can be represented within the PC framework, e.g., Betti curves, persistence landscapes, silhouettes, and entropy curves. See <cit.> for more details.



 

r3in


    
        
    < g r a p h i c s >

        
    
        
    < g r a p h i c s >

        
    
        
    < g r a p h i c s >

        
    Persistence images (PI). PIs are vectorization with 2D (matrix) output. Here, we give examples of PIs for _1() (8-shaped point cloud in <Ref>) with different spread values. The spread parameter σ controls the standard deviation of the Gaussian functions used to smooth the persistence points in the persistence diagram. As illustrated in the figures, as σ increases, the generating Gaussian functions are even more spread out, producing a flatter persistence image.    
        
        





  
*Persistence Images Our next vectorization is Persistence Images, introduced by Adams et al. <cit.>. Unlike most vectorizations, Persistence Images, as the name suggests, produce 2D-arrays (tensors). The idea is to capture the location of the points in the PDs with a multivariable function by using the 2D Gaussian functions centered at these points. For PD()={(b_i,d_i)}, let ϕ_i represent a 2D-Gaussian centered at the point (b_i,d_i)∈^2. Then, one defines a multivariable function, Persistence Surface, μ=∑_iw_iϕ_i where w_i is the weight, mostly a function of the life span d_i-b_i. To represent this multivariable function as a 2D-vector, one defines a k× l grid (resolution size) on the domain of μ, i.e., threshold domain of PD(). Then, one obtains the Persistence Image, a 2D-vector (matrix) μ⃗=[μ_rs]  of size k× l such that 

    μ_rs=∫_Δ_rsμ(x,y) dxdy   where Δ_rs=  pixel with index rs  in the  k× l  grid.



Note that the resolution size k× l is independent of the number of thresholds used in the filtering; the choice of k and l is completely up to the user. There are two other important tuning parameters for persistence images, namely the weight w_i and the variance σ (the width of the Gaussian functions). Like Silhouettes, one can choose w_i=(d_i-b_i)^p to emphasize large or small features in the PD. Similarly, the width parameter σ determines the sharpness of Gaussian, where smaller σ would make the Gaussian functions more like Dirac δ-function, and larger σ would make the Gaussians flat. Depending on the context, σ can be chosen a constant (e.g. σ=0.1) or depending on the point (b_i,d_i), e.g., σ_i=k(d_i-b_i) for some constant k>0.





  
*Kernel Methods  Kernel methods provide an alternative to traditional direct vectorization techniques used to transform PDs into formats suitable for ML algorithms. In ML, these methods utilize a mathematical function known as a kernel to analyze data and discern patterns. Unlike earlier approaches that directly represent PDs using vectors or functions, kernel methods compute a similarity score as an inner product between pairs of PDs in a high-dimensional space without explicitly mapping the data. This approach is advantageous because it adapts well to ML techniques like support vector machines (SVMs) and kernel principal component analysis (KPCA). Therefore, kernel methods are well-suited for tasks such as classification, regression, and principal component analysis.

One common such method is the persistence weighted Gaussian kernel (PWGK) <cit.>. It enhances the Gaussian kernel by incorporating weights based on the significance of topological features. This approach amplifies the influence of important features while reducing noise (i.e., short-lived holes) impact. For instance, it assigns weights to points p=(b,d)∈() proportional to their lifespans w(p)=(d-b). In particular, PWGK is defined as

    𝐊((^+), (^-)) = ∑_p ∈(^+)∑_q ∈(^-) w(p) w(q) k(p, q)


where k(p, q) = exp(-p - q^2/2σ^2) denotes the Gaussian kernel function.

Another significant approach is the sliced Wasserstein kernel (SWK) <cit.>, which computes the Wasserstein distance between PDs and integrates it into a kernel framework. This method employs Optimal Transport Theory to establish a meaningful metric for comparing PDs. Although kernel methods can yield better results in some settings, they can be computationally intensive and impractical for large datasets due to the high computational costs associated with computing the kernel matrix. In particular, computing kernels takes quadratic time in the number of diagrams, while vectorizing PDs takes only linear time.





  §.§.§ Stability
 In most applications, the stability of vectorization is vital for statistical and inferential tasks. Essentially, stability means that a small change in the persistence diagram (PD) should not lead to a significant change in its vectorization. In particular, if two PDs, _k(^+) and _k(^-), are close, their corresponding vectorizations, β⃗_k(^+) and β⃗_k(^-), should also be close. This ensures that the vectorization process preserves the structural properties of the data. Therefore, when two persistence diagrams _k(^±) are similar, it implies that the datasets ^+ and ^- share similar shape characteristics. If these datasets are intuitively expected to belong to the same class, their vectorizations β⃗_k(^±) should likewise remain close.

To formalize this concept, we need to define what constitutes a "small/big change" or what it means for PDs to be close. This requires a distance (metric) in the space of PDs, with the most common being the Wasserstein distance (or matching distance) as defined in <Ref>. Similarly, we need a metric for the space of vectorizations, such as the Euclidean distance β⃗(^+) -β⃗(^-) for β⃗∈^N. Thus, for a given vectorization β⃗, we call it stable (with respect to the Wasserstein-p metric) if it satisfies the stability equation

    β⃗_k(^+) -β⃗_k(^-)≤ C ·_p(_k(^+), _k(^-))

This ensures that if _p(_k(^+), _k(^-)) (the distance between PDs) is small, then the distance β⃗_k(^+) -β⃗_k(^-) between the corresponding vectorizations will also be small.
Among the methods described earlier, persistence landscapes, silhouettes, persistence images, and most kernel methods are stable vectorizations, while Betti functions are generally unstable.




  §.§.§ Choice of Vectorization and Hyperparameters
  
The choice of vectorization method should align with the characteristics of your data and the problem at hand. If your data contains a few prominent topological features that are crucial to the task, Silhouettes with p ≥ 2 or Persistence Images may be the most suitable options. These methods are also effective when dealing with noisy data, allowing you to filter out less significant features. Conversely, if your data generates a high number of small features, where the task hinges on their location and density—in other words, when the noise itself carries important information—Betti curves, Persistence curves, Silhouettes with p≤ 0.5 and kernel methods are likely to yield strong performance. On the other hand, for point cloud data, Turkes et al. <cit.> provide valuable insights on how to effectively utilize PH for shape recognition. Lastly, if interpretability is a priority, Betti Curves stands out as the most interpretable vectorization method. Note that most vectorizations are computationally efficient and require minimal time compared to the computation of PDs.





  
*Hyperparameters:  Each vectorization method comes with its own set of hyperparameters that need to be carefully tuned to maximize performance. In all of them, the choice of thresholds is one of the key hyperparameters we discussed in <Ref>.

In addition to thresholds, many vectorization methods have tuning parameters that are used to adjust sensitivity to topological noise. Typically, topological features with short lifespans (i.e., those close to the diagonal in PDs) are regarded as noise. These features, represented as pairs {(b,d)} with short lifespans (d-b), are generally easy to detect, leading to many tuning parameters being tied to the lifespan.
For example, in persistence landscapes, the number of landscapes (m) specifies how many maxima are included in the representation. A higher number of landscapes provides a richer depiction of topological features but increases computational complexity. The first few landscapes (e.g., 1^st and 2^nd landscapes) emphasize features with the largest lifespans.
In silhouettes, the tuning parameter p in the weight function w_i = (d_i - b_i)^p determines the emphasis of the vectorization. A higher p (e.g., p ≥ 2) de-emphasizes features with short lifespans, while a smaller p (e.g., p ≤ 0.5) increases the importance of topological noise.

For persistence images, the spread (σ) controls the width of the Gaussian kernels applied to each persistence point. A smaller spread results in sharper, more localized features, whereas a larger spread yields smoother images. The resolution parameter sets the number of pixels in the persistence image, thus determining the output dimension of the vectorization. Higher resolution captures finer detail but at a higher computational cost. Similar to silhouettes, a weight function (e.g., linear) can be applied to weigh the contribution of each persistence point based on its significance, such as its lifespan.




 §.§ Integrating PDs to ML tasks 2: Using Neural Networks
 

While vectorization choices provide greater control over the model and data analysis, there are automated methods to avoid dealing with vectorization choices or hyperparameter tuning. These methods optimize the selection for downstream tasks by either finding the best vectorization within a large vectorization space or working directly with persistence diagrams as point clouds, allowing neural networks to handle the data.



  §.§.§ Vectorization with NN: Perslay
  PersLay, introduced by Carrière et al. <cit.>, is an automated vectorization framework for persistence diagrams, which decides the best vectorization method for the downstream task. The method leverages a neural network architecture that processes PDs via a specialized layer called the PersLay layer. This layer incorporates a variety of vectorization strategies, providing a unified and flexible approach to extracting topological features.

The PersLay layer is composed of several key components: 

    
  * Weighting Functions: These assign importance to points in the persistence diagram. Common choices include exponential weighting and Gaussian weighting, which can be learned during training.

  * Transformation Functions: These functions, such as linear transformations or learned neural networks, apply to the coordinates of the persistence points to encode meaningful geometric and topological information.

  * Symmetric Functions: After transformation, symmetric functions like sum, mean, or max aggregate the transformed points into a fixed-size vector, ensuring permutation invariance.


With this architecture, PersLay covers a wide range of known vectorization methods within its framework, i.e., Persistence Landscapes, Persistence Images, Silhouettes, Betti Curves, Entropy Functions and other Persistence Curves <cit.>. By incorporating various vectorization techniques, PersLay can adapt to diverse data types and topological features, making it a robust and versatile tool for integrating PH into ML workflows. For an alternative method involving slightly different generating functions for learnable vectorizations, see <cit.> by Hofer et al.
 
 



  §.§.§ PD Learning with NN: PointNet
  An alternative automated method to leverage PDs in ML tasks without hyperparameter tuning is to directly utilize PointNet <cit.> or similar neural network architectures to process and analyze point cloud data, e.g., DeepSets <cit.>, PointMLP <cit.>. Unlike conventional techniques that transform point clouds into regular grids or structured formats, these networks treat PDs as sets of points in ℝ^2. Perslay and PointNet offer fundamentally different strategies for incorporating PDs into ML models without requiring vectorization or hyperparameter optimization. The choice between these methods depends on the specific task, as there is no universally superior approach.
 




 §.§ Computational Complexity for Persistent Homology
 

 

The computational complexity of persistent homology (PH) is strongly influenced by the choice of filtration complex, as it is directly related to the size, |𝒦|, of the simplicial complex 𝒦, i.e., the number of simplices it contains. Early PH algorithms had a cubic complexity with respect to the number of simplices, i.e., 𝒪(|𝒦|^3) <cit.>. Subsequent advancements have reduced this exponent to w = 2.376, i.e., 𝒪(|𝒦|^w) <cit.>. 

Thus, the computational challenge boils down to the number of simplices in a given simplicial complex 𝒦. For a point cloud with N points, the worst-case scenario for the Rips or Čech complex involves up to 2^N simplices across all dimensions. However, in PH, the computation of d-dimensional topological features requires only simplices up to dimension d+1. Higher-dimensional simplices do not contribute to the calculation of H_d(𝒦). For instance, when focusing on 0- and 1-dimensional features, only simplices up to dimension 2 are needed. 

In a point cloud of size N, for a Rips complex, the number of 0-simplices (vertices) is N, the number of possible 1-simplices (edges) is N2 = N(N-1)/2, and the number of possible 2-simplices (triangles) is N3 = N(N-1)(N-2)/6. Recall that if we include all dimensions, the possible number of simplices would be 2^N.  Hence, increasing the dimension d for topological features thus significantly raises computational costs, which is why most ML applications limit the maximum dimension for simplices in the filtration (d+1)  to 2 or 3, which is enough to calculate up to 1- or 2-holes, respectively.

The discussion above primarily applies to Rips and Čech complexes, but PH computations can be significantly more efficient when using cubical complexes, which are commonly employed in image data. For 2D images, the time complexity of PH is approximately 𝒪(|𝒫|^r), where r ≈ 2.37 and |𝒫| represents the total number of pixels <cit.>. Practically, this implies that the computational cost of PH scales almost quadratically with the image size. For higher-dimensional images, alternative methods for efficient computation exist <cit.>.


In recent years, several works have been published to improve the scalability and computational efficiency of PH. One approach focuses on developing alternative methods for more efficient computation of persistence diagrams <cit.>, while another aims to sparsify datasets while preserving topological information <cit.>.



 §.§ Software Libraries for Persistent Homology
 


Several software libraries provide tools for computing persistent diagrams and vectorizations across different types of data structures, particularly point clouds and graphs (see <Ref>). Notable among these are GUDHI, DIONYSUS, and RIPSER, which have been compared in a benchmarking study <cit.>.  

Most TDA libraries require point cloud data as input, rather than image or network (graph) data. GUDHI can process images, but only after they are converted into a suitable format, such as a point cloud or a cubical complex. Giotto-TDA extends the capability to graphs through its VietorisRipsPersistence and SparseRipsPersistence for undirected graphs and FlagserPersistence for directed ones. However, a common challenge in graph-based TDA is the reliance on shortest path (geodesic) distances, which may not always capture the most relevant topological features. Giotto-TDA also supports the input of images and time series data, enabling the computation of various vectorizations. Beyond PH computations, the Persim library from Scikit-TDA offers a suite of tools for further analysis of persistence diagrams, including metrics like the Bottleneck distance and visualizations like Persistence Landscapes and Persistence Images, enhancing the interpretability of TDA results.

In terms of usage, most libraries (e.g., TDA in R) were language-specific bindings of the underlying GUDHI, Dionysus, and PHAT libraries, which were originally implemented in Matlab or C++ for efficiency and acted as the earliest and most established workhorses for TDA. In recent years, Python implementations like Giotto-TDA have become increasingly popular, and the widely used ML library Scikit-learn has contributed to this trend by creating a new set of tools called Scikit-TDA for TDA integration (<https://github.com/scikit-tda>). 

It is worth mentioning that since TDA has been developed and extended primarily by mathematicians, the popularity of R within this community has made the TDA R library particularly important. Mathematicians often release advanced models and methods that can be directly connected to TDA research in R first, leading to a delay in their implementation in other programming languages.







§ MULTIPARAMETER PERSISTENCE
 

Multiparameter Persistence (often referred to as Multipersistence) introduces a novel concept with the potential to significantly enhance the performance of single-parameter persistence. The concept, introduced in the late 2000s by Carlsson, Zomorodian, and Singh <cit.>, has since been actively investigated for real-world applications <cit.>.



r6cm


Generic Bifiltration 
1.!

_m1     ⊂     _m2     ⊂     …     ⊂     _mn 

∪          ∪          ∪          ∪

…     ⊂     …     ⊂     …     ⊂     …

∪          ∪          ∪          ∪

_21     ⊂     _22     ⊂     …     ⊂     _2n 

∪          ∪          ∪          ∪

_11     ⊂     _12     ⊂     …     ⊂     _1n 






 
For original PH, the term single persistence is applied because we filter the data in just one direction, _1⊂_2⊂…⊂_n. The filtration's construction is pivotal in achieving a detailed analysis of the data to uncover hidden shape patterns. For example, for graph setting, when utilizing a single function f:→ containing crucial domain information (e.g., value for blockchain networks, atomic number for protein networks), it induces a single-parameter filtration as described earlier. However, numerous datasets offer multiple highly relevant domain functions for data analysis. Similarly, in other data types, there are several ways to expand the single filtration into multifiltration to get finer information on the topological patterns hidden in the data. The idea of using multiple filtration functions simultaneously to obtain a finer decomposition of the dataset led to the suggestion of the MP theory as a natural extension of single persistence (SP).

 
 


 §.§ Multifiltrations
 

Multifiltration refers to filtrations constructed using multiple scale parameters. A typical example of a bifiltration, denoted as {_ij}, is illustrated in Table <ref>, where each row and column corresponds to a distinct filtration. This bifiltration is generated by simultaneously varying one scale parameter along the horizontal axis and another along the vertical axis. The specific construction of these filtrations depends on the data type, with additional parameters introduced based on the nature of the data. Below, we summarize common methods for different data types. While the explanation focuses on two parameters for simplicity, this approach generalizes to any number of parameters.


 


  §.§.§ Multifiltrations for Graphs.
 

In graph settings, a single filtration function results in a single-parameter filtration _1 ⊂…⊂_N =. However, employing multiple filtration functions enables a more detailed analysis of the data. For instance, two node functions f: → and g: →, which provide complementary information about the network, can be combined using Multiparameter Persistence to generate a unique topological fingerprint. These functions induce a multivariate filtration function F: →^2 defined by F(v) = (f(v), g(v)).

Next, we define sets of increasing thresholds {α_i}_i=1^m for f and {β_j}_j=1^n for g. Then, we have _ij = {v_r ∈| f(v_r) ≤α_i, g(v_r) ≤β_j}, which can be written as _ij = F^-1((-∞, α_i] × (-∞, β_j]). Let _ij be the subgraph of  induced by _ij, meaning the smallest subgraph of  generated by _ij. This setup induces a bifiltration of complexes {_ij| 1 ≤ i ≤ m, 1 ≤ j ≤ n}, which can be visualized as a rectangular grid of size m × n (see <ref>). For more details on applying multipersistence in graph settings, refer to <cit.>. Additionally, one can combine power filtration with filtration by functions by applying power filtration to each subgraph in the sequence induced by sublevel filtration via functions <cit.>. In <Ref>, we give details of the utilization of MP method on graph setting for computer-aided drug discovery.





  §.§.§ Multifiltrations for Point Clouds.
 

In the context of point clouds, a natural parameter to consider is the radius r, as discussed in <Ref>. However, a point cloud  often has regions of varying density, with some areas being very dense and others quite sparse. Additionally, a few outliers can significantly distort the topological signature due to the way it is constructed. To address this issue, it is common to use a density parameter in conjunction with the radius parameter when constructing the filtration. This parameter reflects the local density of points in the point cloud and helps identify regions of different densities, distinguishing features significant in dense regions from those in sparse ones.
To achieve this, we first need to compute local densities. For each point p ∈, we compute a density measure δ(p), which could be based on the number of points within a fixed radius r_d or kernel density estimation <cit.>.
 

Next, we define our threshold set {(d_i, r_j)}, where d_1 > d_2 > … > d_m corresponds to thresholds for density, and 0 = r_1 < r_2 < … < r_n = diam() are radius thresholds. We then define the multifiltration as follows. Using δ(p), we define a nested sequence of subsets _1 ⊂_2 ⊂…⊂_m =, where _i = {p ∈|δ(p) ≥ d_i}. For each i_0, we treat _i_0 as a separate point cloud and apply the Rips filtration process to it. Specifically, for each i_0, we construct the Rips filtration _i_01⊂_i_02⊂…⊂_i_0n. This gives a bifiltration {_ij} as in <Ref>, with _ij = _ij representing the Rips filtration for point cloud _i with radius parameter r_j. e.g., the first column _i1 = _i corresponds to radius r_1 = 0.

This multifiltration approach provides more detailed information by capturing topological patterns in both dense and sparse regions, offering a richer understanding of the data's underlying shape, especially when varying densities and scales are crucial for the analysis. For further details, see <cit.>.



  §.§.§ Multifiltrations for Images.
 

Similarly, in the image setting, if the image is a color image, one can easily employ all color channels at the same time, similar to the graph setting. In particular, a color image has three color functions, denoted as R (red), G (green), and B (blue). Thus, for every pixel Δ_ij, there exist corresponding color values: R_ij,G_ij,B_ij∈ [0,255]. To proceed, we establish a three-parameter multifiltration with parameters {α_m}_1^N_R, {β_n}_1^N_G, {γ_r}_1^N_B, where α_m,β_n,γ_r∈ [0,255] are threshold values for color channels R,G, and B respectively. By simply defining binary images _m,n,r={Δ_ij⊂|  R_ij≤α_m, G_ij≤β_n, B_ij≤γ_r}, we obtain a three parameter multifiltration of size N_R × N_G × N_B.
Similarly, for grayscale images, one can utilize height, radial, erosion, signed distance, or similar filtration methods <cit.> along with the color filtration to obtain a multifiltration for color images.





 §.§ MP Vectorization Methods
 

By computing the homology groups of the complexes in these multifiltrations, {_k(_ij)}, along with the induced inclusion maps, we obtain the corresponding multipersistence module, which can be visualized as a rectangular grid of size m× n. The goal here is to track k-dimensional topological features through the homology groups {_k(_ij)} within this grid. However, as explained in <cit.>, technical challenges rooted in commutative algebra prevent us from transforming the multipersistence module into a mathematical structure like a "Multipersistence Diagram." The key issue is that for any k-dimensional hole in the multifiltration, we cannot directly assign a birth and death time due to the partial ordering within the module. For example, if the same k-hole σ appears at _1,4 and _2,3, neither (1,4) ⊀ (2,3) nor (2,3) ⊀ (1,4), making it difficult to define a birth or death time for σ unless it exists in a perfectly rectangular region. This issue does not arise in single persistence, where any two indices are always comparable, i.e., i<j or j<i.

As a result, we do not have an effective representation of the MP module <cit.>. While these technical obstacles prevent this promising approach from reaching its full potential in real-life applications, in the past years, several approaches have been proposed to extract key insights from MP modules by skipping MP representations and directly going to the vectorization from multifiltrations, which we detail below. 



The simplest method to derive a vector (or tensor) from a multifiltration involves using the Hilbert function of the MP module, which is a generalization of Betti functions. Let {_ij} denote an m × n multifiltration associated with a given dataset . First, we construct the multipersistence module _k(_ij). Next, by computing the Betti numbers of each simplicial complex in the multifiltration, we obtain an m × n matrix (or tensor) β_k() = [β_k(ij)], i.e., β_k(_ij) = rank(_k(_ij)). These are also called multipersistence (or bigraded) Betti numbers. In graph representation learning, this simple vectorization of the MP module has demonstrated remarkable performance <cit.>.



While bigraded Betti numbers offer a straightforward vectorization method for multiparameter (MP) modules, they fail to capture crucial topological information, especially regarding the significance of dominant features or the presence of a few important topological features. More sophisticated techniques are necessary for effective MP vectorization. One widely used strategy involves the "slicing technique," which focuses on studying one-dimensional fibers within the multiparameter domain. A clearer understanding can be obtained by restricting the multidimensional persistence module to these single directions (slices) and applying single persistence analysis. 

In their work <cit.>, Carriere et al. explored this approach by considering multiple such slices, often referred to as vineyards, and summarizing the resulting persistence diagrams. Another significant advancement is found in multipersistence landscapes <cit.> by Vipond, which extends the concept of persistence landscapes <cit.> from single to multiparameter persistence. 

The vectorization of MP modules is a rapidly evolving area of research. Various techniques have recently been proposed in practical applications, e.g.,  Hilbert decomposition signed measure (MP-HSM-C) <cit.>, Effective MultiPersistence (EMP) <cit.>, Generalized Rank Invariant Landscape (GRIL) <cit.>, and Stable Candidate Decomposition Representations (S-CDR) <cit.>. Although many recent methods are burdened by high computational costs, the MP-HSM-C method introduced by Loiseaux et al. <cit.> stands out as the most efficient in terms of speed and performance among current MP vectorization techniques.
 







§ MAPPER
 

Next to Persistent Homology, another powerful method in TDA is the Mapper, which Singh, Memoli, and Carlsson introduced in the late 2000s <cit.>.  The Mapper stands out as an effective and versatile tool for extracting insights from high-dimensional datasets, which creates a simplified graph representation of data by combining ideas from algebraic topology and data visualization. While PH is mostly used in supervised learning, Mapper is commonly utilized in unsupervised settings. In the past decade, it has made critical contributions to various fields involving point clouds in high dimensional spaces, e.g., biomedicine <cit.>. This method works by projecting data onto a lower-dimensional space, clustering the points within overlapping intervals, and constructing a topological network that captures the essential structure of the dataset. Through this approach, Mapper can reveal hidden patterns, relationships, and the underlying shape of data, making it particularly valuable for unsupervised learning in fields such as genomics, neuroscience, and social network analysis. Therefore, it is also considered a smart dimension reduction technique, too (<Ref>).

 


 §.§ Mapper for Point Clouds
 
For a point cloud  in ^N and a real-valued function f: →, the Mapper algorithm provides a summary of the data by scanning the clusters in  in the direction dictated by the function f, which is commonly referred to as a filter function or lens <cit.>. The output of the Mapper algorithm is the Mapper graph, which is considered a meaningful summary of the data, representing clusters and relations between the clusters in the data. It has been applied in several contexts like diabetes subtype identification <cit.>, ransomware payment detection on blockchains <cit.>, and cancer genotyping from single-cell RNA sequence data (<Ref>).

The Mapper algorithm generates a graph-based summary of a high-dimensional point cloud. In this new Mapper graph, nodes represent clusters within the original point cloud, and edges indicate the interaction between these clusters. In other words, Mapper is a soft clustering approach where a data point may appear in multiple clusters; these clusters are then represented as nodes in the new Mapper graph.

Given a point cloud  in ^N, we first define a lens f:→. In general, the lens function can be a natural function derived from the data domain, or it may be computed using a dimensionality reduction technique. The proper selection of hyperparameters for Mapper is crucial, and we will shortly discuss these in <Ref>.

The lens function then decomposes the data into subregions through a cover of the image f()⊂, allowing for the identification of clusters within each subregion.  
Formally, we cover the image f()⊂ with a collection of open intervals, i.e. = {I_k}_k=1^n where ⋃_k I_k⊃ f(). The intervals I_k=(a_k,b_k) are indexed such that a_k < a_k+1 and typically b_k ≥ a_k+1, allowing for overlap between consecutive intervals.  Figure <ref>b shows 6 overlapping intervals ={I_1,…, I_6} for a toy example. 

Next, we consider the preimage of each interval, U_k=f^-1(I_k)⊂. In particular, the preimage of each interval refers to the set of points in the original high-dimensional space, such as the data points, mapped into a given interval by the lens function. For example,  in Figure <ref>, the union of all points from U_2 and U_3 and lower-part points from U_1 are the preimage of the interval I_2. We apply a clustering algorithm (e.g., k-means, DBSCAN) to each preimage U_k_0, dividing it into m_k_0 clusters {C_k_0l}_l=1^m_k_0. 

For each cluster C_kl, we create a node v_kl in the Mapper graph and let _k={v_k1, v_k2, …, v_km_k} represent the clusters in U_k⊂. Thus, =⋃_k _k forms the node set for the Mapper graph .  In <Ref>, there are 10 nodes in the graph for the 6 preimages: I_2, I_3, I_4, and I_6, each resulting in two clusters. Note that the number of clusters within a preimage can vary. For example, representing U_1 with a single cluster might result from using a density-based clustering algorithm (e.g., DBSCAN), which identifies dense regions as clusters while treating sparse regions as noise, thereby avoiding the creation of clusters in those less populated areas.



 
r3in



    < g r a p h i c s >

    
      Toy Example of Mapper. For a point cloud  (a), we define a lens function f:→ (b), and the induced covering defines a Mapper graph where nodes represent clusters and edges represent related clusters (c).     





After obtaining the nodes, we define the edges based on the intersections of the clusters. Note that as I_k∩ I_k+1≠∅, clusters in subsequent levels might have nontrivial intersections. Then, if the clusters C_i∩ C_j≠∅, we add an edge between the corresponding nodes v_i and v_j. i.e., ={e_ij| C_i∩ C_j≠∅} where e_ij represents the edge between the nodes v_i and v_j. It is also possible to define a weighted graph with weights ω_ij=#(C_i∩ C_j), the count of the points in the intersection. In this case, the weights reflect the strength of the interaction between the corresponding clusters, indicating how many data points are included in the overlapped area. The final Mapper graph gives a rough summary/sketch of the whole point cloud. 

The Mapper graph not only helps in identifying data clusters but also enables the selection of data points that are similar to a given subset by leveraging its structure <cit.>. In this context, nodes that are closer together on the Mapper graph tend to contain more similar data points. This aspect of cluster similarity based on their proximity in the Mapper graph is an area that remains largely understudied and offers significant potential for further exploration. The Mapper graph computations have been detailed in <Ref>. In <Ref>, we outline the utilization of Mapper in a real-life application, i.e., cancer genotyping from RNA sequencing.





  
*Common Filter Functions for Mapper. While we explain the Mapper algorithm for single-valued filter functions to keep the exposition simpler, in practice it is very common to use multivariate filter functions F:→^2. One of the most common filter functions used in applications is the Stochastic Neighborhood Embedding (t-SNE) <cit.>, and a variant of this, called Neighborhood Lens <cit.>. 

Similarly, it is common to use other dimension reduction techniques as filter map F:→^2, e.g., Isomap, PCA, or UMAP <cit.>. Then, for an open covering ={I_ij} of F() in ^2, one can repeat the process of U_ij=F^-1(I_ij), and assign a node w for each cluster in U_ij. Again, assign an edge between the nodes when the corresponding clusters have nontrivial intersections. Note that if the domain of the data offers good filter functions f,g:→,  one can also utilize multivariate filter functions F(x)=(f(x),g(x))∈^2 by combining both functions in the process. 




  §.§.§ Hyperparameters for Mapper
 
For a point cloud , four main parameter choices must be made to obtain its Mapper graph. 





    
  * The lens function. If available, the lens function can be derived directly from the data domain, which in turn greatly enhances the interpretability of the model. In practice, when appropriate lens functions are not readily available from the data domain, dimensionality reduction techniques like t-SNE, UMAP, or PCA are commonly used to create lenses f: 𝒳→ℝ^2. 
    
  * Clustering method. For each subset f^-1(I_k), the clustering method determines the nodes of the Mapper graph. The most common methods employed are DBSCAN and k-means, with their hyperparameters controlling the granularity of the resulting clusters.
    For k-means, it is common to use <10 for clusters (e.g., see <cit.>).
    
    
  * Resolution and Overlap. The resolution and overlap are the parameters for the covering {I_k}_1^n of f()⊂, i.e., ⋃_k=1^n I_k⊃ f(). 
    The resolution, n, is the number of bins (intervals) {I_k} to cover f(). The overlap is the percentage of overlaps of these intervals, i.e. |I_k∩ I_k+1|/|I_k|.  Note that increasing the resolution gives a finer summary by increasing the number of nodes in the Mapper network and making the clusters smaller (See <Ref>). On the other hand, increasing overlap adds more relation (edges) between the nodes (clusters).  In <Ref>, the resolution is 6, and the overlap is the fixed intersection amount (e.g., 20%) between the intervals I_k and I_k+1. Typically, it is common to use 10%-30% overlap and 10-20 intervals (e.g., Figure 10 of GraphPulse <cit.>).


 
 








 §.§ Mapper for Other Data Formats
  
 

While the Mapper approach is most commonly applied to point clouds, it can be adapted for other data formats, such as graphs and images. For graphs, this adaptation can be seen as a form of graph coarsening/skeletonization, where the goal is to summarize clusters of nodes within the graph.

In <cit.>, Bodnar et al. present a natural extension of the Mapper algorithm to the graph setting. Given a graph =(,), we start with a filter function f:→ and define a cover = {I_k}_k=1^n, where ⋃_k I_k ⊃ f(). The set _k = f^-1(I_k) represents a subset of nodes, and _k is the induced subgraph, i.e., _k = (_k,_k) where _k ⊂ are the edges connecting pairs of vertices in _k. Instead of clustering, we directly use the components (connected subgraphs) in _k as the nodes of the Mapper graph.
Specifically, if _k_0 = ⋃_l=1^m_k_0_kl has m_k_0 connected subgraphs, we define the node set _k = {c_k_01, …, c_k_0m_k_0}. For example, if _k has five connected subgraphs, we define five nodes in the Mapper graph, each representing one connected subgraph in _k. Let = {w_1, w_2, …, w_M} be the collection of all such nodes where M = ∑_k=1^n m_k. This defines the node set of the Mapper graph of . 

The edge set is defined similarly: If the subgraphs _kl and _(k+1) l' share a common node, then an edge is added between the corresponding nodes in . In summary, in the original construction, we replace the point cloud with the node set of , and the clustering algorithm is directly derived from the graph structure.  For this approach, one can utilize the common filtering functions f:→ from persistent homology (e.g., degree, closeness, betweenness, eccentricity) or can use other popular functions from graph representation learning, e.g., eigenfunctions of graph Laplacian, pagerank <cit.>. 


 
Alternatively, if the graph = (, , ) includes node attributes, a Mapper graph can be directly derived from the set of node attributes . In particular, we treat the node attribute vectors = {_i}_1^n with _i ∈^N as a point cloud ⊂^N and apply the Mapper algorithm as before. This Mapper graph provides a visual summary of the node attribute space, independent of the graph's structure. In other words, the resulting Mapper graphs summarize the graph based on node attributes alone, without incorporating information about node neighborhoods in the original graph. An effective utilization of this approach can be found in <cit.>.




The Mapper algorithm on images is analogous to its application in graphs, where clustering is derived from the original image <cit.>. For a given r × s image , let f: → represent the color values of the pixels. We define a cover = {I_k}_k=1^n such that ⋃_k I_k ⊃ f(). A node is defined for each connected component in f^-1(_k), and an edge is created between nodes if the corresponding components have a nontrivial intersection. This Mapper graph summarizes the interaction between different color regions in  as a graph.




 §.§ Computational Complexity of Mapper

The computational complexity of TDA Mapper varies depending on the specific choices of filter functions, clustering algorithms, and dimensionality reduction techniques.  


The first step in TDA Mapper involves applying a filter function. If t-SNE is used as the filter, the computational complexity of the Barnes-Hut t-SNE variant is 𝒪(N log N · D) <cit.>, where N is the number of data points (e.g., for graphs, the number of nodes |𝒱|) and D is the dimensionality of the data (e.g., for graphs, the number of node attributes). The complexity of t-SNE is largely driven by the computation of pairwise distances in a D-dimensional space, followed by optimizing the low-dimensional embedding. While dimensionality D influences the time required for these computations, the overall complexity is typically dominated by the number of data points N.

After t-SNE, the data is covered by overlapping intervals, and clustering is performed within each interval. The complexity of this step depends on the clustering algorithm used. For example, if a clustering algorithm with complexity 𝒪(N^2) is used, this step may add significant computational cost. Finally, the Mapper complex is constructed by connecting clusters, which typically has a lower complexity, often 𝒪(N) to 𝒪(N log N), depending on the number of clusters and the method used to connect them.

The overall computational complexity of TDA Mapper when using t-SNE is dominated by the complexity of t-SNE, which is 𝒪(N log N · D).  The subsequent steps in the Mapper pipeline add to this complexity, particularly the clustering step. When using k-means clustering, the complexity is typically 𝒪(k · N · D · T), where k is the number of clusters, N is the number of data points, D is the dimensionality of the data, and T is the number of iterations. Therefore, the overall complexity can be expressed as 𝒪(N log N · D) + 𝒪(k · N · D · T)=𝒪(N log N). 



 §.§ Software Libraries for Mapper
 

Several software libraries provide tools for constructing and analyzing Mapper (see <Ref>). Among these, Kepler-Mapper is a popular choice for Python users, particularly within the Scikit-TDA ecosystem, offering HTML outputs that are ideal for creating shareable, interactive visualizations. This makes it a practical tool for users who need to present their findings to non-technical users. 

Giotto-TDA integrates well with the Python environment and offers parameter visualization capabilities, making it an excellent choice for those needing to iteratively refine their Mapper construction. Its interactive features are especially helpful for real-time exploration of the impact of different filter functions and covering parameters.

For users prioritizing performance over interactivity, tda-mapper-python offers a faster computation engine, making it well-suited for large datasets or scenarios where rapid iteration is necessary. However, the lack of interactive features means that users will need to rely on external tools for visual exploration.

TDAmapper in R is an excellent choice for users already embedded in the R ecosystem. It’s a robust standalone solution that is particularly useful for those who prefer the extensive statistical and data manipulation capabilities available in R. However, the lack of interactivity might require additional effort for visualization and exploration.

Finally, Mapper Interactive is designed for those who value hands-on engagement with the data. It offers a highly interactive Python-based environment, making it ideal for exploratory data analysis where immediate feedback and manipulation of the Mapper graph are crucial.

In practice, the choice of library typically depends on the specific requirements of your project, such as whether you prioritize interactivity, computational speed, or seamless integration with other tools in your data analysis workflow. For example, Kepler-Mapper and Giotto-TDA are ideal for exploratory analysis and visualization, while tda-mapper-python and TDAmapper are more suitable for handling large-scale computations and integrating with ecosystems like Python and R, respectively.

 


 




§ APPLICATIONS
 

In this section, we will provide examples of the applications of TDA methods in machine learning. We provide four examples from published papers and outline their model and how they applied TDA in their project.



 §.§ PH for Point Clouds: Shape Recognition
 
r3in



    < g r a p h i c s >

    Shape Recognition via PH. In <cit.>, the authors obtained 100 point clouds for each animal by subsampling the surfaces above and analyzed their persistence landscapes. On the right, they present the 95% confidence bands for the persistence landscapes of these 400 point clouds, demonstrating that each animal's persistence landscape exhibits significantly distinct characteristics. The figure is adapted from <cit.>.     
    


In this section, we present an illustrative example of the utilization of PH in shape recognition, based on the work <cit.>, where Chazal et al. demonstrated the effectiveness of PH in shape recognition. One specific example involved four different animals: an elephant, flamingo, lion, and camel, each representing a unique shape in ^3, where each shape is normalized to have a diameter of one. They began by selecting 500 random points from the surface of each animal, as shown in <Ref>, resulting in the point clouds _E, _L, _F, and _C for the elephant, lion, flamingo, and camel, respectively. Then, they created 100 subsamples of 300 points each from these point clouds, resulting in 400 different point clouds: E_i ⊂𝒳_E, L_i ⊂𝒳_L, F_i ⊂𝒳_F, and C_i ⊂𝒳_C for 1 ≤ i ≤ 100.

 

Next, they performed Rips filtration on each point cloud {E_i,L_i,F_i,C_i}_i=1^100 (see <Ref>) and computed the corresponding persistence diagrams for dimension one. By vectorizing these persistence diagrams, they obtained the persistence landscapes for each point cloud, resulting in 400 persistence landscape functions: 
{λ(E_i), λ(L_i), λ(F_i),  λ(C_i)}_i=1^100.
In <Ref>, they present the 95% confidence bands for the true average landscapes (bold curves) of each class. e.g., blue confidence band correspond to 100  persistence landscapes coming from lion point clouds {λ(L_1), …, λ(L_100)}. The narrowness of these confidence bands indicates that PH techniques provide robust shape-embedding methods that are minimally affected by noise. Moreover, the distinct confidence bands for each class underscore the effectiveness of PH in shape recognition. We also note that recent work by Türkeş et al. <cit.> studied the effectiveness of PH for shape recognition in different settings.






 §.§ PH for Graphs: Crypto-Token Anomaly Forecasting
 

Three primary types of problems dominate graph representation learning: graph classification, node classification, and link prediction, as outlined by Hamilton et al. <cit.>. These tasks encompass a broad range of real-life applications, including brain connectivity networks, molecular property prediction, recommender systems, fraud detection and transaction networks. In a related study, Li et al. <cit.> employ persistent homology to extract effective feature vectors from weighted and directed Ethereum crypto-token networks, modeled as temporal graphs. The approach is tailored for graph anomaly prediction, a specialized form of graph classification.  
 

 



The research problem is set as a prediction task where the authors aim to predict whether the absolute price return of an Ethereum token will change significantly beyond a predefined threshold |δ| > 0 within the next h days. This involves analyzing the token's transaction network and its price fluctuations over multiple discrete time graph snapshots. The price information is sourced from external ground truth data, as token prices are determined by trading on blockchain exchanges.

For each snapshot, the authors use transferred token amounts on edges to define distances between adjacent vertices. These functions help  establish the similarity between nodes:
ω_uv = [1 + α·A_uv - A_min/A_max - A_min]^-1 
where A_uv represents the amount of tokens transferred between nodes u and v, and A_min and A_max are the minimum and maximum transferred amounts, respectively. The authors set the parameter α = 9 to map these weights to the interval [0.1, 1], thus standardizing the weight values.
 

r3in


    [Tronix ]
        
    < g r a p h i c s >

	[Power Ledger ]
        
    < g r a p h i c s >

    Betti pivots. Comparative functional summaries of the Tronix and Power Ledger token networks over seven days, displaying daily Betti-1 values across various scaling parameters. Each plot captures the evolution of topological features, with the central red lines indicating the Betti pivots. These pivots represent the most stable or 'normal' network structures, providing insights into network behavior and potential anomalies driven by underlying transaction activities. The figure is adapted from <cit.>.
    
 


By treating the edge weights (similarity measures) as distances between nodes, i.e., d(u,v) = ω_uv, the authors construct a filtration of Rips complexes (<Ref>-ii) from the snapshot graphs. For non-adjacent nodes, the distance is defined by the shortest path in the graph. The main insight here is to extract topological patterns developed in the graphs at multiple scales. This process involves forming a simplicial complex where nodes are connected if their distance is within a threshold ϵ, reflecting their similarity based on the edge weights. Essentially, nodes with high similarity (those with large transaction amounts between them) are connected early in the filtration, while more distant nodes appear later. 

The authors then compute persistence diagrams from these filtrations and convert them into Betti functions. Additionally, they introduce new functional summaries of topological descriptors, namely Betti limits and Betti pivots, which track the evolution of topological features as the scale parameter ϵ changes over the snapshots. <Ref> illustrates two token networks and their corresponding functional summaries.



To identify which transaction networks indicate anomalous patterns, the authors use Modified Band Depth to assess how central or peripheral each network's topological descriptors are within the observed data set. Figure <ref> shows two token networks and their snapshot graphs as represented by functions within these figures.  Snapshot graphs with deeper Betti limits are considered more typical or central.

 
The authors integrate these novel topological features with conventional network summaries to predict price anomalies. Daily labels are assigned as anomalous based on significant price changes anticipated in the near future (e.g., in one or two days). To this end, the authors construct a predictive model by utilizing topological vectors they produce for token networks. <Ref> displays the accuracy metrics for their model, specifically for a prediction horizon of h=2 (two days). Achieving an average accuracy of 96% across ten token networks, they demonstrate the effectiveness of topological features in the anomaly forecasting task.

  


 

 



 §.§ PH for Images: Cancer Diagnosis from Histopathological Images
 

Our example will directly apply persistent homology methods to histopathological images. As noted in <Ref>, cubical persistence is the primary method for applying persistent homology in an image context. In Yadav et al. (2023) <cit.>, Yadav et al. successfully applied cubical persistence to analyze histopathological images. Specifically, for each image , the authors generated a topological feature vector β() and utilized standard ML methods on these vectors (image embeddings) for tumor classification.

r3in



    < g r a p h i c s >

    Eight Color Channels. Different color channels are used to create sublevel filtrations. The figure is adapted from <cit.>.    
    


Recall that for a given image X with r × s resolution, the first step is to create a filtration, which is a nested sequence of binary images X_n. A common method to create such a sequence is to use the color values γ_ij of each pixel Δ_ij⊂ X. Specifically, for a sequence of grayscale values (t_1 < t_2 < … < t_N), one obtains a nested sequence of binary images X_1 ⊂ X_2 ⊂…⊂ X_N such that X_n = Δ_ij⊂ X |γ_ij≤ t_n.




In <cit.>, the authors construct filtrations for cubical persistence by first extracting eight color channels γ^k_ij from histopathological images with 1 ≤ k ≤ 8, where each superscript k corresponds to one color channel. The first four channels come from the RGB color space: red, green, blue, and grayscale (the average of R, G, and B). Additionally, they utilize the HSV color space, which includes hue, saturation, value, and their average (see <Ref>). Each color channel defines a different filtration {_n^k}_n=1^N. In the paper, they set N=100.



Next, by using these filtrations, they obtain the corresponding persistence diagrams for dimensions 0 and 1. Then, by applying Betti vectorization to these persistence diagrams, they obtain 100-dimensional Betti vectors β⃗^k_m=[β^k_m(t_1) …β^k_m(t_100)] where k represents the color, and m=0,1 represents the dimension. Hence, β^k_0(t_n) is the number of components (Betti-0 number) in the binary image ^k_n, and β^k_1(t_n) is the number of holes (Betti-1 number) in the binary image ^k_n. Considering that there are eight color channels and two dimensions, one obtains 1600-dimensional vector β⃗() by concatenating these 16 vectors {β⃗^k_m}. By extracting extra features by utilizing local binary patterns (LBP) and Gabor filters, they produced another 800 and 400 dimensional features for each image, respectively. Then, they obtain a 2800-dimensional final vector β(). In particular, one can consider this as an image embedding method, and each image is realized as a point in ^2800. Next, to improve the performance by removing correlated features, they use a feature selection method for the downstream task and reduce the vector size to 500 dimensions. 

They studied five cancer types, namely bone, breast, cervical, prostate, and colon cancer. In <Ref>, they provide the median curves and confidence bands for each class for three cancer types. As discussed earlier, while Betti vectors are considered as a weak vectorization method in general, they can be highly effective when the signature comes from the quantity and distribution of small features like these histopathological images. They utilized benchmark datasets for each cancer type consisting of 20K to 60K histopathological images. In <Ref>, we give their results for various sets of feature vectors. Their results indicate that utilizing filtrations with multiple color channels can significantly improve the classification results in some cancer types, e.g., bone, breast, and colon cancers.














 §.§ Multiparameter Persistence: Computer Aided Drug Discovery
 

 

In this part, we outline an effective application of multiparameter persistence in computer-aided drug discovery (CADD), showcasing its use within the graph setting. For an application in the context of point clouds, see <cit.>.

Virtual screening (VS), a key technique in CADD, is used to identify potential drug candidates from a vast library of compounds that are most likely to bind to a specific molecular target. Demir et al. <cit.> employed a multiparameter persistence approach for virtual screening by framing it as a graph classification task. In this method, compounds are represented as graphs, with atoms as nodes and bonds as edges. They adopted a ligand-based approach, wherein a few positive samples are provided for a given protein target, and the goal is to screen the compound library to find compounds that are most similar to these positive samples. Essentially, this approach can be viewed as a topology-based graph ranking problem. While their approach is more technical, here we outline the key concepts of their method. 




 r0.25 
    

  

    
    < g r a p h i c s >

    Cytosine. Atom types are coded by their color: White=Hydrogen, Gray=Carbon, Blue=Nitrogen, and Red=Oxygen. The decimal numbers next to atoms represent their partial charges. The figure is adapted from <cit.>.
    
  

  

 

Their framework, TODD, generates fingerprints of compounds as multidimensional vectors (tensors), represented as a 2D or 3D array for each compound. The core idea is to simultaneously employ 2 or 3 highly relevant functions or weights (e.g., atomic mass, partial charge, bond type, electron affinity, distance) to obtain a multifiltration, which decomposes the original compound into substructures using these relevant chemical functions. As detailed in <Ref>, for a given compound =(,), they used atomic number A and partial charge P as node functions A:→ and P:→, as well as bond strength B as an edge function B:→, to define multifiltrations. An illustration of such a multifiltration for the compound cytosine <Ref> is provided in <Ref>. In this example, the functions are partial charge and atomic number. 
 



 
After constructing multifiltrations, they obtained compound fingerprints using a slicing technique. In particular, within each multifiltration ij, they took horizontal slices by fixing a row i_0 to get a single persistence filtration i_0j. For each horizontal slice i_0, they then generated a persistence diagram _k(i_0j). These persistence diagrams were vectorized using Betti and Silhouette vectorizations to produce 2D arrays, such as 𝐌β=[β(ij)] (Betti numbers of the clique complex ij).  



Next, using these 2D arrays, they applied two ML classifiers: Random Forest and ConvNeXt Vision Transformers. Both classifiers performed exceptionally well, surpassing all state-of-the-art models on benchmark datasets. In <Ref>, their results are shown for the DUD-E Diverse dataset, which comprises 116K compounds targeting eight proteins. The common metric in virtual screening is enrichment factor, which compares the proportion of active compounds found in the top-ranked subset of a model output to the proportion of active compounds in the entire dataset. EF_1%   represents the enrichment factor for the top 1% of the dataset. 

 
 












 §.§ Mapper: Cancer Genotyping from RNA Sequencing
 

In this section, we will explore a significant application domain for the Mapper algorithm: the analysis of single-cell RNA sequencing data. Recall that Mapper is particularly effective for analyzing high-dimensional point clouds by generating a low-dimensional graph summary that preserves local relationships (<Ref>). In the Mapper summary graph, the nodes represent clusters in the point cloud, and the edges between the nodes indicate that the corresponding clusters are nearby (interacting) in the high-dimensional space.



RNA sequencing data is crucial for cancer genotype analysis, though it presents significant challenges. The expression profile of a cell can be mathematically represented as a point in a high-dimensional expression space, where each dimension corresponds to the RNA level of a gene, and the dimension of the space is the number of expressed genes. Points that are close to each other in this space correspond to cells with similar expression profiles. The set of all possible tumors of a cancer type spans a subspace of the expression space. From an ML perspective, this is a highly sparse point cloud in a high-dimensional space. In general, cells are assigned to some cancer subtype (e.g., malignant, benign, etc.), and the aim is to understand the relation of these types in this high-dimensional space. 






r3in



    < g r a p h i c s >

    Genetic Alterations. Topological representation of the expression space for low-grade gliomas using Mapper. The structure reveals three main clusters corresponding to distinct glioma subtypes. One cluster, labeled IDH2mut, highlights the significance of the IDH2 gene in differentiating this subtype. This visualization emphasizes the distinct expression profiles that separate these glioma subtypes. This figure is adapted from <cit.>.

There are several significant works on cancer genotyping through the use of Mapper techniques <cit.>. In this paper, we will detail the methods in two works. In the first one <cit.>, Wang et al. conducted a comparative analysis of Mapper visualization techniques on RNA-sequencing data. They used the melanoma tumor cells dataset GSE72056, which contains 4,645 cells. Of these, 1,257 are malignant, and 3,256 are non-malignant, with some minority classes included. The dataset includes 23,686 expressed genes. From an ML perspective, this data can be represented as a point cloud = {x^i}_i=1^4645 in very high dimensional space ^23,686, where each dimension corresponds to a gene {γ_j}_j=1^23,686. In particular, for a cell x^i ∈ with x^i = (x^i_1, x^i_2, …, x^i_23,686), x^i_j represents the RNA level of gene γ_j in cell x^i. This is computed as x^i_j = log_2(1 + TPM_ij/10), where TPM_ij is the transcript-per-million (TPM) for gene γ_j in cell x^i. Notice that the number of dimensions far exceeds the number of points. From an ML perspective, this presents a highly challenging setup. In <Ref>, the authors illustrate different dimension reduction techniques on this dataset. In <Ref>, they show how the graph changes when the number of bins (resolution) is varied. 


In the second work <cit.>, Rabadan et al. applied the Mapper algorithm to identify somatic mutations relevant to tumor progression. They analyzed mutation and RNA expression data for 4500 genes across 4476 patients from 12 tumor types, including low-grade glioma (LGG), lung adenocarcinoma, breast invasive carcinoma, and colorectal adenocarcinoma. This study led to the identification of 95 mutated cancer genes, 38 of which were previously unreported. They identified three large expression groups with oligodendrogliomas (enriched for CIC and IDH2 mutations), IDH1-mutant astrocytomas (enriched for TP53 mutations), and IDH1-wild-type astrocytomas (enriched for EGFR mutations) (<Ref>). In their analysis, they used the neighborhood lens (<Ref>) as the filter function and scanned the resolution and gain parameter space to determine the optimal parameters (<Ref>).  














§ FUTURE DIRECTIONS
 
In the preceding sections, we covered the primary topological methods employed in machine learning. In this section, we will explore future directions for advancing these methods, aiming to enhance the practicality of TDA in ML. We will also discuss strategies for expanding the application of topological methods into new and emerging areas.
 



 §.§ Topological Deep Learning


While TDA is already an effective tool for feature extraction from complex data, recent research underscores its potential to augment deep learning models by providing complementary insights. Deep learning algorithms typically emphasize local data relationships, whereas TDA offers a global perspective, delivering supplementary information. 
Furthermore, TDA methods have also been integrated into deep learning algorithms to optimize its process by regulating the loss functions and other layers <cit.>. 
In recent years, there has been significant progress in integrating TDA with deep learning in domains such as image analysis <cit.>, biomedicine <cit.>, graph representation learning <cit.>, genomics <cit.>, cybersecurity <cit.> and time-series forecasting <cit.>. This is a rapidly emerging field, with several recent papers surveying the current state-of-the-art <cit.>.

 

 §.§ TDA and Interpretability


TDA, though often used as a feature extraction technique in machine learning, provides a powerful framework for enhancing model interpretability by capturing the inherent geometric and topological structures of data. Tools like PH and Mapper excel at identifying and quantifying features such as clusters, loops, and voids across multiple scales, unveiling intricate patterns that might escape traditional analytic methods <cit.>. By utilizing these topological insights into ML pipelines, TDA offers a unique perspective on decision boundaries, feature interactions, and model behavior. For instance, Mapper facilitates the visualization and interpretation of complex, high-dimensional datasets by projecting them into lower-dimensional topological spaces<cit.>. This transformation uncovers relationships within the data that are otherwise difficult to grasp, making it easier to interpret the reasoning behind certain predictions or classifications. Ultimately, by incorporating TDA, researchers and practitioners can clarify the decision-making processes of machine learning models, promoting greater transparency and trustworthiness in AI systems.



 §.§ Fully Automated TDA Models


While PH and Mapper have demonstrated effectiveness across several machine learning domains, their successful application still requires considerable expertise, such as hyperparameter tuning and selecting appropriate filtration functions. To make PH and Mapper more accessible to the broader ML community, there is a pressing need for end-to-end algorithms that automate these processes. For PH, this automation can be tailored to specific data formats, such as point clouds, images, or graphs. In each case, the fully automated algorithm would handle the selection of filtration functions, thresholds, vectorization methods, and other hyperparameters for optimal performance in downstream tasks. For instance, an automated PH model designed for graph data could automatically choose the best filtration function (learnable), appropriate threshold values, and vectorization techniques to achieve optimal results in a node classification task. Similarly, a fully automated Mapper algorithm that selects filtration functions, resolution, gain, and other hyperparameters would significantly enhance its usability and accessibility to the ML community. For PH, there are works for learnable filtration functions, which find the best filtration functions in graph setting <cit.>. On the other hand, For Mapper, there is promising progress in this direction with Interactive Mapper <cit.>. The broader adoption of TDA hinges on developing practical software libraries that streamline its use for various applications.



 §.§ Scalability of TDA methods


Another key barrier to the widespread adoption of TDA methods in various application areas is their high computational cost. Although TDA performs well on small to medium datasets, scaling these methods to larger datasets poses significant challenges due to computational demands. Several strategies have been introduced to alleviate these costs, varying by data type. For graph datasets, recent research <cit.> has proposed practical techniques that substantially lower the computational burden of PH. Similarly, for point clouds, recent works <cit.> present algorithms that improve the efficiency of PH computations. While these advances are promising, there remains a critical need to enhance the scalability of PH for large graph and point cloud datasets. In contrast, cubical persistence is already computationally efficient for image datasets and can be applied to large datasets without significant cost concerns. Even in the case of 3D image datasets <cit.>, PH provides a cost-effective alternative to deep learning approaches. Similarly, scalability issues are minimal for Mapper when hyperparameters are appropriately tuned.


 





§ CONCLUSION
 
In this tutorial, we have introduced the key concepts of topological methods, particularly focusing on persistent homology and the mapper algorithm, and demonstrated their practical application in machine learning tasks. By providing a clear and accessible roadmap, we aimed to equip readers with the tools and understanding necessary to integrate TDA techniques into their research workflows. The strength of these methods lies in their ability to capture and quantify intricate, multi-scale topological features that are often missed by traditional ML approaches.

As demonstrated through various case studies, including cancer diagnosis, shape recognition, and drug discovery, TDA offers unique insights and interpretability, which are increasingly valuable in today's era of complex data. The integration of persistent homology for feature extraction and mapper for intuitive data visualization opens new pathways for researchers to explore the underlying structures of their data.

Looking forward, the continued development of software libraries, scalable algorithms, and more interpretable models will play a crucial role in making topological machine learning even more accessible to a broader audience. The future directions we highlighted, such as topological deep learning and automated TDA models, provide exciting opportunities for further advancements in the field.

We hope this tutorial serves as a foundational resource for those new to TDA and inspires further exploration of topological methods in machine learning and beyond. By embracing these techniques, researchers can unlock novel insights and push the boundaries of what is possible in data analysis.


This work was partially supported by the National Science Foundation under grants  DMS-2202584, 2229417, and DMS-2220613 and by the Simons Foundation under grant # 579977.
The authors acknowledge the http://www.tacc.utexas.eduTexas Advanced Computing Center (TACC) at The University of Texas at Austin for providing computational resources that have contributed to the research results reported within this paper. 


ACM-Reference-Format








§ NOTATION TABLE








§ DATASET RESOURCES


As we develop this tutorial to introduce topological methods to the ML community, we also seek to highlight their real-world applications to the mathematics community. For those eager to gain hands-on experience with topological techniques using real-world datasets, Table <ref> lists the commonly utilized datasets. For task-specific datasets, you can explore the collections at <https://paperswithcode.com/datasets> and <https://www.kaggle.com>.







 

