Efektywna konstrukcja SparseArray z LIL (lista list wpisów w kolumnach)

Nov 21 2020

W Pythonie scipy.sparse istnieją metody konwersji między CSR, CSC, LIL, DOK itp. Implementacjami rzadkiej macierzy. Jaki jest najbardziej efektywny sposób w Mathematica na skonstruowanie a mxn SparseArrayna podstawie danych LIL ? (odwrotność tego pytania)

Dokładniej, mam listę ll={l1,...,ln}, gdzie każda lvma postać {{u1,w1},...}, co oznacza, że matryca ma wpis {u,v}->w. Zauważ, że lvmoże być puste (kolumna zerowa). Zauważ, że lvmogą mieć powtarzające się wpisy , które należy zsumować (rozwiązanie tego problemu znajduje się tutaj ). Dla celów testowych moje przypadki są podobne do następującego przykładu (np. Macierz milion X milionów z 10 wpisami na kolumnę, wszystkie z listy R):

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

Moje obecne rozwiązanie to:

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}

Czy jest lepszy sposób? A co z równoległością? Jak mogę skompilować ten kod? (wpisy macierzy są liczbami całkowitymi lub wymiernymi)

PS Wspomnę tylko, że w Pythonie nie znalazłem jeszcze biblioteki dla rzadkich macierzy, która umożliwia wprowadzanie liczb wymiernych (dokładne ułamki). Ponadto, gdy ustawię co drugą kolumnę i co drugi wiersz macierzy na 0, implementacja scipy.sparse jest znacznie wolniejsza niż SparseArray Mathematica (o współczynnik 100). Jestem więc niesamowicie szczęśliwy, że mamy taką strukturę danych zaimplementowaną w Mathematica w tak efektywny sposób.

Odpowiedzi

5 HenrikSchumacher Dec 21 2020 at 05:16

Wydaje się, że robisz coś złego, ponieważ dostarczona przez Ciebie LIL jest bardziej odpowiednia do złożenia transpozycji żądanej macierzy w formacie CRS (lub złożenia żądanej macierzy w formacie CCS). Ponieważ Mathematica używa CRS, pokażę ci, jak złożyć transpozycję.

Pierwsze dwie skompilowane funkcje pomocnicze:

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

Nie jestem z nich zbyt zadowolony, ponieważ oba zadania można właściwie połączyć w jedną pętlę. Ale ponieważ CompiledFunctions nie może zwrócić więcej niż jednej tablicy, a majstrowanie przy rozpakowanych tablicach jest tak kosztowne, na razie zostawiam to w ten sposób.

Oto interfejs; CompiledFunctions nie lubią pustych tablic jako danych wejściowych, więc najpierw muszę je wyczyścić. Niestety wiąże się to z dodatkowymi kosztami.