interp2 returns NaN when given points are within mesh boundaries

I'm writing a scrip to execute gradient descent on a function. I'm generating 150 random starting points, within the bounds of the mesh, to execute the algorithm with. However, when I run this, at some point I get a NaN - I'm pretty sure this comes from interp2. Since all the starting points I generate are within the x,y boundaries, I am very confused as to why. What do you think?

Here's the script:

fun = @(x) cos(x(1)/2)*cos(x(2)) + x(2)/10 - x(1)/5;
sol = fmincon(fun,[0,0],[],[],[],[],[-9,-8],[9,8])
fun(sol)

% Define x and y range
x = -9:0.2:9;
y = -8:0.2:8;

% Compute z values
[X, Y] = meshgrid(x, y);
%Z = comp_z(X, Y);
Z = cos(X/2).*cos(Y) + (Y/10) - (X/5);

% Plot z using imagesc
imagesc(x,y,Z);
%Set the colormap to "hot"
colormap hot;
    
% Add axis labels and a title
xlabel('x');
ylabel('y');
zlabel('z(x,y)');
title('z(x,y) = cos(x/2)cos(y) + (y/10) - (x/5)');
hold on;


for i = 1:150
    disp(i);
    % Set starting points and parameters for gradient descent
    x0_1 = -8 + (rand*16);
    y0_1 = -7 + (rand*14);
    gamma = 0.2;
    epsilon = 0.002;
  
    % Perform gradient descent starting at (x0_1, y0_1)
    descent_matrix = gradient_descent(x, y, Z, x0_1, y0_1, gamma, epsilon);
   
    % Plot end points of gradient decents 
    line(descent_matrix(end,1), descent_matrix(end,2),'Marker','o','MarkerEdgeColor','cyan')
end

And here's the function:

function [descent_matrix] = gradient_descent(x, y, z, x0, y0, alpha, epsilon)

  % Compute gradient of z
  [dzdx, dzdy] = gradient(z); 

  % Compute z at starting point using cubic interpolation
  z0 = interp2(x, y, z, x0, y0, '*cubic');

  % Compute gradient at starting point using cubic interpolation
  gradz0 = [interp2(x, y, dzdx, x0, y0, 'cubic'), interp2(x, y, dzdy, x0, y0, 'cubic')]; 

  % Store starting point in positions matrix
  descent_matrix = [x0, y0, z0]; 
  
  % Initialize steps counter
  steps = 0; 
  
  % Perform gradient descent until step size is sufficiently small
  while true
      % Update x position
      x0 = x0 - alpha .* gradz0(1);

      % Update y position
      y0 = y0 - alpha .* gradz0(2);

      zi = interp2(x, y, z, x0, y0, 'cubic');
     

      % Compute gradient at starting point using cubic interpolation
      gradz0 = [interp2(x, y, dzdx, x0, y0, 'cubic'), interp2(x, y, dzdy, x0, y0, 'cubic')]; 

      % Store starting point in positions matrix
      
      descent_matrix = [descent_matrix; x0, y0, zi]; 
      
      % Increment steps counter
      steps = steps + 1; 
      
      % Check if step size is sufficiently small      
    
      if sqrt((x0 - descent_matrix(end-1,1)).^2 + (y0 - descent_matrix(end-1,2)).^2) < epsilon     
          % Exit loop
          break; 
      end
  end
end

[Torsten]

Looking at this update

x0 = x0 - alpha .* gradz0(1);
% Update y position
y0 = y0 - alpha .* gradz0(2);

I wonder how you guarantee that x0 and y0 remain in the initial bounds

x = -9:0.2:9;
y = -8:0.2:8;

And I ask myself why you apply gradient descent on a matrix of values for z and not on the function z(x,y) itself.

Answer 1

[Naeimeh N]

1- If all the the starting points you are generating are within the range of the x and y arrays, then make sure that you have sufficient data points in the x, y, and z arrays to accurately interpolate the values.

2- Another possible issue could be that the gradient descent algorithm is getting stuck at a local minimum. You may want to try adjusting the step size (alpha) or the convergence criteria (epsilon) to see if this helps.