Fast Evaluation of a series of dot product

Question

Context

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.

Illustration

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!

Answer 1

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

ClearAll[f]
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.

First@AbsoluteTiming[
  ClearAll[f];
  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}];
  ]

0.098479

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:

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

0.038625

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.

First@AbsoluteTiming[
  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;
  ]

0.000677

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]]

1.11022*10^-16

1.77636*10^-15