Parallelize plotting

Question

I'm using Mathematica to make a few 3d plots. However, the machine I'm using has 4 cores and plotting (evaluating the function at different points can be done pretty much independently, so I should be able to parallelize my plots (or even a single plot for that matter). How can I do this?

P.s.: If the question is not clearly posed or more information is needed, let me know.

Edit: @halirutan's answer corrects my mistake, which I've now struck off.

Answer 1

The general answer to your question is: This depends on several factors of your functions and your plot setting. To explain this let me elaborate a bit what happens if you use functions like Parallelize, ParallelTable, ParallelEvaluate or one of their friends. Basically, your main Mathematica kernel starts several sub-kernels which are used as slaves to do the dirty work. If you have a computation which can be broken down into several unrelated tasks and you use this kind of parallelisation, then mainly the following steps will take place

1. The main kernel breaks (with your help) the big task into several sub tasks.
2. The main kernel sends all required data and function definitions to the sub-kernels. Function definitions are often distributed with DistributeDefinitions before the computation, but Mathematica always tries to figure out itself, what definitions are required on the sub-kernels.
3. The computation is executed on each sub-kernel in parallel.
4. The result of the computation is send back to the main kernel.
5. The main kernel composes the final result by combining all sub-results, e.g. it makes a big list from several sub-lists.

This list is not complete, but it gives you a good overview of the heavy parts. If you look at the above list and imagine what happens when you plot something, you'll probably can tell instantly several things which would lead to a slow down instead of a speed up when you do things in parallel. Let me point out some of them for your specific plotting scenario

1. When you try to plot a big array with e.g. ListPlot3D or when you use a large InterpolatingFunction, then the main kernel has to send a lot of data to its slaves.
2. When the plotting, meaning the creation of the Graphics (or Graphics3D) object doesn't take long enough, then the overhead of sending data between the kernels will outbalance what you probably gain from the parallelisation.
3. When you create a very detailed graphics (e.g. with setting PlotPoints or MaxRecursion very high), the graphic may become so big in memory, that sending it back to the main kernel takes too long.

If we take all this into account, we can easily create examples for both, a speed up and a slow down. First, let's look at an example where the heavy part is really the evaluation of the function you try to plot:

LaunchKernels[];

(* Serial evaluation *)
Table[ Plot3D[Arg[BesselK[4,x+I y]],{x,-1,1},{y,-1,1}], {16}];//AbsoluteTiming

(* Parallel evaluation *)
ParallelTable[ Plot3D[Arg[BesselK[4,x+I y]],{x,-1,1},{y,-1,1}], {16}];//AbsoluteTiming

This gives 19.3 seconds for the serial and 4.5 seconds for the parallel evaluation. Next, let's check the pretty much worst example I could think of: Creating an Image3D, when the 3d data was created on the main memory. So what we do basically is sending data back and forth only

data = CellularAutomaton[{224, {2, {{2, 2, 2}, {2, 1, 2}, {2, 2, 2}}}, {1, 1}},
  {SeedRandom[1]; RandomInteger[1, {15, 15}], 0}, 100];

Table[Image3D[data], {16}]; // AbsoluteTiming
ParallelTable[Image3D[data], {16}]; // AbsoluteTiming

This gives 0.007 for the serial and 0.021 for the parallel execution.

As you see we can create both sides of the medal, but this shouldn't prevent you from giving it a try. One last thing of your question is not clear to me

let's say they are unrelated to each other (so I can't imagine making a table out of them)

If you know which graphics you like to create, you can always make a table of them, because you could simply put them in a List and iterate over it. In fact, if you have all your plot commands in a list, then Mathematica will evaluate each plot one after another for you. To make a specific example, here a list of randomly copied documentation examples. I included the $KernelID to make clear that all this is evaluated on your main kernel (0)

{{$KernelID, 
  Plot[{Sin[x], Sin[2 x], Sin[3 x]}, {x, 0, 2 Pi}, 
   PlotLegends -> "Expressions"]},
 {$KernelID, 
  ParametricPlot3D[{(3 + Cos[v]) Cos[u], (3 + Cos[v]) Sin[u], 
    Sin[v]}, {u, 0, 2 Pi}, {v, 0, 2 Pi}, 
   PlotStyle -> Specularity[White, 50], Mesh -> None]},
 {$KernelID, 
  RegionPlot3D[x y z < 1, {x, -5, 5}, {y, -5, 5}, {z, -5, 5}, 
   PlotStyle -> Directive[Yellow, Opacity[0.5]], Mesh -> None]},
 {$KernelID, 
  ParametricPlot3D[ {Sin[u] Sin[v] + 0.05 Cos[20 v], 
    Cos[u] Sin[v] + 0.05 Cos[20 u], 
    Cos[v]}, {u, -π, π}, {v, -π, π}, 
   MaxRecursion -> 4, PlotStyle -> {Orange, Specularity[White, 10]}, 
   Axes -> None, Mesh -> None]}}

Now you could use the exact same list and evaluate it in parallel and you will notice that the kernel id's are different, because each sublist is evaluated by a different slave kernel

Parallelize[{{$KernelID, 
   Plot[{Sin[x], Sin[2 x], Sin[3 x]}, {x, 0, 2 Pi}, 
    PlotLegends -> "Expressions"]},
  {$KernelID, 
   ParametricPlot3D[{(3 + Cos[v]) Cos[u], (3 + Cos[v]) Sin[u], 
     Sin[v]}, {u, 0, 2 Pi}, {v, 0, 2 Pi}, 
    PlotStyle -> Specularity[White, 50], Mesh -> None]},
  {$KernelID, 
   RegionPlot3D[x y z < 1, {x, -5, 5}, {y, -5, 5}, {z, -5, 5}, 
    PlotStyle -> Directive[Yellow, Opacity[0.5]], Mesh -> None]},
  {$KernelID, 
   ParametricPlot3D[ {Sin[u] Sin[v] + 0.05 Cos[20 v], 
     Cos[u] Sin[v] + 0.05 Cos[20 u], 
     Cos[v]}, {u, -π, π}, {v, -π, π}, 
    MaxRecursion -> 4, PlotStyle -> {Orange, Specularity[White, 10]}, 
    Axes -> None, Mesh -> None]}}]

---

Update regarding your comment

When you plot 16 copies of the Bessel function, I suppose that different copies run on different cores, right? What if I was plotting over (say) the domain {x,-10,10} and {y,-10,10} and wanted to reduce the computational time?

Yes to your first question, each plot is run by a different slave kernel. You have to imagine, that the slave kernels are fully equipped Mathematica kernels. Regarding your second question: I believe that this is possible, the questions is with what quality. I'm in no position to tell you exactly what happens inside the plot routines, but it is likely that plotting only small sub-regions of your graphics and putting them back together in the end will maybe do unexpected things if your function is not of good nature.

Some things which need further considerations are the Mesh, the BoundaryStyle and PlotPoints because they are set for the whole graphics and if you subdivide your plot they mess up the final display.

Let me try to construct an example where we can really compare non-parallel vs parallel evaluation of a single plot. We take the BesselK again and will use a fixed setting for PlotPoints and MaxRecursion. We do this because otherwise the smaller sub-graphics will be created with a lot more Polygons than the big serial graphic.

grserial = 
   Plot3D[Arg[BesselK[4, x + I y]], {x, -1, 1}, {y, -1, 1}, 
    PlotPoints -> {128, 128}, MaxRecursion -> 0, Mesh -> None, 
    BoundaryStyle -> None]; // AbsoluteTiming

This requires about 23 seconds on my MacBook. To subdivide the plot, let's subdivide the x-range into 8 equal intervals, create the plot for each interval and combine the results. For this I will use ParallelCombine which lets me specify how I want to combine all the sub-graphics into a big one. The first argument of ParallelCombine is the function which creates the plot for several sub-intervals. The second argument is a list of all the intervals we need to plot and the last argument is the combinator

n = 8;
grparallel = With[{rangeParts = Partition[Range[-1, 1, 2/n], 2, 1]},
    ParallelCombine[
     Function[xrange, 
        Plot3D[Arg[BesselK[4, x + I y]], 
         Evaluate[{x, First[xrange], Last[xrange]}], {y, -1, 1}, 
         Mesh -> None, BoundaryStyle -> None, 
         PlotPoints -> {128/n, 128}, MaxRecursion -> 0]] /@ # &, 
     rangeParts, Show[{##}, PlotRange -> All] &]
    ]; // AbsoluteTiming

This runs fully parallel in about 5 seconds here. A visual inspection shows that both graphics look similar

Let us check the number of 3D points which are used in each graphic

Length[Cases[#, {_?NumericQ, _?NumericQ, _?NumericQ}, Infinity]] & /@ {grserial, grparallel}

(* {66650, 65012} *)

You see the number of created 3d points is almost equal.