拓扑数据分析实战：从点云到机器学习特征提取

1. 拓扑数据分析：从数学直觉到工程实践

如果你处理过点云、图结构或者高维的时间序列数据，可能有过这样的困惑：传统的统计特征（比如均值、方差）或者深度学习的卷积核，似乎总是抓不住数据里那种“形状感”。比如，一个由传感器点构成的环形分布和一个实心球状分布，在特征空间里可能拥有相似的统计分布，但它们的“洞”和“连接性”截然不同。这种对数据“形状”和“结构”的量化需求，正是拓扑数据分析（Topological Data Analysis, TDA）的用武之地。

简单来说，TDA是一套数学工具，它不关心数据点的精确坐标，只关心它们之间的连接关系所构成的“形状”。它的核心武器来自代数拓扑，特别是同调群理论，用来数数据里的“洞”（0维洞是连通分量，1维洞是环，2维洞是空腔等等）。而持久同调则是TDA的“时间机器”，它让我们能看到这些“洞”随着观察尺度（比如点与点之间的连接阈值）变化时，是如何“诞生”和“消亡”的，从而区分出稳定的结构特征和短暂的噪声。

这篇文章不是一篇数学教科书，而是一个实践者的经验分享。我会带你绕开最抽象的代数推导，直接切入核心概念和操作流程。我们将基于像ModelNet（3D形状分类）和UCR（时间序列分类）这样的真实基准数据集，手把手地展示如何用Python工具链（如GUDHI, scikit-tda）从原始数据中提取拓扑特征，并将其融入机器学习管道。你会发现，为你的数据加上一双“拓扑之眼”，往往能发现那些被传统方法忽略的、却至关重要的结构模式。

2. 核心概念拆解：从“形状”到“不变量”

在深入代码之前，我们必须建立清晰的直觉。TDA的核心思想是：将一堆离散的数据点（点云），通过某种规则连接起来，形成一个连续的几何对象（称为单纯复形），然后计算这个对象的拓扑不变量。

2.1 同调群与Betti数：给“洞”编号

想象你有一张渔网（数据点），拓扑学不关心网线有多粗（度量），只关心网眼（洞）的个数和类型。同调群H_k(X)就是用来系统化描述这些“洞”的代数结构。这里的k代表维度：

• H_0(X)：描述0维洞，即连通分量的数量。一个数据集被分成几个互不相连的“孤岛”。
• H_1(X)：描述1维洞，即环状结构的数量。像一个甜甜圈中间的洞，或者一个圆圈。
• H_2(X)：描述2维洞，即封闭空腔的数量。像一个实心球内部挖空的空洞。

Betti数β_k就是这些同调群的“秩”，你可以粗暴地理解为第k维“洞”的个数。β_0=3意味着数据有三个独立的簇；β_1=1意味着数据中有一个显著的环状结构。

实操心得：刚开始很容易把β_0和聚类数量混淆。但聚类算法（如DBSCAN）的结果依赖于参数，且边界模糊。而β_0是一个在给定连接尺度下的、严格的数学定义。它告诉你“在这个特定的连接阈值下，数据恰好有几个连通块”，这个结论是确定无疑的。

2.2 从点云到复形：如何构建“形状”

数据点本身是离散的，没有“形状”。我们需要定义规则把它们连起来。最常用的两种复形是：

1. Vietoris-Rips复形 (Rips Complex)：给定一个距离阈值ϵ。如果一组点中任意两点之间的距离都小于ϵ，那么就用一个单形（点、线段、三角形、四面体…）把它们连起来。

   • 优点：计算相对高效，只需要两两距离。
   • 缺点：可能会“填充”掉本应是空洞的区域。比如一个正方形的四个顶点，在ϵ大于对角线长度时，会被填充成一个实心的四面体（三角形），从而“杀死”了中间本应存在的1维洞（环）。

2. Čech复形 (Čech Complex)：给定一个距离阈值ϵ。以每个数据点为球心，以ϵ/2为半径画球。如果一组点对应的球的交集非空，就用一个单形连接这组点。

   • 优点：拓扑性质更优，能更准确地捕捉空洞。著名的神经定理保证了Čech复形的同伦型与这些球的并集相同。
   • 缺点：计算极其昂贵，需要检查所有高阶交集，在实际中（尤其是高维数据）几乎不可行。

工程中的取舍：由于计算复杂度，Rips复形是绝大多数实际应用的选择。虽然它不完美，但通过后续的持久同调分析，我们依然能有效地识别出稳定的拓扑特征。GUDHI等库对Rips复形的计算有高度优化。

2.3 持续同调：捕捉特征的“生命周期”

单一的阈值ϵ是武断的。一个环在小尺度下可能看不见（点太分散），在中等尺度下出现，在大尺度下又被填充消失。持续同调的核心思想就是：让阈值ϵ从一个最小值连续增长到最大值，观察每个拓扑特征（洞）的“生”与“死”。

这个过程称为过滤 (Filtration){K_i}。我们得到一系列嵌套的复形：K_0 ⊆ K_1 ⊆ ... ⊆ K_n。

• 诞生 (Birth)：当一个新的k维洞在复形K_b中出现时，记录其诞生尺度b。
• 死亡 (Death)：当这个洞在复形K_d中被更高等的单形“填充”时，记录其死亡尺度d。
• 持久性 (Persistence)：d - b。持久性越长的特征，越可能是数据中真实的信号；持久性短的，很可能是噪声。

2.4 可视化：持续图与持续条形码

如何表示这些(birth, death)对？

1. 持续条形码 (Persistence Barcode)：每个特征用一条横线段表示，左端点在birth，右端点在death。所有线段平行排列。一眼就能看出哪些特征“长寿”。

   • 优点：非常直观，比较不同维度特征的持久性很方便。
   • 缺点：当特征很多时，会显得杂乱。

2. 持续图 (Persistence Diagram)：每个特征用二维平面上的一个点(birth, death)表示。由于洞总是在诞生之后才死亡，所有点都位于对角线y=x的上方。点到对角线的垂直距离就是其持久性。

   • 优点：空间利用率高，便于观察点的分布模式。是进行统计分析（如计算Wasserstein距离）的标准形式。
   • 缺点：不同维度的特征需要用不同颜色或形状区分。

注意事项：在持续图中，有些特征可能“永不死亡”（例如，数据整体构成的连通分量在最大尺度下也不会消失）。对于这些特征，其死亡时间记为∞。在实际绘图和计算时，通常用一个远大于所有有限死亡时间的数值（如np.inf或一个很大的数）来替代，并在后续处理中特殊对待。

3. 实战流程：从数据到拓扑特征向量

理论说得再多，不如一行代码。下面我们以经典的ModelNet10数据集（3D点云形状分类）为例，展示一个完整的TDA特征提取流程。我们将使用gudhi和ripser这两个强大的Python库。

3.1 环境准备与数据加载

首先，确保环境就绪。gudhi功能全面但安装稍麻烦（依赖C++库），ripser轻量且接口简单，特别适合快速计算Rips持续同调。

# 安装核心库 pip install numpy scikit-learn pip install gudhi # 可能需要处理系统依赖，如CGAL、Eigen等 # 或者，对于快速上手，更推荐： pip install ripser pip install persim # 用于处理持续图的可视化和距离计算 pip install gtda # scikit-tda，提供了类似scikit-learn的管道接口

我们使用ripser进行演示，因为它最容易跑通。

import numpy as np import matplotlib.pyplot as plt from ripser import Rips from persim import plot_diagrams import requests import io import tarfile # 1. 下载并加载一个ModelNet10样本（这里以离线加载一个.npy文件为例） # 假设我们有一个本地文件 'chair_sample.npy'，形状为 (N, 3) point_cloud = np.load('chair_sample.npy') # 例如，1000个点 print(f"点云形状: {point_cloud.shape}") print(f"前5个点:\n{point_cloud[:5]}")

3.2 计算持续同调

使用ripser计算点云的持续同调。我们通常关注0维和1维特征，因为2维及以上计算量剧增，且在许多数据中意义不明显。

# 初始化Rips计算器 rips = Rips(maxdim=1) # maxdim=1 表示计算H_0和H_1 # 计算持续同调 diagrams = rips.fit_transform(point_cloud) # diagrams 是一个列表，diagrams[0]是H_0的特征，diagrams[1]是H_1的特征 print(f"H_0 特征数量: {len(diagrams[0])}") print(f"H_1 特征数量: {len(diagrams[1])}") # 可视化持续图 plt.figure(figsize=(10, 5)) plt.subplot(121) plot_diagrams(diagrams, show=False) plt.title("Persistence Diagram") # 可视化持续条形码 plt.subplot(122) rips.plot(diagrams, show=False) plt.title("Persistence Barcode") plt.tight_layout() plt.show()

关键参数解析：

• maxdim：计算到第几维同调。maxdim=2会计算H_0, H_1, H_2。
• thresh：距离阈值的最大值。默认是无穷大，但设置为一个合理值（如点云直径的倍数）可以显著加速计算。
• coeff：同调群计算的系数域，默认是素数域2。使用coeff=2是最快且最稳定的，它相当于在模2运算下考虑单形（有/无）。有时为了捕捉更精细的拓扑信息（如定向），会使用coeff=0（整数域），但计算更复杂。

实操心得：对于含有数万个点的点云，直接计算Rips复形是灾难性的。标准预处理步骤是：

1. 下采样：使用最远点采样 (Farthest Point Sampling) 或体素网格下采样，将点云减少到1000-2000个点，同时尽量保持其拓扑结构。
2. 计算距离矩阵：ripser支持传入预计算的距离矩阵D，rips.fit_transform(D, distance_matrix=True)。如果你有自定义的度量（如测地距离），这非常有用。
3. 设置thresh：观察点云的大致范围，设置thresh为最大点间距的1.5-2倍，通常足以捕捉主要特征，并能大幅剪枝不必要的计算。

3.3 从持续图到机器学习特征

持续图（或条形码）本身是点集或线段集，不能直接输入SVM或随机森林。我们需要将其“向量化”。以下是几种主流方法：

方法一：持久性统计量 (Persistence Statistics)最简单直接的方法。对某一维度的所有特征(b, d)，计算一组统计量：

• 各特征的持久性pers = d - b
• 所有持久性的均值、方差、熵、最大值等。
• 所有诞生时间b的统计量。
• 所有死亡时间d的统计量。

def persistence_statistics(diagrams, homology_dim=1): """ 计算持续图的统计特征。 diagrams: ripser返回的图列表 homology_dim: 计算哪一维同调的特征统计 """ diag = diagrams[homology_dim] if len(diag) == 0: return np.zeros(9) # 返回零向量以防无特征 # 移除无限持续的特征（死亡时间为inf） finite_diag = diag[diag[:, 1] != np.inf] if len(finite_diag) == 0: return np.zeros(9) births = finite_diag[:, 0] deaths = finite_diag[:, 1] persistences = deaths - births stats = [] # 持久性统计 stats.append(persistences.mean()) stats.append(np.median(persistences)) stats.append(persistences.max()) stats.append(persistences.std()) # 诞生/死亡时间统计 stats.append(births.mean()) stats.append(deaths.mean()) # 持久性分布的熵（粗糙估计） hist, _ = np.histogram(persistences, bins=10, density=True) hist = hist[hist > 0] entropy = -np.sum(hist * np.log(hist)) stats.append(entropy) # 特征数量 stats.append(len(finite_diag)) # 持久性总和 stats.append(persistences.sum()) return np.array(stats) # 为H_0和H_1分别计算特征向量 feat_h0 = persistence_statistics(diagrams, homology_dim=0) feat_h1 = persistence_statistics(diagrams, homology_dim=1) # 拼接成最终特征向量 topological_feature_vector = np.concatenate([feat_h0, feat_h1]) print(f"拓扑特征向量维度: {topological_feature_vector.shape}") print(f"特征值: {topological_feature_vector}")

方法二：持久性图像 (Persistence Images)将持续图转化为一个二维灰度图像。过程如下：

1. 将每个点(b, d)转换为(b, persistence)，其中persistence = d - b。这相当于把图沿对角线“旋转”了一下。
2. 将每个点用一个高斯核函数（或其他核函数）表示，其权重可以是持久性本身或其他函数。
3. 在指定的像素网格上对所有这些高斯核进行叠加积分，得到一个二维矩阵（图像）。

from persim import PersistenceImager # 创建持久性成像器 pimgr = PersistenceImager(pixel_size=0.1, birth_range=(0, 2), pers_range=(0, 2)) # 设置参数：像素大小，birth轴范围，persistence轴范围 pimgr.fit(diagrams[1]) # 以H_1的图为例进行拟合（确定范围） # 转换持续图为图像 persistence_image = pimgr.transform(diagrams[1]) # 可视化 pimgr.plot_image(persistence_image) plt.title("Persistence Image (H1)") plt.show() # 将图像展平作为特征向量 image_vector = persistence_image.flatten()

方法三：拓扑向量 (Betti曲线/持久性景观)

• Betti曲线：对于每个尺度ϵ，计算该尺度下的Betti数β_k(ϵ)，得到一个关于ϵ的函数。对其进行采样或插值得到向量。
• 持久性景观 (Persistence Landscape)：一种将持续图转化为一系列函数（景观）的方法，具有良好的数学性质（属于希尔伯特空间），便于进行统计分析。

注意事项：选择哪种向量化方法取决于任务和数据。统计量最简单，可解释性强，但可能丢失空间分布信息。持久性图像保留了更多空间信息，适合与CNN结合，但引入了像素大小、带宽等超参数。持久性景观理论性质好，但计算稍复杂。我的经验是，对于初步尝试，从持久性统计量开始，它通常能提供一个不错的基线。

3.4 构建机器学习管道

现在，我们可以将拓扑特征作为传统特征工程的补充，构建分类器。

from sklearn.model_selection import train_test_split from sklearn.ensemble import RandomForestClassifier from sklearn.preprocessing import StandardScaler from sklearn.pipeline import Pipeline from sklearn.metrics import accuracy_score import glob # 假设我们有一个数据加载函数，返回点云列表和标签列表 def load_modelnet_samples(data_path, class_name, num_samples_per_class=50): # 简化示例：加载某个类别的.npy文件 file_pattern = f"{data_path}/{class_name}/*.npy" file_list = glob.glob(file_pattern)[:num_samples_per_class] point_clouds = [np.load(f) for f in file_list] labels = [class_name] * len(point_clouds) return point_clouds, labels # 加载两个类别做二分类示例 chair_pcs, chair_labels = load_modelnet_samples('./ModelNet10', 'chair', 50) table_pcs, table_labels = load_modelnet_samples('./ModelNet10', 'table', 50) all_pcs = chair_pcs + table_pcs all_labels = [0]*50 + [1]*50 # 0 for chair, 1 for table # 提取拓扑特征 def extract_tda_features(point_clouds, maxdim=1): features = [] rips = Rips(maxdim=maxdim, thresh=2.0) # 设置阈值加速 for pc in point_clouds: # 可选：下采样 pc = downsample(pc, num_points=1024) diagrams = rips.fit_transform(pc) feat_h0 = persistence_statistics(diagrams, 0) feat_h1 = persistence_statistics(diagrams, 1) feat_vector = np.concatenate([feat_h0, feat_h1]) features.append(feat_vector) return np.vstack(features) X_tda = extract_tda_features(all_pcs) y = np.array(all_labels) # 划分数据集 X_train, X_test, y_train, y_test = train_test_split(X_tda, y, test_size=0.2, random_state=42) # 构建管道：标准化 + 分类器 pipeline = Pipeline([ ('scaler', StandardScaler()), ('clf', RandomForestClassifier(n_estimators=100, random_state=42)) ]) # 训练与评估 pipeline.fit(X_train, y_train) y_pred = pipeline.predict(X_test) accuracy = accuracy_score(y_test, y_pred) print(f"仅使用TDA特征的分类准确率: {accuracy:.4f}")

与其它特征融合：拓扑特征很少单独使用。通常与几何特征（如点云的PFH、FPFH）、统计特征或深度学习特征拼接，形成混合特征向量，往往能获得最佳效果。

# 伪代码：特征融合示例 def extract_geometric_features(point_cloud): # 例如，计算法向量直方图、曲率统计等 geometric_feat = compute_fpfh(point_cloud) # 假设的函数 return geometric_feat # 对每个样本 tda_feat = extract_tda_features([pc])[0] geom_feat = extract_geometric_features(pc) combined_feat = np.concatenate([tda_feat, geom_feat])

4. 拓展到图数据与时间序列

TDA不仅适用于点云。图和时间序列可以通过巧妙的转换，变成点云或过滤复形。

4.1 图数据的拓扑特征提取

对于一个图G=(V, E)，我们可以用两种主要方式应用TDA：

方法A：将图节点嵌入为点云使用图嵌入算法（如Node2Vec, GraphSAGE）将每个节点映射为一个低维向量。然后，这个向量集合就可以当作点云来处理，计算其持续同调。这能捕捉节点在嵌入空间中的全局拓扑结构。

方法B：构建图本身的过滤复形这是更“原生”的方法。我们基于图的结构定义一个过滤参数（如节点/边的权重、度中心性等），构建一个嵌套的子图序列（过滤），然后计算其持续同调。

• 顶点过滤：根据每个节点的某个属性（如PageRank值）排序，逐步添加节点及其关联边。
• 边过滤：更常用。根据边的权重（如相似度）排序，从空图开始，逐步添加权重从大到小的边。随着边的增加，连通分量（H_0）会合并，环（H_1）会产生和消失。
• 团复形 (Clique Complex)：为了捕捉更高阶的交互，我们不只添加边，当一组节点两两相连（形成一个团）时，添加更高维的单形。这需要计算团复形Cl(G)。

import networkx as nx from gudhi import SimplexTree def graph_persistence_diagram(G, weight_attr='weight'): """ 计算基于边权过滤的图持续同调（H0）。 G: networkx图，边有权重属性。 """ # 1. 获取边及其权重，按权重升序排列（从最弱连接到最强连接） edges = list(G.edges(data=True)) edges_sorted = sorted(edges, key=lambda x: x[2].get(weight_attr, 1.0)) # 2. 初始化一个单纯复形（这里用GUDHI的SimplexTree） st = SimplexTree() # 先插入所有0-单形（顶点） for v in G.nodes(): st.insert([v], filtration=-1e10) # 在最小过滤时间插入顶点 # 3. 逐步添加边，并分配过滤值（权重） for u, v, data in edges_sorted: w = data.get(weight_attr, 1.0) st.insert([u, v], filtration=w) # 注意：这里为了简单只构建了图复形（只有0维和1维单形）。 # 要计算H1，通常需要构建团复形，即当三角形形成时插入2-单形。 # 这需要检测循环，计算量更大。 # 4. 计算持续同调 diag = st.persistence() # diag 是一个列表，元素如 (dim, (birth, death)) # 分离不同维度的特征 persistence_by_dim = {} for dim, (birth, death) in diag: if dim not in persistence_by_dim: persistence_by_dim[dim] = [] persistence_by_dim[dim].append((birth, death)) return persistence_by_dim # 示例：创建一个随机权重图 G = nx.erdos_renyi_graph(n=30, p=0.15) for u, v in G.edges(): G[u][v]['weight'] = np.random.rand() # 随机边权 diagrams = graph_persistence_diagram(G) print(f"H0特征: {len(diagrams.get(0, []))}") # 对于H1，需要更复杂的团复形构建，可使用GUDHI的FlagComplex或RipsComplex从距离矩阵构建。

实操心得：对于图数据，基于边权的H_0持续同调与层次聚类 (Hierarchical Clustering)的树状图有深刻联系。每个H_0特征的死亡时间对应着两个簇合并时的距离。因此，图的持续同调提供了一种多尺度的聚类视角。

4.2 时间序列的拓扑分析

时间序列{x_t}可以通过以下方式转化为拓扑对象：

方法A：时滞嵌入 (Takens Embedding)这是分析动力系统的经典方法。对于一维时间序列，通过时滞坐标重构出一个高维空间中的轨迹。

• 选择嵌入维数m和时滞τ。
• 构造点云：v_t = [x_t, x_{t+τ}, x_{t+2τ}, ..., x_{t+(m-1)τ}]。
• 对这个点云计算持续同调。其拓扑特征反映了原始动力系统的吸引子结构。

方法B：子水平集/超水平集过滤将时间序列视为一个定义在时间轴（1D空间）上的函数f(t)。

• 子水平集过滤：让一个阈值α从最小值扫到最大值。观察函数值低于α的区域（即{t | f(t) < α}）的拓扑变化。这主要捕捉“谷底”的连通性。
• 超水平集过滤：类似，观察函数值高于α的区域（{t | f(t) > α}），捕捉“波峰”的连通性。
• 这种过滤产生的持续同调，可以量化时间序列中“峰”和“谷”的持久性。

UCR数据集实战示例： UCR是时间序列分类的基准库。我们可以对每条时间序列应用时滞嵌入，然后提取拓扑特征。

from sklearn.preprocessing import StandardScaler def time_series_to_point_cloud(series, m=3, tau=1): """使用时滞嵌入将时间序列转化为点云""" n = len(series) if n < m * tau: raise ValueError("序列长度不足以进行时滞嵌入") point_cloud = [] for i in range(n - (m-1)*tau): point = series[i:i + m*tau:tau] point_cloud.append(point) return np.array(point_cloud) # 假设加载了UCR数据集中的一条数据 # series, label = load_ucr_sample(...) series = np.random.randn(100) # 示例随机序列 series = StandardScaler().fit_transform(series.reshape(-1, 1)).flatten() # 标准化 # 参数选择：m和tau的选择是关键，可以用互信息法、虚假最近邻法等，这里简单设置。 m, tau = 3, 5 pc = time_series_to_point_cloud(series, m=m, tau=tau) # 计算拓扑特征 rips = Rips(maxdim=1, thresh=3.0) diagrams = rips.fit_transform(pc) # 后续提取统计特征向量，与3.3节类似

注意事项：时间序列的拓扑分析对噪声非常敏感。预处理（如去噪、平滑）至关重要。此外，时滞嵌入参数(m, tau)的选择极大影响结果，需要通过交叉验证或基于数据本身的方法（如虚假最近邻法确定m，互信息法确定tau）进行优化。

5. 常见陷阱、技巧与高级话题

在实际项目中应用TDA，你会遇到一些典型的坑。这里分享一些我的经验。

5.1 计算效率与可扩展性

问题：Rips复形的计算复杂度随点数和维度呈组合爆炸式增长。对于超过几千个点的数据集，直接计算是不现实的。

解决方案：

• 下采样：如前所述，这是第一步。确保下采样方法能保持拓扑（如最远点采样）。
• 稀疏化/近似算法：
  • Witness复形：从大量点（地标点）中选出一个小子集（见证点）来近似原数据的拓扑。计算量远小于Rips。
  • Graph-induced 复形：先构建一个最近邻图（如k-NN图），然后只考虑这个图中的边和团。这相当于在距离矩阵上施加了一个稀疏约束。
• 使用高效库：ripser库内部使用了高度优化的算法（如稀疏矩阵的持久同调计算）。GUDHI也提供了多种复形和优化选项。
• 并行化：特征提取过程（对每个样本计算持续同调）是独立的，可以轻松并行。

from joblib import Parallel, delayed def process_one_sample(pc): rips = Rips(maxdim=1, thresh=2.0, n_threads=1) # 每个进程一个实例 diagrams = rips.fit_transform(pc) return persistence_statistics(diagrams, 0), persistence_statistics(diagrams, 1) # 并行处理多个点云 results = Parallel(n_jobs=4)(delayed(process_one_sample)(pc) for pc in list_of_point_clouds)

5.2 特征稳定性与参数选择

问题：拓扑特征对噪声、采样密度和算法参数（如Rips的maxdim,thresh）敏感。

解决方案：

• 稳定性理论：持续同调本身具有稳定性定理。小的数据扰动（在豪斯多夫距离或瓶颈距离下）只会导致持续图发生小的移动。这意味着长寿的特征是稳定的。
• 关注持久性长的特征：在分析时，重点关注那些远离对角线（持久性大）的点。它们代表了数据的稳健拓扑信号。
• 多尺度集成：不要只依赖一组参数。可以尝试不同的下采样率、不同的距离度量（如添加密度加权距离），然后集成多组拓扑特征。
• 使用持久性图像或景观：它们对特征的位置微小变化比原始的(b,d)点对更鲁棒。

5.3 与深度学习的结合

这是当前的研究热点。如何让拓扑特征与深度神经网络端到端地协同工作？

1. 作为输入特征：如前所述，将向量化后的拓扑特征与其它特征拼接，输入全连接网络。
2. 作为正则化项：设计一个基于拓扑的损失函数，鼓励网络学习到的隐层表示具有某种期望的拓扑性质。例如，希望某个层的激活值构成的点云具有较少的连通分量（更紧凑）。
3. 拓扑层：研究如何构建可微的拓扑层，使其能够嵌入到神经网络中。这是一个前沿方向，例如“持续同调变换”的尝试，但实现复杂且计算开销大。

一个简单的融合示例：在点云分类网络（如PointNet）中，除了每个点的特征和全局特征外，额外拼接一个由整个点云计算出的拓扑特征向量。

5.4 距离与比较：如何度量两个持续图的差异？

当需要比较两个数据集或两个模型的拓扑特征时，需要度量两个持续图PD1和PD2之间的距离。常用两种：

1. 瓶颈距离 (Bottleneck Distance)W_∞：找到两个点集之间的最佳匹配，使得任意匹配点对之间的距离最大值最小。它对异常值非常敏感。
2. p-Wasserstein距离W_p：是瓶颈距离的推广，考虑了所有匹配点对距离的p次方的和。p=2时就是平方和再开方，更平滑。

from persim import bottleneck, wasserstein # 计算两个持续图（例如都是H1）之间的距离 diagrams1 = ... # 第一个样本的持续图列表 diagrams2 = ... # 第二个样本的持续图列表 # 提取H1的图，并处理无穷大值（通常用np.max替换） dgm1 = diagrams1[1] dgm2 = diagrams2[1] # 处理无穷大死亡时间，用一个大的有限数替代 max_finite = max(np.max(dgm1[dgm1[:,1] != np.inf, 1]), np.max(dgm2[dgm2[:,1] != np.inf, 1])) dgm1[dgm1[:,1] == np.inf, 1] = max_finite * 1.1 dgm2[dgm2[:,1] == np.inf, 1] = max_finite * 1.1 bottleneck_dist = bottleneck(dgm1, dgm2) wasserstein_dist = wasserstein(dgm1, dgm2) print(f"瓶颈距离: {bottleneck_dist}") print(f"2-Wasserstein距离: {wasserstein_dist}")

这些距离可以用于构建拓扑层面的核函数，或者直接用于聚类、检索任务。

拓扑数据分析为理解复杂数据打开了一扇新窗口。它提供的“形状”特征，与传统的数值特征和深度学习学到的表示特征形成了有力的互补。入门的关键在于动手实践：选一个熟悉的数据集（比如MNIST图像转换成点云，或者一个简单的社交网络图），从计算一个持续图开始，观察其特征，思考其含义。开始时可能会被数学符号吓到，但记住，工具库已经为我们封装了大部分复杂性。我们的任务，是成为一个聪明的“调参师”和“特征工程师”，让这些深刻的数学思想为实际的机器学习问题服务。在实践中，你可能会发现，对于某些问题，拓扑特征带来了显著的提升；对于另一些问题，它可能只是提供了一个有趣但辅助的视角。但这正是探索的乐趣所在。