我想要一个稀疏的n x n 矩阵A，其中B 是一个m x m 稀疏矩阵和n = m * m：
我知道怎么弄B：
data = np.vstack([-np.ones(n), 4 * np.ones(n), -np.ones(n)])
diag = np.array([-1, 0, 1])
B = scipy.sparse.spdiags(data, diag, m, m).tocsr()
当然还有如何制作身份I。不知何故，我不能做同样的事情来制作A。
我知道如何以愚蠢的方式做到这一点：将 B 中的所有数字替换为 A 并以与 B 相同的方式使 A 具有 5 个对角线。但我确实相信scipy.sparse 可以更简单地做到这一点。
我已经花了两个小时思考它，也许找到一个更简单的方法并不值得。
【问题讨论】：
标签： python scipy sparse-matrix
我找到了一个sparse.block_diag 函数。
In [383]: n,m=3,3
In [385]: data = np.vstack([-np.ones(n), 4 * np.ones(n), -np.ones(n)]).astype(int)
...: diag = np.array([-1, 0, 1])
...: B = sparse.spdiags(data, diag, m, m).tocsr()
In [386]: B
Out[386]:
<3x3 sparse matrix of type '<class 'numpy.int64'>'
with 7 stored elements in Compressed Sparse Row format>
In [387]: B.A
Out[387]:
array([[ 4, -1, 0],
[-1, 4, -1],
[ 0, -1, 4]], dtype=int64)
In [388]: I = -sparse.eye(m).astype(int)
block_diag 没有偏移功能，所以我先将B 和I 组合成一个更大的块。
In [389]: M1 = sparse.bmat([[B,I],[I,B]])
In [390]: M1
Out[390]:
<6x6 sparse matrix of type '<class 'numpy.int64'>'
with 20 stored elements in COOrdinate format>
In [391]: M1.A
Out[391]:
array([[ 4, -1, 0, -1, 0, 0],
[-1, 4, -1, 0, -1, 0],
[ 0, -1, 4, 0, 0, -1],
[-1, 0, 0, 4, -1, 0],
[ 0, -1, 0, -1, 4, -1],
[ 0, 0, -1, 0, -1, 4]], dtype=int64)
In [392]: M2 =sparse.block_diag((M1,M1))
In [393]: M2
Out[393]:
<12x12 sparse matrix of type '<class 'numpy.int64'>'
with 40 stored elements in COOrdinate format>
In [394]: M2.A
Out[394]:
array([[ 4, -1, 0, -1, 0, 0, 0, 0, 0, 0, 0, 0],
[-1, 4, -1, 0, -1, 0, 0, 0, 0, 0, 0, 0],
[ 0, -1, 4, 0, 0, -1, 0, 0, 0, 0, 0, 0],
[-1, 0, 0, 4, -1, 0, 0, 0, 0, 0, 0, 0],
[ 0, -1, 0, -1, 4, -1, 0, 0, 0, 0, 0, 0],
[ 0, 0, -1, 0, -1, 4, 0, 0, 0, 0, 0, 0],
[ 0, 0, 0, 0, 0, 0, 4, -1, 0, -1, 0, 0],
[ 0, 0, 0, 0, 0, 0, -1, 4, -1, 0, -1, 0],
[ 0, 0, 0, 0, 0, 0, 0, -1, 4, 0, 0, -1],
[ 0, 0, 0, 0, 0, 0, -1, 0, 0, 4, -1, 0],
[ 0, 0, 0, 0, 0, 0, 0, -1, 0, -1, 4, -1],
[ 0, 0, 0, 0, 0, 0, 0, 0, -1, 0, -1, 4]], dtype=int64)
看block_mat代码；它构造一个嵌套的输入列表和None，并将其传递给bmat。而bmat 的代码将所有输入的coo 属性与适当的偏移量结合起来，为结果矩阵制作coo 样式的输入。 sparse.hstack 和 sparse.vstack 是使用 bmat 的其他示例。
因此，按照这些模型，您可以直接从 B 和 I 构建复合材料。
===
直接5条对角线：
In [405]: data = -np.ones((5,12),int)
In [406]: data[0,:] *= -4
In [407]: offsets=np.array([0,-1,-3,1,3]) # offsets dont have to be in order
In [408]: M4 = sparse.spdiags(data, offsets,12,12)
In [409]: M4
Out[409]:
<12x12 sparse matrix of type '<class 'numpy.int64'>'
with 52 stored elements (5 diagonals) in DIAgonal format>
In [410]: M4.A
Out[410]:
array([[ 4, -1, 0, -1, 0, 0, 0, 0, 0, 0, 0, 0],
[-1, 4, -1, 0, -1, 0, 0, 0, 0, 0, 0, 0],
[ 0, -1, 4, -1, 0, -1, 0, 0, 0, 0, 0, 0],
[-1, 0, -1, 4, -1, 0, -1, 0, 0, 0, 0, 0],
[ 0, -1, 0, -1, 4, -1, 0, -1, 0, 0, 0, 0],
[ 0, 0, -1, 0, -1, 4, -1, 0, -1, 0, 0, 0],
[ 0, 0, 0, -1, 0, -1, 4, -1, 0, -1, 0, 0],
[ 0, 0, 0, 0, -1, 0, -1, 4, -1, 0, -1, 0],
[ 0, 0, 0, 0, 0, -1, 0, -1, 4, -1, 0, -1],
[ 0, 0, 0, 0, 0, 0, -1, 0, -1, 4, -1, 0],
[ 0, 0, 0, 0, 0, 0, 0, -1, 0, -1, 4, -1],
[ 0, 0, 0, 0, 0, 0, 0, 0, -1, 0, -1, 4]])
直接到coo 输入：
In [411]: data,rows,cols = [],[],[]
In [413]: for v,z in zip([-1,-1,4,-1,-1],[-3,-1,0,1,3]):
...: len = 12-abs(z)
...: d = np.ones(len, int)*v
...: r = np.arange(len)+max(0,z)
...: c = np.arange(len)+max(0,-z)
...: data.append(d); rows.append(r); cols.append(c)
...:
In [414]: data
Out[414]:
[array([-1, -1, -1, -1, -1, -1, -1, -1, -1]),
array([-1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1]),
array([4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4]),
...]
In [415]: rows
Out[415]:
[array([0, 1, 2, 3, 4, 5, 6, 7, 8]),
array([ 0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10]),
....]
In [416]: cols
Out[416]:
[array([ 3, 4, 5, 6, 7, 8, 9, 10, 11]),
array([ 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11]),
...]
In [417]: M5 = sparse.coo_matrix((np.hstack(data), (np.hstack(rows), np.hstack(cols))))
In [418]: M5
Out[418]:
<12x12 sparse matrix of type '<class 'numpy.int64'>'
with 52 stored elements in COOrdinate format>
或者只是使用data 列表（和偏移量）：
In [427]: M6=sparse.diags(data, [-3,-1,0,1,3],dtype=int,format='csr')
In [428]: M6
Out[428]:
<12x12 sparse matrix of type '<class 'numpy.int64'>'
with 52 stored elements in Compressed Sparse Row format>
【讨论】：
B，最终矩阵 A（或您的符号中的 M2）的大小必须为 9 x 9。我还不确定如何正确执行此操作，但无论如何我会尝试使用你的方法。感谢B.A，太酷了。我一直用B.toarray()
coo 样式数组（或列表）中填充 5 个对角线并不难（或更慢）。