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

Posted 2023-02-19

【中文标题】计算 scipy 稀疏矩阵的稀疏传递闭包

【英文标题】：Compute sparse transitive closure of scipy sparse matrix

【发布时间】：2018-06-04 00:54:11

【问题描述】：

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

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

我在scipy.sparse.csgraph 中找到了Floyd-Warshall algorithm（德文页面有更好的伪代码），这比它应该做的要多一点：没有仅适用于Warshall 算法的函数 - 这是一回事。

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

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

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

---

编辑：由于存在许多混淆，我将指出理论背景：

我正在寻找传递闭包（不是自反或对称的）。

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

我有两种关系：

自反 自反和对称

我想对这两个关系应用传递闭包。这对矩阵功率非常有效（只是在某些情况下它太昂贵了）：

>>> 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]])

所以在第一种情况下，我们得到一个节点的所有后代（包括它自己），而在第二种情况下，我们得到所有组件，即在同一个组件中的所有节点。

【问题讨论】：

关于 PS，我想看一些例子（在另一个问题中？）。稀疏矩阵与密集数组的乘法返回密集数组。为提高效率，sparse.csr 不会在每次更改值时更改稀疏索引。有时你必须运行eliminate_zeros 来清理它。

过去的例子：***.com/questions/37674435/…, ***.com/questions/43146968/…

如果乘法返回密集矩阵，请先尝试将other 数组转换为稀疏矩阵。 sparse*sparse 产生 sparse。

floyd_warshall 在 sparse.csgraph.shortest_path.so 中，一个已编译的模块。

你说得对，我刚刚在“对 scipy 的失望”下收集了这个......我会为此做一个新问题。

【参考方案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 方法，但它的文档说

此算法最适合密集图。运行时间为O(n^3)，运行空间为O(n^2)，其中n为G中的节点数。

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

传递性和自反性

如果你想要 transitive and reflexive 闭包（包含给定的最小传递和反身关系）而不是 传递闭包（这是包含给定的最小传递关系） one) ，代码简化了，因为我们不再担心 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]]

【讨论】：

我有点担心创建 python 字典作为输出的开销，但也许真的没有基于 C 的实现返回一个稀疏的 scipy 矩阵......

也许我误解了传递闭包的含义。我实际上想要你所说的“自反传递闭包”。

node_connected_component(G,n) 也是一种选择，但我希望有一个解决方案可以处理通常的依赖项 numy 或 scipy

连通分量也执行对称闭合，因为“a 和 b 在同一个分量中”的关系是对称的。请向您自己和其他人说明，无论您是寻求传递闭包、传递自反闭包、传递自反对称闭包（因此，由给定闭包产生的等价关系）还是其他什么.

添加了使用 connected_components 的等价关系的完整示例。