Efficient memory usage while building a large sparse matrix

Question

I am sure the following problem has been solved already, but I am unable to find any solution... Any help appreciated!

So I am building a pretty huge matrix (or tensor, actually) using

mat=SparseArray[ParallelTable[UnitStep[f[a,b,c,d]],{a,1,D},{b,1,D},{c,1,D},{d,1,D}]]]

In my specific case I have D=180 and f is a simple function (basically addition or subtraction of certain elements stored in a short list), such that the majority of the values UnitStep[f[...]] are zeros and the resulting matrix is indeed sparse.

Now, constructing this array takes quite some time, but more problematically, it uses around 12 GB of memory during the computation. I want to use less memory, but I do not know how. However, I think it should be possible because it is a very sparse array containing only zeros and a few ones. In particular, using

Export["mat.mx",max]

the final file on my computer is only 12 MB (instead of GB). Is there a better way?

Answer 1

Even though you have shared no details about f, as others have mentioned, I think I understand your problem.

Consider this (note that D is a protected symbol in Mathematica):

dim = 40;
x = ParallelTable[
  UnitStep[a - b - c - d], 
  {a, 1, dim}, {b, 1, dim}, {c, 1, dim}, {d, 1, dim}];
y = SparseArray[x];

On my machine, x occupies 65 MB of memory while y just 3 MB. So when you do

z = SparseArray[ParallelTable[
    UnitStep[a - b - c - d], 
    {a, 1, dim}, {b, 1, dim}, {c, 1, dim}, {d, 1, dim}]]; // AbsoluteTiming
(* {0.600104, Null} *)

Mathematica uses 65 MB of memory to store the intermediate calculation result.

Instead, you can generate the list of rules for SparceArray directly:

zz = SparseArray[Reap[Do[
       If[UnitStep[a - b - c - d] == 1, Sow[{a, b, c, d} -> 1]],
       {a, 1, dim}, {b, 1, dim}, {c, 1, dim}, {d, 1, dim}
     ]][[2, 1]], {dim, dim, dim, dim}]; // AbsoluteTiming
(* {3.55455, Null} *)

z == zz (* True *)

The code above is not parallelized and is 5.6x slower than the original on my machine, which has 8 cores. You cannot naively parallelize Sow, but you can use Mr.Wizard's trick with minimal performance loss over the original solution (and significant memory gain!).

ParallelEvaluate[foo = {}];
sow[x_] := (foo = {foo, x};)
ParallelDo[
   If[UnitStep[a - b - c - d] == 1, sow[{a, b, c, d} -> 1]],
   {a, 1, dim}, {b, 1, dim}, {c, 1, dim}, {d, 1, dim}]; // AbsoluteTiming
(zzz = SparseArray[
     Flatten@ParallelEvaluate[Flatten@foo], 
     {dim, dim, dim, dim}];) // AbsoluteTiming

{0.444935, Null}
{0.299598, Null}

zzz == z (* True *)