Matrice symétrique "intelligente" Numpy

Question

Existe-t-il une matrice symétrique intelligente et compacte dans numpy qui remplit automatiquement (et de manière transparente) la position à [j][i] lorsque [i][j] est écrit ?

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) # pour toute matrice symétrique

Un Hermitien automatique serait également agréable, bien que je n'en aurai pas besoin au moment de l'écriture.

Answer 1

Si vous pouvez vous permettre de symétriser la matrice juste avant de faire des calculs, ce qui suit devrait être raisonnablement rapide :

def symmetrize(a):
    """
    Renvoie une version symétrisée du tableau NumPy a.

    Les valeurs 0 sont remplacées par la valeur du tableau à la position symétrique
    (par rapport à la diagonale), c'est-à-dire si a_ij = 0,
    alors le tableau renvoyé a' est tel que a'_ij = a_ji.

    Les valeurs diagonales sont laissées intactes.

    a -- tableau NumPy carré, tel que a_ij = 0 ou a_ji = 0, 
    pour i != j.
    """
    return a + a.T - numpy.diag(a.diagonal())

Cela fonctionne sous des hypothèses raisonnables (comme ne pas faire à la fois a[0, 1] = 42 et le contraire a[1, 0] = 123 avant d'exécuter symmetrize).

Si vous avez vraiment besoin d'une symétrisation transparente, vous pouvez envisager de sous-classer numpy.ndarray et simplement redéfinir __setitem__ :

class SymNDArray(numpy.ndarray):
    """
    Sous-classe de tableau NumPy pour les matrices symétriques.

    Un tableau SymNDArray arr est tel que lorsque arr[i, j] = valeur
    fait automatiquement arr[j, i] = valeur, de sorte que les mises à jour
    du tableau restent symétriques.
    """

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

def symarray(tableau):
    """
    Renvoie une version symétrisée du tableau-like input_array.

    Le tableau renvoyé a la classe SymNDArray. Les futures affectations au tableau
    sont donc automatiquement symétrisées.
    """
    return symmetrize(numpy.asarray(tableau)).view(SymNDArray)

# Exemple :
a = symarray(numpy.zeros((3, 3)))
a[0, 1] = 42
print a  # a[1, 0] == 42 aussi !

(ou l'équivalent avec des matrices au lieu d'arrays, selon vos besoins). Cette approche gère même des affectations plus complexes, comme a[:, 1] = -1, qui définit correctement les éléments de a[1, :].

Notez que Python 3 a supprimé la possibilité d'écrire def …(…, (i, j),…), alors le code doit être légèrement adapté avant de s'exécuter avec Python 3 : def __setitem__(self, indexes, valeur): (i, j) = indexes…

Answer 2

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

Après y avoir réfléchi, je pense que la réponse est probablement que numpy est quelque peu limité par la disposition en mémoire supportée 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 les matrices symétriques, elles utilisent toujours la même structure mémoire qu'une matrice complète, c'est-à-dire un espace de n^2 au lieu de n(n+1)/2. Elles utilisent simplement les valeurs situées dans le triangle supérieur ou inférieur de la matrice.

Certaines des routines de scipy.linalg acceptent des indicateurs (comme sym_pos=True sur linalg.solve) qui sont transmis aux routines BLAS, bien qu'un support plus important pour cela dans numpy soit souhaitable, en particulier des wrappers pour des routines comme DSYRK (mise à jour de rang k symétrique), qui permettrait de calculer une matrice de Gram assez rapidement par rapport à dot(M.T, M).

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

Answer 3

Il existe plusieurs façons bien connues de stocker des matrices symétriques de sorte 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 courantes pour accéder à ces moyens de stockage révisés. L'œuvre définitive est Golub et Van Loan, Matrix Computations, 3ème édition 1996, Johns Hopkins University Press, sections 1.27-1.2.9. Par exemple, en les citant de la forme (1.2.2), dans une matrice symétrique, il suffit de stocker A = [a_{i,j} ] pour i >= j. Ensuite, en supposant que le vector contenant la matrice soit désigné V, et que A soit de taille n par n, placez a_{i,j} dans

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

Ceci suppose un indexage à partir de 1.

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

Golub et Van Loan offrent également une manière de stocker une matrice sous forme de diagonale dominante. Cela ne permet pas d'économiser de l'espace de stockage, mais facilite l'accès à certaines autres opérations.

Answer 4

Il s'agit simplement de Python pur et non de numpy, mais j'ai simplement assemblé une routine pour remplir une matrice symétrique (et un programme de test pour m'assurer que tout est correct) :

import random

# remplir une matrice symétrique avec des coûts (c'est-à-dire m[x][y] == m[y][x]
# À des fins de démonstration, cette routine connecte chaque nœud à tous les autres
# Comme une matrice stocke les coûts, des nombres sont utilisés pour représenter les nœuds
# de sorte que les indices de ligne et de colonne peuvent représenter des nœuds

def fillCostMatrix(dim):        # tableau carré de tableaux
    # Créer une matrice nulle
    new_square = [[0 for row in range(dim)] for col in range(dim)]
    # Remplir la diagonale principale
    for v in range(0,dim):
        new_square[v][v] = random.randrange(1,10)

    # Remplir les triangles supérieur et inférieur de manière symétrique en répliquant en diagonale
    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

# test de cohérence
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 = 'PAS'
    print "La matrice est", isSymmetric, "symétrique"

def showSquare(square):
    # Afficher la matrice carrée
    columnHeader = ' '
    for i in range(len(square)):
        columnHeader += '  ' + str(i)
    print columnHeader

    i = 0;
    for col in square:
        print i, col    # imprimer le numéro de ligne et les données
        i += 1

def myMain(argv):
    if len(argv) == 1:
        nodeCount = 6
    else:
        try:
            nodeCount = int(argv[1])
        except:
            print  "l'argument doit être numérique"
            quit()

    # garder nodeCount <= 9 pour que la matrice de coûts soit jolie
    costMatrix = fillCostMatrix(nodeCount)
    print  "Matrice de coûts"
    showSquare(costMatrix)
    test_symmetry(costMatrix)   # test de cohérence
if __name__ == "__main__":
    import sys
    myMain(sys.argv)

# vim:tabstop=8:shiftwidth=4:expandtab

Answer 5

Pour construire une matrice NxN qui est symétrique le long de la diagonale principale, et avec des 0 sur la diagonale principale, vous pouvez faire :

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

C'est un cas un peu spécial, mais récemment j'ai utilisé ce genre de matrice pour représenter l'adjacence d'un réseau.

J'espère que cela aide. Santé.