Effiziente Erstellung eines SparseArray aus LIL (Liste der Listen der Spalteneinträge)
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
SparseArray
aus den LIL- Daten eine zu erstellen? (Umkehrung dieser Frage)
Genauer gesagt habe ich eine Liste ll={l1,...,ln}
, in der jede lv
Form vorliegt {{u1,w1},...}
, was bedeutet, dass die Matrix einen Eintrag hat {u,v}->w
. Beachten Sie, dass dies lv
möglicherweise leer ist (Spalte Null). Beachten Sie, dass lv
mö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
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 CompiledFunction
s 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; CompiledFunction
s 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