It is not very difficult to face a function for which ContourPlot works too slow. And it seems natural that this function can be parallelized well. Anyway, naive Parallelize@ContourPlot produces "ContourPlot[...] cannot be parallelized; "...

So, is it possible to parallelize ContourPlot?

  • 5
    $\begingroup$ Have you tried generating a discrete mesh and using ListContourPlot instead? The mesh does not have to be regular. $\endgroup$
    – Verbeia
    Jan 30, 2012 at 20:53
  • $\begingroup$ @Verbeia: actually I have not. You suggest generate points using ParallelMap and then plot them? $\endgroup$
    – faleichik
    Jan 30, 2012 at 20:59
  • $\begingroup$ Yes, that is what I had in mind. I wasn't sure if that was worth posting as an answer. $\endgroup$
    – Verbeia
    Jan 30, 2012 at 21:01
  • 1
    $\begingroup$ @verbeia I assume that would work, but it would mean you work on a regular grid. It doesn't use adaptive gridding that ContourPlot applies to increase sample density in rapidly changing spots. So the result can be rather different. $\endgroup$ Jan 30, 2012 at 21:02
  • $\begingroup$ @Sjoerd one could use Szabolcs's approach here: mathematica.stackexchange.com/questions/216/… $\endgroup$
    – acl
    Jan 30, 2012 at 21:04

5 Answers 5


I second @Verbeia's suggestion: compute the function on a mesh of points and use ListContourPlot. The disadvantage is that ListContourPlot has no adaptive sampling, so it'd be preferable if we could do our own adaptive sampling somehow. Adaptive sampling can give you a much better result while needing to compute the function in far less points---and the problem here is indeed computation time. So ContourPlot with its adaptive sampling might give a better result in less time on a single CPU than ListContourPlot will with a high resolution mesh computed on many CPUs.

Adaptive sampling is what I asked about (and solved) here: Adaptive sampling for slow to compute functions in 2D

The method I implemented there is usable (I am using it for something very similar to what you describe) but it is not nearly as good as ContourPlot's own. So one might still try to somehow make use of it. I'm quoting one suggestion I received from Leonid Shifrin there (in a comment):

You probably can control the DensityPlot, although not directly. Since it calls your function, you can simply Sow the values until some criteria (which you define) is violated (or satisfied). Then, you stop via throwing an exception, and catching it in the outer function, but still inside Reap. Alternatively, you could just start fooling DensityPlot by supplying faked values (perhaps, interpolated, or whatever), and it will stop by itself, I guess. Not sure this will work for you, but it may be worth trying.

I have not tried to implement this before, but I think it could work if your function is sufficiently smooth (which mine is definitely not, but yours may be).

Here's a quick sample implementation of how it could work:

First, let's define a sample function to plot:

fun[{x_, y_}] := 1/(1 + Exp[10 (Norm[{x, y}] - 3)])

Let's divide both the $x$ and $y$ axes into 5 parts on the interval $[0,5]$ and generate a mesh of points:

initialDivision = Range[0, 5];

points = N@Tuples[initialDivision, {2}];

Calculate function values on the intial mesh. This can be parallelized (just use ParallelMap)

values = fun /@ points;

This counter i will be used to control the maximal subdivisions in ContourPlot:

i = 0;

Now put the following code into a single cell, and evaluate it several times. Each time a finer and finer approximation will be computed. The points where function values have been computed will also be visualized. Note that I fixed the plot points in ContourPlot to force it to use the same initial mesh that I used, and I also fixed the number of contours.

if = Interpolation@ArrayFlatten[{{points, List /@ values}}]

{plot, {newpoints}} = Reap[
   ContourPlot[if[x, y], {x, 0, 5}, {y, 0, 5}, 
    Contours -> Range[0, 1, .1], MaxRecursion -> (++i), 
    PlotPoints -> Length[initialDivision], 
    EvaluationMonitor :> Sow[{x, y}]]

newpoints = Complement[newpoints, points];
newvalues = fun /@ newpoints;  (* <-- this can be parallelized *)
points = Join[points, newpoints];
values = Join[values, newvalues];


After a few iterations the contour plot and the point mesh will look like this (note that the code above only plots the contours for the previous step, not the current results):

Mathematica graphics Mathematica graphics

After 3 iterations, this method has computed the function value in 3809 points for this particular function.

Let's compare this with a plain ContourPlot using the same parameters:

ContourPlot[fun[{x, y}], {x, 0, 5}, {y, 0, 5}, 
    PlotPoints -> 6, MaxRecursion -> 3]

Mathematica graphics

The quality of the plot is about the same with a plain ContourPlot as well.

How many points did the plain CountoutPlot use?

Reap[ContourPlot[fun[{x, y}], {x, 0, 5}, {y, 0, 5}, PlotPoints -> 6, 
    MaxRecursion -> 3, EvaluationMonitor :> Sow[{x, y}]]][[2, 1]] // Length

(* ==> 3790 *)

It uses almost the same number of points, so if the bottleneck is computing f, the method I described is going to be almost as fast as ContourPlot on a single core, with the advantage that it is parallelizable for multiple cores.

The next step would be packaging this up into a self-contained function, but seeing how the quality improves step by step is also valuable as you can make decisions about when to stop calculating (and avoid excessive computation times).

I find it quite disappointing that all those nice and fast algorithms that plotting functions use (fast Voronoi cells, Delaunay trinagulation, adaptive sampling) are not directly accessible by users. We either have to use hacks to access these algorithms or reimplement them.

  • $\begingroup$ Very nice! I thought it could work, but now you've implemented this, which is much better than just an idea. +1. $\endgroup$ Jan 30, 2012 at 22:21
  • $\begingroup$ I liked the idea. I'll try to make it neater and post it in a little while $\endgroup$
    – Rojo
    Jan 31, 2012 at 1:30
  • 1
    $\begingroup$ The first time I read this, I missed the fact that you were pulling the point info from ContourPlot instead of calculating it on your own. Clever. $\endgroup$
    – rcollyer
    Jan 31, 2012 at 16:51

One idea, which was originally appealing to me, was to simply split the plot domain in four quarter size pieces and have these pieces calculated in parallel using ParallelTable and then combined using Show. It appears though that the overhead of getting the graphics data back to Show is pretty large, so this only yields some extra speed if the function is computationally slow.

a = 36;
(parPlot = Show[
       Mod[Sqrt[x^2 + y^2] - 7/2 ArcTan[x, y] +Sin[x] +Cos[y], \[Pi]] < \[Pi]/2, 
       {x, -a + i a, a + i a}, {y, -a + i a, a + i a}, 
       PlotRange -> {{-a, a}, {-a, a}}, PlotPoints -> 50], 
       {i, 0, 1}, {j, 0, 1}]
     ]); // AbsoluteTiming

{25.5344605, Null}


Mathematica graphics

The original serial calculation:

a = 36;
g = RegionPlot[
    Mod[Sqrt[x^2 + y^2] - 7/2 ArcTan[x, y] + Sin[x] + 
       Cos[y], \[Pi]] < \[Pi]/2, {x, -a, a}, {y, -a, a}, 
    PlotPoints -> 100]; // AbsoluteTiming

==> {51.2839332, Null}


Mathematica graphics

It's twice as slow, but it looks slightly better in places. RegionPlot is worse than ContourPlot in this sense.


Try generating a mesh of points using ParallelMap and then plot them using ListContourPlot.

You might find the techniques in this question on adaptive sampling for hard-to-compute functions useful.


The Presentations Application (which I am the author of) has a few routines to aid in the parallel processing of graphics.

Here is a serial contour plot of a complex function. (Presentations has routines for drawing complex functions directly with complex variables.)


Module[{f = z \[Function] Sin[z/2], z, zmin = 0, zmax = 4 (1 + I)},
   serial0 =
     {ComplexCartesianContour[f[z], {z, zmin, zmax}, Abs,
       Contours -> Range[0, 4, 0.1],
       ColorFunctionScaling -> False,
       ColorFunction -> (ColorData["Rainbow"][Rescale[#, {0, 4}]] &),
       PlotPoints -> 5,
       MaxRecursion -> 4]},
     Frame -> True,
     ImageSize -> 300]
   ]; // AbsoluteTiming

{0.298017, Null}

enter image description here

The following statement parallel loads all the Presentations package routines to the processors.


The following routine will subdivide the complex plane into rectangular regions.

ComplexParallelPartitionDomain[{var, nRe, nIm, start, end, overlap}] will generate a list of nRe x nIm complex subiterators in the form {{var, start1, end1},...} for partitioning complex drawing object domains. The parameter overlap has a default value of 0 but otherwise may be used to overlap each subdomain except on the right and top boundaries.

Here, since I have six processors, we will use

ComplexParallelPartitionDomain[{z, 2, 3, 0, 4 (1 + I), 0}]
{{z, 0, 2 + (4 I)/3}, {z, (4 I)/3, 2 + (8 I)/3},
 {z, (8 I)/3, 2 + 4 I}, {z, 2, 4 + (4 I)/3},
 {z, 2 + (4 I)/3, 4 + (8 I)/3}, {z, 2 + (8 I)/3, 4 + 4 I}}

I find that a convenient way to do parallel processing is to use the ParallelSubmit and WaitAll mechanism. Here is the code for the parallel processing version.

Module[{f = z \[Function] Sin[z/2], z, zmin = 0, zmax = 4 (1 + I)},
   par1 =
          ComplexCartesianContour[f[z], #, Abs,
           Contours -> Range[0, 4, 0.1],
           ColorFunctionScaling -> False,

           ColorFunction -> (ColorData["Rainbow"][
               Rescale[#, {0, 3.7}]] &),
           PlotPoints -> 6,
           MaxRecursion -> 3]] & /@ 
        ComplexParallelPartitionDomain[{z, 2, 3, 0, 4 (1 + I), 0}]
     Frame -> True,
     ImageSize -> 300]
   ]; // AbsoluteTiming
{0.101006, Null}

enter image description here

Next I will try Sjoerd's example above.

region = {x, y} \[Function] (Mod[Sqrt[x^2 + y^2] - 7/2 ArcTan[x, y] + Sin[x] + 
      Cos[y], \[Pi]] < \[Pi]/2)

(serial1 = RegionPlot[region[x, y], {x, -35, 35}, {y, -35, 35},
     PlotPoints -> 20,
     MaxRecursion -> 5]); // AbsoluteTiming
{44.546548, Null}

enter image description here

There is a similar partitioning function for x-y domains.

ParallelPartitionDomain[{x, 2, -35, 35, 0}, {y, 3, -35, 35, 0}]
{{{x, -35, 0}, {y, -35, -(35/3)}}, {{x, -35, 0}, {y, -(35/3), 35/3}},
 {{x, -35, 0}, {y, 35/3, 35}}, {{x, 0, 35}, {y, -35, -(35/3)}},
 {{x, 0, 35}, {y, -(35/3), 35/3}}, {{x, 0, 35}, {y, 35/3, 35}}}

There is one problem with partitioning RegionPlots in this way. If we want the automatic outlining of the region, then with partitioning we obtain some boundary lines on the edges of the partition. These are also on the edges of the serial plot but much less noticable. Drawing the boundary lines definitely improves the image. One way to handle this is to use a darker color for the region and the same color for the boundary lines.

par2 =
         RegionDraw[region[x, y], Evaluate[Sequence @@ #], 
          MaxRecursion -> 3, PlotPoints -> 20,
          PlotStyle -> Darker@Orange,
          BoundaryStyle -> Darker@Orange]] & /@ 
       ParallelPartitionDomain[{x, 2, -35, 35, 0}, {y, 3, -35, 35, 0}]
    Frame -> True,
    ImageSize -> 350]; // AbsoluteTiming
{2.196126, Null}

enter image description here

Notice that we obtained a speedup of a factor of 20, despite using only six processors, and the Parallel Kernel Status showing a speedup of 5.45. This is because we were able to use much less recursion to obtain the same quality of plot. These comparisions are a bit tricky because we have to compare the two plots by eye and we have two variables, PlotPoints and MaxRecursion to adjust. Nevertheless, a smaller piece of a plot will generally be topologically simpler and we can often achieve speedups that exceed the number of processors.


The following version is a mixture between Sjoerd's answer and Szabolcs answer. I'm using the CountorPlot example, and not the RegionPlot example, since I believe there's a bug on the EvaluationMonitor option of the RegionPlot (and also, because there's no such thing as a ListRegionPlot).

First, I start by making the function a little slower, just to be on a "real" scenario.

fun[{x_, y_}] := (Pause[0.001]; 1/(1 + Exp[10 (Norm[{x, y}] - 3)]))

The method is straightforward (and if someone is willing to make my code simpler, please go ahead). We just divide the 2D space.

Let's start by the serial case for speed benchmarking:

AbsoluteTiming@ContourPlot[fun[{x, y}], {x, 0, 5}, {y, 0, 5},
 Contours -> Range[0, 1, .1], ColorFunctionScaling -> False, 
 MaxRecursion -> 3, PlotPoints -> {6, 6}, Frame -> False, 
 ImageSize -> 200, ClippingStyle -> Automatic

enter image description here

Now, we parallelize it:

coreMesh = {10, 10};
xRange = {0, 5};
yRange = {0, 5};
totalPlotPoints = {6, 6};
GraphicsGrid[allPlots = ParallelTable[
  fun[{x, y}],
   {x, (xI - 1) (xRange[[2]] - xRange[[1]])/coreMesh[[1]],
    xI (xRange[[2]] - xRange[[1]])/coreMesh[[1]]},
   {y, (yI - 1) (yRange[[2]] - yRange[[1]])/coreMesh[[2]],
    yI (yRange[[2]] - yRange[[1]])/coreMesh[[2]]}, 
 Contours -> Range[0, 1, .1], ColorFunctionScaling -> False, 
 MaxRecursion -> 3, 
 PlotPoints -> {Max[2, Floor[totalPlotPoints[[1]]/coreMesh[[1]]]], 
  Max[2, Floor[totalPlotPoints[[2]]/coreMesh[[2]]]]}, 
 Frame -> False, ImageSize -> Floor[200/coreMesh[[1]]], 
 ClippingStyle -> Automatic], {yI, coreMesh[[2]], 1, -1}, {xI, 1, 
 coreMesh[[1]]}], Spacings -> 0]

enter image description here

Show[allPlots, PlotRange -> All, ImageSize -> 200]

enter image description here

I get similar artefacts as the ones identified by Sjoerd for the RegionPlot.

On this speed testing case, I'm going to divide on a coarser mesh (and I slightly improved the code readability).

coreMesh = {4, 4};
xRange = {0, 5};
yRange = {0, 5};
{stepX, stepY} = {(xRange[[2]] - xRange[[1]])/coreMesh[[1]],
 (yRange[[2]] - yRange[[1]])/coreMesh[[2]]};
totalPlotPoints = {6, 6};
{time, allPoints} =AbsoluteTiming@Flatten[ParallelTable[Reap[
 ContourPlot[z = fun[{x, y}],
  {x, (xI - 1) stepX, xI*stepX}, {y, (yI - 1) stepY, yI*stepY}, 
  Contours -> Range[0, 1, .1], MaxRecursion -> 3, 
  PlotPoints -> {Max[2,Floor[totalPlotPoints[[1]]/coreMesh[[1]]]], 
    Max[2, Floor[totalPlotPoints[[2]]/coreMesh[[2]]]]}, 
  EvaluationMonitor :> Sow[{x, y, z}]]][[2, 1]],
{yI, coreMesh[[2]], 1, -1}, {xI, 1, coreMesh[[1]]}], 2];

{time, Length[allPoints]}
(*{1.82858, 4472}*)

A little more points, and yet, less than half the time (on 8 threads...). I got the best performance on a 3 by 3 mesh, but since the amount of points was smaller, the comparison wouldn't be honest.

Graphics[Point[allPoints[[All, {1, 2}]]], ImageSize -> Small]

enter image description here

Putting all points together:

AbsoluteTiming@ListContourPlot[allPoints, Contours -> Range[0, 1, .1], 
 ColorFunctionScaling -> False, Frame -> False, ImageSize -> 200

enter image description here

So, on total, 2.1 seconds, instead of 5.4 seconds. Not huge... but this method is very susceptible to non optimized work balances, since all the complexity may end up being on just one kernel (which, on this example, is not very far from the truth).


Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.