Matrice symétrique "intelligente" Numpy

La question plus générale du traitement optimal des matrices symétriques dans numpy m'a aussi dérangé.

Après l'avoir examiné, je pense que la réponse est probablement que numpy est quelque peu contraint par le support de la disposition de la mémoire par les routines Blas sous-jacentes pour les matrices symétriques.

Alors que certaines routines BLAS exploitent la symétrie pour accélérer les calculs sur des matrices symétriques, elles utilisent toujours la même structure de mémoire qu'une matrice complète, c'est-à-dire n^2 espace plutôt que n(n+1)/2. Juste qu'ils on vous dit que la matrice est symétrique et d'utiliser uniquement les valeurs dans le triangle supérieur ou inférieur.

Certaines des routines scipy.linalg acceptent les drapeaux (comme sym_pos=True sur linalg.solve) qui sont transmis aux routines BLAS, bien que plus de support pour cela dans numpy serait bien, en particulier les wrappers pour les routines comme DSYRK (symetric rank K update), ce qui permettrait de calculer une matrice Gram un peu plus rapidement que dot(M. T, M).

(peut sembler nitpicky à s'inquiéter de l'optimisation pour un 2x facteur constant sur le temps et / ou l'espace, mais il peut faire une différence à ce seuil de la taille d'un problème que vous pouvez gérer sur une seule machine...)

Il existe un certain nombre de façons bien connues de stocker des matrices symétriques afin qu'elles n'aient pas besoin d'occuper n^2 éléments de stockage. De plus, il est possible de réécrire des opérations communes pour accéder à ces moyens de stockage révisés. Le travail définitif est Golub et Van Loan, Matrix Computations , 3e édition 1996, Johns Hopkins University Press, sections 1.27-1.2.9. Par exemple, les citant de la forme (1.2.2), dans une matrice symétrique seulement besoin de stocker A = [a_{i,j} ] pouri >= j. Ensuite, en supposant que le vecteur la tenue de la matrice est notée V, et que A est n-par-n, mettre a_{i,j} dans

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

Cela suppose 1-indexation.

Golub et Van Loan offrent un algorithme 1.2.3 qui montre comment accéder à un tel V stocké pour calculer y = V x + y.

Golub et Van Loan fournissent également un moyen de stocker une matrice sous forme dominante diagonale. Cela ne sauve pas le stockage, mais prend en charge l'accès facile pour certains autres types d'opérations.

C'est Python simple et pas numpy, mais je viens de jeter ensemble une routine à remplir une matrice symétrique (et un programme de test pour s'assurer qu'il est correct):

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

Il est trivial de remplir Pythoniquement [i][j] si [j][i] est rempli. La question du stockage est un peu plus intéressante. On peut augmenter la classe numpy array avec un attribut packed qui est utile à la fois pour enregistrer le stockage et pour lire plus tard les données.

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)

```