Python-Scipy 稀疏矩阵 - A[i, j] 在做什么?
Python-Scipy sparse Matrices - what is A[i, j] doing?
根据我之前的问题(Python - Multiply sparse matrix row with non-sparse vector by index),稀疏矩阵的直接索引是不可能的(至少如果你不想使用 sparse.csr
定义了矩阵,data
、indices
、indptr
)。
但我刚刚发现,给定一个 csr-sparse 矩阵 A
,这个操作工作正常并产生正确的结果:A[i, j]
。
我还注意到:它非常慢,甚至比使用密集矩阵还要慢。
我找不到有关此索引方法的任何信息,所以我想知道:A[i, j]
到底在做什么?
如果你想让我猜一猜,我会说它正在生成一个密集矩阵,然后像往常一样对其进行索引。
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>
再次尝试索引会给出另一个稀疏矩阵或错误。那是因为它仍然是一个二维对象 (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.matrix
和 sparse
return 是二维矩阵,而不是一维矩阵。这是我们不推荐 np.matrix
的原因之一。随着 sparse
矩阵性质变得更深,所以我们不能简单地描述它。
我们可以通过触发错误了解 select 从稀疏矩阵中获取元素所涉及的代码:
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 索引转换为平面索引,并 select 在 1d 平面数据缓冲区中输入所需的字节 - 其中大部分是在快速编译代码中完成的。但是正如您从回溯中看到的那样,select从稀疏矩阵中提取一个项目更为复杂,而且很多 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
到 select 来自其他部分的切片:
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))
然后要选择列 1,我们必须检查该值是否在 indices[slc]
中。如果是,return data[slc]
中对应的元素。如果不是 return 0.
对于更复杂的索引,稀疏实际上使用矩阵乘法,从索引创建了一个 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]
根据我之前的问题(Python - Multiply sparse matrix row with non-sparse vector by index),稀疏矩阵的直接索引是不可能的(至少如果你不想使用 sparse.csr
定义了矩阵,data
、indices
、indptr
)。
但我刚刚发现,给定一个 csr-sparse 矩阵 A
,这个操作工作正常并产生正确的结果:A[i, j]
。
我还注意到:它非常慢,甚至比使用密集矩阵还要慢。
我找不到有关此索引方法的任何信息,所以我想知道:A[i, j]
到底在做什么?
如果你想让我猜一猜,我会说它正在生成一个密集矩阵,然后像往常一样对其进行索引。
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>
再次尝试索引会给出另一个稀疏矩阵或错误。那是因为它仍然是一个二维对象 (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.matrix
和 sparse
return 是二维矩阵,而不是一维矩阵。这是我们不推荐 np.matrix
的原因之一。随着 sparse
矩阵性质变得更深,所以我们不能简单地描述它。
我们可以通过触发错误了解 select 从稀疏矩阵中获取元素所涉及的代码:
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 索引转换为平面索引,并 select 在 1d 平面数据缓冲区中输入所需的字节 - 其中大部分是在快速编译代码中完成的。但是正如您从回溯中看到的那样,select从稀疏矩阵中提取一个项目更为复杂,而且很多 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
到 select 来自其他部分的切片:
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))
然后要选择列 1,我们必须检查该值是否在 indices[slc]
中。如果是,return data[slc]
中对应的元素。如果不是 return 0.
对于更复杂的索引,稀疏实际上使用矩阵乘法,从索引创建了一个 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]