問題描述

處理矩陣結構分析中的一個問題.我正在用 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

這篇關于在numpy中將n個子矩陣編譯成NxN矩陣的文章就介紹到這了，希望我們推薦的答案對大家有所幫助，也希望大家多多支持html5模板網！