# Efficient construction of a SparseArray from LIL (list of lists of column entries)

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.

• I have been wondering why there exist so many (CSR, CSC, LIL, DOK, etc.) ways like in Python to construct sparse matrices. Nov 21, 2020 at 9:16
– Leo
Nov 21, 2020 at 10:29
• I suspect a typo in the line SparseArray[Transpose[{uu, vv}] -> we]. we is supposed to be ww, correct? Nov 21, 2020 at 11:11
• @SjoerdSmit I have corrected the typo, thank you! I suggest we delete our comments.
– Leo
Nov 21, 2020 at 15:42
• 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... Dec 27, 2020 at 19:22

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 = CompileGetElement[p, pointer];
newlabel = CompileGetElement[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 = CompileGetElement[p, pointer];
newlabel = CompileGetElement[a, pos, 1];
boolean = Unitize[label - newlabel];
counter += boolean;
label += boolean (newlabel - label);
b[[counter]] += CompileGetElement[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

• 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?
– Leo
Dec 27, 2020 at 18:08
• 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. Dec 27, 2020 at 19:13
• Can your functions be used on rational entries? If I understand correctly, Compile only allows _Integer, _Real, _Complex, and booleans.
– Leo
Jan 31, 2021 at 17:13
• Right, exact rational numbers cannot be used in CompiledFunction`s. Jan 31, 2021 at 17:17