【问题标题】:Python-Scipy sparse Matrices - what is A[i, j] doing?Python-Scipy 稀疏矩阵 - A[i, j] 在做什么?
【发布时间】:2025-12-31 22:20:03
【问题描述】:

根据我之前的问题(Python - Multiply sparse matrix row with non-sparse vector by index),稀疏矩阵的直接索引是不可能的(至少如果您不想使用定义sparse.csr 矩阵的三个数组,@987654324 @、indicesindptr)。 但我刚刚发现,给定一个 csr 稀疏矩阵 A,这个操作可以正常工作并产生正确的结果:A[i, j]。 我还注意到:它非常慢,甚至比处理密集矩阵还要慢。

我找不到有关此索引方法的任何信息,所以我想知道:A[i, j] 到底在做什么?

如果你喜欢我猜测一下,我会说它会产生一个密集矩阵,然后像往常一样对其进行索引。

【问题讨论】:

  • 是什么让您产生无法索引稀疏矩阵的想法?这样做通常是个坏主意,但很有可能。
  • @user2357112 因为典型的A[i][j] 给了我一个错误。但是A[i,j] 有效——你知道这是在做什么吗?
  • A[i][j] 不应被视为索引任何 NumPy 或 SciPy 数据结构的“典型”方式。它甚至不适用于正则矩阵,更不用说稀疏矩阵了。
  • 给定一个二维 numpy 数组 AA[i][j] 确实是可能的(至少在我的系统上)。你介意告诉我为什么A[i, j]这么慢吗?
  • 对于A[i, j],它为您提供i 行和j 列的值,就像常规矩阵或数组索引一样。它的作用因 SciPy 版本和稀疏矩阵格式而异,但它并没有具体化整个矩阵。那将是一种荒谬的奢侈浪费。您可以从 source code 中的相关 __getitem__ 方法跟踪代码路径,但由于涉及到所有继承,您需要在文件之间跳转很多。

标签: python matrix scipy sparse-matrix


【解决方案1】:
In [211]: M = sparse.csr_matrix(np.eye(3))                                   
In [212]: M                                                                  
Out[212]: 
<3x3 sparse matrix of type '<class 'numpy.float64'>'
    with 3 stored elements in Compressed Sparse Row format>

使用 [0] 进行索引会产生一个新的稀疏矩阵,(1,3) 形状:

In [213]: M[0]                                                               
Out[213]: 
<1x3 sparse matrix of type '<class 'numpy.float64'>'
    with 1 stored elements in Compressed Sparse Row format>

尝试索引,再次给出另一个稀疏矩阵,或错误。那是因为它仍然是 2d 对象 (1,3) 形状。

In [214]: M[0][1]                                                            
---------------------------------------------------------------------------
IndexError                                Traceback (most recent call last)
<ipython-input-214-0661a1f27e52> in <module>
----> 1 M[0][1]

/usr/local/lib/python3.6/dist-packages/scipy/sparse/csr.py in __getitem__(self, key)
    290             # [i, 1:2]
    291             elif isinstance(col, slice):
--> 292                 return self._get_row_slice(row, col)
    293             # [i, [1, 2]]
    294             elif issequence(col):

/usr/local/lib/python3.6/dist-packages/scipy/sparse/csr.py in _get_row_slice(self, i, cslice)
    397 
    398         if i < 0 or i >= M:
--> 399             raise IndexError('index (%d) out of range' % i)
    400 
    401         start, stop, stride = cslice.indices(N)

IndexError: index (1) out of range

使用 [0,1] 语法的索引确实有效,两个数字适用于两个不同的维度:

In [215]: M[0,1]                                                             
Out[215]: 0.0

A[0][1] 确实可以与 np.ndarray 一起使用,但这是因为第一个 [0] 会生成一个维度少 1 的数组。但是np.matrixsparse 返回一个二维矩阵,而不是一维矩阵。这是我们不推荐np.matrix 的原因之一。有了sparse,矩阵的性质就更深了,所以我们不能简单地贬低它。

我们可以通过触发错误来了解从稀疏矩阵中选择元素所涉及的代码:

In [216]: M[0,4]                                                             
---------------------------------------------------------------------------
IndexError                                Traceback (most recent call last)
<ipython-input-216-4919ae565782> in <module>
----> 1 M[0,4]

/usr/local/lib/python3.6/dist-packages/scipy/sparse/csr.py in __getitem__(self, key)
    287             # [i, j]
    288             if isintlike(col):
--> 289                 return self._get_single_element(row, col)
    290             # [i, 1:2]
    291             elif isinstance(col, slice):

/usr/local/lib/python3.6/dist-packages/scipy/sparse/compressed.py in _get_single_element(self, row, col)
    868         if not (0 <= row < M) or not (0 <= col < N):
    869             raise IndexError("index out of bounds: 0<=%d<%d, 0<=%d<%d" %
--> 870                              (row, M, col, N))
    871 
    872         major_index, minor_index = self._swap((row, col))

IndexError: index out of bounds: 0<=0<3, 0<=4<3

===

是的,在稀疏矩阵中索引项目比在密集数组中索引要慢。这不是因为它首先转换为稠密。使用密集数组索引项目只需将 n-d 索引转换为平面索引,并在 1d 平面数据缓冲区中选择所需的字节 - 其中大部分是在快速编译的代码中完成的。但正如你从回溯中看到的那样,从稀疏矩阵中选择一个项目更加复杂,其中很多是 Python。

稀疏lil 格式旨在更快地进行索引(尤其是设置)。但即使这样也比索引密集数组要慢得多。如果您需要迭代或重复访问单个元素,请勿使用稀疏矩阵。

===

要了解索引 M 所涉及的内容,请查看其关键属性:

In [224]: M.data,M.indices,M.indptr                                          
Out[224]: 
(array([1., 1., 1.]),
 array([0, 1, 2], dtype=int32),
 array([0, 1, 2, 3], dtype=int32))

要选择第 0 行,我们必须使用 indptr 从其他切片中选择一个切片:

In [225]: slc = slice(M.indptr[0],M.indptr[1])                               
In [226]: M.data[slc], M.indices[slc]                                        
Out[226]: (array([1.]), array([0], dtype=int32))

然后要选择 col 1,我们必须检查该值是否在 indices[slc] 中。如果是,则返回data[slc] 中的对应元素。如果不返回 0。

对于更复杂的索引,sparse 实际上使用矩阵乘法,从索引创建了一个extractor 矩阵。它还使用乘法来执行行或列求和。

矩阵乘法是一种稀疏矩阵强度 - 只要矩阵足够稀疏。稀疏格式,尤其是csr 的数学根源在于稀疏线性方程问题,例如有限差分和有限元 PDES。

===

这是lil 矩阵的基本属性

In [227]: ml=M.tolil()                                                       
In [228]: ml.data                                                            
Out[228]: array([list([1.0]), list([1.0]), list([1.0])], dtype=object)
In [229]: ml.rows                                                            
Out[229]: array([list([0]), list([1]), list([2])], dtype=object)
In [230]: ml.data[0],ml.rows[0]                                              
Out[230]: ([1.0], [0])          # cf Out[226]

【讨论】: