# Parallelize evaluation of function with memoization

I have a complicated function that I need multiple times, so I want to memoize it and have the first evaluation done in parallel. Unlike in my example below it's not a continuous function, so interpolation is not an option. (In fact, it's values are also functions.) The naive approach clearly does not work, because the memoized value is then only known on the kernel it was evaluated on:

LaunchKernels;
f[x_] := f[x] = (Pause; N[Sin[x]]); (*Expensive calculation*)
ParallelDo[f[x], {x, 3}];

ParallelEvaluate[AbsoluteTiming[f]]

(* ==> {{3.000632, 0.841471}, {0.000024, 0.841471}} *)


I believe I found a workaround by doing something like this:

f[x_] := (Pause; N[Sin[x]]); (*Expensive calculation - NO memoization*)
t = ParallelTable[f[x], {x, 3}];
Do[f[x] = t[[x]], {x,3}];


Using SetSharedFunction[f] before ParallelDo also yields a non-optimal result:

{{0.012051, 0.841471}, {0.012202, 0.841471}}


0.01 s is still a long time to look up a value (see above, it should be < 1 ms). Is there something more elegant or do I have to keep it like this?

Edit: Just to be clear, the workaround above without Shared Functions works, it runs in parallel and the main kernel knows the values afterwards, but it strikes me as an ugly hack. I was wondering if there was an "official" solution.

• The problem with SetSharedFunction is that it forces f to be evaluated on the main kernel: this means that you will lose parallelization and that there is extra communication overhead. Would a solution where long evaluations are done on the subkernels, but retrieval of memoized values is done on the main kernel be useful for you? If the function is expensive enough, this should be useful. If it is not extremely expensive to calculate, the communication overhead might hurt more (as you noticed above where you got ~12 ms timings) Feb 3, 2012 at 18:17
• @Szabolcs Sounds like a very meaningful suggestion. Why don't you make an answer out of it? Feb 3, 2012 at 18:25
• @Szabolcs I for one could definitely use something like that. Feb 3, 2012 at 18:30
• @Leonid and David I wrote an answer now. I got confused by the fact that undefined shared functions return Null on parallel kernels, and couldn't figure out why my initial code did not work. Feb 3, 2012 at 18:49
• I updated my answer a bit. Feb 4, 2012 at 11:56

The problem with SetSharedFunction is that it forces f to be evaluated on the main kernel: this means that if you simply do

SetSharedFunction[f]


then you will lose parallelization (a timing of ParallelTable[f[x], {x, 3}] will give about 9 seconds).

This property of SetSharedFunction is not clear from the documentation in my opinion. I learned about it from this answer. It is also not clear if the behaviour is the same in version 7 (can someone test? I tested my answer on version 8 only).

We can however store the memoized vales on the main kernel, while evaluating the expensive computations on the paralell kernels. Here's one way to do it:

f[x_] :=
With[{result = g[x]},
If[result === Null,
g[x] = (Pause; N[Sin[x]]),
result
]
]

SetSharedFunction[g]


Here I used the special property of shared functions that they return Null on the paralell kernels when they have no value for a given argument.

The first time we run this, we get a 6 s timing, as expected:

AbsoluteTiming@ParallelTable[f[x], {x, 3}]

(* ==> {6.0533462, {0.841471, 0.909297, 0.14112}} *)


The second time it will be very fast:

AbsoluteTiming@ParallelTable[f[x], {x, 3}]

(* ==> {0.0260015, {0.841471, 0.909297, 0.14112}} *)


However, as you noticed, evaluating f on the parallel kernels is a bit slow. On the main kernel it's much faster. This is due to the communication overhead: every time f is evaluated or changed on a subkernel, it needs to communicate with the main kernel. The slowdown does not really matter if f is really expensive (like the 3 seconds in your toy example), but it can be significant if f is very fast to execute (comparable in time to the apparently ~10 ms communication overhead)

ParallelTable[AbsoluteTiming@f[x], {x, 3}]

(* ==> {{0.0100006, 0.841471}, {0.0110006, 0.909297}, {0.0110007,  0.14112}} *)

Table[AbsoluteTiming@f[x], {x, 3}]

(* ==> {{0., 0.841471}, {0., 0.909297}, {0., 0.14112}} *)


Finally a note about benchmarking: in general, measuring very short times like the 10 ms here should be done with care. On older versions of Windows, the timer resolution is only 15 ms. On Windows 7, the resolution is much better. These timings are from Windows 7.

Update:

Based on @Volker's @Leonid's suggestion in the comments, and @Volker's original solution, we can combine subkerel and main kernel caching like this:

f[x_] := With[{result = g[x]}, (f[x] = result) /; result =!= Null];
f[x_] := g[x] = f[x] = (Pause; N[Sin[x]])


## Packaged up solution

We can bundle all these ideas into a single memoization function (see the code at the end of the post).

Here is an example use:

Clear[f]
f[x_] := (Pause; Sin[x])

AbsoluteTiming@ParallelTable[AbsoluteTiming@pm[f][Mod[x, 3]], {x, 15}]

(* ==>
{6.0683471, {{3.0181726, Sin}, {0.0110007, Sin},
{3.0181726, 0}, {0., Sin}, {3.0191727, Sin},
{3.0181726, 0}, {0.0110007, Sin}, {0., Sin},
{0., 0}, {0., Sin}, {0., Sin},
{0., 0}, {0., Sin}, {0., Sin}, {0., 0}}}
*)


The function simply needs to be called using pm[f][x] instead of f[x]. Memoized values are associated with f as UpValues, so I thought automatic distribution of definitions should take care of both synchronizing memoized values and clearing them when necessary. Unfortunately this mechanism doesn't seem to be reliable (sometimes it works, sometimes it doesn't), so I provided a function clearParallelCache[f] that will clear all memoized values on all kernels.

Caching happens at two levels: on the main kernel level and on subkernels. Computed or main-kernel-cached values are copied to the subkernels as soon as possible. This is visible in the timings of the example above. Sometimes retrieving cached values takes 10 ms, but eventually it becomes very fast. Note that it might happen that the two kernels will each compute the same value (if they start computing it at the same time). This can sometimes be avoided by using a different setting for the Method option of Parallelize (depending on the structure of the input data).

For simplicity, I restricted pm to only work with functions that take a single numerical argument (and return anything). This is to avoid having to deal with more complex conditional definitions (especially cases when the function won't evaluate for certain argument types). It could safely be changed to accept e.g. a vector or matrix of values instead.

## The code

pm[f_Symbol][x_?NumericQ] :=
With[{result = memo[f, x]},
If[result === Null,
With[{value = valueOnMain[f, x]},
If[value === Null,
f /: memo[f, x] = setValueOnMain[f, x, f[x]],
f /: memo[f, x] = value]
],
result]
]

memo[___] := Null
DistributeDefinitions[memo];

valueOnMain[f_, x_] := memo[f, x]

setValueOnMain[f_, x_, fx_] := f /: memo[f, x] = fx

SetSharedFunction[valueOnMain, setValueOnMain]

clearParallelCache[f_Symbol] :=
(UpValues[f] = {}; ParallelEvaluate[UpValues[f] = {}];)

• Can't you then just distribute the definition of g instead of making it a shared function? But actually it might be sufficient if the main kernel knows the values after the hard part of the computation. Feb 3, 2012 at 18:58
• @Volker Distributing definitions is not enough (even if DistributeDefinitions were directly available on subkernels---it is not), we'd need to merge them, merge all results from all subkernels. Generally, we have a tradeoff here: we can either synchronize the memoized values often, like I did, and accept the communication overhead, or synchronize them only occasionally, and accept that sometimes two different kernels will compute the slow-to-evaluate function for the same argument. Feb 3, 2012 at 19:03
• Maybe the best solution is to mix my approach and yours: after the ParalellTable has finished, distribute all memoized values to all kernels, and in subsequent computations only access the main kernel from subkernels when really needed. This is what you actually mean, right? Feb 3, 2012 at 19:05
• @Szabolcs Here is a version which solves both problems: ClearAll[f, g]; f[x_] := With[{result = g[x]}, f[x] = result /; result =!= Null];f[x_] := g[x] = f[x] = (Pause; N[Sin[x]]);SetSharedFunction[g]. The f on the sub-kernel will use its memoized value if that exists. If not, it will "ask" g (so only in this case there will be communication overhead). This has an advantage that there is no need to explicitly distribute definitions after the fact (for the cost of being a tiny bit slower until all kernels "learn" a given value of f, but this is really a small cost IMO). Feb 3, 2012 at 19:43
• @Leonid Your version, as you posted it here simply does not work correctly on this machine I'm testing on. It does a lot of extra calculations and takes 9 seconds (on this dual-core). However, if I add a few debug Prints, then it does work. It's just driving me crazy. It's like quantum mechanics: if I look at it, it behaves differently than when you don't look at it. I need a break now, and if I still can't debug it this afternoon, I'll mail you. Feb 4, 2012 at 12:44

The problem with a normal shared memoized function definition

f[x_] := f[x] = (Pause; N[Sin[x]])


is indeed that any evaluation of f[n] on a parallel subkernel causes the rhs f[x] = (Pause; N[Sin[x]]) to be sent back and evaluated on the master kernel.

There is an undocumented feature that if such a definition is made on a subkernel, the rhs will also be flagged for evaluation on the subkernel. So you could do instead

SetSharedFunction[f]; (* no definitions on master *)
ParallelEvaluate[f[x_] := f[x] = (Pause; N[Sin[x]]), First[Kernels[]]]


If you look at the result, you see

?f
f[x_]:=ParallelDeveloperSendBack[f[x]=(Pause;N[Sin[x]])]


Knowing this you can get the desired effect with a definition on the master kernel:

g[x_] := ParallelDeveloperSendBack[g[x] = (Pause; N[Sin[x]])]
SetSharedFunction[g]


and then

In:= AbsoluteTiming[ParallelTable[g[Mod[i, 4]], {i, 40}]]

Out= {3.13635, {0.841471, 0.909297, 0.14112, 0., 0.841471,
0.909297, 0.14112, 0., 0.841471, 0.909297, 0.14112, 0., 0.841471,
0.909297, 0.14112, 0., 0.841471, 0.909297, 0.14112, 0., 0.841471,
0.909297, 0.14112, 0., 0.841471, 0.909297, 0.14112, 0., 0.841471,
0.909297, 0.14112, 0., 0.841471, 0.909297, 0.14112, 0., 0.841471,
0.909297, 0.14112, 0.}}

?g
g=0.

g=0.841471

g=0.909297

g=0.14112

g[x_]:=ParallelDeveloperSendBack[g[x]=(Pause;N[Sin[x]])]


Each evaluation of g[n] on the subkernel will then first send g[n] to the master. If there is already a specific value for that n it is sent back; otherwise the expression g[x] = (Pause; N[Sin[x]]) is sent back. The subkernel then evaluates its rhs and makes the definition g[x]=val, which expression is now sent back to the master where it is performed, for all kernels to see.

• Welcome to Mathematica.SE and thank you for sharing this information. Very interesting and useful feature! Sep 1, 2016 at 14:57
• The timing 3.13635 above corresponds to a 4 subkernels machine (or more). Sep 2, 2016 at 23:03

I realize this is an old question, but I recently had the same issue and have come across (link to google groups question) what I think is a cleaner solution. I don't want to take credit for coming up with that solution, but I thought it would be helpful to add it to this site.

I'll use a simple example function to demonstrate.

 f[x_] := f[x] = x
ParallelDo[f[i], {i, 1, 3}]


The problem is that, while f, f, and f have now been computed, they are stored on the subkernels and not on the main one. But we can access the values using ParallelEvaluate:

 DownValues[f] = DeleteDuplicates[Flatten[ParallelEvaluate[DownValues[f]]]]


DeleteDuplicates is needed because each subkernel includes the original function definition as a down-value.

• Welcome to Mathematica.SE! I suggest the following: 1) As you receive help, try to give it too, by answering questions in your area of expertise. 2) Take the tour! 3) When you see good questions and answers, vote them up by clicking the gray triangles, because the credibility of the system is based on the reputation gained by users sharing their knowledge. Also, please remember to accept the answer, if any, that solves your problem, by clicking the checkmark sign! Jan 13, 2016 at 23:17

Here is an awkward but potentially useful method:

SetAttributes[newSet, HoldFirst]
newSet[lhs_, rhs : Except[_EvaluationObject]] := Set[lhs, rhs]

LaunchKernels;

f[x_] := f[x] ~newSet~ ParallelSubmit@Sum[N@Sin[i], {i, 1000000 x}]

WaitAll@Table[f[x], {x, {1, 2, 3}}] // AbsoluteTiming
WaitAll@Table[f[x], {x, {3, 2, 1}}] // AbsoluteTiming


{2.7351565, {-0.11711, -0.103631, 0.0387313}}

{0., {0.0387313, -0.103631, -0.11711}}

Depending on your application this may be completely useless, since nothing not already memoized will evaluate until you use WaitAll on explicit output. On the other hand it performs well, and you can do without the Parallel* version of functions, which at times may be useful.

• See my comment above on what's wrong with this. In version 8 you lose all parallelization (just time the ParallelDo --- you'll get 9 seconds). EDIT I just realized you have v7, is it different there? Feb 3, 2012 at 18:19
• I was lucky enough to have witnessed the moment when newSet was not yet infixed. It was, however, hopelessly brief. +1. Feb 3, 2012 at 19:58
• @Leonid LOL -- I actually started with it infixed, then thought it would just stir up trouble, then rethought that it looked more like Set if it was infix. Maybe I should change it back again... Feb 3, 2012 at 20:04