N個の部分行列をnumpyでNxN行列にコンパイルする

Question

行列構造解析の問題に取り組んでいます。私は、Pythonで（Anaconda 3を使って）トラスを分析するプログラムを書いています。それぞれの個々のトラス部材は、合計4x4行列のために1つの4x4行列を生成する。次いで、これらの4×4行列は、N×Nのマトリックスにコンパイルされるので、行列A、B、Cのためのように配置された：あなたが見ることができるように

、各連続する部分行列を上1行1行下に配置されています前のものから。さらに、トラスのサイズおよびトラスジョイント（ノード）の数がユーザによって指定されるので、N×Nマトリックスのサイズは動的に決定されなければならない（サブマトリックスは常に4×4である）。

私はNxNゼロ行列を得ました。サブ行列を正しくコンパイルする方法を理解しようとしています。

私はいくつか似たような疑問を見つけましたが、大きな行列を動的にスケーリングしたものはありませんでした。

皆様のご支援をいただき、ありがとうございます。

Answer 1

シンプルfor -loopな方法は、大きなゼロ行列の適切なスライスに各4×4行列を追加することです：

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

あなたはまたのストライドビューを作成することにより、無Pythonのループでこれを行うことができますバッファなし加算にはnp.add.atを使用します。あなたの4×4行列を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) 

Answer 2

n潜在的に大きいので、結果は、に沿って集中非ゼロ値を持つ大規模なスパース行列であります対角線？スパース行列は、この種の行列を念頭に置いて設計されています（FDおよびFE PDEの問題から）。私はこれをMATLABで多く行い、一部はscipyスパースモジュールで行いました。

このモジュールには動作するブロック定義モードがありますが、私がよく知っているのはcooからcsrまでです。 coo形式で

は、ゼロ以外の要素が3つのベクター、i、j、及びdataによって定義されます。重複を心配することなく、A、Bなどのすべての値をこれらの配列に集めることができます（Bの値に適切なオフセットを適用するなど）。その後、その形式がcsr（行列計算のため）に変換されると、重複する値が合計されます。これは正確にあなたが望むものです。

私はsparseドキュメントには、この簡単な例がいくつかあると思います。概念的には、最も簡単なことは、n部分行列を繰り返し実行し、それら3つの配列の値を集めることです。しかし、私はまた、1つの大きな配列操作として、またはより小さい次元を反復することによって、より複雑なシステムを作成しました。例えば、各サブマトリックスは16個の値を有する。現実的なケースでは、nは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が1つの大きな配列に集められている場合はループしません。

Iはcoo行列を返すようにfoo1を変更する場合、所与として、私は重複の加算せず、（配列として）i,j,dataリストを参照してください。 5つのマトリックスを有する例では、私はループのないもの、特にrowを生成しcolできるはず

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]]) 

ように整形することができる80の素子アレイを得ます。

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

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

：

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

放送

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

介しオフセットでそれらを複製する配列の要素のための

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

座標が一つの大きなアレイにデータを収集しますまたはcompacとしてt関数

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番目のバージョンは2倍高速です。最初のスケールはリストの長さに比例して変化します。第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