9

Suppose we have a very big very sparse matrix, e.g.:

m=10^5; n=10^5; r=10^6; a =SparseArray[Transpose[{RandomInteger[{1,m},r],RandomInteger[{1,n},r]}]->RandomChoice[{-1,1,2,-3},r],{m,n}]

How can I efficiently transform a SparseArray (which is in the CRS format) into the list of lists format, i.e. a table, which contains in the i-th place all entries of the i-th column? One solution is:

Table[{#[[1,1]],#[[2]]}& /@ ArrayRules[a[[All,v]]][[;;-2]],{v,Dimensions[a][[2]]}]

This is way too slow. Another solution (which I copied from here and don't understand), is: 

Module[{b=Transpose[a],ci,rp,v,l1,l2}, ci=b["ColumnIndices"]; rp=b["RowPointers"]; v=b["NonzeroValues"]; 
 l1 = Internal`PartitionRagged[Flatten[ci], Differences[rp]];   
 l2 = Internal`PartitionRagged[v, Differences[rp]];   
 Table[ Transpose[{l1[[j]],l2[[j]]}], {j,Length@l1}]]

This is slow as well. A third attempt:

Module[{l=GatherBy[Sort@Transpose@Join[Reverse@Transpose@a["NonzeroPositions"],{a["NonzeroValues"]}],First]},   
  l=AssociationThread[Map[First,l,2],Map[Rest,l,{2}]];  Table[Lookup[l,j,{}],{j,Dimensions[a][[2]]}]];

This is very fast, but it constructs an Association, which eats up a lot of RAM. Is there a better way?

Leo
  • 1,155
  • 5
  • 14
  • Why do you need the ListOfLists format? Is it for export? I ask because perhaps there might be more efficient methods than converting to an explicit list. – MarcoB May 24 '20 at 22:39
  • 1
    @MarcoB Not for export, I need LOL in an efficient-to-be algorithm, in which I search for paths in a bipartite acyclic directed graph with weighted edges. The mentioned LOL represents the edges of the graph. – Leo May 24 '20 at 23:48
  • 1
    I modified the third one to avoid associations: Block[{pos, id, l, nu}, pos = a["NonzeroPositions"]; l = GatherBy[Sort[MapThread[Append, {Reverse[pos, 2], a["NonzeroValues"]}]], First]; nu = Transpose[{Complement[Range[Last[Dimensions[a]]], a["NonzeroPositions"][[All, 2]]]}]; Insert[Drop[#, None, 1] & /@ l, {}, nu]]. Please test it against your actual application. – J. M.'s missing motivation May 25 '20 at 02:47
  • @J.M. Thank you! For a $10^6\times10^6$ matrix with $10^7$ nonzero entries, it needs 20seconds and 3GB RAM, an improvement! – Leo May 25 '20 at 08:59

1 Answers1

11

Try

result = With[{b = Transpose[a]},
     Internal`PartitionRagged[
      Transpose[{Flatten[b["ColumnIndices"]], b["NonzeroValues"]}],
      Differences[b["RowPointers"]]
      ]
     ]; // AbsoluteTiming // First

0.074976

This is only difference from the second code snippet that you posted: Internal`PartitionRagged is performed only once and in the very end. This allows us to skip Table (slow!) and to reduce the amount of operations on ragged lists.

How does this work? First of all, LOL is a column-major format while CRS is row-major. That is why we transpose in the beginning to obtain a matrix b. The CRS is such that b["ColumnIndices"] stores all the $j$-indices from the array rules {i,j} -> v of b in row major order in a single flat list. (That corresponds to the i-indices of a.) The according nonzero values are stored in b["NonzeroValues"] in the same order. So Transpose[{Flatten[b["ColumnIndices"]], b["NonzeroValues"]}] contains the LOL data, but flattened out by one level. The row pointers in b["RowPointers"] tell one where to break this list apart to obtain the lists for each row: The $k$-th row of b has Differences[b["RowPointers"]][[k]] entries. And the actual chopping can be performed with Internal`PartitionRagged; it is a specialized and thus faster variant of TakeList. (Actually, it is vice versa: TakeList is a generalized but slower version of Internal`PartitionRagged. The former was introduced in version 11.2 while the latter has been around in the system for quite some time before.)

Henrik Schumacher
  • 106,770
  • 7
  • 179
  • 309