LILからのSparseArrayの効率的な構築（列エントリのリストのリスト）

Nov 21 2020

Python scipy.sparseには、スパース行列のCSR、CSC、LIL、DOKなどの実装間で変換するメソッドがあります。LILデータmxn SparseArrayからを構築するためのMathematicaで最も効率的な方法は何ですか？（この質問の逆）

より具体的には、リストll={l1,...,ln}があり、それぞれlvがの形式{{u1,w1},...}であり、マトリックスにエントリがあることを意味し{u,v}->wます。空（列がゼロ）のlv場合があることに注意してください。エントリが繰り返されている可能性があることに注意してください。これは合計する必要があります（これに対する解決策はここにあります）。テストの目的で、私のケースは次の例に似ています（たとえば、すべてリストRからの列ごとに10個のエントリを持つmillionXmillion行列）。lv

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}

もっと良い方法はありますか？並列化はどうですか？このコードをどのようにコンパイルできますか？（行列のエントリは整数または有理数です）

PS Pythonでは、有理数の入力（正確な分数）を可能にするスパース行列のライブラリをまだ見つけていません。また、行列の1つおきの列と1つおきの行を0に設定すると、scipy.sparseの実装はMathematicaのSparseArrayよりもかなり遅くなります（100倍）。ですから、このデータ構造がMathematicaにこのように効率的な方法で実装されていることを非常に嬉しく思います。

回答

5 HenrikSchumacher Dec 21 2020 at 05:16

提供するLILは、目的の行列の転置をCRS形式でアセンブルする（または目的の行列をCCS形式でアセンブルする）のに適しているため、何か問題があるようです。MathematicaはCRSを使っているので、転置を組み立てる方法を紹介します。

最初の2つのコンパイル済みヘルパー関数：

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

両方のタスクが実際に1つのループに融合される可能性があるため、私はそれらに本当に満足していません。ただし、CompiledFunctionsは複数の配列を返すことはできず、解凍された配列をいじくり回すのは非常にコストがかかるため、今のところはこのままにしておきます。

これがインターフェースです。CompiledFunctions入力として空の配列が好きではないので、最初にクリーンアップする必要があります。残念ながら、これには追加のコストがかかります。