Ripser.py学习 (4)：近似空间滤形 (Approximate Sparse Filtrations)

1 概述

近似空间滤形用于降低持续同伦算法的运行时间，关于该算法的更多信息可以参考：
https://www.youtube.com/watch?v=3WxuSwQhAgU

2 一些必须的库

import numpy as np
import matplotlib.pyplot as plt
import time
import tadasets
from ripser import ripser
from persim import plot_diagrams
from sklearn.metrics.pairwise import pairwise_distances
from scipy import sparse

3 测试数据

绘制测试数据，并输出持续同伦算法的运行时间：

np.random.seed(1)
X = tadasets.infty_sign(n=2000, noise=0.1)
plt.scatter(X[:, 0], X[:, 1])
plt.show()

tic = time.time()
resultfull = ripser(X)
toc = time.time()
timefull = toc-tic
print("Elapsed Time: %.3g seconds, %i Edges added" % (timefull, resultfull['num_edges']))

输出如下：

以及：

Elapsed Time: 4.73 seconds, 1285560 Edges added

4 近似空间滤形

4.1 贪心序列 (Greedy Permutation)

首先定义贪心序列，其用于执行最远点的采用，以确定插入半径${\lambda }_i$：

def getGreedyPerm(D):
    """
    一个指定最远点采样的简单方法O(N^2)

    Parameters
    ----------
    D : ndarray (N, N)
        N x N的距离句子

    Return
    ------
    lamdas: list
        每一个点的插入半径
    """

    N = D.shape[0]
    # 默认将序列中的第一个点作为点云中的第一个点，当然也可以采用随机操作
    perm = np.zeros(N, dtype=np.int64)
    lambdas = np.zeros(N)
    ds = D[0, :]
    for i in range(1, N):
        idx = np.argmax(ds)
        perm[i] = idx
        lambdas[i] = ds[idx]
        ds = np.minimum(ds, D[idx, :])

    return lambdas[perm]

4.2 稀疏距离矩阵

现在定义一个函数，其使用一个表示点云的距离矩阵、一个插入半径的集合，以及近似因子$ϵ$作为输入，并返回重权衡边界的稀疏距离矩阵。每一个持续图都保证是真实持续图的$(1+ϵ)$乘法近似：

def getApproxSparseDM(lambdas, eps, D):
    """
    目的： 返回通过权重排序后的稀疏边界列表

    Parameters
    ----------
    lambdas: list
        插入半径
    eps: float
        近似因子
    D: ndarray
        NxN距离矩阵

    Return
    ------
    DSparse: scipy.sparse
        NxN稀疏距离矩阵
    """
    N = D.shape[0]
    E0 = (1+eps)/eps
    E1 = (1+eps)**2/eps

    # 创建稀疏列表候选
    nBounds = ((eps**2+3*eps+2)/eps) * lambdas

    # 设置搜索领域外的所有距离为正无穷
    D[D > nBounds[:, None]] = np.inf
    [I, J] = np.meshgrid(np.arange(N), np.arange(N))
    idx = I < J
    I = I[(D < np.inf)*(idx == 1)]
    J = J[(D < np.inf)*(idx == 1)]
    D = D[(D < np.inf)*(idx == 1)]

    # 修剪稀疏列表并更新边界长度 (14页Algorithm 3)
    minlam = np.minimum(lambdas[I], lambdas[J])
    maxlam = np.maximum(lambdas[I], lambdas[J])

    # 排除顶点之间的边缘，这些顶点的球在接触之前停止增长，或者其中一个顶点将被删除
    # M存储先发生的那一个
    M = np.minimum((E0 + E1)*minlam, E0*(minlam + maxlam))

    t = np.arange(len(I))
    t = t[D <= M]
    (I, J, D) = (I[t], J[t], D[t])
    minlam = minlam[t]

    # 计算出实际添加的边界的度量
    t = np.ones(len(I))

    # 如果圆锥体没有变成圆柱体，度量不变
    t[D <= 2*minlam*E0] = 0

    # 否则，如果它们在上面的M条件之前满足，则度量被扭曲
    D[t == 1] = 2.0*(D[t == 1] - minlam[t == 1]*E0) # Multiply by 2 convention
    return sparse.coo_matrix((D, (I, J)), shape=(N, N)).tocsr()

4.3 获取近似解

tic = time.time()
D = pairwise_distances(X, metric='euclidean')
# 贪婪采样
lambdas = getGreedyPerm(D)
# 近似稀疏矩阵
DSparse = getApproxSparseDM(lambdas, 0.1, D)

# 稀疏解
resultsparse = ripser(DSparse, distance_matrix=True)
toc = time.time()
timesparse = toc - tic
percent_added = 100.0 * float(resultsparse['num_edges']) / resultfull['num_edges']
print("Elapsed Time: %.3g seconds, %i Edges added" % (timesparse, resultsparse['num_edges']))

输出如下：

Elapsed Time: 0.963 seconds, 354585 Edges added

可以发现计算时间大幅降低，且连接的边的数量减少。

4.4 绘制原始解和稀疏解

plt.figure(figsize=(10, 5))
plt.subplot(121)
D = pairwise_distances(X, metric='euclidean')
plt.imshow(D)
plt.title("Original Distance Matrix: %i Edges"%resultfull['num_edges'])
plt.subplot(122)
DSparse = DSparse.toarray()
DSparse = DSparse + DSparse.T
plt.imshow(DSparse)
plt.title("Sparse Distance Matrix: %i Edges"%resultsparse['num_edges'])

输出如下：

4.5 绘制持续图

plt.figure(figsize=(10, 5))
plt.subplot(121)
plot_diagrams(resultfull['dgms'], show=False)
plt.title("Full Filtration: Elapsed Time %g Seconds"%timefull)
plt.subplot(122)
plt.title("Sparse Filtration (%.3g%% Added)\nElapsed Time %g Seconds"%(percent_added, timesparse))
plot_diagrams(resultsparse['dgms'], show=True)

输出如下：

4.6 调整稀疏矩阵的稀疏程度

增大近似矩阵输入的eps即可：

DSparse = getApproxSparseDM(lambdas, 0.4, D)

部分输出如下：

• 时间：Elapsed Time: 0.188 seconds, 86524 Edges added
• 距离矩阵：