将n个子矩阵编译成NxN矩阵

Question

研究矩阵结构分析中的一个问题。我正在用Python (使用Anaconda 3)编写一个程序来分析桁架。每个单独的桁架构件生成一个4x4矩阵，总共生成n个4x4矩阵。然后，对于矩阵A、B、C，这些4x4矩阵被编译成NxN矩阵，排列如下：

​

​

正如您所看到的，每个相继的子矩阵都比前一个子矩阵上一行，下一行。此外，由于桁架的大小和桁架节点(节点)的数量由用户指定，因此NxN矩阵的大小必须动态确定(子矩阵始终为4x4)。

我得到了一个NxN零矩阵；我正在尝试找出如何正确编译子矩阵。

我发现了一些类似的问题，但它们都没有动态地缩放较大的矩阵。

我很感谢你们能提供的任何帮助。

Answer 1

n可能很大，所以结果是一个沿对角线集中了非零值的大型稀疏矩阵？稀疏矩阵的设计考虑到了这种矩阵(来自FD和FE PDE问题)。我在MATLAB中做了很多这样的事情，还有一些是用scipy稀疏模块做的。

该模块有一个可能工作的块定义模式，但我更熟悉的是coo到csr的路由。

在coo格式中，非零元素由3个向量定义：i、j和data。您可以在这些数组中收集A、B等的所有值(为B等中的值应用适当的偏移量)，而无需担心重叠。然后，当该格式转换为csr (用于矩阵计算)时，将对重叠的值求和-这正是您想要的。

我认为sparse文档中有一些简单的例子。从概念上讲，要做的最简单的事情是迭代n子矩阵，并收集这3个数组中的值。但我还设计了一个更复杂的系统，它可以作为一个大的数组操作来完成，或者通过在较小的维度上迭代来完成。例如，每个子矩阵有16个值。在现实情况下，16将远远小于n。

为了给出一个更具体的例子，我会玩弄代码。

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

下面是一个包含3个模块的简单示例-功能强大，但不是最高效的

定义3个块：

In [620]: A=np.ones((4,4),int)    
In [621]: B=np.ones((4,4),int)*2
In [622]: C=np.ones((4,4),int)*3

要在其中收集值的列表；可以是数组，但添加或扩展列表很容易，并且相对有效：

In [623]: i, j, dat = [], [], []

In [629]: def foo(A,n):
   # turn A into a sparse, and add it's points to the arrays
   # with an offset of 'n'
   ac = sparse.coo_matrix(A)
   i.extend(ac.row+n)
   j.extend(ac.col+n)
   dat.extend(ac.data)


In [630]: foo(A,0)

In [631]: i
Out[631]: [0, 0, 0, 0, 1, 1, 1, 1, 2, 2, 2, 2, 3, 3, 3, 3]    
In [632]: j
Out[632]: [0, 1, 2, 3, 0, 1, 2, 3, 0, 1, 2, 3, 0, 1, 2, 3]

In [633]: foo(B,1)
In [634]: foo(C,2)  # do this in a loop in the real world

In [636]: M = sparse.csr_matrix((dat,(i,j)))

In [637]: M
Out[637]: 
<6x6 sparse matrix of type '<class 'numpy.int32'>'
    with 30 stored elements in Compressed Sparse Row format>

In [638]: M.A
Out[638]: 
array([[1, 1, 1, 1, 0, 0],
       [1, 3, 3, 3, 2, 0],
       [1, 3, 6, 6, 5, 3],
       [1, 3, 6, 6, 5, 3],
       [0, 2, 5, 5, 5, 3],
       [0, 0, 3, 3, 3, 3]], dtype=int32)

如果我做得对的话，A，B，C的重叠值是相加的。

更一般地说：

In [21]: def foo1(mats):
      i,j,dat = [],[],[]
      for n,mat in enumerate(mats):
          A = sparse.coo_matrix(mat)
          i.extend(A.row+n)
          j.extend(A.col+n)
          dat.extend(A.data)
      M = sparse.csr_matrix((dat,(i,j)))
      return M
   ....:   

In [22]: foo1((A,B,C,B,A)).A
Out[22]: 
array([[1, 1, 1, 1, 0, 0, 0, 0],
       [1, 3, 3, 3, 2, 0, 0, 0],
       [1, 3, 6, 6, 5, 3, 0, 0],
       [1, 3, 6, 8, 7, 5, 2, 0],
       [0, 2, 5, 7, 8, 6, 3, 1],
       [0, 0, 3, 5, 6, 6, 3, 1],
       [0, 0, 0, 2, 3, 3, 3, 1],
       [0, 0, 0, 0, 1, 1, 1, 1]], dtype=int32)

想出一种更有效地做到这一点的方法可能取决于单个子矩阵是如何生成的。如果它们是迭代创建的，那么您也可以迭代地收集i，j，数据值。

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

由于子矩阵是密集的，因此我们可以直接获得适当的i,j,data值，而不需要经过coo中介。如果A,B,C被收集到一个更大的数组中，则不会出现循环。

如果我修改foo1以返回coo矩阵，我会看到给定的i,j,data列表(作为数组)，而不是重复项的求和。在包含5个矩阵的示例中，我得到了80个元素数组，可以将其重塑为

In [110]: f.col.reshape(-1,16)
Out[110]: 
array([[0, 1, 2, 3, 0, 1, 2, 3, 0, 1, 2, 3, 0, 1, 2, 3],
       [1, 2, 3, 4, 1, 2, 3, 4, 1, 2, 3, 4, 1, 2, 3, 4],
       [2, 3, 4, 5, 2, 3, 4, 5, 2, 3, 4, 5, 2, 3, 4, 5],
       [3, 4, 5, 6, 3, 4, 5, 6, 3, 4, 5, 6, 3, 4, 5, 6],
       [4, 5, 6, 7, 4, 5, 6, 7, 4, 5, 6, 7, 4, 5, 6, 7]], dtype=int32)

In [111]: f.row.reshape(-1,16)
Out[111]: 
array([[0, 0, 0, 0, 1, 1, 1, 1, 2, 2, 2, 2, 3, 3, 3, 3],
       [1, 1, 1, 1, 2, 2, 2, 2, 3, 3, 3, 3, 4, 4, 4, 4],
       [2, 2, 2, 2, 3, 3, 3, 3, 4, 4, 4, 4, 5, 5, 5, 5],
       [3, 3, 3, 3, 4, 4, 4, 4, 5, 5, 5, 5, 6, 6, 6, 6],
       [4, 4, 4, 4, 5, 5, 5, 5, 6, 6, 6, 6, 7, 7, 7, 7]], dtype=int32)

In [112]: f.data.reshape(-1,16)
Out[112]: 
array([[1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1],
       [2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2],
       [3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3],
       [2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2],
       [1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1]])

我应该能够在没有循环的情况下生成它们，特别是row和col。

In [143]: mats=[A,B,C,B,A]

数组元素的坐标

In [144]: I,J=[i.ravel() for i in np.mgrid[range(A.shape[0]),range(A.shape[1])]] 

通过广播使用offset复制它们

In [145]: x=np.arange(len(mats))[:,None]
In [146]: I=I+x    
In [147]: J=J+x

将数据收集到一个大型数组中：

In [148]: D=np.concatenate(mats,axis=0)

In [149]: f=sparse.csr_matrix((D.ravel(),(I.ravel(),J.ravel())))

或者作为紧凑函数

def foo3(mats):
    A = mats[0]
    n,m = A.shape
    I,J = np.mgrid[range(n), range(m)]
    x = np.arange(len(mats))[:,None]
    I = I.ravel()+x
    J = J.ravel()+x
    D=np.concatenate(mats,axis=0)
    f=sparse.csr_matrix((D.ravel(),(I.ravel(),J.ravel())))
    return f

在这个简单的例子中，第二个版本快了2倍；第一个版本与列表的长度成线性关系；第二个版本几乎与它的长度无关。

In [158]: timeit foo1(mats)
1000 loops, best of 3: 1.3 ms per loop

In [160]: timeit foo3(mats)
1000 loops, best of 3: 653 µs per loop

Answer 2

简单的for-loop方法是将每个4x4矩阵添加到大零矩阵的适当切片中：

for i, small_m in enumerate(small_matrices):
    big_m[i:i+4, i:i+4] += small_m

您还可以在不使用Python循环的情况下完成此操作，方法是创建零矩阵的步进式视图，并使用np.add.at进行无缓冲加法。如果您的4x4矩阵被打包到一个k×4×4的数组中，这将特别有效：

import numpy as np
from numpy.lib.stride_tricks import as_strided

# Create a view of big_m with the shape of small_matrices.
# strided_view[i] is a view of big_m[i:i+4, i:i+4]
strides = (sum(big_m.strides),) + big_m.strides
strided_view = as_strided(big_m, shape=small_matrices.shape, strides=strides)

np.add.at(strided_view, np.arange(small_matrices.shape[0]), small_matrices)