Kompilieren von n Submatrizen in eine NxN-Matrix in numpy

Question

Arbeiten an einem Problem in der Matrixstrukturanalyse. Ich schreibe ein Programm in Python (mit Anaconda 3), um ein Fachwerk zu analysieren. Jedes einzelne Fachwerkelement erzeugt eine 4x4-Matrix für insgesamt n 4x4-Matrizen. Dann werden diese 4x4-Matrizen in eine NxN Matrix zusammengestellt, wie so angeordnet sind, für Matrices A, B, C:

Wie man sehen kann, die jeweils aufeinanderfolgende Submatrix wird eine Zeile über platziert und eine Reihe nach unten von dem vorhergehenden. Da die Größe des Fachwerks und die Anzahl der Fachwerkverbindungen (Knoten) vom Benutzer festgelegt wird, muss die Größe der NxN-Matrix dynamisch bestimmt werden (die Untermatrizen sind immer 4x4).

Ich habe eine NxN Nullmatrix; Ich versuche herauszufinden, wie man die Submatrizen korrekt kompiliert.

Ich fand ein paar ähnliche Fragen, aber keiner von ihnen skalierte die größere Matrix dynamisch.

Ich schätze jede Hilfe, die Sie Leute bieten können.

Answer 1

Die einfache for -loop Weise wäre jedes 4x4-Matrix zu einer entsprechenden Scheibe der großen Nullmatrix hinzuzufügen:

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

Sie könnten auch ohne Python Schleifen tun dies durch eine strided Blick auf die Schaffung Nullmatrix und Verwendung von np.add.at für ungepufferte Addition. Dies sollte besonders effizient sein, wenn Ihre 4x4 Matrizen zu einer k-by-4-by-4-Anordnung gepackt sind:

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

Ist n potenziell groß, so ist das Ergebnis eine große Sparse Matrix mit Nicht-Null-Werten entlang der konzentrierten Diagonale? Sparse-Matrizen werden mit dieser Art von Matrix entworfen (aus FD- und FE-PDE-Problemen). Ich habe dies viel in MATLAB getan, und einige mit dem spärlichen Modul scipy.

Dieses Modul hat einen Blockdefinitionsmodus, der funktionieren könnte, aber was mir vertrauter ist, ist die Route coo zu csr.

Im coo Format werden von Null verschiedene Elemente, die durch 3-Vektoren definiert sind, i, j und data. Sie können alle Werte für A, B usw. in diesen Arrays sammeln (indem Sie den entsprechenden Offset für die Werte in B usw. anwenden), ohne sich über Überlappungen Gedanken machen zu müssen. Wenn dann dieses Format in csr (für Matrixberechnungen) konvertiert wird, werden die überlappenden Werte summiert - was genau das ist, was Sie wollen.

Ich denke, die sparse Dokumentation hat einige einfache Beispiele dafür. Vom Konzept her ist es am einfachsten, über die n Submatrizen zu iterieren und die Werte in diesen 3 Arrays zu sammeln. Aber ich habe auch ein komplexeres System ausgearbeitet, bei dem es sich um eine große Array-Operation oder um eine kleinere Dimension handelt. Zum Beispiel hat jede Untermatrix 16 Werte. In einem realistischen Fall wird 16 viel kleiner als n sein.

Ich hätte mit Code herumspielen, um ein konkreteres Beispiel zu geben.

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

Hier ist ein einfaches Beispiel mit 3 Blöcken - funktional, aber nicht die effizienteste

definieren 3 Blöcke:

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 

Listen in Werte zu sammeln; sein Arrays könnte, aber es ist einfach und relativ effizient anhängen oder Listen erweitern:

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) 

Wenn ich das richtig gemacht haben, überlappende Werte von A, B, C werden summiert.

Allgemeiner:

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) 

kommend mit einem Weg, dies zu tun effizienter abhängen, wie sich die einzelnen Untermatrizen erzeugt werden. Wenn sie iterativ erstellt werden, können Sie die i, j, Datenwerte auch iterativ sammeln.

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

Da die Untermatrizen dicht sind, können wir die entsprechenden i,j,data Werte direkt erhalten, ohne durch einen coo Vermittler gehen. Und ohne Schleife, wenn die A,B,C in einem größeren Array gesammelt werden.

Wenn I foo1 modifizieren, um eine Matrix coo zurückzukehren, sehe ich die i,j,data Listen (als Arrays), wie angegeben, ohne Summierung von Duplikaten. In dem Beispiel mit der 5-Matrices, erhalten I 80 Elementanordnungen, die

neu geformt werden können, wie

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

sollte ich die ohne eine Schleife erzeugen können, insbesondere die row und col.

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

die Koordinaten für die Elemente eines Arrays

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

sie über Rundfunk

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

sammeln, die Daten in eine große Anordnung mit versetzten Replizieren:

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

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

oder als Kompak t Funktion

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 

In diesem bescheidenen Beispiel ist die 2. Version 2x schneller; der erste skaliert linear mit der Länge der Liste; der 2. ist fast unabhängig von seiner Länge.

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