Numpy 'intelligente' symmetrische Matrix

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

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 == a.T) # for any symmetric matrix

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

Wenn Sie es sich leisten können, die Matrix unmittelbar vor den Berechnungen zu symmetrisieren, sollte Folgendes relativ schnell gehen:

def symmetrize(a):
    """
    Return a symmetrized version of NumPy array a.

    Values 0 are replaced by the array value at the symmetric
    position (with respect to the diagonal), i.e. if a_ij = 0,
    then the returned array a' is such that a'_ij = a_ji.

    Diagonal values are left untouched.

    a -- square NumPy array, such that a_ij = 0 or a_ji = 0, 
    for i != j.
    """
    return a + a.T - numpy.diag(a.diagonal())

Dies funktioniert unter vernünftigen Voraussetzungen (z. B. nicht beides a[0, 1] = 42und der Widerspruch a[1, 0] = 123vor dem Ausführen symmetrize).

Wenn Sie wirklich eine transparente Symmetrisierung benötigen, können Sie die Unterklasse numpy.ndarray in Betracht ziehen und einfach neu definieren __setitem__:

class SymNDArray(numpy.ndarray):
    """
    NumPy array subclass for symmetric matrices.

    A SymNDArray arr is such that doing arr[i,j] = value
    automatically does arr[j,i] = value, so that array
    updates remain symmetrical.
    """

    def __setitem__(self, (i, j), value):
        super(SymNDArray, self).__setitem__((i, j), value)                    
        super(SymNDArray, self).__setitem__((j, i), value)                    

def symarray(input_array):
    """
    Return a symmetrized version of the array-like input_array.

    The returned array has class SymNDArray. Further assignments to the array
    are thus 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 anstelle von Arrays, je nach Ihren Anforderungen). Dieser Ansatz behandelt sogar kompliziertere Zuweisungen, wie z. B. a[:, 1] = -1das korrekte Festlegen von a[1, :]Elementen.

Beachten Sie, dass Python 3 die Möglichkeit des Schreibens entfernt def …(…, (i, j),…)hat. Daher muss der Code vor der Ausführung mit Python 3 leicht angepasst werden: def __setitem__(self, indexes, value): (i, j) = indexes…

Das allgemeinere Problem der optimalen Behandlung symmetrischer Matrizen in Numpy nervte mich ebenfalls.

Nach eingehender Prüfung denke ich, dass die Antwort wahrscheinlich darin besteht, dass Numpy durch das Speicherlayout, das von den zugrunde liegenden BLAS-Routinen für symmetrische Matrizen unterstützt wird, etwas eingeschränkt wird.

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

Einige der scipy.linalgRoutinen akzeptieren Flags (wie sym_pos=Trueon linalg.solve), die an BLAS-Routinen weitergegeben werden, obwohl eine stärkere Unterstützung in numpy sinnvoll wäre, insbesondere Wrapper für Routinen wie DSYRK (Symmetric Rank K Update), die eine Gram-Matrix ermöglichen würden ein gutes Stück schneller als Punkt (MT, M) berechnet werden.

(Es mag nicht schwierig erscheinen, sich Gedanken über die Optimierung auf einen 2x konstanten Faktor für Zeit und / oder Raum zu machen, aber es kann einen Unterschied zu dieser Schwelle machen, wie groß ein Problem ist, das Sie auf einer einzelnen Maschine bewältigen können ...)

Es gibt eine Reihe bekannter Möglichkeiten zum Speichern symmetrischer Matrizen, sodass sie nicht n ^ 2 Speicherelemente belegen müssen. Darüber hinaus ist es möglich, allgemeine Vorgänge neu zu schreiben, um auf diese überarbeiteten Speichermittel 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 zitieren sie von der Form (1.2.2) in einer symmetrischen Matrix müssen Speicher nur A = [a_{i,j} ]für i >= j. Dann wird die unter der Annahme Vektor der Matrix hält , wird V bezeichnet, und daß ein n-von-n, PUT a_{i,j}in

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

Dies setzt eine 1-Indizierung voraus.

Golub und Van Loan bieten einen Algorithmus 1.2.3 an, der zeigt, wie auf ein solches gespeichertes V zugegriffen werden kann, um zu berechnen y = V x + y.

Golub und Van Loan bieten auch die Möglichkeit, eine Matrix in diagonal dominanter Form zu speichern. Dies spart keinen Speicherplatz, unterstützt jedoch den sofortigen Zugriff für bestimmte andere Arten von Vorgängen.

Dies ist einfaches Python und nicht numpy, aber ich habe gerade eine Routine zusammengestellt, um eine symmetrische Matrix zu füllen (und ein Testprogramm, um sicherzustellen, dass es korrekt ist):

import random

# fill a symmetric matrix with costs (i.e. 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]es ausgefüllt ist. Die Speicherfrage ist etwas interessanter. Man kann die numpy-Array-Klasse mit einem packedAttribut 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[:(m-i)]
                    A[i:, i] = l[:(m-i)]
                    l = l[(m-i):]
                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)

`` `

So erstellen Sie eine NxN-Matrix, die entlang der Hauptdiagonale symmetrisch ist und auf der Hauptdiagonale Nullen enthält:

a = np.array([1, 2, 3, 4, 5])
b = np.zeros(shape=(a.shape[0], a.shape[0]))
upper = np.triu(b + a)
lower = np.tril(np.transpose(b + a))
D = (upper + lower) * (np.full(a.shape[0], fill_value=1) - np.eye(a.shape[0]))

Dies ist eine Art Sonderfall, aber kürzlich habe ich diese Art von Matrix für die Darstellung der Netzwerkadjazenz verwendet.

Hoffentlich hilft das. Prost.