Numpy 'smart' symmetrische Matrix

Gibt es eine intelligente und platzsparende symmetrische Matrix in numpy, die automatisch (und transparent) die Position bei [j][i] füllt, wenn [i][j] geschrieben wird?

 import numpy a = numpy.symmetric((3, 3)) a[0][1] = 1 a[1][0] == a[0][1] # True print(a) # [[0 1 0], [1 0 0], [0 0 0]] assert numpy.all(a == aT) # for any symmetric matrix 

Ein automatischer Hermitianer wäre auch nett, obwohl ich das nicht zum Zeitpunkt des Schreibens brauche.

  • Erzeugen eine gemusterte numpy Matrix
  • Wie spalte ich eine Matrix in 4 Blöcke mit numpy?
  • Eine Loopless 3D Matrix Multiplikation in Python
  • So finden Sie degenerierte Zeilen / Spalten in einer Kovarianzmatrix
  • Sehr große Matrizen mit Python und NumPy
  • Tensor Punkt Betrieb in Python
  • 5 Solutions collect form web for “Numpy 'smart' symmetrische Matrix”

    Wenn Sie es sich leisten können, die Matrix zu simulieren, bevor Sie Berechnungen durchführen, sollte das Folgende relativ schnell sein:

     def symmetrize(a): return a + aT - numpy.diag(a.diagonal()) 

    Dies funktioniert unter vernünftigen Annahmen (z. B. nicht beides a[0, 1] = 42 und das widersprüchliche a[1, 0] = 123 vor dem symmetrize ).

    Wenn Sie wirklich eine transparente Symmetrisierung benötigen, können Sie die Unterklassen von numpy.ndarray betrachten und einfach neu definieren __setitem__ :

     class SymNDArray(numpy.ndarray): def __setitem__(self, (i, j), value): super(SymNDArray, self).__setitem__((i, j), value) super(SymNDArray, self).__setitem__((j, i), value) def symarray(input_array): """ Returns a symmetrized version of the array-like input_array. Further assignments to the array are automatically symmetrized. """ return symmetrize(numpy.asarray(input_array)).view(SymNDArray) # Example: a = symarray(numpy.zeros((3, 3))) a[0, 1] = 42 print a # a[1, 0] == 42 too! 

    (Oder das Äquivalent mit Matrizen statt Arrays, je nach Ihren Bedürfnissen). Dieser Ansatz behandelt sogar kompliziertere Zuordnungen wie a[:, 1] = -1 , das korrekt a[1, :] Element setzt.

    Beachten Sie, dass Python 3 die Möglichkeit des Schreibens von def …(…, (i, j),…) , also muss der Code etwas angepasst werden, bevor er mit Python 3 läuft: def __setitem__(self, indexes, value): (i, j) = indexes

    Die allgemeinere Frage der optimalen Behandlung von symmetrischen Matrizen in numpy bugged mich auch.

    Nachdem ich hingegangen bin, denke ich, dass die Antwort wahrscheinlich ist, dass das Numpy etwas durch das Gedächtnis-Layout, das von den zugrunde liegenden BLAS-Routinen für symmetrische Matrizen unterstützt wird, beschränkt ist.

    Während einige BLAS-Routinen Symmetrie ausnutzen, um Berechnungen auf symmetrischen Matrizen zu beschleunigen, verwenden sie immer noch dieselbe Speicherstruktur wie eine vollständige Matrix, dh n^2 Raum statt n(n+1)/2 . Nur man sagt ihnen, dass die Matrix symmetrisch ist und nur die Werte im oberen oder unteren Dreieck verwenden darf.

    Einige der scipy.linalg Routinen akzeptieren Flaggen (wie sym_pos=True auf linalg.solve ), die an BLAS-Routinen weitergegeben werden, obwohl mehr Unterstützung für diese in numpy wäre schön, insbesondere Wrapper für Routinen wie DSYRK (symmetrischer Rang k Update), die eine Grammatrix erlauben würde, ein fairer kleiner als dot (MT, M) zu berechnen.

    (Vielleicht scheinen nitpicky Sorgen über die Optimierung für einen 2x konstanten Faktor auf Zeit und / oder Raum, aber es kann einen Unterschied zu dieser Schwelle, wie groß ein Problem, das Sie auf einer einzigen Maschine zu verwalten …)

    Es gibt eine Reihe von bekannten Arten der Speicherung von symmetrischen Matrizen, so dass sie nicht brauchen, um n ^ 2 Speicherelemente zu besetzen. Darüber hinaus ist es möglich, gemeinsame Operationen umzuschreiben, um auf diese überarbeiteten Lagerräume zuzugreifen. Die endgültige Arbeit ist Golub und Van Loan, Matrix Computations , 3. Auflage 1996, Johns Hopkins University Press, Abschnitte 1.27-1.2.9. Zum Beispiel, zitiert sie aus Form (1.2.2), in einer symmetrischen Matrix müssen nur A = [a_{i,j} ] für i >= j speichern. Dann wird angenommen, daß der Vektor , der die Matrix hält, mit V bezeichnet wird und daß A n-by-n ist, setzen a_{i,j} in

     V[(j-1)n - j(j-1)/2 + i] 

    Dies setzt eine 1-indexierung voraus.

    Golub und Van Loan bieten einen Algorithmus 1.2.3 an, der zeigt, wie man auf einen solchen gespeicherten V zugreift, um y = V x + y zu berechnen.

    Golub und Van Loan bieten auch eine Möglichkeit, eine Matrix in diagonal dominanter Form zu speichern. Dies spart nicht Speicherplatz, sondern unterstützt den beruflichen Zugriff für bestimmte andere Arten von Operationen.

    Das ist einfach python und nicht numpy, aber ich warf nur eine Routine zusammen, um eine symmetrische Matrix zu füllen (und ein Testprogramm, um sicherzustellen, dass es richtig ist):

     import random # fill a symmetric matrix with costs (ie m[x][y] == m[y][x] # For demonstration purposes, this routine connect each node to all the others # Since a matrix stores the costs, numbers are used to represent the nodes # so the row and column indices can represent nodes def fillCostMatrix(dim): # square array of arrays # Create zero matrix new_square = [[0 for row in range(dim)] for col in range(dim)] # fill in main diagonal for v in range(0,dim): new_square[v][v] = random.randrange(1,10) # fill upper and lower triangles symmetrically by replicating diagonally for v in range(1,dim): iterations = dim - v x = v y = 0 while iterations > 0: new_square[x][y] = new_square[y][x] = random.randrange(1,10) x += 1 y += 1 iterations -= 1 return new_square # sanity test def test_symmetry(square): dim = len(square[0]) isSymmetric = '' for x in range(0, dim): for y in range(0, dim): if square[x][y] != square[y][x]: isSymmetric = 'NOT' print "Matrix is", isSymmetric, "symmetric" def showSquare(square): # Print out square matrix columnHeader = ' ' for i in range(len(square)): columnHeader += ' ' + str(i) print columnHeader i = 0; for col in square: print i, col # print row number and data i += 1 def myMain(argv): if len(argv) == 1: nodeCount = 6 else: try: nodeCount = int(argv[1]) except: print "argument must be numeric" quit() # keep nodeCount <= 9 to keep the cost matrix pretty costMatrix = fillCostMatrix(nodeCount) print "Cost Matrix" showSquare(costMatrix) test_symmetry(costMatrix) # sanity test if __name__ == "__main__": import sys myMain(sys.argv) # vim:tabstop=8:shiftwidth=4:expandtab 

    Es ist trivial, pythonisch auszufüllen [i][j] wenn [j][i] ausgefüllt ist. Die Speicherfrage ist ein wenig interessanter. Man kann die numpy Array-Klasse mit einem packed Attribut erweitern, das sowohl zum Speichern von Speicher als auch zum späteren Lesen der Daten nützlich ist.

     class Sym(np.ndarray): # wrapper class for numpy array for symmetric matrices. New attribute can pack matrix to optimize storage. # Usage: # If you have a symmetric matrix A as a shape (n,n) numpy ndarray, Sym(A).packed is a shape (n(n+1)/2,) numpy array # that is a packed version of A. To convert it back, just wrap the flat list in Sym(). Note that Sym(Sym(A).packed) def __new__(cls, input_array): obj = np.asarray(input_array).view(cls) if len(obj.shape) == 1: l = obj.copy() p = obj.copy() m = int((np.sqrt(8 * len(obj) + 1) - 1) / 2) sqrt_m = np.sqrt(m) if np.isclose(sqrt_m, np.round(sqrt_m)): A = np.zeros((m, m)) for i in range(m): A[i, i:] = l[:(mi)] A[i:, i] = l[:(mi)] l = l[(mi):] obj = np.asarray(A).view(cls) obj.packed = p else: raise ValueError('One dimensional input length must be a triangular number.') elif len(obj.shape) == 2: if obj.shape[0] != obj.shape[1]: raise ValueError('Two dimensional input must be a square matrix.') packed_out = [] for i in range(obj.shape[0]): packed_out.append(obj[i, i:]) obj.packed = np.concatenate(packed_out) else: raise ValueError('Input array must be 1 or 2 dimensional.') return obj def __array_finalize__(self, obj): if obj is None: return self.packed = getattr(obj, 'packed', None) 

    “ `

    Python ist die beste Programmiersprache der Welt.