Maximum of 2D data interpolation

Hi, I have 2D data for which I know there is a single maximum, but the points spacing may not be ideal so the maximum may lie in between. Therefore I would like to use griddedInterpolant to find the maximum point. X and Y contain the points, Z the corresponding function value.

It should look something as:

F = griddedInterpolant(x_mat,y_mat,z_mat);
xmaxs = arrayfun(@(xx,yy)fminsearch(@(x,y)-F(x,y),xx,yy),[x_mat(ii) y_mat(ii)]);

Where I provide an initial guess ii. But this is not the right way of calling it for a 2D function I'm afraid...

Answer 1

[John D'Errico]

You will NEED to use a spline based method of interpolation here, specifically spline or cubic. NOT makima, or pchip. Those interpolants will not generate functions with an extremum between the data points. And of course you certainly cannot use the default linear interpolant, as that will NEVER generate a solution that is not at a data point.

[X,Y] = ndgrid(0:0.5:6);
Z = cos(X+Y).*sin(X-Y);
G = griddedInterpolant(X,Y,Z,'spline');
fsurf(@(x,y) G(x,y),[0 6 0 6])

Warning: Function behaves unexpectedly on array inputs. To improve performance, properly vectorize your function to return an output with the same size and shape as the input arguments.

We can clearly see several peaks. Of course, since I used a trig function here, they will lie at transcendental locations, so never exactly at a grid point.

[xymin,fval] = fminsearch(@(xy) -G(xy(1),xy(2)),[1,1])

xymin = 1×2

0.7930 2.3606

fval = -0.9945

To within the tolerances used by fminsearch, this should be a solution.

xymin/pi

ans = 1×2

0.2524 0.7514

And that would clearly be the point at [pi/4,3*pi/4]. If you want a better solution you need a finer grid spacing, as the spline has limited approximative capabilities with that course grid.

[Albert Zurita]

Thanks John, that's very good. However, I don't seem to find an 'interpolated' maximum. The fval value that I get seems to be always equal to the maximum of the original function, not an interpolated maximum ?

[G_val_max,ixmax] = max(G_val(:));
G_int = griddedInterpolant(X,Y,G_val,'cubic');
[xymin,fval] = fminsearch(@(xy) -G_int(xy(1),xy(2)),[X(ixmax) Y(ixmax)]);

It seems G_val_max is equal to fval. I am however re-assessing this, as perhaps it's a condition of my data only...

Thanks!

[John D'Errico]

I did show the use of 'spline' there, did I not? Perhaps cubic might have an issue that I have not thought about. I'd need to look carefully at the implementation of the cubic interpolant to know, and it has been some time since I looked at the codes for the cubic option. But I will assert that spline should not have that problem. And I very explicitly used spline in my example, so you can see that it did work as desired.

[Albert Zurita]

I will investigate further, but I tested for a simple paraboloid and it seems to do things right, both for cubic and spline interpolations, as it should be. I will check what it's about my data then, perhaps there is some special pattern.

Brilliant answer, thank you!

[Albert Zurita]

I am thinking that given the fact that I'm looking at very precise maximum determination the problem is very likely related with the tolerance options to the fminsearch function. Changing the default value of 1e-4 to a lower value seems to provide better results!

Answer 2

[Torsten]

That doesn't make sense.

F is only an approximating function between the grid points - you don't know if the "real" underlying function behaves like this.

If you cannot supply the "real" function for F, better just choose the z with maximum value and the corresponding x and y.

[Albert Zurita]

Sorry. I rely on the grid points and the corresponding z values, so that f(x,y) = z is known and well sampled. However even if I can sample denser and do max(z) I would like to find (xi,yi) interplated points which would correspond to a better approximation/interpolation of the maximum. I found some answers in 1D and thought it could be extrapolated to 2D.