max value of poynomial surface fit

I'm trying to fit a 2D surface and extract the max value from the fit. The fit seems to be working:

x=[1 2 3 1 2 3 1 2 3]'
y=[1 1 1 2 2 2 3 3 3]'
z=[1 3 2 6 10 5 2 4 3]'
f=fit([x,y],z, 'poly22','Normalize','off')
figure
subplot(1,2,1);plot(f, [x,y],z)

So one way I thought I could find the max value was to use the fit at finer sampling, differentiate and find the smallest values in the differentials:

%Differentiate fit with higher x,y sampling:

[xx, yy] = meshgrid( 1:0.1:3, 1:0.1:3 )   %Finer sampling
[fx, fy] = differentiate( f, xx, yy )

then look at each fx, fy and find indicies of lowest values, then get the coordinates using these indices

X=0;
Y=0;
xm=min(abs(fx(:)))
ym=min(abs(fy(:)))
for i=1:numel(fx)
    for j=1:numel(fy)        
        if (fx(i)==xm)&&(fy(j)==ym)
            X=xx(i);
            Y=yy(j);
            disp('found')
            [fx(i) fy(j)]
        end    
    end        
end    

subplot(1,2,2); plot(fx,fy)
subplot(1,2,1); hold on; plot(X,Y,'ro')

But X & Y always come out as 0. I've notices that condition:

if (fx(i)==xm)&&(fy(j)==ym)

is never met and I can't see why. Thanks for any help

Answer 1

[Image Analyst]

Since you're doing it numerically anyway, why not simply use the max() function:

[maxOfF, linearIndexOfMax] = max(f(:));
[row, col] = ind2sub(size(f), linearIndexOfMax)

[Jason]

Hi I.A.

So I get this error message:

To evaluate an SFIT object you must specify either two arguments, or one argument with two columns.

[Torsten]

Use

[xx, yy] = meshgrid( 1:0.1:3, 1:0.1:3 );
ff=f(xx,yy);
[maxOfF, linearIndexOfMax] = max(ff(:));
[row, col] = ind2sub(size(ff), linearIndexOfMax)

Best wishes

Torsten.

[Jason]

Thanks, Im still not able to plot this max location, I have tried:

plot([xx(row),yy(col)],maxOfF,'r*','MarkerSize',10)

but xx(row) , yy(col) just come out as 1 & 2 respectively?

[Image Analyst]

If you want to plot a star in 3-D space, you need to use plot3(), not plot().

v = 1 : 0.1 : 3;
plot3(v(row), v(col), maxOfF, 'r*', 'MarkerSize', 10);

[Jason]

Thanks I.A

Answer 2

[John D'Errico]

Simple enough. Using your data, I'll use my polyfitn tools, found on the file exchange for download. That toolbox has a utility in it that can translate the model into a sym.

P = polyfitn([x,y],z,2)
P = 
  struct with fields:

        ModelTerms: [6×2 double]
      Coefficients: [-2.5 1.8821e-15 10.167 -4.5 18.5 -20.667]
      ParameterVar: [0.88889 0.44444 16.296 0.88889 16.296 29.432]
      ParameterStd: [0.94281 0.66667 4.0369 0.94281 4.0369 5.4251]
               DoF: 3
                 p: [0.076885 1 0.086294 0.017474 0.019509 0.0318]
                R2: 0.91111
        AdjustedR2: 0.76296
              RMSE: 0.7698
          VarNames: {'X1'  'X2'}

The same coefficients as you got, but in a different order.

Next, convert it to a symbolic polynomial.

Psym = polyn2sym(P)
Psym =
(61*X1)/6 + (37*X2)/2 + (4771727564350351*X1*X2)/2535301200456458802993406410752 - (5*X1^2)/2 - (9*X2^2)/2 - 62/3

Then differentiate with respect to each variable, and solve.

maxloc = solve(diff(Psym,X1) == 0, diff(Psym,X2) == 0)
maxloc = 
  struct with fields:

      X1: [1×1 sym]
      X2: [1×1 sym]

vpa(maxloc.X1)
ans =
2.0333333333333341070915836534138

vpa(maxloc.X2)
ans =
2.0555555555555559807740534792035

The value at that location is

vpa(subs(Psym,{X1,X2},[maxloc.X1,maxloc.X2]))
ans =
8.6833333333333411998755449208187

You should see that the max of the surface is lower than the max datapoint. if you plot the surface and the data (as you already did) you will see that one point stands quite a bit higher than the surface.

No need to do any search or approximation with meshgrid. I'm not sure if the curvefitting TB has a similar utility in it to convert to a sym. If not, I should probably write one someday.

[Jason]

Thanks for your answer, it seems a bit more complex than the other method, but I will try your suggestion anyway, thanks again.

[John D'Errico]

Because this gives you the correct answer, instead of an approximation?

[Florian Marco Lipp]

I converted my model (higher dimension) to a symbolic polynomial. It's not possible to solve it because Operator '==' is not supported for operands of type 'sympoly'. Is there another opportunity to calculate the maximum value of polynomial surface fit?

[John D'Errico]

No, of course you cannot use solve with the sympoly toolbox. I did not write a solve function. Nor does it use the same syntax. HOWEVER, if you read the code I wrote, I did show how to solve the problem.

In what I wrote, I converted the polynomial to a sym, NOT to a sympoly.

Psym = polyn2sym(P)

If you do something sompletely different you must expect it to fail.