高效切片三角稀疏矩阵

Question

我有一个稀疏的三角矩阵（例如距离矩阵）。实际上，这将是一个 > 1M x 1M 的距离矩阵，具有高稀疏性。

from scipy.sparse import csr_matrix
X = csr_matrix([
      [1, 2, 3, 3, 1],
      [0, 1, 3, 3, 2],
      [0, 0, 1, 1, 3],
      [0, 0, 0, 1, 3],
      [0, 0, 0, 0, 1],
])

我想将此矩阵子集化为另一个三角距离矩阵。 索引的顺序可能不同和/或重复。

idx = np.matrix([1,2,4,2])
X2 = X[idx.T, idx]

这可能会导致生成的矩阵不是三角形的，其中一些值缺失 上三角形，一些值在下三角形中重复。

>>> X2.toarray()
array([[1, 3, 2, 3],
       [0, 1, 3, 1],
       [0, 0, 1, 0],
       [0, 1, 3, 1]])

如何才能尽可能高效地得到正确的上三角矩阵？ 目前，我在子集之前镜像矩阵，然后将其子集到三角形，但这感觉不是特别有效，因为它至少需要复制所有条目。

# use transpose method, see https://stackoverflow.com/a/58806735/2340703
X = X + X.T - scipy.sparse.diags(X.diagonal())
X2 = X[idx.T, idx]
X2 = scipy.sparse.triu(X2, k=0, format="csr")

>>> X2.toarray()
array([[1., 3., 2., 3.],
       [0., 1., 3., 1.],
       [0., 0., 1., 3.],
       [0., 0., 0., 1.]])

Answer 1

总结所有优秀的贡献，这个问题的简短回答是：

不要使用三角矩阵。与使用方阵相比，在速度或内存方面没有任何优势。

原因在@hpaulj's answer中有解释：

• 对稀疏矩阵进行切片使用矩阵乘法，这是非常高效的。将矩阵重新排列成三角形会很慢。
• 使用triu 是一项昂贵的操作，因为它实现了密集掩码。

将@jakevdp's solution 与仅使用方阵进行比较时，这一点变得很明显。使用方形会更快，并且使用更少的内存。

测试使用一个稀疏的三角形 800k x 800k 距离矩阵，具有高稀疏度 (%nnz here.

# Running benchmark: Converting to square matrix
./benchmark.py squarify   6.29s  user 1.59s system 80% cpu 9.738 total
max memory:                4409 MB

# Running benchmark: @jakevdp's solution
./benchmark.py sparse_triangular   67.03s  user 3.01s system 99% cpu 1:10.15 total
max memory:                5209 MB

如果有人迫切希望在使用方阵之外对其进行优化，@CJR's comment 是一个很好的起点：

我会考虑将其实现为与 pdist 相同的样式的压缩距离矩阵，但作为 1xN CSR 矩阵，然后在需要获取特定值时使用坐标数学对其进行重新索引。

Answer 2

这不是一个改进的工作答案，而是探索稀疏索引和triu 的作用。它可能会为您提供进行更直接计算的想法。您从 tri 开始并期望 tri 的事实意味着这不是一项简单的任务，即使是密集数组（索引速度要快得多）也不行。

sparse.csr 索引使用矩阵乘法。我将用密集数组来说明这一点：

In [304]: X = np.array([
     ...:       [1, 2, 3, 3, 1],
     ...:       [0, 1, 3, 3, 2],
     ...:       [0, 0, 1, 1, 3],
     ...:       [0, 0, 0, 1, 3],
     ...:       [0, 0, 0, 0, 1],
     ...: ])
In [305]: idx = np.array([1,2,4,2])
In [306]: X[idx[:,None],idx]
Out[306]: 
array([[1, 3, 2, 3],
       [0, 1, 3, 1],
       [0, 0, 1, 0],
       [0, 1, 3, 1]])
In [307]: m = np.array([[0,1,0,0,0],[0,0,1,0,0],[0,0,0,0,1],[0,0,1,0,0]])
In [308]: m@X@m.T
Out[308]: 
array([[1, 3, 2, 3],
       [0, 1, 3, 1],
       [0, 0, 1, 0],
       [0, 1, 3, 1]])

以及全距离数组：

In [309]: X2 = X+X.T-np.diag(np.diag(X))
In [311]: X2[idx[:,None],idx]
Out[311]: 
array([[1, 3, 2, 3],
       [3, 1, 3, 1],
       [2, 3, 1, 3],
       [3, 1, 3, 1]])
In [312]: m@X2@m.T
Out[312]: 
array([[1, 3, 2, 3],
       [3, 1, 3, 1],
       [2, 3, 1, 3],
       [3, 1, 3, 1]])

我不知道是否可以直接从X（或X2）构造一个提供所需结果的m，上三或不直接

In [316]: sparse.triu(Out[312])
Out[316]: 
<4x4 sparse matrix of type '<class 'numpy.int64'>'
    with 10 stored elements in COOrdinate format>
In [317]: _.A
Out[317]: 
array([[1, 3, 2, 3],
       [0, 1, 3, 1],
       [0, 0, 1, 3],
       [0, 0, 0, 1]])

sparse.triu 会：

In [331]: A = sparse.coo_matrix(_312)
     ...: mask = A.row <= A.col 
In [332]: A
Out[332]: 
<4x4 sparse matrix of type '<class 'numpy.int64'>'
    with 16 stored elements in COOrdinate format>
In [333]: mask
Out[333]: 
array([ True,  True,  True,  True, False,  True,  True,  True, False,
       False,  True,  True, False, False, False,  True])

这个mask数组有16个词，A.nnz。

然后它创建一个新的coo 矩阵，其中包含从A 属性中选择的数据/行/列数组：

In [334]: d=A.data[mask]
In [335]: r=A.row[mask]
In [336]: c=A.col[mask]
In [337]: d
Out[337]: array([1, 3, 2, 3, 1, 3, 1, 1, 3, 1])
In [338]: sparse.coo_matrix((d, (r,c)))
Out[338]: 
<4x4 sparse matrix of type '<class 'numpy.int64'>'
    with 10 stored elements in COOrdinate format>
In [339]: _.A
Out[339]: 
array([[1, 3, 2, 3],
       [0, 1, 3, 1],
       [0, 0, 1, 3],
       [0, 0, 0, 1]])

np.triu 使用 mask 之类的：

In [349]: np.tri(4,4,-1)
Out[349]: 
array([[0., 0., 0., 0.],
       [1., 0., 0., 0.],
       [1., 1., 0., 0.],
       [1., 1., 1., 0.]])

Answer 3

这是一种不涉及镜像数据的方法，而是对稀疏索引进行操作以达到所需的结果：

import scipy.sparse as sp

X2 = X[idx.T, idx]

# Extract indices and data (this is essentially COO format)
i, j, data = sp.find(X2)

# Generate indices with elements moved to upper triangle
ij = np.vstack([
  np.where(i > j, j, i),
  np.where(i > j, i, j)
])

# Remove duplicate elements
ij, ind = np.unique(ij, axis=1, return_index=True)

# Re-build the matrix
X2 = sp.coo_matrix((data[ind], ij)).tocsr()

Answer 4

好吧，我无法将其发送至 triu，但这应该会更快：

idx = np.array([1,2,4,2])
i = np.stack(np.meshgrid(idx, idx))
X2 = X[i.min(0), i.max(0)]
 
array([[1, 3, 2, 3],
       [3, 1, 3, 1],
       [2, 3, 1, 3],
       [3, 1, 3, 1]])

所以整个过程是：

idx = np.array([1,2,4,2])
i = np.stack(np.meshgrid(idx, idx))
X2 = scipy.sparse.triu(X[i.min(0), i.max(0)], k=0, format="csr")

但我无法摆脱必须有更优化的方法的感觉。