python - 计算 scipy 稀疏矩阵的稀疏传递闭包

Question

我想计算 transitive closure Python 中的稀疏矩阵。目前我正在使用 scipy 稀疏矩阵。

矩阵幂(**12 在我的例子中)适用于非常稀疏的矩阵，无论它们有多大，但对于定向的不那么稀疏的情况我想使用更智能的算法。

我找到了 Floyd-Warshall algorithm (德语页面有更好的伪代码)在scipy.sparse.csgraph ，它做的比它应该做的多一点:没有仅针对 Warshall 算法 的函数 - 这是一回事。

主要问题是我可以将稀疏矩阵传递给函数，但这完全没有意义，因为函数将始终返回稠密矩阵，因为传递闭包中应该为 0 的现在是 的路径inf 长度，有人认为这需要显式存储。

所以我的问题是:是否有任何 python 模块可以计算稀疏矩阵的传递闭包并使其保持稀疏？

我不是 100% 确定他使用相同的矩阵，但是 Gerald Penn 在 his comparison paper 中展示了令人印象深刻的加速, 这表明有可能解决该问题。

---

编辑: 由于存在一些混淆，我将指出理论背景:

我正在寻找传递闭包(非自反或对称)。

我将确保我在 bool 矩阵中编码的关系具有所需的属性，即对称性或自反性。

我有两种关系:

1. 反身
2. 自反 和对称

我想对这两个关系应用传递闭包。这与矩阵幂一起工作得很好(只是在某些情况下它太昂贵了):

>>> reflexive
matrix([[ True,  True, False,  True],
        [False,  True,  True, False],
        [False, False,  True, False],
        [False, False, False,  True]])
>>> reflexive**4
matrix([[ True,  True,  True,  True],
        [False,  True,  True, False],
        [False, False,  True, False],
        [False, False, False,  True]])
>>> reflexive_symmetric
matrix([[ True,  True, False,  True],
        [ True,  True,  True, False],
        [False,  True,  True, False],
        [ True, False, False,  True]])
>>> reflexive_symmetric**4
matrix([[ True,  True,  True,  True],
        [ True,  True,  True,  True],
        [ True,  True,  True,  True],
        [ True,  True,  True,  True]])

因此，在第一种情况下，我们获取节点的所有后代(包括其自身)，在第二种情况下，我们获取所有组件，即同一组件中的所有节点。

Answer 1

这是提出来的on SciPy issue tracker .问题不在于输出格式； Floyd-Warshall 的实现是从充满无穷大的矩阵开始，然后在找到路径时插入有限值。稀疏性立即消失。

networkx 库提供了一个替代方案，其 all_pairs_shortest_path_length .它的输出是一个迭代器，它返回

形式的元组

(source, dictionary of reachable targets) 

需要一些工作才能转换为 SciPy 稀疏矩阵(这里很自然地使用 csr 格式)。一个完整的例子:

import numpy as np
import networkx as nx
import scipy.stats as stats
import scipy.sparse as sparse

A = sparse.random(6, 6, density=0.2, format='csr', data_rvs=stats.randint(1, 2).rvs).astype(np.uint8)
G = nx.DiGraph(A)       # directed because A need not be symmetric
paths = nx.all_pairs_shortest_path_length(G)
indices = []
indptr = [0]
for row in paths:
  reachable = [v for v in row[1] if row[1][v] > 0]
  indices.extend(reachable)
  indptr.append(len(indices))
data = np.ones((len(indices),), dtype=np.uint8)
A_trans = A + sparse.csr_matrix((data, indices, indptr), shape=A.shape)
print(A, "\n\n", A_trans)

加回A的原因如下。 Networkx 输出包括长度为 0 的路径，它会立即填充对角线。我们不希望这种情况发生(你想要传递闭包，而不是自反传递闭包)。因此 reachable = [v for v in row[1] if row[1][v] > 0] 行。但是我们根本没有得到任何对角线条目，即使 A 有它们(0 长度空路径击败由自循环形成的 1 长度路径)。所以我将 A 添加回结果。它现在有条目 1 或 2，但只有它们非零的事实才有意义。

运行上面的示例(为了输出的可读性，我选择 6 x 6 大小)。原始矩阵:

  (0, 3)    1
  (3, 2)    1
  (4, 3)    1
  (5, 1)    1
  (5, 3)    1
  (5, 4)    1
  (5, 5)    1 

传递闭包:

  (0, 2)    1
  (0, 3)    2
  (3, 2)    2
  (4, 2)    1
  (4, 3)    2
  (5, 1)    2
  (5, 2)    1
  (5, 3)    2
  (5, 4)    2
  (5, 5)    1

您可以看到它工作正常:添加的条目是 (0, 2)、(4, 2) 和 (5, 2)，它们都是通过路径 (3, 2) 获取的。

---

顺便说一下，networkx 也有 floyd_warshall方法，但它的文档说

This algorithm is most appropriate for dense graphs. The running time is O(n^3), and running space is O(n^2) where n is the number of nodes in G.

输出再次密集。我的印象是这个算法本质上被认为是密集的。似乎 all_pairs_shortest_path_length 是一种 Dijkstra's algorithm .

传递和自反

如果您想要传递闭包(包含给定闭包的最小传递关系)而不是传递闭包闭包(包含给定闭包的最小传递闭包和自反关系)一)，代码简化了，因为我们不再担心 0 长度路径。

for row in paths:
  indices.extend(row[1])
  indptr.append(len(indices))
data = np.ones((len(indices),), dtype=np.uint8)
A_trans = sparse.csr_matrix((data, indices, indptr), shape=A.shape)

传递、自反和对称

这意味着找到包含给定关系的最小等价关系。等价地，将顶点划分为连通分量。为此你不需要去networkx，有connected_components SciPy 的方法。在那里设置 directed=False。示例:

import numpy as np
import scipy.stats as stats
import scipy.sparse as sparse
import itertools

A = sparse.random(20, 20, density=0.02, format='csr', data_rvs=stats.randint(1, 2).rvs).astype(np.uint8)
components = sparse.csgraph.connected_components(A, directed=False)
nonzeros = []
for k in range(components[0]):
  idx = np.where(components[1] == k)[0]
  nonzeros.extend(itertools.product(idx, idx))
  row = tuple(r for r, c in nonzeros)
  col = tuple(c for r, c in nonzeros)
  data = np.ones_like(row)
B = sparse.coo_matrix((data, (row, col)), shape=A.shape)

这是一个随机示例的输出 print(B.toarray()) 的样子，20 x 20:

[[1 0 1 0 0 0 1 0 0 0 0 0 1 0 0 0 1 0 0 0]
 [0 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0]
 [1 0 1 0 0 0 1 0 0 0 0 0 1 0 0 0 1 0 0 0]
 [0 0 0 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0]
 [0 0 0 0 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0]
 [0 0 0 0 0 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0]
 [1 0 1 0 0 0 1 0 0 0 0 0 1 0 0 0 1 0 0 0]
 [0 0 0 0 0 0 0 1 0 0 0 0 0 0 1 0 0 0 0 0]
 [0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0 0 0]
 [0 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0 0]
 [0 0 0 0 0 0 0 0 0 0 1 1 0 1 0 0 0 0 0 0]
 [0 0 0 0 0 0 0 0 0 0 1 1 0 1 0 0 0 0 0 0]
 [1 0 1 0 0 0 1 0 0 0 0 0 1 0 0 0 1 0 0 0]
 [0 0 0 0 0 0 0 0 0 0 1 1 0 1 0 0 0 0 0 0]
 [0 0 0 0 0 0 0 1 0 0 0 0 0 0 1 0 0 0 0 0]
 [0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 0 0 0 0]
 [1 0 1 0 0 0 1 0 0 0 0 0 1 0 0 0 1 0 0 0]
 [0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 1 0]
 [0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 1 0]
 [0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1]]