Permutez une matrice en place dans numpy

27

Je veux modifier une matrice de transition carrée dense en place en changeant l'ordre de plusieurs de ses lignes et colonnes, en utilisant la bibliothèque numpy de python. Mathématiquement, cela correspond à la pré-multiplication de la matrice par la matrice de permutation P et à la post-multiplication par P ^ -1 = P ^ T, mais ce n'est pas une solution calculable.

En ce moment, j'échange manuellement des lignes et des colonnes, mais je m'attendais à ce que numpy ait une belle fonction f (M, v) où M a n lignes et colonnes, et v a n entrées, de sorte que f (M, v) se mette à jour M selon la permutation d'index v. Peut-être que je ne réussis pas à chercher sur Internet.

Quelque chose comme cela pourrait être possible avec "l'indexation avancée" de numpy mais je crois comprendre qu'une telle solution ne serait pas en place. De plus, pour certaines situations simples, il peut être suffisant de suivre séparément une permutation d'index, mais ce n'est pas pratique dans mon cas.

Ajouté:
Parfois, lorsque les gens parlent de permutations, ils ne signifient que l'échantillonnage de permutations aléatoires, par exemple dans le cadre d'une procédure pour obtenir des valeurs de p dans les statistiques. Ou ils signifient compter ou énumérer toutes les permutations possibles. Je ne parle pas de ces choses.

Ajouté:
La matrice est assez petite pour tenir dans la RAM du bureau mais assez grande pour que je ne veuille pas la copier sans réfléchir. En fait, je voudrais utiliser des matrices aussi grandes que possible, mais je ne veux pas gérer l'inconvénient de ne pas pouvoir les conserver dans la RAM, et je fais des opérations O (N ^ 3) LAPACK sur la matrice qui serait également limiter la taille de la matrice pratique. Je copie actuellement des matrices de cette taille inutilement, mais j'espère que cela pourrait être facilement évité pour la permutation.

aucun
la source
3
Ce serait bien si vous pouviez mettre à jour la question pour donner la taille de vos matrices. "Gigantesque" ne signifie pas la même chose pour tout le monde.
Bill Barth
2
Vous avez raison que l'indexation avancée (ou soi-disant fantaisie) crée une copie. Mais si vous acceptez de vivre avec ce fait, votre code consiste simplement M[v]à permuter les lignes.
Daniel Velkov
@daniel: Et ce serait M [v,:] [:, v] de faire toute la permutation? Serait-ce le meilleur moyen d'obtenir la permutation en utilisant une indexation sophistiquée? Et utiliserait-il 3 fois la mémoire de la matrice, y compris la taille de la matrice d'origine, la matrice permutée ligne + colonne et la matrice permutée ligne temporaire?
aucun
C'est exact, vous auriez votre matrice originale et 2 copies. Btw pourquoi avez-vous besoin de permuter les lignes et les colonnes en même temps?
Daniel Velkov
4
Qu'allez-vous faire avec la matrice permutée? Il peut être préférable de simplement permuter le vecteur lors de l'application de l'opérateur.
Jed Brown

Réponses:

9

Selon les documents, il n'y a pas de méthode de permutation sur place dans numpy, quelque chose comme ndarray.sort .

Donc, vos options sont (en supposant que Mc'est une matrice et le vecteur de permutation)N×Np

  1. implémenter votre propre algorithme en C comme module d'extension (mais les algorithmes sur place sont difficiles, du moins pour moi!)
  2. N mémoire

    for i in range(N):
        M[:,i] = M[p,i]
    for i in range(N):
        M[i,:] = M[i,p]
    
  3. N2 mémoire

    M[:,:] = M[p,:]
    M[:,:] = M[:,p]
    

J'espère que ces hacks sous-optimaux sont utiles.

Stefano M
la source
@none is hack 2. ce que vous appelez «permuter manuellement des lignes et des colonnes»?
Stefano M
1
Je combinerais les options 1 et 2: écrire du code C qui utilise un tampon d'ordre N pour écrire chaque colonne permutée, puis la réécrire d'où elle vient; faites de même pour les lignes. Comme l'écrit @Stefano, cela ne prend que de mémoire supplémentaire, que vous dépensez déjà pour stocker la permutation en premier lieu. pO(N)p
Erik P.
@ErikP. pour une implémentation C, la mémoire supplémentaire est raisonnable et votre approche d'écriture en dispersion sur temp et de copie en retour est sûre. La question intéressante est cependant de savoir s'il existe des algorithmes plus efficaces, étant donné de mémoire supplémentaire. La réponse est difficile, je pense, car nous devons prendre en compte l'architecture du processeur, les modèles d'accès à la mémoire, les accès au cache, ... Cela dit, je suivrais vos conseils et opterais pour un algorithme simple et facile à implémenter. O ( N )O(N)O(N)
Stefano M
2
C'est un très bon canidate pour une fonction cython. Ne doit pas dépasser 10 lignes. . . voulez-vous que je lui donne une fissure?
meawoppl
Lol. J'ai commencé à Cython ceci, puis j'ai trouvé la bonne réponse dans une fonction que j'utilise tout le temps. Doh. Voir ma réponse publiée.
meawoppl
6

Avertissement: L'exemple ci-dessous fonctionne correctement, mais l' utilisation de l'ensemble complet des paramètres suggérés à la fin du post expose un bogue , ou au moins une "fonctionnalité non documentée" dans la fonction numpy.take (). Voir les commentaires ci-dessous pour plus de détails. Rapport de bogue déposé .

Vous pouvez le faire sur place avec la fonction take () de numpy , mais cela nécessite un peu de saut de cerceau.

Voici un exemple de réalisation d'une permutation aléatoire des lignes d'une matrice d'identité:

import numpy as np
i = np.identity(10)
rr = range(10)
np.random.shuffle(rr)
np.take(i, rr, axis=0)
array([[ 0.,  1.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.],
       [ 0.,  0.,  0.,  0.,  0.,  1.,  0.,  0.,  0.,  0.],
       [ 0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  1.,  0.],
       [ 0.,  0.,  0.,  0.,  0.,  0.,  1.,  0.,  0.,  0.],
       [ 0.,  0.,  0.,  1.,  0.,  0.,  0.,  0.,  0.,  0.],
       [ 0.,  0.,  0.,  0.,  1.,  0.,  0.,  0.,  0.,  0.],
       [ 0.,  0.,  0.,  0.,  0.,  0.,  0.,  1.,  0.,  0.],
       [ 1.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.],
       [ 0.,  0.,  1.,  0.,  0.,  0.,  0.,  0.,  0.,  0.],
       [ 0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  1.]])

Pour le faire sur place, tout ce que vous devez faire est de spécifier que le paramètre "out" doit être le même que le tableau d'entrée ET vous devez définir le mode = "clip" ou mode = "wrap". Si vous ne définissez pas le mode, il fera une copie pour restaurer l'état du tableau sur une exception Python (voir ici) .

Sur une note finale, take semble être une méthode de tableau, donc au lieu de

np.take(i, rr, axis=0)

tu pourrais appeler

i.take(rr, axis=0)

si c'est plus à votre goût. Donc, au total, vous appelez devrait ressembler à ceci:

#Inplace Rearrange
arr = makeMyBixMatrix()
pVec0, pVec1 = calcMyPermutationVectors()
arr.take(pVec0, axis=0, out=arr, mode="clip")
arr.take(pVec1, axis=1, out=arr, mode="clip")

Pour permuter les lignes et les colonnes, je pense que vous devez soit l'exécuter deux fois, soit tirer de vilains manigances avec numpy.unravel_index qui me font penser à la tête.

meawoppl
la source
Comme dit, les algorithmes en place sont difficiles. Votre solution ne fonctionne pas avec numpy 1.6.2. et 1.7.1 (lignes / colonnes en double). N'a pas eu le temps de vérifier si la version 1.8.x résout ce problème
Stefano M
Hmmm. Pouvez-vous publier du code de test quelque part? Dans ma tête, j'ai l'impression qu'il doit y avoir une opération de tri sur les indices qui se produit d'abord avant le plumage. Je vais enquêter davantage sur ce PM.
meawoppl
1
Quand je lance ce code que je reçois 1.6.2, test take, not overwriting: True, test not-in-place take: True, test in-place take: False, rr [3, 7, 8, 1, 4, 5, 9, 0, 2, 6], arr [30 70 80 70 40 50 90 30 80 90], ref [30 70 80 10 40 50 90 0 20 60]. Donc np.takeau moins pour numpy 1.6.2 n'est pas conscient de faire une permutation sur place et gâche les choses.
Stefano M
Yeouch. Bien démontré. Cela se qualifie probablement comme un bug à mon humble avis. Au minimum, les documents devraient dire que l'entrée et la sortie ne peuvent pas être le même tableau, vérifiez probablement pour voir et sauf si c'est le cas.
meawoppl
D'accord sur le bogue: vous devriez peut-être ajouter une note à votre message pour avertir les lecteurs que votre solution peut produire de mauvais résultats.
Stefano M
2

Si vous avez une matrice clairsemée stockée au COOformat, les éléments suivants peuvent être utiles

    A.row = perm[A.row];
    A.col = perm[A.col];

en supposant que Acontient la COOmatrice, et permest un numpy.arraycontenant la permutation. Il n'y aura que surcharge mémoire, où est le nombre d'éléments non nuls de la matrice.mmm

Vincent Traag
la source
mais quelle est la surcharge de mémoire pour stocker une matrice dense complète comme C00matrice clairsemée en premier lieu?
Federico Poloni
Comme le nombre d'éléments est égal à la fois dans une représentation dense et dans une représentation dense (complète), la différence de mémoire est simplement une constante (2 ints et 1 floatdans une représentation clairsemée par élément par rapport à un seul floatdans la représentation dense). Mais la surcharge de mémoire de cette méthode seran2 dans un cas dense, donc on pourrait probablement mieux s'en tenir aux numpy.ndarrays réguliers .
Vincent Traag
1

Je n'ai pas assez de réputation pour commenter, mais je pense que la question SO suivante pourrait être utile: /programming/4370745/view-onto-a-numpy-array

Les points de base sont que vous pouvez utiliser le découpage de base et créerez une vue sur le tableau sans copier, mais si vous le faites tranchage / indexation avancée , alors il va créer une copie.

hadsed
la source
L'OP demande une permutation, ce qui n'est pas possible avec un découpage de base.
Stefano M
Vous avez bien sûr raison. J'ai pensé qu'il serait utile pour le PO de comprendre ce qui se passait avec le découpage (au cas où ils ne le sauraient pas) car ils se demandaient quand les copies auraient lieu. S'il a utilisé quelque chose de votre réponse, je pense que ce serait bien de le savoir puisque vous les utilisez dans vos boucles.
hadsed
-1

Qu'en est-il de

my_array [:, [0, 1]] = my_array [:, [1, 0]]

johnsankey
la source
1
Cela construit un temporaire, c'est exactement ce qu'il veut éviter.
Michael Grant