Quickly creating a sparse array

Question

The motivation for my question is that I'm creating very large, very sparse symmetric matrices and finding the largest eigenvalue. (Think upwards of $10^5\times 10^5$ size matrices, with at most $10^2$ nonzero entries on each row. However, to give a workable example, below I'll use $10^4\times 10^4$ size matrices, with $30$ nonzero entries.) The entries are mostly independent of one another, although there is some time savings in working with one row at a time. (In the example below, I model this type of behavior by having the constant $c$ in the computation depend only on the row index. In the actual work I do, the dependence isn't quite that simple.) By symmetry, we may focus on the part of the matrix above the main diagonal. (For simplicity, I'm assuming the main diagonal is zero in the example below.)

So, to give a toy example consider the following code:

M = Table[0.0, {i, 1, 10000}, {j, 1, 10000}];
Do[
   c = RandomInteger[{1, 10}];
   PossibleColumns = RandomSample[Range[i + 1, 10000], Min[10000 - i - 1, 30]];
   Do[
      M[[PossibleColumns[[j]], i]] = c/RandomInteger[{1, 10}];
      M[[i, PossibleColumns[[j]]]] = c/RandomInteger[{1, 10}],
   {j, 1, Length[PossibleColumns]}
   ],
{i, 1, 9999}
];
Eigenvalues[M, 1];

The creation of the large array takes a couple seconds, and seems to fill up about 1 GB of RAM. (For the larger matrices I want to build, I quickly run out of RAM.) It takes about 3 seconds to update the nonzero entries in the matrix $M$. Finally, the eigenvalue function takes about 3 minutes.

If I replace the top line with the new line

M = SparseArray[{}, {10000, 10000}, 0.0];

and rerun the entire code, then there are no RAM problems (even for much, much larger matrices). The eigenvalue step is almost instantaneous. However, the updating takes over an hour! I was very surprised that the updating takes so long, so I made the following intermediate version, which rather than updating each entry of $M$ independently, it updates each row of $M$ once.

M = SparseArray[{}, {10000, 10000}, 0.0];
Do[
   c = RandomInteger[{1, 10}];
   jRow = SparseArray[{}, {10000}, 0.0];
   PossibleColumns = RandomSample[Range[i + 1, 10000], Min[10000 - i - 1, 30]];
   Do[
      jRow[[PossibleColumns[[j]]]] = c/RandomInteger[{1, 10}],
   {j, 1, Length[PossibleColumns]}
   ];
   M[[i]] = jRow,
{i, 1, 9999}
];
M = M + Transpose[M];
Eigenvalues[M, 1];

This has the RAM savings, and the eigenvalue savings, and finishes the updating of $M$ in a little over 4.5 minutes.

My question is this: Is there a way to keep (most of) the RAM savings of using a sparse array, along with the time saving in the eigenvalue computation, but also reduce the time used in updating $M$ down to the timing in the first example (which only used 3 seconds)?

Answer 1

You can use the first documented usage for SparseArray:

So what you want to do is collect all of the rules {i,j}->val before passing them to SparseArray. There is already a built-in method for efficiently collecting a set of rules with no duplicate keys, namely AssociateTo.

So with a slight modification of the OP's code:

data=<||>;
Do[
   c = RandomInteger[{1, 10}];
   PossibleColumns = RandomSample[Range[i + 1, 10000], Min[10000 - i - 1, 30]];
   Do[
      AssociateTo[data, {PossibleColumns[[j]],i} -> c/RandomInteger[{1, 10}]];
      AssociateTo[data, {i,PossibleColumns[[j]]} -> c/RandomInteger[{1, 10}]];
      ,
   {j, 1, Length[PossibleColumns]}
   ],
{i, 1, 9999}
];

The above runs in ~3.5 seconds on my machine, and constructing the SparseArray is very fast as well:

SparseArray[Normal[data]] // AbsoluteTiming

Answer 2

Let's assume we are given the following data:

n = 100000;
k = 100;

g[i_] := RandomSample[i + 1 ;; n, Min[n - i - 1, k]];
f[i_, jlist_] := RandomReal[{1, 10}, Length[jlist]];

rowcolumnindices = g /@ Range[1, n - 1];
rowlengths = Length /@ rowcolumnindices;
rowvalues = MapIndexed[
   {jlist, i} \[Function] f[i[[1]], jlist], 
   rowcolumnindices
   ];

Probably the most efficient way to populate the sparse array with this data is by preparing the CRS data by hand and to use the following undocumented constructor:

First @ AbsoluteTiming[
 orderings = Ordering /@ rowcolumnindices;
 columnindices = Partition[Join @@ MapThread[Part, {rowcolumnindices, orderings}], 1];
 rowpointers = Accumulate[Join[{0}, Length /@ rowcolumnindices, {0}]];
 values = Join @@ MapThread[Part, {rowvalues, orderings}];
 M = # + Transpose[#] &[
   SparseArray @@ {Automatic, {n, n}, 0., {1, {rowpointers, columnindices}, values}}
   ];
 ]

1.53083

I used machine real values here (not rational numbers) because they will be processed much faster, in particular in Eigenvalues.

The CRS format assumes that the column indices for each row are in ascending order. This is why we have to sort the column indices and the corresponding values first (this explains the extra fuzz around orderings).