Calculer la fermeture transitive épars de la matrice creuse scipy

Question

Je veux calculer la fermeture transitive en Python. Actuellement, j'utilise des matrices creuses de scipy.

La puissance de la matrice (`12` dans mon cas) fonctionne bien sur des matrices très clairsemées, peu importe leur taille, mais pour des cas dirigés pas très clairsemés, j'aimerais utiliser un algorithme plus intelligent.

J'ai trouvé l'_algorithme de Floyd-Warshall_ (la page allemande a un meilleur pseudo-code) dans scipy.sparse.csgraph, qui fait un peu plus que ce qu'il devrait: il n'y a pas de fonction uniquement pour l'algorithme de Warshall - c'est une chose.

Le problème principal est que je peux passer une matrice creuse à la fonction, mais cela n'a aucun sens car la fonction retournera toujours une matrice dense, car ce qui devrait être 0 dans la fermeture transitive est maintenant un chemin de longueur inf et quelqu'un a pensé que cela devait être stocké explicitement.

Ma question est donc : Y a-t-il un module python qui permet de calculer la fermeture transitive d'une matrice creuse et de la garder clairsemée?

Je ne suis pas sûr à 100% qu'il travaille avec les mêmes matrices, mais Gerald Penn montre des accélérations impressionnantes dans son article de comparaison, ce qui suggère qu'il est possible de résoudre le problème.

---

EDIT: Comme il y a eu un certain nombre de confusions, je vais expliquer le contexte théorique:

Je cherche la fermeture transitive (pas réflexive ou symétrique).

Je m'assurerai que ma relation encodée dans une matrice booléenne a les propriétés requises, c'est-à-dire la symétrie ou la réflexivité.

J'ai deux cas de relation:

1. réflexive
2. réflexive et symétrique

Je veux appliquer la fermeture transitive sur ces deux relations. Cela fonctionne parfaitement bien avec la puissance de la matrice (seulement que dans certains cas cela est trop coûteux):

>>> réflexive
matrix([[ True,  True, False,  True],
        [False,  True,  True, False],
        [False, False,  True, False],
        [False, False, False,  True]])
>>> reflexive**4
matrix([[ True,  True,  True,  True],
        [False,  True,  True, False],
        [False, False,  True, False],
        [False, False, False,  True]])
>>> reflexive_symmetric
matrix([[ True,  True, False,  True],
        [ True,  True,  True, False],
        [False,  True,  True, False],
        [ True, False, False,  True]])
>>> reflexive_symmetric**4
matrix([[ True,  True,  True,  True],
        [ True,  True,  True,  True],
        [ True,  True,  True,  True],
        [ True,  True,  True,  True]])

Ainsi, dans le premier cas, nous obtenons tous les descendants d'un nœud (y compris lui-même) et dans le deuxième cas, nous obtenons tous les composants, c'est-à-dire tous les nœuds qui sont dans le même composant.

**

Answer 1

Cela a été abordé sur le tracker de problèmes de SciPy. Le problème n'est pas tant le format de sortie; l'implémentation de Floyd-Warshall consiste à commencer avec la matrice pleine d'infinis et ensuite insérer des valeurs finies lorsqu'un chemin est trouvé. La parcimonie est immédiatement perdue.

La bibliothèque networkx offre une alternative avec sa longueur de chemin la plus courte pour toutes les paires. Sa sortie est un itérateur qui renvoie des tuples de la forme

(source, dictionnaire des cibles atteignables) 

qui prend un peu de travail pour convertir en une matrice creuse SciPy (le format csr est naturel ici). Un exemple complet:

import numpy as np
import networkx as nx
import scipy.stats as stats
import scipy.sparse as sparse

A = sparse.random(6, 6, density=0.2, format='csr', data_rvs=stats.randint(1, 2).rvs).astype(np.uint8)
G = nx.DiGraph(A)       # dirigé car A n'a pas besoin d'être symétrique
chemins = nx.all_pairs_shortest_path_length(G)
indices = []
indptr = [0]
pour ligne in chemins:
  atteignable = [v for v in ligne[1] if ligne[1][v] > 0]
  indices.extend(atteignable)
  indptr.append(len(indices))
data = np.ones((len(indices),), dtype=np.uint8)
A_trans = A + sparse.csr_matrix((data, indices, indptr), shape=A.shape)
print(A, "\n\n", A_trans)

La raison d'ajouter à nouveau A est la suivante. La sortie de Networkx comprend des chemins de longueur 0, qui rempliraient immédiatement la diagonale. Nous ne voulons pas que cela se produise (vous vouliez la fermeture transitive, pas la fermeture réflexive et transitive). D'où la ligne atteignable = [v for v in ligne[1] if ligne[1][v] > 0]. Mais alors, nous n'obtenons aucune entrée diagonale du tout, même là où A en avait (le chemin vide de longueur 0 bat le chemin de longueur 1 formé par une boucle sur soi). Alors j'ajoute A au résultat. Il a maintenant des entrées de 1 ou 2 mais seule leur non-nullité est significative.

Un exemple de l'exécution de ce qui précède (je choisis une taille de 6 par 6 pour la lisibilité de la sortie). Matrice d'origine:

  (0, 3)    1
  (3, 2)    1
  (4, 3)    1
  (5, 1)    1
  (5, 3)    1
  (5, 4)    1
  (5, 5)    1 

Fermeture transitive:

  (0, 2)    1
  (0, 3)    2
  (3, 2)    2
  (4, 2)    1
  (4, 3)    2
  (5, 1)    2
  (5, 2)    1
  (5, 3)    2
  (5, 4)    2
  (5, 5)    1

Vous pouvez voir que cela a fonctionné correctement : les entrées ajoutées sont (0, 2), (4, 2) et (5, 2), toutes acquises via le chemin (3, 2).

---

Par ailleurs, networkx a aussi la méthode floyd_warshall mais sa documentation indique

Cet algorithme est le plus approprié pour les graphes denses. Le temps d'exécution est en O(n^3), et l'espace d'exécution est en O(n^2) où n est le nombre de nœuds dans G.

La sortie est à nouveau dense. J'ai l'impression que cet algorithme est simplement considéré dense par nature. Il semble que la longueur de chemin la plus courte pour toutes les paires soit une sorte d'algorithme de Dijkstra.

Transitive et Réflexif

Si au lieu de la fermeture transitive (qui est la plus petite relation transitive contenant celle donnée) vous vouliez la fermeture transitive et réflexive (la plus petite relation transitive et réflexive contenant celle donnée), le code se simplifie car nous ne nous préoccupons plus des chemins de longueur 0.

pour ligne in chemins:
  indices.extend(ligne[1])
  indptr.append(len(indices))
data = np.ones((len(indices),), dtype=np.uint8)
A_trans = sparse.csr_matrix((data, indices, indptr), shape=A.shape)

Transitive, Réflexif et Symétrique

Cela signifie trouver la plus petite relation d'équivalence contenant celle donnée. De manière équivalente, diviser les sommets en composants connectés. Pour cela, vous n'avez pas besoin de passer par networkx, il y a la méthode connected_components de SciPy. Définissez directed=False là-bas. Exemple:

import numpy as np
import scipy.stats as stats
import scipy.sparse as sparse
import itertools

A = sparse.random(20, 20, density=0.02, format='csr', data_rvs=stats.randint(1, 2).rvs).astype(np.uint8)
composants = sparse.csgraph.connected_components(A, directed=False)
non_nuls = []
for k in range(composants[0]):
  idx = np.where(composants[1] == k)[0]
  non_nuls.extend(itertools.product(idx, idx))
  ligne = tuple(l for l, c in non_nuls)
  colonne = tuple(c for l, c in non_nuls)
  data = np.ones_like(ligne)
B = sparse.coo_matrix((data, (ligne, colonne)), shape=A.shape)

Voici à quoi ressemble la sortie de print(B.toarray()) pour un exemple aléatoire, 20 par 20 :

[[1 0 1 0 0 0 1 0 0 0 0 0 1 0 0 0 1 0 0 0]
 [0 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0]
 [1 0 1 0 0 0 1 0 0 0 0 0 1 0 0 0 1 0 0 0]
 [0 0 0 1 0 0 0 0 0 0 0 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 0 0 0]
 [0 0 0 0 0 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0]
 [1 0 1 0 0 0 1 0 0 0 0 0 1 0 0 0 1 0 0 0]
 [0 0 0 0 0 0 0 1 0 0