Adaptive sampling for slow to compute functions in 2D

Question

EDIT: Although I have posted an answer based on my current progress, this in incomplete. Please see the "open issues" section in the answer.

---

Most plotting functions in Mathematica adjust the sampling density dynamically, based on the slope of the function:

DensityPlot[1/(
 1 + Exp[10 (Sqrt[x^2 + y^2] - 3)]), {x, 0, 5}, {y, 0, 5}, 
 Mesh -> All, PlotRange -> All]

Unfortunately the internal algorithm that does this is not directly callable.

---

Question:

I have a slow to compute 2D function (takes up to 10-40 seconds for a single point even though it's in optimized C++ called through LibraryLink) that I need to sample. How can I sample it adaptively, in a convenient and controlled way?

Since the function is so slow to compute,

• I'd like to be able to take the existing points, and refine them more if needed (i.e. continue the computation using the existing results)

• I can't use DensityPlot directly because I can't control how many points it is going to compute, and I need an upper limit on that (i.e. an upper limit on computation time). It also can't be interrupted and continued later.

I am looking for

1. effective methods to do this

2. implementations in Mathematica (either as an answer or pointers to libraries)

The messy details:

I am trying to compute a phase diagram and map the boundaries between the phases precisely. So I don't need the function value everywhere, only where it very suddenly drops. The function is either of magnitude 1 (say, between $0.1 \div 1$), or very small (close to zero).

The function is computed using Monte Carlo methods, so at a small scale it doesn't appear smooth, and I might get inconsistent results close to the phase boundaries on subsequent runs of the function.

This should give you an idea of what sort of function I'd like to apply this to, which might be important when choosing a method.

Answer 1

Update: I described an alternative approach based on built in plotting functions in this answer. That approach is not very practical here though because I need to be able to handle points at arbitrary positions while built in functions work with a rectangle-based mesh. I am still looking for improvements.

---

I came up with this very naive approach and implementation (I know that the implementation is not optimal at all):

First let's define a test function (same one as in the question):

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

These functions will subdivide lines in the Delaunay triangulation of the points if 1. the points are further apart than a threshold (i.e. the resolution is controlled) and 2. the function values in the two points also differ by more than another threshold.

<< ComputationalGeometry`

makeLines[tri_] := Union[Sort /@ Flatten[Thread /@ tri, 1]]

subdivision[points_, values_, valueThreshold_, distanceThreshold_] :=
 Module[
  {tri, lines, linesToDivide},
  tri = DelaunayTriangulation[points];
  lines = makeLines[tri];
  linesToDivide = 
   Pick[lines, (Abs[values[[#1]] - values[[#2]]] > valueThreshold && 
        Norm[points[[#1]] - points[[#2]]] > distanceThreshold ) & @@@ lines];
  Mean /@ (linesToDivide /. n_Integer :> points[[n]])
  ]

Let's define an initial point grid to compute the function in:

points = Tuples[Range[0, 5, 1], 2];

We can iterate this function to add more and more points and recursively subdivide the grid (evaluate the following commands together repreatedly):

values = fun /@ N[points];
newpoints = subdivision[points, values, 0.1, 0.1];

ListDensityPlot[Flatten /@ Thread[{points, values}], 
 InterpolationOrder -> 0, Mesh -> All, ColorFunction -> "MintColors", 
 Epilog -> {PointSize[Large], Point[points], Red, Point[newpoints]}]

points = Join[points, newpoints];

The result after several iterations:

values = fun /@ N[points];
ListDensityPlot[Flatten /@ Thread[{points, values}], 
 InterpolationOrder -> 0, Mesh -> All, ColorFunction -> "MintColors"]

---

Open question: My aim is to minimize the number of points I need to compute while getting a precise approximation. This is probably not the best subdivision method for it. What are some easy-to-implement better methods?

I think ideally the decision for refining the grid should be made based on some sort of curvature. Take for example the following function:

ContourPlot[Erf[1/(1 + 20 x^2) - y], {x, -3, 3}, {y, -3, 3}]

Using a valueThreshold of 0.3 and distanceThreshold of 0.1, and a starting grid with a spacing of 0.5 produces this:

Let's turn on interpolation (because I can't turn interpolation off in DensityPlot) and compare it with a DensityPlot made using similar options (PlotPoints -> 12, MaxRecursion -> 15):

The curvature-based DensityPlot (right) is clearly much better. Furthermore, my method won't properly detect "fjord-like" structures (similar to the one in this example). It tend to "jump" over them, this is why some artefacts are visible in the middle of the plot.

---

Thanks to @ruebenko for the hints and ideas he sent me!

Answer 2

I have cooked up a routine that performs adaptive sampling in parallel for integration. With a little tweaking, it should work for multidimensional plotting routines such as this--I think they are related problems as "parallel adaptive sampling" which Mathematica uses in several areas such as plotting, polynomial fitting, and integration.

My use case is very expensive integrand. The idea is to trick NIntegrate into giving yielding where it would like to integrate, but not actually doing the integration right away, but rather deferring it to be done in parallel batch. These values are then memoized, and provided to a second invocation of NIntegrate, which merely looks up the (pre-computed) values that it needs.

(It was made in desperation and probably is a horrible implementation it but I hope the idea/solution are sound.)

ParallelNInt[f_, {xmin_, xmax_}, maxr_] := Module[
   {memoize, mtemp, res, fake, mem, xnew, known, minr, integrand},
   
   (* keep a list of known x values *)
   known = {};
   
   (* ensure numeric evaluation inside NIntegrate *)
   
   mem[x_?NumericQ] := memoize[x];
   fake[x_?NumericQ] := 0;
   
   (* increment MaxRecursion from 0 to maxr *)
   Do[
    
    (* First Invocation: "trick" NIntegrate 
    so that it either looks up values from previous
    lower recursive levels, or quickly yields a fake value *)
    
    (*
    Print["Recursion r = " <> ToString[r]];
    Print["---------"];
    *)
    
    xnew = Quiet[
      Reap[
        NIntegrate[
         If[MemberQ[known, x], mem[x], fake[x]],
         {x, xmin, xmax},
         Method -> {"LocalAdaptive", "SymbolicProcessing" -> 0},
         MinRecursion -> r,
         MaxRecursion -> r,
         EvaluationMonitor :> (
           (*
           Print[{"1st: ",x,If[MemberQ[known,x],mem[x],
           fake[x]],MemberQ[known,x]}];
           *)
           Sow[x])
         ]
        ][[2, 1]]
      ];
    
    (* return the previous Second Invocation result if no new data *)

        If[Length[xnew] == 0, Return[res]];
    
    (* for values that we have not seen, evaluate in parallel *)
    
    xnew = Complement[xnew, known];
    mtemp = ParallelMap[
      f[#] &,
      xnew
      , Method -> "FinestGrained"];
    
    (* memoize those values that were computed in parallel *)
    Do[
     memoize[xnew[[i]]] = mtemp[[i]]
     , {i, Length[mtemp]}];
    
    (* append those x-values to the known list *)
    
    known = Join[known, xnew];
    
    (* Second Invocation: NIntegrate with memoized values *)
    
    res = Quiet[
      NIntegrate[
       mem[x],
       {x, xmin, xmax},
       Method -> {"LocalAdaptive", "SymbolicProcessing" -> 0},
       MinRecursion -> r,
       MaxRecursion -> r
       (*,EvaluationMonitor:>Print[{"2nd: ", x,mem[x]}]*)
       ]
      ];
    
    (* increment up to maxr *)
    , {r, 0, maxr}];
   res
   ];

Define a slow function:

f[x_?NumericQ] := (Pause[5]; 1/x Cos[Log[x]/x]);

Integrate the normal way:

AbsoluteTiming[
 NIntegrate[f[x], {x, 0.1, 0.2}, 
  Method -> {"LocalAdaptive", "SymbolicProcessing" -> 0}, 
  MinRecursion -> 1, MaxRecursion -> 1]
 ]

{95.0225, -0.100505}

Integrate with this parallel function:

AbsoluteTiming[
 ParallelNInt[f, {0.1, 0.2}, 1]
 ]

{35.0655, -0.100505}

Addressing a few points in your question:

• I use "LocalAdaptive" because I believe that better-reuses data from previous recursions.
• One may use an integration technique to generate data and then use a List___Plot command. Or, adapt this code to make graphics instead of integrations. Anything that can Reap/Sow is amenable.
• Uses native Mathematica routines and no 3rd party packages (such as Cubature)
• Perhaps can be used with NIntegrate and Monte Carlo methods, I am not sure if they obey MaxRecursion or just how terribly sequential they are.
• Speed-up is merely through parallelization and it is easy to check against serial versions since this solution is supposed to be a order-of-operations difference tops.