如何解决高效切片三角稀疏矩阵
我有一个稀疏的三角矩阵(例如距离矩阵)。实际上,这将是一个 > 1M x 1M 的具有高稀疏性的距离矩阵。
from scipy.sparse import csr_matrix
X = csr_matrix([
[1,2,3,1],[0,1,2],3],])
我想将这个矩阵子集到另一个三角距离矩阵。 索引的顺序可能不同和/或重复。
idx = np.matrix([1,4,2])
X2 = X[idx.T,idx]
这可能会导致生成的矩阵不是三角形的,一些值从 上三角,一些值在下三角中重复。
>>> X2.toarray()
array([[1,0],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.,1.],0.,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,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()
,
好吧,我本机无法将其转至 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,3],[3,1,1],[2,1]])
所以整个过程是:
idx = np.array([1,idx))
X2 = scipy.sparse.triu(X[i.min(0),i.max(0)],k=0,format="csr")
但我无法摆脱这种感觉,必须有更优化的方法。
,这不是改进的工作答案,而是探索稀疏索引和 triu
的作用。它可能会为您提供进行更直接计算的想法。事实上,你从一个 tri 开始,并期待一个 tri,这意味着这不是一项微不足道的任务,即使是密集数组(索引速度要快得多)也是如此。
sparse.csr
索引使用矩阵乘法。我将用密集数组来说明这一点:
In [304]: X = np.array([
...: [1,...: [0,2],...: ])
In [305]: idx = np.array([1,2])
In [306]: X[idx[:,None],idx]
Out[306]:
array([[1,[0,0],1]])
In [307]: m = np.array([[0,0]])
In [308]: m@X@m.T
Out[308]:
array([[1,1]])
和全距离数组:
In [309]: X2 = X+X.T-np.diag(np.diag(X))
In [311]: X2[idx[:,idx]
Out[311]:
array([[1,1]])
In [312]: m@X2@m.T
Out[312]:
array([[1,1]])
我不知道是否可以直接从 m
(或 X
)构建一个提供所需结果的 X2
或不提供所需的结果
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,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,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,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,1]])
np.triu
使用 mask
像:
In [349]: np.tri(4,-1)
Out[349]:
array([[0.,0.,0.],[1.,1.,0.]])
,
总结所有出色的贡献,对这个问题的简短回答是:
不要使用三角矩阵。与使用方阵相比,在速度或内存方面没有任何好处。
原因在@hpaulj's answer中解释:
- 对稀疏矩阵进行切片使用非常高效的矩阵乘法。将矩阵重新排列成三角形会很慢。
- 使用
triu
是一项代价高昂的操作,因为它实现了一个密集的掩码。
将 @jakevdp's solution 与仅使用方阵进行比较时,这一点变得明显。使用正方形的形式更快,占用的内存更少。
该测试使用具有高稀疏度(%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 矩阵,然后在需要获取特定值时使用坐标数学对其重新索引。
>版权声明:本文内容由互联网用户自发贡献,该文观点与技术仅代表作者本人。本站仅提供信息存储空间服务,不拥有所有权,不承担相关法律责任。如发现本站有涉嫌侵权/违法违规的内容, 请发送邮件至 dio@foxmail.com 举报,一经查实,本站将立刻删除。