Timing for generating long matrix

Question

I need a code for thermodynamics research that generates a matrix of 32768 by 32768 elements.

The code is the followng:

s0[MicroNumber_, k_] := If[RealDigits[MicroNumber, 2, 15][[1]][[k]] == 0, -1, 1]
s0[MicroNumber_, 0] = 1;
s0[MicroNumber_, 16] = 1;
s = Table[s0[i - 1, j], {i, 1, 32768}, {j, 0, 16}];

AbsoluteTiming[With[{MicroElements = s}, Table[If[Abs[Sum[(MicroElements[[i]]- MicroElements[[j]])[[d]], {d, 1, 
    16}]] == 2, 1, 0], {i, 1, 32768}, {j, 1, 32768}]]]

This code takes a very long time to run and I am trying to make it run faster. Any ideas on how to make it more efficent? I tried using ParallelTable instead of Table but it didn't help.

Answer 1

I suspect you only want to optimize the last line, i.e. I haven't bothered speeding up s. In what follows I'll use the first 1000 x 1000 submatrix to show some comparative timings.

Your current implementation:

With[{MicroElements = s}, 
   Table[If[Abs[Sum[(MicroElements[[i]] - MicroElements[[j]])[[d]], {d, 1,16}]] == 2, 1, 0],
      {i, 1, 1000}, {j, 1, 1000}]]; // AbsoluteTiming
(*{5.44532,Null}*)

The main reason your implementation is slow is the use of Sum, you should instead use Total since you're trying to sum over the entire list.

With[{MicroElements = s}, 
   Table[If[Abs[Total[(MicroElements[[i]] - MicroElements[[j]])]] == 2, 1,0],
   {i, 1, 1000}, {j, 1, 1000}]]; // AbsoluteTiming
(*{0.681037, Null}*)

We can go a bit further by compiling the function:

cf = Compile[{{x, _Integer, 1}, {y, _Integer, 1}}, 
  Boole[Abs[Total[x - y]] == 2], CompilationTarget -> "C", 
  RuntimeOptions -> "Speed", RuntimeAttributes -> {Listable}, 
  Parallelization -> True]

This seems to help very slightly, but can be used with DistanceMatrix and DistanceFunction to give fairly clean code:

AbsoluteTiming[DistanceMatrix[s[[;;1000]], s[[;;1000]], DistanceFunction -> cf];]
(*{0.469302,Null}*)

Since we specified Listable in our compiled function above, we could also consider using the following:

AbsoluteTiming[Map[cf[s[[;; 1000]], #] &, s[[;; 1000]]];]
(*{0.160225,Null}*)

Since we came this far, we might as-well compile all the way:

cf2 = Compile[{{mat, _Integer, 2}}, 
  Table[Boole[Abs[Total[i - j]] == 2], {i, mat}, {j, mat}], 
  CompilationTarget -> "C", RuntimeOptions -> "Speed"]

cf2[s[[;;1000]]];//AbsoluteTiming
(*{0.105187,Null}*)

As for the actual calculation, the methods above appear to exhaust the memory on my machine (32768 x 32768 is a large array!), but some simple scaling suggests it'll take ~2 mins to complete

With[{timings = {2^#, First[AbsoluteTiming[cf2[s[[;; 2^#]]];]]} & /@ Range[9, 14]},
 Echo[ListLogLogPlot[timings]];
 2^(Fit[Log[2, timings], {1, x}, x] /. x -> Range[9, 15])
 ]

(*{0.0275781, 0.109413, 0.434082, 1.72217, 6.83248, 27.107, 107.544}*)

Finally, note that your matrix is actually symmetric - so you might even consider using Subsets to only compute the upper-triangular part of the matrix, as demonstrated here.

Answer 2

This seems to produce the correct result and is much faster for small values of n.

n = 32768;

x = Total[2 RealDigits[Range[0, n - 1], 2, 15][[All, 1]] - 1, {2}] + 1;
nf = Nearest[x -> "Index"];
B = SparseArray[
   Join @@ MapThread[
      Thread[{#1, Complement[#2, #3]}] &,
      {Range[Length[x]], nf[x, {\[Infinity], 2}], 
       nf[x, {\[Infinity], 1}]}
      ] -> 1,
   {n, n}
   ];

Note that x is effectively the same as Total[s[[All,1;;16]],2]. So you only have to compute only differences of scalars instead of sums of differences of vectors.