Numpy matrice symétrique "intelligente"

Question

Existe-t-il une matrice symétrique intelligente et peu encombrante 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 bienvenu, 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):
    """
    Retourne 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 que si a_ij = 0,
    alors le tableau retourné a' est tel que a'_ij = a_ji.

    Les valeurs diagonales restent inchangées.

    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 (telles que ne pas faire à la fois a[0, 1] = 42 et le contradictoire a[1, 0] = 123 avant d'exécuter symmetrize).

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

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

    Un SymNDArray arr est tel que faire arr[i,j] = valeur
    fait automatiquement arr[j,i] = valeur, de sorte que les mises à jour de 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_d_entree):
    """
    Retourne une version symétrisée du tableau d'entrée de type tableau-like.

    Le tableau retourné a la classe SymNDArray. Les affectations ultérieures au tableau
    sont donc automatiquement symétrisées.
    """
    return symmetrize(numpy.asarray(tableau_d_entree)).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 de tableaux, en fonction de vos besoins). Cette approche gère même des affectations plus compliquées, 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),…), donc le code doit être légèrement adapté avant de l'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 avoir examiné la question, je pense que la réponse est probablement que numpy est quelque peu contraint par la disposition 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 n^2 plutôt que n(n+1)/2. Ils sont simplement informés que la matrice est symétrique et doivent utiliser uniquement les valeurs dans le triangle supérieur ou inférieur.

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 soutien supplémentaire de la part de numpy serait appréciable, en particulier des wrappers pour les routines comme DSYRK (mise à jour symétrique du rang k), ce qui permettrait de calculer une matrice de Gram beaucoup plus rapidement que dot(M.T, M).

(Il peut sembler pointilleux de s'inquiéter de l'optimisation pour un facteur constant de 2 fois sur le temps et/ou l'espace, mais cela peut faire une différence dans 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 telle 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 revisités. L'œuvre définitive est de Golub et Van Loan, Matrix Computations, 3e édition 1996, Johns Hopkins University Press, sections 1.27-1.2.9. Par exemple, en les citant sous 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 est désigné par V, et que A est de taille n-par-n, mettre a_{i,j} dans

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

Cela suppose un indexage à partir de 1.

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

Golub et Van Loan offrent également une méthode de stockage d'une matrice sous une forme dominant la diagonale. Cela ne permet pas d'économiser du stockage, mais facilite l'accès pour certains autres types d'opérations.

Answer 4

Ceci est du python pur et non de numpy, mais j'ai juste créé une routine pour remplir une matrice symétrique (et un programme de test pour vérifier sa correction) :

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 de zéros
    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 les répliquant en diagonale
    for v in range(1,dim):
        itérations = dim - v
        x = v
        y = 0
        while itérations > 0:
            new_square[x][y] = new_square[y][x] = random.randrange(1,10)
            x += 1
            y += 1
            itérations -= 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    # afficher le numéro de la 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 garder la matrice des coûts jolie
    costMatrix = fillCostMatrix(nodeCount)
    print  "Matrice des 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 type de matrice pour la représentation de l'adjacence du réseau.

J'espère que cela vous aidera. À bientôt.