What simple methods are there for adaptively sampling a 2D function?

Question

I have a two-dimensional function $f(x,y)$ whose values I would like to sample. The function is very expensive to compute and it has a complex shape, so I need to find a way to get the most information about its shape using the least number of sample points.

What good methods are there to do this?

What I have so far

• I start from an existing set of points where I have already computed the function value (this could be a square lattice of points or something else).

• Then I compute a Delaunay triangulation of these points.

• If two neighbouring points in the Delaunay triangulation are far enough ($ > \Delta X$) and the function value differs sufficiently in them ($> \Delta f$), then I insert a new point midway inbetween them. I do this for each neighbouring point-pair.

What's wrong with this method?

Well, it works relatively well, but on functions similar to this one it's not ideal because sample points tend to "jump over" the ridge and not even notice it's there.

It produces results like this (if the resolution of the initial point grid is sufficiently rough):

This plot above shows the points where the function value is calculated (actually Voronoi cells around them).

This plot above shows the linear interpolation generated from the same points, and compares it to Mathematica's built-in sampling method (for about the same starting resolution).

How to improve it?

I think the main issue here is that my method decides whether to add a refinement point or not based on the gradient.

It would be better to take into account the curvature or at least the second derivative when adding refinement points.

Question

What is a very simple to implement way to take into account the second derivative or curvature when the locations of my points are not constrained at all? (I don't necessarily have a square lattice of starting points, this should ideally be general.)

Or what other simple ways are there to calculate the position of refinement points in an optimal way?

I am going to implement this in Mathematica, but this question is mainly about the method. For the "easy to implement" bit it does count that I'm using Mathematica though (i.e. this was easy to do so far because it has a package for doing Delaunay triangulation)

What practical problem I am applying this to

I am calculating a phase diagram. It has a complex shape. In one region its value is 0, in another region it's between 0 and 1. There's a sharp jump between the two regions (it's discontinuous). In the region where the function is greater than zero there is both some smooth variation and a couple of discontinuities.

The function value is calculated based on a Monte Carlo simulation, so occasionally an incorrect function value or noise is to be expected (this very is rare, but for a large number of points it happens, e.g. when the steady state is not reached to due some random factor)

_{I have asked this already on Mathematica.SE but I can't link to it because it's still in private beta. This question here is about the method, not the implementation.}

---

Reply to @suki

Is this the type of division you suggest, i.e. putting a new point in the middle of the triangles?

My concern here is that it seems to require special handling at the edges of the region, otherwise it will give very long and very thin triangles, as shown above. Did you correct for this?

UPDATE

A problem which appear both with the method I describe and with @suki's suggestion to put subdivide based on triangles and put the subdivision points inside the triangle is that when there are discontinuities (as in my problem), recomputing the Delaunay triangulation after a step may cause triangles to change and perhaps some big triangles to appear which have different function values in the three vertices.

Here are two examples:

The first shows the end result when sampling around a straight discontinuity. The second shows the sampling point distribution for a similar case.

What simple ways are there to avoid this? Currently I am simply subdividing those egdes that disappear after a retriangulation, but this feels like a hack and needs to be done with care as in the case of symmetrical meshes (like a square grid) there are several valid Delaunay triangulations, hence edges might change randomly after retriangulation.

Answer 1

I worked on a problem similar to this a while back.

I think the main difference between between our implementations is that I was choosing where to add points based on the triangles, not the edges. I also choose new points inside the triangles instead of on the edges.

I have the feeling that adding points inside the triangles would make it more efficient by giving a small increase in the average distance from old points to the new one.

Anyway another nice thing about using triangles instead of edges is that it gives an estimate of the gradient vector, instead of the slope along this particular edge.

In my matlab code I used a base class to take care of most of the machinery, with a few abstract methods:

• weight(self) to decided the priority for which triangles to subdivide next.
• choosePoints(self,npoints = "auto") to decide new points to evaluate based on the weight of each triangle.

I found this setup very flexible:

• setting a sub-class's weight() function to the area of the triangle produces a constant mesh density.
• setting weight() to calculate the average function value times the area of the triangle gives a sort of quasi-random probability sampling.
• using var(triangle.zs) could do, for functions that have binary output, what I feel is a generalization of a bisection search to more than 1 dimension.
• using area + var(triangle.zs) was pretty effective at putting a constant density everywhere, and an increased density along any slope (almost what you have now).

I used the variance of the z values, to approximate the importance of first order effects (slope) because the variance will never go to infinity like the slope can.

For the last example the background density was good because I was searching for discontinuous blobs of high value in a space of low value. So it would slowly fill in the entire mesh and when it would find a blob it would concentrate on following the edge of the blob all the way around because of the high weight I put on the gradient (and that it only filled in the top n triangles on each iteration). At the end I could know there were no (reasonably shaped) blobs (or holes in my blobs) of a size greater than the resulting background mesh density.

Like you I did get some bad points in my results, they weren't a problem for me because the error was such that if you re-ran nearby points they would probably give the correct answer. I'd just end up with blips of increased mesh density around my bad points.

Whatever you do I always recommend making the weights related to the triangle size so that, all else being equal, large triangles are broken up first.

Maybe a solution for you is to take my approach one step further and instead of evaluating triangles based on the contents of that triangular cell, evaluate based on that one and all three adjacent triangles.

That will contain enough information to get an estimate of the complete Hessian matrix. you can get it by doing a least squares fit of z = c1*x + C2*y c11*x^2+c12*x*y+c22*y^2 over all the vertexes in the triangles of interest (center the coordinate system on the triangle first).

I wouldn't use the gradient or Hessian (those constants) directly because they'll go to infinity at a discontinuity.

Maybe the sum-square-error of the z values relative to a planar approximation of those points would be a useful measure of how interesting second order effects are.

---

Updated:

That looks reasonable to me.

I never actually got around to special casing the edges. It bothered me a little but for what I was doing, it was enough to just start with a lot of points around the edges.

more elegant would be to combine our two approaches, weighting edges and triangles. Then if an edge is too long, cut it in half... I like the way that concept generalizes to higher dimensions (but the numbers get big fast)...

But since you don't expect the main body of the mesh to have high aspect ratio triangles, so you could use a function like Matlab's freeboundary function to find the boundary, then run the the same algorithm in one less dimension on the boundary. If done right, on a cube for example, you could get the same mesh density on the edges, on the faces, and inside the cube. Interesting.

One thing that I never found a good solution for was the fact that my version would never explore outside the convex hull of the initial point set.

Answer 2

I think the main problem in your heuristic is that you're considering the gradient in one dimension only and thus, in regions where dfdx is small but dfdy is large (as happens in the middle of your example), you will miss points when looking in the "wrong" dimension.

One quick fix would be to consider sets of four points, taking their center of gravity, and approximating |dfdx|+|dfdy| using those four points. Another alternative is to take three points (i.e. a triangle) and take the maximum gradient of the surface over them.