Fast evaluation of a function in many points

Question

I need to feed to an external program a number of points (in Complex128 format) generated from the numerical evaluation of some function, e.g. $e^{i \vec{k}\cdot\vec{x}}$. The function is evaluated at points of some evenly spaced grid. I currently generate the values with Table (possibly multidimensional, and then Flatten'd) and send it to a file with BinaryWrite. However, this is rather slow, because I use a number of points in the order of $2^{24}$ or more. I tried with ParallelTable, but I inspected with perf and it appears that when transferring the data to the main kernel there's a lot of UTF string encoding going on for whatever reason, which defeats the speedup in the parallelization.

Is there any smarter way to accomplish what I need?

Answer 1

In cases like this where the function to be evaluated a) has the Listable attribute, and b) is not very compicated, we can get a huge benefit from exploiting this listability. In my experience, exploiting vectorized operations, listability and packed arrays, along with occational compilation to C, usually has a much larger impact on speed than parallellizing.

One dimension

For k = 2.5 and $2^{24}$ points:

First@AbsoluteTiming[ans1 = Exp[I*2.5*Range[2^24]]];
(* 1.270682 *)

First@AbsoluteTiming[ans2 = Table[Exp[I*2.5*x], {x, 2^24}]];
(* 19.620243 *)

First@AbsoluteTiming[ans3 = Exp /@ (I*2.5*Range[2^24])];
(* 3.358687 *)

ans1 == ans2 == ans3
(* True *)

Note that the last one uses the Listable attribute of Times, but not of Exp. On my machine it's still way faster than Table though.

Three dimensions

Since $e^{a+b} = e^a e^b$ we can still do most of the calculation vectorized. Note that we don't need a square grid. For $\vec{k} = \{k_x, k_y, k_z\} = \{2.5, 3.5, 4.5\}$:

(* using Table *)
First@AbsoluteTiming[
  ans1 = Table[Exp[I (k1 x + k2 y + k3 z)]
  , {x, 100}, {y, 101}, {z, 102}];
]
(* 1.979192 *)

(* using Table, but vectorizing over one coordinate *)
First@AbsoluteTiming[
  ans2 = Table[Exp[I (k1 x + k2 y + k3 Range[102])]
  , {x, 100}, {y, 101}];
]
(* 0.225951 *)

(* My favorite *)
First@AbsoluteTiming[
  x = Exp[I k1 Range[100]];
  y = Exp[I k2 Range[101]]; 
  z = Exp[I k3 Range[102]];
  ans3 = Outer[Times, x, y, z];
]
(* 0.028202 *)

They are equal up to roundoff:

Chop[ans1 - ans2] == Chop[ans1 - ans3] == ConstantArray[0, {100, 101, 102}]
(* True *)

More complicated example

When the function to be evaluated is not easily split into parts that depend only on one component, I usually check this list to see if it can be compiled. If for instance the function is

f[x_, y_, z_] := Sin[2x + Exp[I(x + y*z)]]

we find:

First@AbsoluteTiming[ans1 = Table[f[x, y, z], {x, 100}, {y, 101}, {z, 102}];]
(* 4.599261 *)

First@AbsoluteTiming[ans2 = Outer[f, Range[100], Range[101], Range[102]];]
(* 3.056824 *)

(* We can still vectorize over one coordinate *)
First@AbsoluteTiming[ans3 = Table[f[x, y, Range[102]], {x, 100}, {y, 101}];]
(* 1.463473 *)

But here, a compiled function is much faster:

comp = Compile[{},
  Table[Sin[2 x + Exp[I (x + y*z)]]
  ,{x, 100}, {y, 101}, {z, 102}], 
  CompilationTarget -> "C", RuntimeOptions -> "Speed"
];
First@AbsoluteTiming[ans4 = comp[];]
(* 0.145347 *)