LIL (열 항목 목록)에서 SparseArray의 효율적인 구성

Nov 21 2020

Python scipy.sparse 에는 희소 행렬의 CSR, CSC, LIL, DOK 등 구현간에 변환하는 메서드가 있습니다. Mathematica mxn SparseArray에서 LIL 데이터 를 구성하는 가장 효율적인 방법은 무엇입니까 ? ( 이 질문의 반대 )

좀 더 구체적으로 말하면, 목록이 ll={l1,...,ln}있는데, 각각 lv의 형식 {{u1,w1},...}은 행렬에 항목이 있음을 의미합니다 {u,v}->w. 참고 lv있을 수 있습니다 (제로 열). 반복되는 항목lv있을 수 있으며 ,이 항목 은 합산되어야합니다 (이에 대한 해결책은 여기에 있음 ). 테스트 목적으로 내 사례는 다음 예와 유사합니다 (예 : 열당 10 개의 항목이있는 millionXmillion 행렬, 모두 목록 R에서).

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

내 현재 솔루션은 다음과 같습니다.

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}

더 좋은 방법이 있습니까? 병렬화는 어떻습니까? 이 코드를 어떻게 컴파일 할 수 있습니까? (행렬 항목은 정수 또는 합리적입니다)

추신 : 파이썬에서 유리수 항목 (정확한 분수)을 허용하는 희소 행렬을위한 라이브러리를 아직 찾지 못했다고 말씀 드리겠습니다. 또한 행렬의 모든 두 번째 열과 두 번째 행을 0으로 설정하면 scipy.sparse 구현이 Mathematica의 SparseArray (100 배)보다 느립니다. 그래서 저는이 데이터 구조가 Mathematica에서 효율적인 방식으로 구현되어 매우 기쁩니다.

답변

5 HenrikSchumacher Dec 21 2020 at 05:16

제공하는 LIL이 CRS 형식으로 원하는 행렬의 전치를 어셈블하거나 CCS 형식으로 원하는 행렬을 어셈블하는 데 더 적합하기 때문에 뭔가 잘못된 것 같습니다. Mathematica 는 CRS를 사용 하므로 조옮김을 어셈블하는 방법을 보여 드리겠습니다.

처음 두 개의 컴파일 된 도우미 함수 :

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

두 작업이 실제로 하나의 루프로 통합 될 수 있기 때문에 나는 그들에별로 만족하지 않습니다. 그러나 CompiledFunctions는 하나 이상의 배열을 반환 할 수없고 압축을 푼 배열을 엉망으로 만드는 것이 너무 비싸기 때문에 지금은 이렇게 둡니다.

다음은 인터페이스입니다. CompiledFunctions는 빈 배열을 입력으로 좋아하지 않으므로 먼저 정리해야합니다. 불행히도 이것은 약간의 추가 비용이 있습니다.

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}}
  ]

두 가지 방법을 비교하는 방법은 다음과 같습니다.

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

진실