Construction efficace d'un SparseArray à partir de LIL (liste des listes d'entrées de colonnes)

Nov 21 2020

Dans Python scipy.sparse , il existe des méthodes pour convertir entre les implémentations CSR, CSC, LIL, DOK, etc. d'une matrice creuse. Quelle est la manière la plus efficace dans Mathematica de construire un à mxn SparseArraypartir des données LIL ? (inverse de cette question)

Plus précisément, j'ai une liste ll={l1,...,ln}, où chacun lvest de la forme {{u1,w1},...}, ce qui signifie que la matrice a une entrée {u,v}->w. Notez que lvpeut être vide (colonne zéro). Notez qu'il lvpeut y avoir des entrées répétées , qui doivent être additionnées (la solution est ici ). À des fins de test, mes cas sont similaires à l'exemple suivant (par exemple, matrice millionXmillion avec 10 entrées par colonne, toutes issues de la liste R):

m=n=10^6; r=10; R={-1,1}; 
ll=Table[Transpose@{RandomInteger[{1,m},r],RandomChoice[R,r]},n];

Ma solution actuelle est:

SetSystemOptions["SparseArrayOptions"->{"TreatRepeatedEntries"->1}]; 
LIL[ll_,m_,n_] := Module[{l,uu,vv,ww}, l=Length/@ll; 
  If[Plus@@l==0,Return@SparseArray[{},{m,n}]]; 
  vv=Flatten[Table[ConstantArray[v,l[[v]]],{v,n}],1]; 
  {uu,ww}=Transpose@Flatten[ll,1];   SparseArray[Transpose[{uu,vv}]->ww] ];
AbsoluteTiming[LIL[ll,m,n];]

{5.07803, Null}

Y a-t-il un meilleur moyen? Qu'en est-il de la parallélisation? Comment pourrais-je compiler ce code? (les entrées de la matrice sont des entiers ou des rationnels)

PS Permettez-moi juste de mentionner qu'en Python, je n'ai pas encore trouvé de bibliothèque pour les matrices clairsemées permettant des entrées de nombres rationnels (fractions exactes). De plus, lorsque je mets à 0 chaque deuxième colonne et chaque deuxième ligne d'une matrice, l'implémentation scipy.sparse est waaay plus lente que SparseArray de Mathematica (d'un facteur 100). Je suis donc extrêmement heureux que cette structure de données soit implémentée dans Mathematica de manière aussi efficace.

Réponses

5 HenrikSchumacher Dec 21 2020 at 05:16

Vous semblez faire quelque chose de mal car le LIL que vous fournissez est plus adapté pour assembler la transposition de la matrice souhaitée au format CRS (ou pour assembler la matrice souhaitée au format CCS). Puisque Mathematica utilise CRS, je vous montre comment assembler la transposition.

Deux premières fonctions d'assistance compilées:

getColumnIndices = Compile[{{p, _Integer, 1}, {a, _Integer, 2}},
   Block[{b, label, newlabel, counter, pointer, n, pos, boolean},
    n = Min[Length[p], Length[a]];
    b = Table[0, {n}];
    counter = 0;
    pointer = 0;
    label = 0;
    pos = 0;
    While[pointer < n,
     pointer++;
     pos = Compile`GetElement[p, pointer];
     newlabel = Compile`GetElement[a, pos, 1];
     boolean = Unitize[label - newlabel];
     counter += boolean;
     label += boolean (newlabel - label);
     b[[counter]] = label;
     ];
    b[[1 ;; counter]]
    ],
   CompilationTarget -> "C",
   RuntimeAttributes -> {Listable},
   Parallelization -> True,
   RuntimeOptions -> "Speed"
   ];

getNonzeroValues = Compile[{{p, _Integer, 1}, {a, _Integer, 2}},
   Block[{b, label, newlabel, counter, pointer, n, pos, boolean},
    n = Min[Length[p], Length[a]];
    b = Table[0, {n}];
    counter = 0;
    pointer = 0;
    label = 0;
    pos = 0;
    While[pointer < n,
     pointer++;
     pos = Compile`GetElement[p, pointer];
     newlabel = Compile`GetElement[a, pos, 1];
     boolean = Unitize[label - newlabel];
     counter += boolean;
     label += boolean (newlabel - label);
     b[[counter]] += Compile`GetElement[a, pos, 2];
     ];
    b[[1 ;; counter]]
    ],
   CompilationTarget -> "C",
   RuntimeAttributes -> {Listable},
   Parallelization -> True,
   RuntimeOptions -> "Speed"
   ];

Je ne suis pas vraiment satisfait car les deux tâches peuvent en fait être fusionnées en une seule boucle. Mais comme CompiledFunctions ne peut pas renvoyer plus d'un tableau et parce que jouer avec des tableaux décompressés coûte si cher, je laisse les choses comme ça pour le moment.

Voici l'interface; CompiledFunctions n'aiment pas les tableaux vides en entrée, donc je dois d'abord nettoyer. Malheureusement, cela a un coût supplémentaire.