文章目录
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)
首先定义贪心序列,其用于执行最远点的采用,以确定插入半径 λ i \lambda_i λ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 稀疏距离矩阵
现在定义一个函数,其使用一个表示点云的距离矩阵、一个插入半径的集合,以及近似因子 ϵ \epsilon ϵ作为输入,并返回重权衡边界的稀疏距离矩阵。每一个持续图都保证是真实持续图的 ( 1 + ϵ ) (1 + \epsilon) (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
- 距离矩阵: