I have a function that depends on 3 real variables x,y and z and that is defined by a series of matrix products. The evaluation of f for a specific (x,y,z) is fast ~0.02 sec but I want to evaluate the function on a huge number of points (a regularly spaced grid of x,y and z values) which in the end makes the evaluation really slow if not unmanageable. I have already tried what was proposed in this answer, but my function is not compilable, and ParalellTable is faster that vectorizing on my laptop.


For the sake of simplicity let me illustrate this with

weight = RandomReal[1, 200];
pts = RandomReal[1, 200];
M = RandomReal[1, {200, 200}];
f[x_,y_] = (weight*Exp[-pts*x]).Exp[M].(weight*Exp[-pts*y])//N

How would one make the evaluation of f on multiple couples of (x,y) faster than relying on ParallelTable ?

ans = ParallelTable[f[x,y],{x,Range[100]},{y,Range[100]}]

Thanks a lot for your help!


First of all, it helps a great deal to define f as follows:

A = Exp[M];
f[x_, y_] := (weight*Exp[-pts*x]).A.(weight*Exp[-pts*y]);

Here, using SetDelayed (:=) is better than Set (=), because the latter creates a humongous symbolic expression that is very expensive to numericise. Defining A first ensures that Exp[M] is computed only once.

  A = Exp[M];
  f[x_, y_] := (weight*Exp[-pts*x]).A.(weight*Exp[-pts*y]);
  ans = ParallelTable[f[x, y], {x, 1., 100}, {y, 1., 100}];


Still, for each x the vector (weight*Exp[-pts*x]) is computed 100 times -- unless the JIT compiler does something clever. This can be mended as follows:

  A = Exp[M];
  X = Table[weight Exp[-x pts], {x, 1., 100.}];
  Y = Table[weight Exp[-y pts], {y, 1., 100.}];
  ans2 = Block[{u},
     u = X[[i]].A;
     Table[u.Y[[j]], {j, 1, 100}]
     , {i, 1, 100}]


Moreover, we can fuse all the $2 \times 100 \times 100 = 20000$ Dot operations in the ParallelTable into two matrix-matrix Dots. This does not reduce the number of flops, but it reduces the calling overhead and allows the Dot routine to exploit horizontal vectorization and cache friendly task decomposition.

  A = Exp[M];
  U = Exp[-KroneckerProduct[Range[1., 100.], pts]].DiagonalMatrix[
     SparseArray[weight ]];
  V = weight Exp[-KroneckerProduct[pts, Range[1., 100.]]];
  ans3 = U.A.V;


This is about $200000$ times faster than the original code.

Just to check the relative errors:

Max[Abs[ans/ans2 - 1]]
Max[Abs[ans/ans3 - 1]]



  • $\begingroup$ Thanks a lot ! This is indeed much smarter than what I was doing :) $\endgroup$ – Hermanter Mar 1 at 11:35
  • $\begingroup$ You're welcome. $\endgroup$ – Henrik Schumacher Mar 1 at 11:36

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.