6
$\begingroup$

In Python scipy.sparse, there are methods to convert between CSR, CSC, LIL, DOK, etc. implementations of a sparse matrix. What is the most efficient way in Mathematica to construct a mxn SparseArray from the LIL data? (inverse of this question)

More specifically, I have a list ll={l1,...,ln}, where each lv is of the form {{u1,w1},...}, which means the matrix has an entry {u,v}->w. Note that lv may be empty (zero column). Note that lv may have repeated entries, which should be summed (solution for this is here). For testing purposes, my cases are similar to the following example (e.g. millionXmillion matrix with 10 entries per column, all from the list R):

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

My current solution is:

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}

Is there a better way? What about parallelization? How could I compile this code? (the matrix entries are integers or rationals)

P.S. Let me just mention that in Python, I haven't yet found a library for sparse matrices that allows rational number entries (exact fractions). Also, when I set every second column and every second row in a matrix to 0, the scipy.sparse implementation is waaay slower than Mathematica's SparseArray (by a factor of 100). So I'm incredibly happy we have this data structure implemented in Mathematica in such an efficient way.

$\endgroup$
5
  • 1
    $\begingroup$ I have been wondering why there exist so many (CSR, CSC, LIL, DOK, etc.) ways like in Python to construct sparse matrices. $\endgroup$ – Αλέξανδρος Ζεγγ Nov 21 '20 at 9:16
  • $\begingroup$ @ΑλέξανδροςΖεγγ Because different data types have different (dis)advantages. You can read about those in the scipy link. $\endgroup$ – Leo Nov 21 '20 at 10:29
  • 1
    $\begingroup$ I suspect a typo in the line SparseArray[Transpose[{uu, vv}] -> we]. we is supposed to be ww, correct? $\endgroup$ – Sjoerd Smit Nov 21 '20 at 11:11
  • $\begingroup$ @SjoerdSmit I have corrected the typo, thank you! I suggest we delete our comments. $\endgroup$ – Leo Nov 21 '20 at 15:42
  • 3
    $\begingroup$ As a side remark: You probably won't find an efficient sparse matrix library that supports rational numbers because rational numbers are not well covered by the native machine types (i.e., (long) integers, and single and double floats). That means that the arithmetic for rational numbers has to be emulated in software. And apparently, there is not enough demand for rational arithmetic to encourage CPU developers to design dedicated circuits for that... $\endgroup$ – Henrik Schumacher Dec 27 '20 at 19:22
5
$\begingroup$

You seem to do something wrong because the LIL you provide is more suitable to assemble the transpose of the desired matrix in CRS format (or to assemble the desired matrix in CCS format). Since Mathematica uses CRS, I show you how to assemble the transpose.

First two compiled helper functions:

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

I am not really happy with them because both tasks can actually be fused into one loop. But since CompiledFunctions cannot return more than one array and because messing around with unpacked arrays is so expensive, I leave it like this for now.

Here is the interface; CompiledFunctions don't like empty arrays as input, so I have to clean up first. Unfortunately, this has some extra cost.

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

Here is how the two methods compare:

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

True

$\endgroup$
4
  • $\begingroup$ You continue to impress, Henrik : ). I compared my and your method for m=n=10^6 and r=300 with AbsoluteTiming@MaxMemoryUsed[...], the results were {92.2022, 45630934968} and {70.6525, 16846589608} an impressive decrease of RAM consumption, thank you! Do you happen to know if it is planned in future releases of Mathematica to make compiling functions more flexible? $\endgroup$ – Leo Dec 27 '20 at 18:08
  • 2
    $\begingroup$ Thank you for the flowers! "Do you happen to know if it is planned in future releases of Mathematica to make compiling functions more flexible?" Well, they are working on a new compiler (see FunctionCompile). But the results so far have been quite underwhelming. So far, best performance can only be obtained with LibraryLink, IMHO. It is not really user-friendly, though. But Szabolcs Horvat's LTemplate package helps quite a lot in automating a lot of the process. $\endgroup$ – Henrik Schumacher Dec 27 '20 at 19:13
  • $\begingroup$ Can your functions be used on rational entries? If I understand correctly, Compile only allows _Integer, _Real, _Complex, and booleans. $\endgroup$ – Leo Jan 31 at 17:13
  • $\begingroup$ Right, exact rational numbers cannot be used in CompiledFunctions. $\endgroup$ – Henrik Schumacher Jan 31 at 17:17

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.