在Python中访问矩阵特定行中非零元素

Question

我有一个稀疏矩阵在Python中非零元素的位置索引，形式如下：

(array([0, 1, 2], dtype=int32), array([2, 0, 0], dtype=int32), array([2, 1, 3]))

或者以矩阵形式表示

[[0 2]
 [1 0]
 [2 0]]

我想使用这个（或者如果有其他方法）来仅对其他矩阵中相应的非零元素进行逐行操作，例如：

for r in range(rows):
    A[r,:] = np.dot(B[r,:],C.T)

基本上，我需要一种指定行的方法，并仅选择与矩阵B中非零元素对应的该行元素。

我无法理解的部分是由于每行/列的条目数量可能会有所不同。

- Deirdre Meehan

---

仅使用非零元素会如何改变最终答案？ - hpaulj

你设置了 scipy 和 sparse-matrix 标签。这是否意味着你打算使用 scipy.sparse 矩阵？它们可以在编译代码中执行矩阵乘法。 - hpaulj

np.nonzero会给你非零元素的索引。但即使你知道这些元素的索引，在A和B中对它们进行索引比乘以零元素更耗费效率。 - hpaulj

我希望在Python中实现类似的东西，那么保留零元素可能比尝试复制Matlab版本更快吗？ - Deirdre Meehan

这种索引方式在MATLAB中是否真的能提高速度？如果密度小于0.01，可能会有所帮助。但是多年前我在有限元计算中使用了MATLAB稀疏矩阵。但稀疏矩阵更多地是用来节省内存而不是加速工具。 - hpaulj

Answer 1

我发现在Python中可以使用布尔数组索引，因此以下代码可以实现我想要的功能：

for r in range(rows):
    A[r,B[r,:]!=0] = np.dot(B[r , B[r,:]!=0], C[: , B[r,:]!=0].T)

看起来有点复杂，但它可以得到正确的元素进行计算。唯一的问题是当B的尺寸大于它所索引的尺寸时，会抛出索引越界错误。

- Deirdre Meehan

Answer 2

我有点困惑第一个元组代表什么。它是稀疏数组的索引和值吗？例如：

In [4]: arrays=(np.array([0, 1, 2], dtype=int), np.array([2, 0, 0], dtype=int), np.array([2, 1, 3], dtype=float))
...
In [6]: from scipy import sparse
In [7]: M=sparse.csr_matrix((arrays[2],(arrays[0],arrays[1])))
In [8]: M
Out[8]: 
<3x3 sparse matrix of type '<class 'numpy.float64'>'
    with 3 stored elements in Compressed Sparse Row format>
In [9]: M.A
Out[9]: 
array([[ 0.,  0.,  2.],
       [ 1.,  0.,  0.],
       [ 3.,  0.,  0.]])

In [10]: print(M)
  (0, 2)    2.0
  (1, 0)    1.0
  (2, 0)    3.0

矩阵乘法是针对这样的数组定义的：

In [12]: M*M.T
Out[12]: 
<3x3 sparse matrix of type '<class 'numpy.float64'>'
    with 5 stored elements in Compressed Sparse Row format>
In [13]: (M*M.T).A
Out[13]: 
array([[ 4.,  0.,  0.],
       [ 0.,  1.,  3.],
       [ 0.,  3.,  9.]])
In [14]: M.A.dot(M.A.T)  # dense equivalent
Out[14]: 
array([[ 4.,  0.,  0.],
       [ 0.,  1.,  3.],
       [ 0.,  3.,  9.]])

我可以使用这个数组实现逐行相乘：

In [21]: A=M.A      # dense array
In [22]: for r in range(3):
    print(np.dot(A[r,:],A[r,:]))
4.0
1.0
9.0
# actually this is just the diagonal

In [23]: for r in range(3):   # or with the nonzero elements
    I=np.nonzero(A[r,:])
    dot = np.dot(A[r,I[0]],A[r,I[0]])
    print(dot)
4.0
1.0
9.0

就算不值得，nonzero 返回我从您的帖子开头复制的数组：

In [24]: np.nonzero(A)
Out[24]: (array([0, 1, 2], dtype=int32), array([2, 0, 0], dtype=int32))
In [25]: A[np.nonzero(A)]
Out[25]: array([ 2.,  1.,  3.])

稀疏矩阵还具有nonzero方法：

In [26]: M.nonzero()
Out[26]: (array([0, 1, 2], dtype=int32), array([2, 0, 0], dtype=int32))

我感觉自己正在挣扎中，试图理解问题和示例。

- hpaulj