Parallelization of the stochastic Euler scheme

I wrote a simulation to approximate the law of a stochastic differential equation via a Monte Carlo method using the stochastic Euler scheme. Then I thought, it would be a good idea to speed up things by means of parallelization, but although I think it is an easy example (and similar problems were discussed here lots of times), none of my attempts was successful.

First the code (you can find some explanations below):

(* Parameters of the SDE and radius of the Ball *)
T = 10; sigma = 3.0; z = 1; R = 5*10^6;

(* w denotes the noise and h the step size *)
Euler[x_, w_, h_] := x - x^3*h + sigma*w;

(* Noise is a List of Lists of random numbers *)
MCEuler[h_, StepsMC_, Noise_] := Quiet[Check[
Total[Table[
Fold[ Euler[#1, #2, h] &, z, Noise[[i]]]
, {i, 1, StepsMC}]
]/StepsMC,
Indeterminate
]];

MCEulerMT[h_, StepsMC_, Noise_] :=
Total[Table[Catch[Fold[If[Abs[Euler[#1, #2, h]] >= R, Throw[0],
Euler[#1, #2, h]] &, z, Noise[[i]]]], {i, 1, StepsMC}]]/StepsMC;

Comparison[h_, StepsMC_] := Module[{Noise},
Noise = Table[RandomReal[NormalDistribution[0, Sqrt[h]], T/h], {StepsMC}];
Return[{StepsMC, h, MCEuler[h, StepsMC, Noise], MCEulerMT[h, StepsMC, Noise]}];
];

Result =
Table[Table[Comparison[h, StepsMC], {h, {2^0, 2^-2, 2^-4, 2^-6}}], {StepsMC, {10^3}}]


MCEuler implements the Monte Carlo Simulation using the standard stochastic Euler scheme. Since the coefficients of the stochastic differential equations doesn't fulfill standard assumptions like Lipschitz-continuity, it is possible, that some "discrete trajectories" explode. Therefore, I added the Quiet and Check commands.

MCEulerMT uses again the Euler scheme, but if some "discrete trajectory" leaves a ball of radius R, we set the value to zero.

The Comparison-method should test both procedures on the same noise.

Then I replaced the Table-Commands in MCEuler, MCEulerMT and the last line of the code (I didn't use it creating the noise, but I think there are no side-effects!?) by "ParallelTable" and the runtime was awful. It was much slower than the code above, and I tried other implementations and commands (e.g. Parallelize) but I didn't make any progress.

Perhaps here are some guys, who could explain me, what's wrong about my idea, and how I could fix it.

-
I gave so many answers on parallel issues but I'm not able at the moment to list them. Please see my profile, there you may find your answers and please be aware that you need to specify to share your definitions among the kernels.... – Stefan Jul 6 '13 at 21:21
I looked over your answers, but unfortunately at first sight I didn't see a solution, or an explanation of my mistake. I didn't use "DistributeDefinitions", since in the documentation of "ParallelEvaluation" in Mathematica they say: "Parallel commands such as ParallelTable will automatically distribute the values and functions needed, using effectively DistributeDefinitions." But even if I add some "DistributeDefinitions[Euler]", it doesn't seem to change anything... – Fisher Jul 6 '13 at 23:07
you have several Table commands. if you just exchange them with ParallelTabe you will loose speed, since a parallelisation is really expensive and since you don't specify how the work-load shall be distributed, they work like a regular Table plus the additional time it takes to copy the definitions to all the available kernels on your machine. – Stefan Jul 6 '13 at 23:18

The problem you have and what I've shown very often here is, that you can't assume, just that you're using a parallelised version of a function, to gain speed at least double.

The problem with parallelised algorithms is always that you've to package your work-load (the stuff you want to calculate on) on the available processors (cores).

If you start to go parallel you have to start to measure and to find out how to best package your work-load in such sub-units, as that your processors are working on them and don't mutual prevent themselves.

Something that you don't want is that processor A needs to wait on some data, which processor B is actually working on.

It is very often the case, that people expect through the prefix Parallel a parallelisation which will make their algorithms faster and are astonished that it took even longer than before.

You have to keep in mind, that all the definitions are needed to be copied to all the available kernels etc. this is really expensive. So in order to gain something from parallelisation you need to rethink your data layout and how it should be distributed amongst the available kernels.

Every parallel algorithm that I've seen so far does look very different from its serial counterpart. The beautiful thing about Mathematica is that you could leave most of your algorithm parts untouched, but you need to take care about how your data is distributed amongst the kernels.

You could see the parallelisation of Mathematica as a data driven parallelisation. This is because most of Mathematica's internals aren't thread-safe, that's why data gets distributed.

Edit 1: Gaining speed

1. Speeding up Euler:

First thing I've done is running a Profile on your implementation:

It comes at no surprise (reading your code), that the function Euler is called very often. So the first thing that came to my mind was: "if that function is eating up all the processor cycles, why not try to reimplement it to gain some speed?"

The function Euler is rather simple and there is no point to parallelise it, but why don't Compile it?

cEuler = Compile[{{x, _Real}, {w, _Real}, {h, _Real}},
Module[{sigma = 3.0}, x - x^3*h + sigma*w],
CompilationTarget -> "C",
CompilationOptions -> {
"ExpressionOptimization" -> True,
"InlineCompiledFunctions" -> True,
"InlineExternalDefinitions" -> True}
];


Let me say something about "InlineExternalDefinitions". You've declared some global variables like sigma. The problem with Compile is the following, if we leave the CompilationOption "InlineExternalDefinitions" undefined the code will be compiled as follows:

This means, that every time the function cEuler is called it'll call the main evaluator in order got get the value of sigma. If something like that happens, all the compilation does not make sense anymore and it is even slower compared to a non-compiled version.

If you define "InlineExternalDefinitions" the variable sigma will get inlined and no call to the main evaluator is needed anymore.

2. Creating a global Noise table:

Although Mathematica has quite different parameter passing semantics when compared to C/C++ it gives me allways a bad feeling if arrays/array indices are passed around. That's why I rewrote the Noise creation functionality.

MyNoise = {};

SetAttributes[createNoise, HoldFirst];
createNoise[noise_, h_, StepsMC_] := Module[{T = 10},
noise = ParallelTable[RandomReal[NormalDistribution[0, Sqrt[h]], T/h], {StepsMC}];
]


We'll call this function with an external list per reference, so that the list is filled up with the list of random numbers.

Actually the Noise creation is the only part of your code that is parallelisable, if you leave the code mostly as it is.

If you execute createNoise and observe the parallel kernel status at the same time, you'll realize that the random number creation is well distributed amongst the kernels, but they aren't really busy.

We should take notice, that we have to weigh the cost of parallelisation with the actual range of the data. It looks that way that a range of StepsMC (1000 in your case) is even not a big enough to really benefit from parallelisation.

But let's see, maybe later we'll increase StepsMC, so it is good to have it and there is no pain in executing it.

3. The rest of the code:

MCEuler[h_, StepsMC_] :=
Quiet[Check[
Total[Table[
Fold[cEuler[#1, #2, h] &, z, MyNoise[[i]]], {i, 1, StepsMC}]]/
StepsMC, Indeterminate]];

MCEulerMT[h_, StepsMC_] :=
Total[Table[
Catch[Fold[
If[Abs[cEuler[#1, #2, h]] >= R, Throw[0], cEuler[#1, #2, h]] &,
z, MyNoise[[i]]]], {i, 1, StepsMC}]]/StepsMC

Comparison[h_, StepsMC_] :=
createNoise[MyNoise, h, StepsMC];
Return[{StepsMC, h, MCEuler[h, StepsMC], MCEulerMT[h, StepsMC]}];]


Here is nothing to parallelize. Although Fold per se can get parallelised, but it is not if you use complex slot patterns. Again, I doubt that with the current range of StepsMC it worth's the effort to rewrite the algorithms to become parallelised on a large degree.

So now let's execute the code:

LaunchKernels[] (*let's save some time *)

Comparison[#, 10^3] & /@ {2^0, 2^-2, 2^-4, 2^-6} // AbsoluteTiming


This yields the following result:

{4.803885, {{1000, 1, Indeterminate, -563.98}, {1000, 1/4, Indeterminate, -2051.07}, {1000, 1/16, 0.0108648, 0.0108648}, {1000, 1/64, -0.0300972, -0.0300972}}}

Comparing that with your original implementation this is a speedup of roughly factor 2:

{10.374996, {...}

-
My (naive) imagination of parallelization looks like (e.g. in MCEuler): We have a List Noise of subnoises and kernel 1 performs Fold with subnoise[[1]], kernel 2 uses subnoise[[2]] and so on... But it seems that this idea doesn't match reality. On the basis of your comments I think the problem is the distribution of the Noise, perhaps it's copied, even if it's not necessary, but at the moment I suffer from a lack of experience in this topic. I would be deeply grateful, if you would help me with the parallelization, but I don't want you to spend too much time on it. – Fisher Jul 7 '13 at 0:01
you are a real gentleman but the copying is needed, since the kernels need to get the data... I don't know about your time zone, but here it is late at night and i'm not able to help you out for now. If you're fit please have a look on my latest answer about how to distribute data and consult please the documentation for ParallelTable... but i'll return and try to do my best to help you out there... ok? – Stefan Jul 7 '13 at 0:23
According to your profile we are in the same time zone, it's likely we could even switch the language. Update: Well, I made some attempts in the past few hours. First of all I used "CoarsestGrained" as Method in ParallelTable, then I tried to structure the program in a different way, generating the noise in the MCEuler-method, so that it doesn't have to be copied to the kernels (at least that's my hope), but all ideas failed. I would suspect that the kernels interact using the random number generator, but I have no idea how to fix it. – Fisher Jul 7 '13 at 14:21
@Fisher passt das? ;) – Stefan Jul 8 '13 at 23:53
Ich schaue mir das heute genau an, und werde dann nochmal antworten. Auf alle Fälle schonmal vielen Dank. :) – Fisher Jul 9 '13 at 10:08