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?

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

1 Answer 1

Reset to default
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}} *)
Improve this answer
$\endgroup$
5
  • $\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, 2016 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, 2016 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, 2016 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, 2016 at 20:49
  • $\begingroup$ Ah. Nice detective work. $\endgroup$
    – Daniel Lichtblau
    May 20, 2016 at 21:12

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

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