为什么在创建这个块矩阵/数组时出现错误?(SciPy和NumPy)

Question

import scipy as sp
import numpy as np

def H1(M1, M2, N):
row = [0, 1]
col = [0, 0]
data = [M1, M2]
return sp.sparse.bsr_matrix((data, (row,col)), blocksize=(2,2), shape=(N,N)).toarray()

M1 = np.array([[1,1], [1,1]])
M2 = np.array([[2,2], [2,2]])
print(H1(M1, M2, 6))

line 248, in getnnz
raise ValueError('row, column, and data arrays must be 1-D')
ValueError: row, column, and data arrays must be 1-D

我不太确定为什么这个错误会发生，因为bsr_matrix允许数据数组为2-D。知道为什么吗?

文档显示，只要定义了块大小，它就应该工作。知道为什么会发生这个错误吗?

对于一个300 x 300的矩阵，我的主要任务实际上是沿着对角线迭代M1和M1以下的M2。有人有更好的方法吗?

Answer 1

In [70]: data
Out[70]: 
[array([[1, 1],
[1, 1]]),
array([[2, 2],
[2, 2]])] 
In [71]: row, col
Out[71]: ([0, 1], [0, 0])
In [72]: N=6

带完整回溯的函数调用:

In [73]: sparse.bsr_matrix((data, (row,col)), blocksize=(2,2), shape=(N,N))
---------------------------------------------------------------------------
ValueError                                Traceback (most recent call last)
Cell In[73], line 1
----> 1 sparse.bsr_matrix((data, (row,col)), blocksize=(2,2), shape=(N,N))
File ~anaconda3libsite-packagesscipysparse_bsr.py:157, in bsr_matrix.__init__(self, arg1, shape, dtype, copy, blocksize)
152     self.indptr = np.zeros(M//R + 1, dtype=idx_dtype)
154 elif len(arg1) == 2:
155     # (data,(row,col)) format
156     self._set_self(
--> 157         self._coo_container(arg1, dtype=dtype, shape=shape).tobsr(
158             blocksize=blocksize
159         )
160     )
162 elif len(arg1) == 3:
163     # (data,indices,indptr) format
164     (data, indices, indptr) = arg1
File ~anaconda3libsite-packagesscipysparse_coo.py:196, in coo_matrix.__init__(self, arg1, shape, dtype, copy)
193 if dtype is not None:
194     self.data = self.data.astype(dtype, copy=False)
--> 196 self._check()
File ~anaconda3libsite-packagesscipysparse_coo.py:281, in coo_matrix._check(self)
278 self.col = np.asarray(self.col, dtype=idx_dtype)
279 self.data = to_native(self.data)
--> 281 if self.nnz > 0:
282     if self.row.max() >= self.shape[0]:
283         raise ValueError('row index exceeds matrix dimensions')
File ~anaconda3libsite-packagesscipysparse_base.py:299, in spmatrix.nnz(self)
291 @property
292 def nnz(self):
293     """Number of stored values, including explicit zeros.
294 
295     See also
296     --------
297     count_nonzero : Number of non-zero entries
298     """
--> 299     return self.getnnz()
File ~anaconda3libsite-packagesscipysparse_coo.py:248, in coo_matrix.getnnz(self, axis)
243         raise ValueError('row, column, and data array must all be the '
244                          'same length')
246     if self.data.ndim != 1 or self.row.ndim != 1 or 
247             self.col.ndim != 1:
--> 248         raise ValueError('row, column, and data arrays must be 1-D')
250     return int(nnz)
252 if axis < 0:
ValueError: row, column, and data arrays must be 1-D

关键在于，对于2元素的arg1，它尝试制作coo矩阵，然后将其转换为bsr。它不以特殊的bsr方式处理输入。文档应该更清晰。

self._coo_container(arg1, dtype=dtype, shape=shape).tobsr(
158             blocksize=blocksize

换句话说，当它说a[ij[0, k], ij[1, k]] = data[k]时，data是一个一维的值数组，而不是一个n块的3d数组。所以ij的值是普通的coo行和col.

如果你想提供块和块'坐标'，你必须使用csr样式的输入。

coo_matrix有类似的描述，A[i[k], j[k]] = data[k]。这支持了bsr_matrix((data, ij), [blocksize=(R,C), shape=(M, N)])就是coo_matrix((data,ij), shape=(M,N)).tobsr(blocksize=(R,C))的想法

---

你对在对角线上迭代M1等的描述不明显。这听起来就像你想说:

In [74]: M=sparse.bmat([[M1,None,None],[M2,M1,None],[None,M2,M1]])
In [75]: M
Out[75]: 
<6x6 sparse matrix of type '<class 'numpy.intc'>'
with 20 stored elements in COOrdinate format>
In [76]: M.A
Out[76]: 
array([[1, 1, 0, 0, 0, 0],
[1, 1, 0, 0, 0, 0],
[2, 2, 1, 1, 0, 0],
[2, 2, 1, 1, 0, 0],
[0, 0, 2, 2, 1, 1],
[0, 0, 2, 2, 1, 1]], dtype=int32)

转换为bsr，我们可以看到bsr_matrix输入应该是什么样子:

In [78]: Mb = M.tobsr(blocksize=(2,2))    
In [79]: Mb
Out[79]: 
<6x6 sparse matrix of type '<class 'numpy.intc'>'
with 20 stored elements (blocksize = 2x2) in Block Sparse Row format>
In [80]: Mb.data
Out[80]: 
array([[[1, 1],
[1, 1]],
[[2, 2],
[2, 2]],
[[1, 1],
[1, 1]],
[[2, 2],
[2, 2]],
[[1, 1],
[1, 1]]], dtype=int32)
In [81]: Mb.indptr
Out[81]: array([0, 1, 3, 5], dtype=int32)
In [82]: Mb.indices
Out[82]: array([0, 0, 1, 1, 2], dtype=int32)

还有一个对角线输入:

In [83]: M1 = sparse。block_diag (M1, M1, M1,格式= bsr)

In [84]: M1

Out[84]: 
<6x6 sparse matrix of type '<class 'numpy.intc'>'
with 12 stored elements (blocksize = 2x2) in Block Sparse Row format>
In [85]: M1.data, M1.indptr, M1.indices
Out[85]: 
(array([[[1, 1],
[1, 1]],

[[1, 1],
[1, 1]],

[[1, 1],
[1, 1]]], dtype=int32),
array([0, 1, 2, 3], dtype=int32),
array([0, 1, 2], dtype=int32))
In [86]: M1.A
Out[86]: 
array([[1, 1, 0, 0, 0, 0],
[1, 1, 0, 0, 0, 0],
[0, 0, 1, 1, 0, 0],
[0, 0, 1, 1, 0, 0],
[0, 0, 0, 0, 1, 1],
[0, 0, 0, 0, 1, 1]], dtype=int32)

您可能会发现sparse.bmat的代码很有指导意义。它构建所有块的coo版本，然后将它们的data,row,col连接到3个数组中，然后使用它们来创建一个新的大coo。block_diag也使用了这个集合coo方法。

'raw'coo输入是稀疏矩阵的原始输入形式，并且仍然是最基本的样式。几年前，当我在MATLAB中使用FEM时，有效地构建这三个coo输入是制作刚度矩阵的最大部分。

编辑

csr样式的输入将是:

In [100]: M=sparse.bsr_matrix((data, [0,0],[0,1,2,2]),shape=(6,6))    
In [101]: M
Out[101]: 
<6x6 sparse matrix of type '<class 'numpy.int32'>'
with 8 stored elements (blocksize = 2x2) in Block Sparse Row format>
In [102]: M.A
Out[102]: 
array([[1, 1, 0, 0, 0, 0],
[1, 1, 0, 0, 0, 0],
[2, 2, 0, 0, 0, 0],
[2, 2, 0, 0, 0, 0],
[0, 0, 0, 0, 0, 0],
[0, 0, 0, 0, 0, 0]])