Effiziente Erstellung eines SparseArray aus LIL (Liste der Listen der Spalteneinträge)

Nov 21 2020

In Python scipy.sparse gibt es Methoden zum Konvertieren zwischen CSR- , CSC- , LIL- , DOK- usw. Implementierungen einer dünn besetzten Matrix. Was ist der effizienteste Weg in Mathematica, um mxn SparseArrayaus den LIL- Daten eine zu erstellen? (Umkehrung dieser Frage)

Genauer gesagt habe ich eine Liste ll={l1,...,ln}, in der jede lvForm vorliegt {{u1,w1},...}, was bedeutet, dass die Matrix einen Eintrag hat {u,v}->w. Beachten Sie, dass dies lvmöglicherweise leer ist (Spalte Null). Beachten Sie, dass lvmöglicherweise wiederholte Einträge vorhanden sind , die summiert werden sollten (Lösung hierfür finden Sie hier ). Zu Testzwecken ähneln meine Fälle dem folgenden Beispiel (z. B. Millionen-Millionen-Matrix mit 10 Einträgen pro Spalte, alle aus der Liste R):

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

Meine aktuelle Lösung ist:

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}

Gibt es einen besseren Weg? Was ist mit Parallelisierung? Wie könnte ich diesen Code kompilieren? (Die Matrixeinträge sind ganze Zahlen oder Rationalen.)

PS Lassen Sie mich nur erwähnen, dass ich in Python noch keine Bibliothek für spärliche Matrizen gefunden habe, die rationale Zahleneinträge (exakte Brüche) ermöglicht. Wenn ich jede zweite Spalte und jede zweite Zeile in einer Matrix auf 0 setze, ist die Implementierung von scipy.sparse möglicherweise langsamer als das SparseArray von Mathematica (um den Faktor 100). Ich bin unglaublich froh, dass wir diese Datenstruktur so effizient in Mathematica implementiert haben.

Antworten

5 HenrikSchumacher Dec 21 2020 at 05:16

Sie scheinen etwas falsch zu machen, da die von Ihnen bereitgestellte LIL besser geeignet ist, die Transponierte der gewünschten Matrix im CRS-Format (oder die gewünschte Matrix im CCS-Format) zusammenzusetzen. Da Mathematica CRS verwendet, zeige ich Ihnen, wie Sie die Transponierte zusammensetzen.

Die ersten beiden kompilierten Hilfsfunktionen:

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"
   ];

Ich bin nicht wirklich glücklich mit ihnen, weil beide Aufgaben tatsächlich zu einer Schleife verschmolzen werden können. Aber da CompiledFunctions nicht mehr als ein Array zurückgeben kann und das Herumspielen mit entpackten Arrays so teuer ist, lasse ich es vorerst so.

Hier ist die Schnittstelle; CompiledFunctions mögen keine leeren Arrays als Eingabe, deshalb muss ich zuerst aufräumen. Dies hat leider einige zusätzliche Kosten.

LIL2[ll_, m_, n_] := Module[{idx, llclean, orderings, vals, rp, ci},
  idx = Pick[Range[Length[ll]], Unitize[Length /@ ll], 1];
  llclean = ll[[idx]];
  rp = ConstantArray[0, Length[ll] + 1];
  orderings = Ordering /@ llclean;
  vals = Join @@ getNonzeroValues[orderings, llclean];
  With[{data = getColumnIndices[orderings, llclean]},
   ci = Partition[Join @@ data, 1];
   rp[[idx + 1]] = Length /@ data;
   ];
  rp = Accumulate[rp];
  SparseArray @@ {Automatic, {n, m}, 0, {1, {rp, ci}, vals}}
  ]

So vergleichen sich die beiden Methoden:

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

A = LIL[ll, m, n]; // AbsoluteTiming // First
B = LIL2[ll, m, n]; // AbsoluteTiming // First
A == Transpose[B]

4,02563

1,81523

Wahr