2

I have an enormous $m\!\times\!n$ sparse matrix $d$ with $l$ nonzero entries, stored as d=SparseArray[{{i1,j1}->w1,...,{il,jl}->wl},{m,n}];.

  1. How do I access the nonzero entries in the $i$-th row or $j$-th column? If I use d[[i,1;;-1]] or d[[1;;-1,j]], I get all entries in the $i$-th row or $j$-th column.

  2. How does 1. work? Do they search through all $l$ entries? Does $d$ have pointers to all $i$-row or $j$-column entries? What is the time complexity of these operations?

  3. Is there a better data structure, to store all the information of $d$, such that I can do 1. very fast, and that the memory requirements are not much larger than with SparseArray?

I'm a beginner in Mathematica, so simple (but effective) answers are welcome :).

EDIT: basically what I need are functions or data structures row and col (with low CPU and low RAM) that return the first and last entry (=position and value) in each row or column.

I will be calling this function many times, so I'm wondering if it is better to construct fast functions, or to store this data in 2 new lists. This depends how SparseArray works...

Leo
  • 1,155
  • 5
  • 14

1 Answers1

3

$i$-th row:

d[[i]]["NonzeroValues"]

$j$-th column:

d[[All,j]]["NonzeroValues"]

or (if you want to access many colums):

dt = Transpose[d]; 
dt[[j]]["NonzeroValues"]

Here a concrete example

n = 1000000;
m = 6000000;
d = SparseArray[
   Sort@RandomInteger[{1, n}, {m, 2}] -> RandomReal[{-1, 1}, m],
   {n, n}
   ];
k = 100000;
ilist = RandomInteger[{1, Length[d]}, k];
jlist = RandomInteger[{1, Dimensions[d][[2]]}, k];

nonzerorowvalues = Table[d[[i]]["NonzeroValues"], {i, ilist}]; // 
  RepeatedTiming // First

dt = Transpose[d];
nonzerocolvalues = Table[dt[[j]]["NonzeroValues"], {j, jlist}]; // 
  RepeatedTiming // First

0.265

0.28

Edit

Here is a significantly faster method to get the values and column indices of each row:

nonzerorowvalues3 = Internal`PartitionRagged[vals, Differences[rp]]; //
   RepeatedTiming // First

nonzerorowpositions3 = Internal`PartitionRagged[Flatten[ci], Differences[rp]]; // RepeatedTiming // First

0.248

0.26

A compiled function that uses nonzero values and row pointers to read the rows:

getRowValues = 
  Compile[{{vals, _Real, 1}, {rp, _Integer, 1}, {i, _Integer}},
   Block[{a, b},
    a = Compile`GetElement[rp, i];
    b = Compile`GetElement[rp, i + 1];
    If[a < b,
     Table[Compile`GetElement[vals, j], {j, a + 1, b}],
     {}
     ]
    ],
   CompilationTarget -> "C",
   RuntimeAttributes -> {Listable},
   Parallelization -> True,
   RuntimeOptions -> "Speed"
   ];

And the analogon for the column inidces:

getRowPositions = 
  Compile[{{ci, _Integer, 2}, {rp, _Integer, 1}, {i, _Integer}},
   Block[{a, b},
    a = Compile`GetElement[rp, i];
    b = Compile`GetElement[rp, i + 1];
    If[a < b,
     Table[Compile`GetElement[ci, j, 1], {j, a + 1, b}],
     Most[{0}]
     ]
    ],
   CompilationTarget -> "C",
   RuntimeAttributes -> {Listable},
   Parallelization -> True,
   RuntimeOptions -> "Speed"
   ];

This is how this gets applied:

ci = d["ColumnIndices"];
rp = d["RowPointers"];
vals = d["NonzeroValues"];

nonzerorowvalues = Table[d[[i]]["NonzeroValues"], {i, 1, Length[d]}]; // RepeatedTiming // First

nonzerorowvalues2 = getRowValues[vals, rp, Range[1, Length[d]]]; // 
  RepeatedTiming // First

nonzerorowvalues == nonzerorowvalues2

2.26

0.33

True

You can obtain the column indices of nonzero entries of rows by

nonzerorowpositions2 = getRowPositions[ci, rp, Range[1, Length[d]]]; // RepeatedTiming // First

Note the sparse matrices are stored in the CRS format that is optimized for row access in order to get maximum performance for matrix-vector multiplication. So the fasted method to get the columns would be to perform these tasks on the transposed matrix. Since you need only to store the nonzerovalues of and the rowpointers of the transposed matrix, the additional memory requirement will be roughly half the memory that you need to store d. So, if you first compute the nonzero values and row pointers of d and Transpose[d] and then erase d and dt, you won't need significantly more than the storage space for d.

Edit 2:

An even faster method to get the values and column indices of each row:

nonzerorowvalues3 = Internal`PartitionRagged[vals, Differences[rp]]; // RepeatedTiming // First

nonzerorowpositions3 = Internal`PartitionRagged[Flatten[ci], Differences[rp]]; // RepeatedTiming // First

0.248

0.26

See also this post for understanding what is going on.

Henrik Schumacher
  • 106,770
  • 7
  • 179
  • 309
  • Hmm, I will be accessing the nonzero entries in the $i$th row and $j$th column for every $i$ and $j$, so Transpose is not practical (going through all entries of $d$ every time). See the edited opening post. – Leo Jul 21 '18 at 13:43
  • That's why I would suggest to transpose only once. Accessing rows is significantly faster than accessing columns. – Henrik Schumacher Jul 21 '18 at 13:46
  • So that's like having another copy of $d$ in my RAM, right? In practice, $d$ will take up a large amount of memory, so having additional copies seems wasteful. Do you know which internal pointers SparseArray contains? If you want the $i$-th row, does it go through all entries, or directly to the ones of the form {i,...}->...? – Leo Jul 21 '18 at 13:58
  • I am curious: What are you going to do afterwards you have retrieved the nonzero values of the rows? – Henrik Schumacher Jul 21 '18 at 14:06
  • I'm writing a new multicore algorithm for computing the homology of a chain complex. On $d$, I must be able to find entries $d_{i,j}$, such that left and below it are only zeros (i.e. $d_{i',j}=0=d_{i,j'}$ for $i'!<!i$ and $j'!>!j$); collect these entries into $M$. Then I must be able to find all columns $i$ / rows $j$, that contain no entry from $M$. For each such $i$, I will be computing some recursive functions that require $d$ and $M$. I'm not sure whether RAM or CPU will be a bigger problem... – Leo Jul 21 '18 at 14:51
  • "I'm not sure whether RAM or CPU will be a bigger problem..." Maybe you should figure that out before bothering people with that. What is the source of your chain complexes? How size are the sparse matrices? How dense are they? What is the average number of nonzeros per row/per column? Please don't let us guess what you really need. – Henrik Schumacher Jul 21 '18 at 15:09
  • The algorithm will be used on all sorts of chain complexes, coming from simplicial complexes, groups, associative and Lie algebras, knots, etc. Up to now (with SNF algorithms), I always ran out of memory. I didn't mean to bother anyone, sorry, but it helps knowing what SparseArray is capable of and how it works, before writing my code... Thank you for your answer, I'm off trying to implement it, later I might have some question about your method. – Leo Jul 21 '18 at 15:26
  • Sorry, I don't understand what's going on in Compile[..], but could you replace your function with a similar fast function row[d_,i_] that returns a list {{j,v},...} of positions $j$ and nonzero values $v$ of the $i$-th row in $d$? – Leo Jul 22 '18 at 00:20
  • 1
    Thinking of it: You will mostly run out of memory because of $O(n^3)$ algorithms. A further copy of the matrix for the transposed matrix should be your least problem. – Henrik Schumacher Jul 22 '18 at 10:10