高效地切割三角形稀疏矩阵

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://dev59.com/tWQn5IYBdhLWcg3w9K73#58806735
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.]])

- Gregor Sturm

---

请澄清一下 - 您是将样本采样回到与原始距离矩阵相同的大小，还是将其子集化为较小的大小？ - CJR

由于重复元素的存在，翻译后的文本大小可能会更大，甚至没有显著缩小。 - Gregor Sturm

是的，使用 triu 来节省内存是个好主意。但是我开始有这样的印象，即这并不值得。创建 triu 似乎会消耗大量内存。它是否实现了布尔掩码？ - Gregor Sturm

从实际例子来看：X.nnz / X.shape[0]**2 = 8.28e-05，翻译为中文是：远小于1%。 - Gregor Sturm

我会考虑以压缩距离矩阵的方式来实现，就像pdist一样，但是作为一个1xN CSR矩阵，并且在需要获取特定值时使用坐标数学重新索引它。不过这有点像XY解决方案；我只是认为没有好的方法来做你要求做的具体事情。 - CJR

Answer 1

这里有一种方法不涉及数据镜像，而是使用稀疏索引进行操作，以得到所需结果：

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

- jakevdp

Answer 2

好的，我无法直接将其转换为 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")

但我总感觉一定有更优化的方式。

- Daniel F

---

有趣的方法！但在当前形式下，它不会扩展，因为meshgrid是一个密集矩阵，其维度与“X”相同。 - Gregor Sturm

Answer 3

这不是一个改进的工作答案，而是对于稀疏索引和 "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]])

的作用是：

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。
然后它从A的属性中选择数据/行/列数组，并生成一个新的coo矩阵：

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

- hpaulj

Answer 4

总结所有出色的贡献，对于这个问题的简短回答是：

不要使用三角形矩阵。与使用方形矩阵相比，无论是速度还是内存方面都没有任何优势。

原因在于@hpaulj's answer中已经解释了：

稀疏矩阵上的切片使用矩阵乘法非常高效。重新排列矩阵成三角形形状将会很慢。

使用triu是一个昂贵的操作，因为它会实现一个密集的掩码。

当比较@jakevdp's solution和仅使用方形矩阵时，这一点变得明显。使用方形形式更快，使用的内存更少。 测试使用稀疏三角形800k x 800k距离矩阵，具有高稀疏性(%nnz << 1%)。数据和代码可在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矩阵，并在需要获取特定值时使用坐标计算进行重新索引。

- Gregor Sturm