如何重复使用block_diag

Question

我的问题很简单，但仍然无法解决。我想要一个块对角线 n^2*n^2 矩阵。这些块是稀疏的 n*n 矩阵，只有对角线，首先是对角线，然后是对角线。对于 n=4 的简单情况，这很容易完成

datanew = ones((5,n1))
datanew[2] = -2*datanew[2]
diagsn = [-4,-1,0,1,4]
DD2 = sparse.spdiags(datanew,diagsn,n,n)
new = sparse.block_diag([DD2,DD2,DD2,DD2])

因为这只对小 n 有用，有没有更好的方法来使用 block_diag？想到 n -> 1000

Answer 1

构建一长串 DD2 矩阵的简单方法是使用列表理解：

In [128]: sparse.block_diag([DD2 for _ in range(20)]).A
Out[128]: 
array([[-2,  1,  0, ...,  0,  0,  0],
       [ 1, -2,  1, ...,  0,  0,  0],
       [ 0,  1, -2, ...,  0,  0,  0],
       ..., 
       [ 0,  0,  0, ..., -2,  1,  0],
       [ 0,  0,  0, ...,  1, -2,  1],
       [ 0,  0,  0, ...,  0,  1, -2]])

In [129]: _.shape
Out[129]: (80, 80)

至少在我的版本中，block_diag 想要一个数组列表，而不是 *args:

In [133]: sparse.block_diag(DD2,DD2,DD2,DD2)
...
TypeError: block_diag() takes at most 3 arguments (4 given)

In [134]: sparse.block_diag([DD2,DD2,DD2,DD2])
Out[134]: 
<16x16 sparse matrix of type '<type 'numpy.int32'>'
    with 40 stored elements in COOrdinate format>

这可能不是构建这种块对角线阵列的最快方法，但它是一个开始。

================

查看 sparse.block_mat 的代码，我推断它确实如此：

In [145]: rows=[]
In [146]: for i in range(4):
    arow=[None]*4
    arow[i]=DD2
    rows.append(arow)
   .....:     

In [147]: rows
Out[147]: 
[[<4x4 sparse matrix of type '<type 'numpy.int32'>'
    with 10 stored elements (5 diagonals) in DIAgonal format>,
  None,
  None,
  None],
 [None,
  <4x4 sparse matrix of type '<type 'numpy.int32'>'
  ...
  None,
  <4x4 sparse matrix of type '<type 'numpy.int32'>'
    with 10 stored elements (5 diagonals) in DIAgonal format>]]

换句话说，rows 是 None 的 'matrix'，对角线上有 DD2。然后它将这些传递给 sparse.bmat.

In [148]: sparse.bmat(rows)
Out[148]: 
<16x16 sparse matrix of type '<type 'numpy.int32'>'
    with 40 stored elements in COOrdinate format>

bmat 依次从所有输入矩阵的 coo 格式中收集 data,rows,cols，将它们连接到主数组中，并从中构建一个新的 coo 矩阵他们。

因此，另一种方法是直接构建这 3 个数组。

如何重复使用block_diag

How to use block_diag repeatedly

scipy

sparse-matrix

python-2.7