I am interested in the most efficient way to transform a huge rectangular array to a sparse array structure.
Consider a numeric rectangular table that is rather sparse, with zero positions being set to Null rather than 0, which saves time and memory during the fill. The problem is that SparseArray function can not be directly applied to that kind of table. First, every Null should be replaced by 0, as can be seen in the following example:
Q = 100;(*the matrix dimension is QxQ.*)
f[p1_,p2_]:=If[PrimeQ[p1+p2],1,Null];(*the specific form of the matrix elements chosen for illustrative purposes only*)
mat=ParallelTable[f[p1,p2],{p1,Q},{p2,Q}];(*the matrix to be later compressed to SparseArray*)
mat0=mat/.(Null)->0;(*Here we replace Nulls with zeroes*)
sparse=SparseArray[mat0];(*Now the SparseArray may be applied*)
However, during such replacement procedure, the ByteCount[mat0] is dramatically increased, and at sufficiently large Q one runs out of memory. Besides, it does not seem logical to artificially inflate the matrix in order to be able to compress it afterwards.
So, how can I transform a huge rectangular array full of Nulls to a sparse array structure, without replacing Nulls with zeroes?
Thank you!
P.S. Other ways to form a sparse array that I am aware of are getting too speed- and memory inefficient as Q is taken large (say, Q > 50000), e.g. this way
SparseArray @ Flatten[ParallelTable[{p1, p2} -> f[p1, p2], {p1,Q}, {p2,Q}] /. (_ -> Null) -> Sequence[]];
takes 100Gb of memory to form a sparse array of 10 Gb.
P.P.S. One could ask, what's the need to compress a matrix that is already quite sparse? The SparseArray version takes less memory, anyway: ByteCount[sparse] is less than ByteCount[mat].
P.P.P.S. Above was the toy problem. The more realistic code is as follows:
<< CompiledFunctionTools`
Compiler`$CCompilerOptions = {"SystemCompileOptions" ->
"-fPIC -Ofast -march=native"};
On[ Compile::noinfo]
n = 4;(*the number of electrons*)
Upp = 6;(*Upp is the upper occupied
quantum state*)
a = Subsets[Range[Upp], {n}];(*a is the set of all many-particle states*)
Q = Binomial[Upp, n];(*The number of many-particle states in a*)
(*Fast compiled function that compares two vectors and returns the positions of different elements.*)
VectorCompare =
Compile[{{v1, _Integer, 1}, {v2, _Integer, 1}}, Block[{i1 = 1, i2 = 1, d1 = Internal`Bag@Most[{0}],
d2 = Internal`Bag@Most[{0}]},
(*Run along the lists,recording differences as we go*)
While[i1 <= Length[v1] && i2 <= Length[v2],
Which[v1[[i1]] < v2[[i2]], Internal`StuffBag[d1, i1]; i1++,
v1[[i1]] > v2[[i2]], Internal`StuffBag[d2, i2]; i2++, True, i1++;
i2++]];
(*Fix up in case we ran off the end of one of the lists*)
While[i1 <= Length[v1], Internal`StuffBag[d1, i1]; i1++];
While[i2 <= Length[v2], Internal`StuffBag[d2, i2]; i2++];
{Internal`BagPart[d1, All], Internal`BagPart[d2, All]}],
"CompilationTarget" -> "C",
CompilationOptions -> {"ExpressionOptimization" -> True,
"InlineExternalDefinitions" -> True}, RuntimeOptions -> "Speed",
RuntimeAttributes -> {Listable}, Parallelization -> True];
CompareMatrix = SparseArray@
Developer`ToPackedArray[
ParallelTable[
If[p1 < p2, vc = VectorCompare[a[[p1]], a[[p2]]];
diff = Length@Flatten@vc;
Which[diff == 2, {vc[[1, 1]], vc[[2, 1]], 0, 0},
diff == 4, {vc[[1, 1]], vc[[2, 1]], vc[[1, 2]],
vc[[2, 2]]}]], {p1, Q}, {p2, Q}] /. (Null) -> {0, 0, 0, 0}];
The array I need could be built without any unnecessary zeros by running the following code:
CMB2 = Internal`Bag@Most[{0}];
CMB4 = Internal`Bag@Most[{0}];
Do[vc = VectorCompare[a[[p1]], a[[p2]]]; l = Length@Flatten@vc;
Which[l == 2,Internal`StuffBag[CMB2, {p1, p2, vc[[1, 1]], vc[[2, 1]]}],
l == 4,Internal`StuffBag[CMB4, {p1, p2, vc[[1, 1]], vc[[2, 1]], vc[[1, 2]], vc[[2, 2]]}]], {p1, 1, Q - 1}, {p2, p1 + 1, Q}];
If only it were parallelizable...
Edit: The fastest and memory-efficient solution is as follows:
CM=Join @@ ParallelMap[Developer`ToPackedArray,
Table[vc=VectorCompare[a[[p1]],a[[p2]]];diff = Length@Flatten@vc;
Which[diff > 4, empty,
diff == 4, {p1,p2,vc[[1,1]],vc[[1,2]],vc[[2,1]],vc[[2,2]]},
True, {p1,p2,vc[[1,1]],vc[[2,1]],0,0}],
{p1,1,Q-1},{p2,p1+1,Q}]];
In essence, this closes the question.
Null
here instead of0
(or0.
for a matrix of reals) is much more likely to increase memory usage ... $\endgroup$ByteCount
is accurate here. It's known to double count in certain cases when several parts of an expression are really a reference to the same thing in memory.ByteCount@Outer[If[PrimeQ[#1 + #2], 1] &, Range[100], Range[100]]
is indeed about half the size ofByteCount@Outer[If[PrimeQ[#1 + #2], 1, 0] &, Range[100], Range[100]]
. But if I useMaxMemoryUsed
(introduced in v9 I think, at least for use this way) instead ofByteCount
, I get the same result for both. $\endgroup$