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

Question

我想在Python中计算稀疏矩阵的传递闭包.目前我正在使用scipy稀疏矩阵.

矩阵功率(**12在我的情况下)在非常稀疏的矩阵上运行良好,无论它们有多大,但对于有针对性的不那么稀疏的情况,我想使用更智能的算法.

我发现了Floyd-Warshall算法(德语页面有更好的伪代码)scipy.sparse.csgraph,它比它应该做的多一点:只有Warshall的算法没有功能- 这是一回事.

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

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

我并不是100%确定他使用相同的矩阵,但Gerald Penn在他的比较文件中显示了令人印象深刻的加速,这表明有可能解决问题.

---

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

我正在寻找传递闭包(不反身或对称).

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

我有两种关系:

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

这是在SciPy问题跟踪器上提出的.问题不在于输出格式; 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乘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的算法.

传递性和反身性

如果不是传递闭包(它是包含给定的传递关系的最小传递关系),那么你需要传递和自反闭包(包含给定关系的最小传递和自反关系),代码简化了,因为我们不再担心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_componentsSciPy的方法.在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乘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]]