06 距离矩阵、稀疏与贪心置换

第 6 章 距离矩阵、稀疏过滤与贪心置换

当 n 很大或距离很稀疏时，预计算距离矩阵、稀疏存储、贪心置换子采样 是三类关键工程手段。本章讲清原理、API 用法与精度—速度权衡。

6.1 为什么需要这三招

VR 过滤的边数随 thresh 增长，最坏接近 O(n²) 边，更高维单纯形更爆炸。

6.2 预计算距离矩阵

from sklearn.metrics import pairwise_distances
from ripser import ripser

# 任意 sklearn 支持的 metric
D = pairwise_distances(X, metric='manhattan')
result = ripser(D, distance_matrix=True, maxdim=1, thresh=3.0)
​

注意：

• D 必须对称、非负、对角线为 0
• 不满足度量公理的距离有时仍可用，但理论保证减弱

6.3 稀疏距离矩阵

只保留 距离 ≤ thresh 或 k 近邻 的边，用 scipy.sparse 存储：

import numpy as np
from scipy import sparse
from sklearn.metrics import pairwise_distances
from ripser import ripser

X = np.random.randn(500, 10)
D_dense = pairwise_distances(X)
thresh = 2.0

# 只保留 <= thresh 的边（含对角线 0）
D_sparse = sparse.csr_matrix(D_dense * (D_dense <= thresh))

result = ripser(
    D_sparse,
    distance_matrix=True,
    maxdim=1,
    thresh=thresh,  # 与稀疏结构一致
)
print("edges in filtration:", result['num_edges'])
​

稀疏路径原理（直觉）

稠密 D  → C++ doRipsFiltrationDM（全矩阵）
稀疏 D  → doRipsFiltrationDMSparse（只遍历非零边）
​

仅存在于稀疏矩阵中的边才会进入过滤——等价于在 VR 中禁止「超长边」，与设置 thresh 一致时最自然。

⚠️ 禁止组合：n_perm 与稀疏距离矩阵 不能同时使用（会 ValueError）。大规模稀疏场景用 thresh 控制，或先子采样再稠密计算。

6.4 近似稀疏过滤

若只有点云、想要「局部 VR」：

1. 建 k-NN 图，边权为距离
2. 非边设为 0 或足够大（稀疏矩阵中不出现）
3. distance_matrix=True 传入

这得到 近似 VR——远距离点永不连边，可能漏掉应出现的跨簇短桥。需结合领域判断。

6.5 贪心置换（Greedy Permutation）

算法：最远点采样（Farthest Point Sampling）

1. 选起点 p0（默认第一个点）
2. 维护每个点到已选集合的最小距离 ds
3. 选 ds 最大的点作为下一个 p_i
4. 重复直到选满 n_perm 个点
​

性质：子集形成 ε-net，覆盖半径 r_cover 有理论界。

API

result = ripser(large_data, maxdim=1, n_perm=500)

idx = result['idx_perm']           # 500 个索引
subsample = large_data[idx]        # 子集点云
cover = result['r_cover']          # 覆盖半径
D_sub = result['dperm2all']        # (500, n) 子集到全集距离
​

精度保证（直观）

子采样持久图与全点云持久图在 瓶颈距离 意义下误差有界，界与 r_cover 相关。n_perm 越大越准、越慢。

全点云 (n=10000)          贪心 n_perm=500
        │                         │
        └──── 瓶颈距离 ≤ f(r_cover) ────┘
​

6.6 三招组合策略

                    开始
                      │
          n > 5000 ? ─┴─ 否 ──► 稠密 ripser(X)
                      │
                     是
                      │
          需要全图距离 ? ─ 否 ──► n_perm 子采样
                      │
                     是
                      │
          边很稀疏 ? ─ 是 ──► 稀疏 D + thresh
                      │
                     否
                      │
              thresh 截断 + 可选 n_perm
​

6.7 性能观测

import time
import numpy as np
from ripser import ripser

X = np.random.randn(2000, 5)

t0 = time.perf_counter()
r1 = ripser(X, maxdim=1, thresh=2.0)
t1 = time.perf_counter()

t2 = time.perf_counter()
r2 = ripser(X, maxdim=1, thresh=2.0, n_perm=400)
t3 = time.perf_counter()

print(f"full: {t1-t0:.2f}s edges={r1['num_edges']}")
print(f"perm: {t3-t2:.2f}s edges={r2['num_edges']} r_cover={r2['r_cover']:.3f}")
​

6.8 动手实验

1. 对 n=3000 随机 3D 点云，比较 n_perm=300/600/1200 与全量的 H1 瓶颈距离。
2. 构造 k-NN 稀疏矩阵（k=15），与稠密 thresh 结果对比 H1 显著点数量。
3. 验证 n_perm + sparse 触发 ValueError。
4. 记录 num_edges 随 thresh 变化，估计你的数据「可承受」边数上限。

本章小结

• 预计算 D + distance_matrix=True 支持任意距离。
• 稀疏矩阵走专用 C++ 路径，适合 k-NN / 阈值图；不可与 n_perm 同用。
• n_perm 用最远点采样子集，以有界误差换速度；看 r_cover 评估近似。
• 大规模实战常组合 thresh + 稀疏或 n_perm。

下一章：Lower-Star 过滤——图像与时序上的非 VR 过滤。