3
$\begingroup$

Compare two ways of computing a table of values of a multivariable Gaussian:

First@AbsoluteTiming[
  ta = Table[Exp[-Total[Array[r, 5]^2]],
  {r[1], 1., 2.},
  {r[2], 1., 3.},
  {r[3], 1., 3.},
  {r[4], 1., 23.},
  {r[5], 1., 421.}
  ];
]
(* 2.1069 *)

comp = Compile[{},
  Table[Exp[-(r[1]^2 + r[2]^2 + r[3]^2 + r[4]^2 + r[5]^2)],
  {r[1], 1., 2.},
  {r[2], 1., 3.},
  {r[3], 1., 3.},
  {r[4], 1., 23.},
  {r[5], 1., 421.}
  ], CompilationTarget->"C"
];
First@AbsoluteTiming[
  tb = comp[];
]
(* 0.002853 *)

I would like to get such a speedup in a case where I specify the number of iterators/variables programmatically, and where the iterator values are stored in a list itvals. Non-compiled is easy, e.g.:

n = 5; (* number of iterators *)
itvals = RandomReal[{0,1}, {n, 8}];
tc = Fold[Table[#1,{r[#2],itvals[[#2]]}]&, Exp[-Total[Array[r,n]^2]], Range[n]];

Or:

td = Table[Exp[-Total[Array[r,n]^2]],##] & @@ Table[{r[i], itvals[[i]]}, {i,n,1,-1}];

tc == td
(* True *)

But I can't get this to work inside Compile! If I try

c1 = Compile[{{n, _Integer}},
  Block[{itvals},
   itvals = RandomReal[{0, 1}, {n, 5}];
   Fold[Table[#1, {x[#2], itvals[[#2]]}] &, 1, Range[n]]
   ]
  ]

where I use the "function" 1 instead of the Gaussian for simplicity. CompilePrint[c1] shows that the whole Fold is wrapped in a MainEvaluate. I can't do it similarly to td either because Compile only supports @@ with the functions Plus, Times or List...

Anybody know a way to do it?

$\endgroup$
  • $\begingroup$ Could you show the code you attempted as your compiled version? $\endgroup$ – MarcoB May 19 '16 at 19:30
  • $\begingroup$ @MarcoB, I've added it. $\endgroup$ – Marius Ladegård Meyer May 19 '16 at 20:52
6
$\begingroup$

Verily, this is a headache, since Table wants its iterators as Sequence rather than nested list of lists. Here is a method I used quite recently to get sufficiently fast code for this.

c1c[n_] := 
 With[{itvals = RandomReal[{0, 1}, {n, 5}]}, 
  With[{iters = Apply[Sequence, Table[{x[j], itvals[[j]]}, {j, n}]]},
   c1[n] = cCompile[{},
      tTable[1, iters]] /. {cCompile -> Compile, tTable -> Table}]]

In[55]:= c1c[2][]

(* Out[55]= {{1, 1, 1, 1, 1}, {1, 1, 1, 1, 1}, {1, 1, 1, 1, 1}, {1, 1, 1,
   1, 1}, {1, 1, 1, 1, 1}} *)
$\endgroup$
  • $\begingroup$ This looks promising! I can see that each list on level 1 of itvals gets copied at the start of the compiled code, which seems unnecessary. But since it only happends once at the start and each element of itvals is short, it should not matter, right? $\endgroup$ – Marius Ladegård Meyer May 20 '16 at 5:34
  • $\begingroup$ As a side question, do you know why tc with Fold evaluates much faster than td with iterators spliced in? I had not expected that. $\endgroup$ – Marius Ladegård Meyer May 20 '16 at 8:28
  • $\begingroup$ I do not know why the Fold multi-iterator is faster. I had not done a time test and I confess I would not have expected that. Maybe related to your observation about the list copying? $\endgroup$ – Daniel Lichtblau May 20 '16 at 15:31
  • 2
    $\begingroup$ I figured it out, it was simply a matter of how long the body of the Table was held. In the Fold case, the Gaussian is evaluated as soon as the iterator in the innermost Table has a value, and the inner Table has elements like Exp[-0.527 - r[2]^2 - r[3]^2 - ...]. This gets passed as the body to the next step in the Fold. OTOH when splicing in, the body is held all the way, so for each set of iterator values it has to turn Exp[-Total[Array[r,n]^2]] into Exp[-r[1]^2-r[2]^2-...] before inserting values. Simply wrapping the body of td in Evaluate made it faster than tc. $\endgroup$ – Marius Ladegård Meyer May 20 '16 at 20:49
  • $\begingroup$ Ah. Nice detective work. $\endgroup$ – Daniel Lichtblau May 20 '16 at 21:12

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.