为什么创建此块矩阵/数组时出错?(SciPy和NumPy)

umuewwlo  于 2023-01-02  发布在  其他
关注(0)|答案(1)|浏览(258)
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允许数据数组是二维的。
文档显示,只要定义了块大小,它就应该工作。知道为什么会出现这个错误吗?
实际上我的主要任务是沿着对角线迭代M1,然后在M1下面迭代M2,得到一个300 x 300的矩阵,有人有更好的方法吗?

eivnm1vs

eivnm1vs1#

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

使用FULL回溯的函数调用:

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 ~\anaconda3\lib\site-packages\scipy\sparse\_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 ~\anaconda3\lib\site-packages\scipy\sparse\_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 ~\anaconda3\lib\site-packages\scipy\sparse\_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 ~\anaconda3\lib\site-packages\scipy\sparse\_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 ~\anaconda3\lib\site-packages\scipy\sparse\_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是一个一维数组,而不是一个三维数组,所以ij的值是普通的coo的行和列。
如果你想提供块和块"坐标",你必须使用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)

还有一个块对角输入:
在[83]中:M1 =稀疏块诊断([M1,M1,M1],格式="bsr")
在[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个数组中,然后使用这些数组创建一个新的大cooblock_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]])

相关问题