在numpy中将n个子矩阵编译成NxN矩阵
问题描述
处理矩阵结构分析中的一个问题.我正在用 Python(使用 Anaconda 3)编写一个程序来分析桁架.每个单独的桁架成员生成一个 4x4 矩阵，总共有 n 个 4x4 矩阵.然后，这些 4x4 矩阵被编译成一个 NxN 矩阵，排列如下，对于矩阵 A、B、C:
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())))
或者作为一个紧凑的函数
or as a compact function
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
相关文章