将N个子矩阵编译为numpy中的NxN矩阵 [英] Compiling n submatrices into an NxN matrix in numpy

问题描述

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

Working on a problem in matrix structural analysis. I am writing a program in Python (using Anaconda 3) to analyze a truss. Each individual truss member generates one 4x4 matrix, for a total of n 4x4 matrices. Then, these 4x4 matrices are compiled into an NxN matrix, arranged like so, for matrices A, B, C:

如您所见，每个连续的子矩阵都比前一个子矩阵上一行，下一行.此外，由于桁架的大小和桁架接头(节点)的数量是由用户指定的，因此必须动态确定NxN矩阵的大小(子矩阵始终为4x4).

As you can see, each successive submatrix is placed one row over and one row down from the preceding one. Further, because the size of the truss and the number of truss joints (nodes) is specified by the user, the size of the NxN matrix has to be determined dynamically (the submatrices are always 4x4).

我有一个NxN零矩阵；我试图弄清楚如何正确地编译子矩阵.

I've got an NxN zero matrix; I am trying to figure out how to compile the submatrices correctly.

我发现了一些类似的问题，但是没有一个问题可以动态缩放较大的矩阵.

I found a few similar questions but none of them scaled the larger matrix dynamically.

感谢您提供的任何帮助.

I appreciate any assistance you folks can provide.

推荐答案

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

Is n potentially large, so the result is a large sparse matrix with nonzero values concentrated along the diagonal? Sparse matrices are designed with this kind of matrix in mind (from FD and FE PDE problems). I did this a lot in MATLAB, and some with the scipy sparse module.

该模块具有可能会起作用的块定义模式，但是我更熟悉的是从coo到csr的路由.

That module has a block definition mode that might work, but what I'm more familiar with is the coo to csr route.

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

In the coo format, nonzero elements are defined by 3 vectors, i, j, and data. You can collect all the values for A, B, etc in these arrays (applying the appropriate offset for the values in B etc), without worrying about overlaps. Then when that format is converted to csr (for matrix calculations) the overlapping values are summed - which is exactly what you want.

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

I think the sparse documentation has some simple examples of this. Conceptually the simplest thing to do is iterate over the n submatrices, and collect the values in those 3 arrays. But I also worked out a more complex system whereby it can be done as one big array operation, or by iterating over a smaller dimension. For example each submatrix has 16 values. In a realistic case 16 will be much smaller than n.

我会尝试使用代码来给出更具体的示例.

I'd have play around with code to give a more concrete example.

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

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

这是一个包含3个块的简单示例-功能性，但不是最有效的

Here's a simple example with 3 blocks - functional, but not the most efficient

定义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

列表，用于收集值；可以是数组，但是添加或扩展列表很容易且相对有效:

lists to collect values in; could be arrays, but it is easy, and relatively efficient to append or extend lists:

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的重叠值相加.

If I've done this right, overlapping values of A,B,C are summed.

更一般地:

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，data值.

Coming up with a way of doing this more efficiently may depend on how the individual submatrices are generated. If they are created iteratively, you might as well collect the i,j,data values iteratively as well.

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

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

由于子矩阵是密集的，因此我们可以直接获取适当的i,j,data值，而无需通过coo中介.而且，如果A,B,C被收集到一个更大的数组中，也不会循环.

Since the submatrices are dense, we can get the appropriate i,j,data values directly, without going through a coo intermediary. And without looping if the A,B,C are collected into one larger array.

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

If I modify foo1 to return a coo matrix, I see the i,j,data lists (as arrays) as given, without summation of duplicates. In the example with 5 matrices, I get 80 element arrays, which can be reshaped as

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.

I should be able generate those without a loop, especially the row and 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])]] 

通过广播用偏移量复制它们

replicate them with offset via broadcasting

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

将数据收集到一个大数组中:

Collect the data into one large array:

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 this modest example the 2nd version is 2x faster; the first scales linearly with the length of the list; the 2nd is almost independent of its length.

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

这篇关于将N个子矩阵编译为numpy中的NxN矩阵的文章就介绍到这了，希望我们推荐的答案对大家有所帮助，也希望大家多多支持IT屋！