Topology Data Analysis(Coursework of Mathematical Frontier)

1. Introduction

This report is related to persistent topology and Mapper in TDA and it is divided into 3 parts, book report, project of persistent topology and my slides of Mapper.

2. Book Report

2.1 背景介绍

在这篇读书报告中，我主要对拓扑数据处理进行汇报，其中包括了两类算法：Mapper以及Persistent Homology，在这一部分中，我将介绍为什么要使用Mapper以及Persistent Homology的使用背景。拓扑数据处理是数据处理以及可视化中的一个研究领域，这个领域主要对复杂的高维数据进行分析，而其基础思想是：数据具有一定的形状，而研究数据的形状以及从数据中提取不变的几何特征有助于分析数据中的关系以及模式。与传统的方法，比如主成分分析以及聚类算法，不同在于，拓扑数据分析能够更好的兼具高维性以及对数据形状，尤其是拓扑不变量的分析。

拓扑数据处理的关键要点在于：

（1）坐标不变性：由于拓扑数据处理主要对数据的拓扑性质进行分析，所以分析结果不依赖与对于数据的的坐标轴的选取。

（2）形变不变性：对于拓扑数据分析研究结构的拉伸等变形不会对结果产生影响。比如如果最后的结果是一个圆环，然后对这个圆环进行拉伸变成一个椭圆的话，这里拓扑中的环状结构并不会被改变。

（3）压缩不变性：这里对数据的压缩表示，使用更少的点表示原先的数据或者通过离散的点表示连续的数据，在这一个过程中，拓扑数据分析的研究对象不会发生改变。一个典型的例子是：一个圆通过离散的六个点以及六条边构成的六边形进行表示，而拓扑中的环的结构不会发生该表。

在后面几个部分中，我将一次对拓扑数据处理的基础知识，拓扑数据处理中的Persistent Homology以及拓扑数据处理中的Mapper进行总结。

2.2 拓扑数据处理的基础知识

2.2.1 Simplex 与 Simplicial Complex

2.2.1.1 Simplex

在$R^n$中，一个k维的Simplex由$R^n$上k+1个线性无关点组成，设点$v_0,…,v_k\in R^n$线性无关，则由$v_0,…,v_k$构成的Simplex定义为下面的拓扑空间：
$$
\begin{align}
{\sum_{i=0}^k t_i v_i |t_0+t_1+…+t_k=1,t_i\ge0}
\end{align}
$$
这个Simplex可以表示为$[v_0,…,v_k]$,而$t_i$是Simplex中的点$x=\sum_i t_i v_i\in[v_0,…,v_k]$的坐标。

在这里有一个概念“点列的线性性”，点列的线性无关性定义为：
$$
\begin{aligned}
&~~~~~~~~\text{points }v_0,…,v_k\in R^n\text{ is linearly independent }\triangleq\
&~~~~~~~~\text{vectors } v_0v_1,…,v_0v_k\in R^n\text{ are linearly independent }
\end{aligned}
$$
在下图中，分别表示了维数分别为0,1,2,3的Simplex

对于Simplex可以定义从属关系，需要注意的是，Simplex的从属关系不单是其对应的拓扑空间的从属关系，而是其构成定点的从属关系，即：
$$
\sigma=[x_1,…,x_n] \subset \tau=[y_1,…,y_m]\Leftrightarrow {x_1,…,x_n} \subset {y_1,…,y_m}
$$

在上面的定义中$\sigma$被称为$\tau$的面。

2.2.1.2 Simplicial Complex

Simplicial Complex是若干个单形的集合，并且这些单形需要满足一定的性质，Simplicial Complex严格的定义为：设X是一个有限集合，则定义在X上的Simplicial Complex是有X的子集构成的若干Simplex，这些Simplex满足要求

（1）对于任意的$x\in X$,$[x]\in K$

（2）如果$\sigma \in K$且$\tau \subset \sigma$,则$\tau \in K$

在Simplex与Simplicial Complex上定义了一系列的常用符号，包括i-chain $K_i$以及 $d_i \sigma$

（1）i-chain: $K_i={\sigma\in K | dim \sigma =i}$，即i-chain $K_i$表示所有由K中维数为i的单形组成的集合

（2） $d_i \sigma$：For any n-dimensional, simplex $\sigma={x_0<x_1<…<x_n}\in K_n$ and $0\le i \le n$, define $d_i\sigma$ to be the (n−1)-dimensional simplex in K given by omitting $x_i$ from $\sigma$.

由于在Simplicial Complex上定义运算需要用到Free Module的概念，下面我们先对Free Module的概念进行回顾，然后再引入Simplicial Complex上的运算以及Homology。

2.2.2 Free Module

2.2.2.1 R-Module

首先我们对R-Module进行定义：Suppose that R is a ring, and 1 is its multiplicative identity. A left R-module M consists of an abelian group (M, +) and an operation · : R × M → M such that for all r, s in R and x, y in M, we have :
$$
\begin{align}
&r\cdot (x+y)=r\cdot x+r\cdot y\
&( r + s ) ⋅ x = r ⋅ x + s ⋅ x (r+s)\cdot x=r\cdot x+s\cdot x\
&( r s ) ⋅ x = r ⋅ ( s ⋅ x )\ &(rs)\cdot x=r\cdot (s\cdot x)
1 ⋅ x = x . {\displaystyle 1\cdot x=x.}
\end{align}
$$
One may write $_{R}M$ to emphasize that M is a left R-module. A right R-module $M_R$ is defined similarly in terms of an operation $·$ : $M \times R \rightarrow M$.

2.2.2.2 Free Module

Free module是一个具有基的Module，其变现得和向量空间类似。首先，我们对R-Module上面的基进行定义：For a ring R and an R-module M, the set $E\subseteq M$ is a basis for M if:
(1) E is a generating set for M; that is to say, every element of M is a finite sum of elements of E multiplied by coefficients in R;
(2) E is linearly independent, that is, for every subset ${e_1,e_2,…,e_n }$ of distinct elements of E, $r_{1}e_{1}+r_{2}e_{2}+\cdots +r_{n}e_{n}=0_{M}$ implies that $r_{1}=r_{2}=\cdots =r_{n}=0_{R}$ (where $0_{M}$ is the zero element of M and $0_{R}$ is the zero element of R).

给定R-Module的基E后，这个Free R-module可以定义为通过E索引的若干R的直和：
$$
\begin{align}
{\displaystyle R^{(E)}=\bigoplus _{e\in E}R}.
\end{align}
$$

2.2.3 Homology

在定义Homology前，我们首先来看Simplicial Complex上的运算，然后基于该运算，对Homology Group进行讨论：

2.2.3.1 Simplicial Complex上的运算

1. 边界映射$\partial_n$
   Let R be a commutative ring, and X be a set and K be a simplicial complex on X. Define $C_n(K) := \bigoplus_{\sigma \in K_n} R$. $\partial_n$ is a map between free R-module $C_n(K)$ and $C_{n-1}(K)$:
   $$
   \begin{align}

   \partial_n(e_\sigma)=e_{d_0\sigma}-e_{d_1\sigma}+...+(-1)^n e_{d_n\sigma}
   

   \end{align}
   $$
   where $e_\sigma$ is the generator of the components of R in $\bigoplus_{\sigma \in K_n} R$ indexed by $\sigma \in K_n$

   简单地说，运算$\partial_n$将Simplex集合中的一个Simplex映射成了另一个Simplex集合中的另一个Simplex，只不过，Simplex集合上定义运算需要通过R-Module来进行定义（R-Module类型与线性空间中的系数域）。在数据分析中，通常设置R Module为域 $F_2$

   不同的$\partial_n$组合在一起形成了一个映射链
   $$
   \begin{align}

      \cdots \stackrel{\partial_{n+1}}{\longrightarrow} C_n(K) \stackrel{\partial_n}{\rightarrow} C_{n-1}(K) \rightarrow \cdots \rightarrow C_1(K) \stackrel{\partial_1}{\rightarrow} C_0(K) \rightarrow 0
   

   \end{align}
   $$

2. $\partial_n$的性质

The composition $\partial_n \circ \partial_{n+1}$ is trivial, which means that $Im \partial_{n+1} \subseteq Ker \partial_n$. And the elements of $Im \partial_{n+1}$ are called the boundaries and the elements of $Ker \partial_n$ are called cycles.

举一个例子：
$$
\begin{aligned}
&[v_1,v_2,v_3]\in C_2(K), \partial_2([v_1,v_2,v_3])=[v_1,v_2]+[v_2,v_3]-[v_1,v_3]\Rightarrow\ &\partial_1([v_1,v_2]+[v_2,v_3]-[v_1,v_3])=(v_1-v_2)+(v_2-v_3)-(v_1-v_3)=0
\end{aligned}
$$

3. $\partial_n$的计算：边界矩阵

   在实际的计算中，尤其是在计算机的计算中，通常将 $\partial_n$表示成矩阵的形式：考虑将Simplicial Complex的n-Chain以及(n-1)-chain的Free Module之间的边界映射(即线性映射)表示为一个矩阵D，D的每一列表示n-Chain中的每一个n维单形，D的每一行表示(n-1)-Chain中的每一个n-1维单形。$D[i,j]$等于1 ，如果第i个n-1维单形是第j个n维单形列的面，否则$D[i,j]$等于0。这里的矩阵D被称为边界矩阵。

   在计算的过程中，经常用到的两个量为$Im \partial_n$ 以及$Ker \partial_n$，为了计算这个量，通常需要将边界矩阵华为Smith标准型，Smith标准型如下所示，实际上即为非方阵的对角矩阵：

   

   设边界矩阵的Smith标准型中，设有p行，即$# K_{n-1}=p$;有q列，即$# K_{n}=q$;非零元素共有k列(行)，则$Im \partial_n = k$，而 $Ker \partial_n = q-k$

4. 化矩阵为Smith标准型的算法

   我们默认$K$的R-Module中取的R为$F_2$，即下面的运算都在$F_2$上进行。在将矩阵华为Smith标准型时，可以在矩阵上使用下面的操作：

   （1）在矩阵中交换任意两列或任意两行。

   （2）向另一个列添加一个列，或者将一行添加到另一行。

   该算法的具体流程为：

   def _reduce(x):
   	nonzero=False
   	for i in range(x,m):
   		for j in range(x,n):
   			if matrix[i,j]==1:
   				matrix[[x,i],:]=matrix[[i,x],:]
   				matrix[:,[x,j]]=matrix[:,[j,x]]
   				nonzero=True
   				break
   		if nonzero:
   			break
   	if nonzero:
   		for i in range(x+1,m):
   			if matrix[i,x]==1:
   				matrix[i,:] = np.logical_xor(matrix[x,:], matrix[i,:])
   		for i in range(x+1,n):
   			if matrix[x,i]==1:
   				matrix[:,i] = np.logical_xor(matrix[:,x], matrix[:,i])
   		_reduce(x+1)
   		return;
   	else:
   		return;

2.2.3.2 Homology Group

1. Homology Group $H_n(K, R)$

   The quotient $Ker \partial_n/Im \partial_{n+1}$ is called the homology group of K with coefficients in R and is denoted by $H_n(K, R)$. The isomorphism type of the R-module $H_n(K, R)$ does not depend on the ordering we have chosen on X. For our purposes we are mainly interested in the case R being a finite field $F_p$.

2. n-th betti数 $\beta_n(K,R)$

   The dimension of the R-vector space denoted by $H_n(K, R)$ is called the n-th Betti number of K with respect to R and is denoted by $\beta_n(K, R)$

2.2.3.3 Q&A

Q: 为什么Homology Group以及betti数如此重要。

A: The homology $H_n(K, R)$ is an algebraic representation of certain topological features that X possesses. Its dimension(which is called Betti numbers) represents a count of those features.
For example, $\beta_0(K, R)$ is the number of connected components of K, $\beta_1(K, R)$ is the number of certain loops and $\beta_2(K, R)$ of certain voids.

2.3 Persistent Homology

在上一节最后的Q&A中我们知道，在Homology Group $H_n(K, R)$ 表示了定义Complex K的集合X的某些拓扑性质，而这些拓扑性质也正好是我们所需要的。但是，在现实世界中，我们一般只能够得到一些相互无关的离散的点，而不能够得到这些点构成的Simplicial Complex等拓扑结构。而这些点通常是在一个图形上面采样得到的，而我们又期望获得这个图形的拓扑性质。这个时候，我们期望恢复原来图形上的拓扑结构。但是问题就在于，在采样的过程中，可能有较大的噪声，我们期望在我们恢复拓扑结构的过程中不受这些噪声的影响，这个时候我们可以采用一个persistent的想法。

persistent的想法的含义是，我们从最开始的离散点开始，一步一步的构造若干个更加复杂的Simplicial Complex，然后考虑这些Simplicial Complex的拓扑性质，如果某个拓扑性质在这个过程中能够长时间保留，那么这个拓扑性质就可以被认为是原来被采样图形的拓扑性质。这个Simplicial Complex需要满足从属结构，即早生成的Simplicial Complex 从属新生成的Simplicial Complex。首先，我们需要了解如何在一个离散点集上构造Simplicial Complex。其次，我们探寻通过persistent homology能够得到什么样的拓扑性质。最后我们讨论如何计算。

2.3.1 离散点集上构造Simplicial Complex

在离散点集上构造Simplicial Complex有常用的两种复型： Vietoris-Rips complex与 Cech complex（下面提到的subsets的含义是通过这个subset中的点构造Simplex）

1. Vietoris-Rips Complex: Define $V(X, \epsilon)$ to be a simplicial complex on the set X consisting of these subsets $\sigma \subseteq X$ where $d(x, y) <\epsilon$ for any $x,y ∈ \sigma$
2. Cech Complex: Define $C(X, \epsilon)$ to be a simplicial complex on the set X consisting of these subsets
   $\sigma \subseteq X$ for which there is $y \in X$ such that $d(x, y) < \epsilon$ for any $x \in \sigma$

在实际计算中，由于上面的构造需要生成较多的复型，所以，我们一般使用他们的witness版本，例如VR complex的witness版本为： Choose a subset $X_w \subseteq X$ of witnesses of X. Define $V_w (X, \epsilon)$, the witness version of the Vietoris-Rips complex, to be a complex on $X_w$ that consists of these $\sigma \subseteq X_w$ where for any $x, y \in \sigma$, there is $z \in X$ such that $d(x, z) < \epsilon$ and $d(y, z) < \epsilon$.

2.3.2 Persistent Homology的计算

1. Complex列与边界矩阵

   设 $X_\epsilon$ 为上一节定义的Vietoris-Rips Comple 或者 Cech Complex. 然后考虑递增实数列$\epsilon_1<\epsilon_2<…<\epsilon_i,…$, 可以Complex序列 $X_{\epsilon_1},…,X_{\epsilon_n}$,对于某一特定的维数n,可以得到这些Simplicial Complex上的Homology Group 序列
   $$
   \begin{align}

   H_n\left(X_{\varepsilon_1}, R\right) \rightarrow H_n\left(X_{\varepsilon_2}, R\right) \rightarrow \cdots \rightarrow H_n\left(X_{\varepsilon_i}, R\right) \rightarrow \cdots
   

   \end{align}
   $$
   考虑$\epsilon_1,…,\epsilon_m$,我们考虑$X_{\epsilon_1},…,X_{\epsilon_m}$上所有的Simplex，将这些Simplex进行排序，排序结果$\sigma_1,…,\sigma_k$需要满足

   1. 如果 $\sigma_i \subset \sigma_j$, 则$i<j$
   2. 如果$dim\sigma _i<dim\sigma _j$,则$i<j$

   排序完成后，考虑矩阵m维方阵D，D的每一行每一列均和Simplex $\sigma_1,…,\sigma_k$进行对应，如果$\sigma_i \subset \sigma_j$, $D[i,j]=1$。这个矩阵D也被我们称作边界矩阵

2. Persisent Homology

   （1）对边界矩阵D使用Reduction 算法

   （2）对Reduction 算法的结果使用ReadInterval 算法

   （3）对ReadInterval算法的结果进行可视化

Reduction 算法，ReadInterval 算法以及可视化我们在下面进行介绍

2.3.3 Reduction 算法 与 ReadInterval 算法

1. Reduction 算法

   在Reduction 算法中需要使用Low(j): 对于矩阵D, Low(j) 表示矩阵D的第j列最后一个非零元素的行数，如果该列全部为0，则Low(j)=-1。

   Reduction 算法具体为：

   for j = 1 to n do
   	while ∃j‘ < j with low(j‘ ) = low(j) ！= -1 do
   		add column j 0 to column j
   	endwhile
   endfor.

   Reduction 算法将边界矩阵化为上三角矩阵。Reduction 算法的C++实现为：

   static SparseMatrixLIL<int> ReduceBoundaryMatrix_Echelon(
       const SparseMatrixLIL<int>& boundary_m,
       SparseMatrixLIL<int>& memory_m
   ) {
       SparseMatrixLIL<int> result_boundary_m = boundary_m;
       memory_m = SparseMatrixLIL<int>::IndentityMatrix(result_boundary_m.row_num, result_boundary_m.col_num);
       for (int i = 0; i < result_boundary_m.row_num; i++) {
           bool end_loop = false;
           while (!end_loop) {
               end_loop = true;
               for (int j = 0; j < i; j++) {
                   int low_i = Low(result_boundary_m, i);
                   int low_j = Low(result_boundary_m, j);
                   if (low_i == low_j && low_i != -1) {
                       result_boundary_m.SetSparseRow(
                           result_boundary_m.GetSparseRow(i) ^ result_boundary_m.GetSparseRow(j)
                           , i
                       );
                       memory_m.SetSparseRow(
                           memory_m.GetSparseRow(i) ^ memory_m.GetSparseRow(j)
                           , i
                       );
                       end_loop = false;
                   }
               }
   
           }
       }
       return result_boundary_m;
   }

2. ReadInterval 算法

   在ReadInterval 算法中，我们需要从reduction matrix中读取出一个特征开始的时间以及结束的时间（这里的时间指的是对应的$\epsilon$)。该算法的流程为：

   从左只有读取reduction matrix的每一列，在读取的过程中需要记录三维组(k,s,e)的集合,k表示特征的维数。

   设现在读取到了第j列，第j列对应的Simplex维数为t：

   1. 如果Low(j)==-1, 即这一列全为0，则为集合添加元素（t, j,-1）
   2. 如果Low(j)！=-1，则如果在集合中存在元素的开始时间为Low(j),则把该元素修改为（t, Low(j),j）。否则，向集合中添加元素（t, Low(j),j）

   读取完成后，现在需要将集合中的元素进行替换：

   对于(k,s,e)，设第s个simplex第一次出现在Complex $X_{\epsilon_s}$,第e个simplex第一次出现在Complex $X_{\epsilon_e}$,将(k,s,e)替换为(k,$\epsilon_s$,$\epsilon_e$)， 现在(k,$\epsilon_s$,$\epsilon_e$)表示有一个k维特征，最开始出现在$\epsilon_s$，在$\epsilon_e$结束。

   该算法通过C++可以表示为：

   static void ReadIntervals(
       const SparseMatrixLIL<int>& boundary_m,
       std::vector<std::pair<int, int>>& interval
   ) {
       std::unordered_map<int, int> start_map;
       for (int i = 0; i < boundary_m.row_num; i++) {
           int low_i = Low(boundary_m, i);
           if (low_i == -1) {
               interval.push_back(std::pair<int, int>(i, -1));
               start_map.insert(std::pair<int, int>(i, interval.size() - 1));
           }
           else {
               auto iter = start_map.find(low_i);
               if (iter == start_map.end()) {
                   interval.push_back(std::pair<int, int>(low_i, i));
                   start_map.insert(std::pair<int, int>(low_i, interval.size() - 1));
               }
               else {
                   interval[iter->second].second = i;
               }
           }
       }
   }

2.3.4 可视化表示

对于ReadInterval 算法的结果${(t,\epsilon_s,\epsilon_t)}$，我们采用两种方法进行表示：

1. barcode

   在barcode中，使用x轴代表特征开始的时间和结束的时间，对于特征$(t,\epsilon_s,\epsilon_t)$,添加一个从$\epsilon_s$开始到$\epsilon_t$结束的条形码。不同维数的特征通过不同的颜色进行表示。在barcode中，越长的特征是越可能在原始形状上出现的。

2. persistent diagram

   横轴表示开始时间，纵轴表示结束时间。

black bars represent the 0-dimensional cycles or vertices, a red bar represents a one-dimensional cycle(loop), and the blue bars represent cycles in H2(trapped voids) We like the feature that last for a long time.

2.3.5 实例

在这一个例子中，我们将观察VR complex构造的过程以及从persistent diagram和persistent barcode中读取出的拓扑特征。

1. 离散图形如下所示

   

2. VR Complex构造如下所示

   

   

   

3. 可视化

   

Reference

[1] Illustrations of Data Analysis Using the Mapper Algorithm and Persistent Homology

[2] Persistent Homology — a Survey, Herbert Edelsbrunner and John Harer

3. Project

Following is the report for our porject.

3. Mapper

Following is the report for our porject.