Two dimensional interpolation polynomial

Hello everyone.

I have a file which size is 681*441, it contains values measured at 300321 points. I want to find the interpolation polynomial. I have looked for functions which return what I am seeking; but I have found none. griddedinterpolant does not return the polynomial in question, but only values in a more precise grid. On the other hand, the function polyfit does provide with the nodes of a polynomial of any desired order, but it seems to work only in one dimension.

Can someone please tell me if there is a way to do as in the function polyfit but returning a two dimensional polynomial?

Thanks.

Jaime.

[David K.]

I do not have any certainties on how to do this well, but if no one else has a way to do this I have some ideas. First, fitting a polynomial is at it's core an optimization problem so you could probably create your own 2d fitting with some of Matlab's optimization with fmincon. Where it is searching for the coefficients of the 2d polynomial and the cost that it is minimizing is the error off the real set of data. You would need to decide what degree of polynomial you want to search for. If you want to attach the actual file you have I could probably check it out and see how it goes.

[Jaime De La Mota Sanchis]

Thanks for the offer; this is the file in question. In it, I have generated a function of X and Y and I have it's value in a rectangle.

I hope it is clear; wind is the matrix with the values that I want to interpolate.

close all
clear
clc

x = 0:440;
y = 0:680;

[X,Y] = meshgrid(x,y);
wind=X+Y.^2;

surf(X,Y,wind,'EdgeColor','none','LineStyle','none','FaceLighting','phong')

[David K.]

Hm, actually I do not know if I understand your question. I thought you had data and was looking for the polynomial equation that fit it, but you already have said equation. Is this just a test example and you plan to apply to actual data to find the equations for?

You are also mentioning interpolation and I am not sure why gridded interpolant ot interp2 would not work with what you need it for.

[Jaime De La Mota Sanchis]

Yes, it is a test example; I know that the data is generated from the polynomial X+Y.^2, therefore the result also has to look similar. I am doing this to make sure that the code works before going to my data, where I have no means to check the results.

About the functions that you mention, they do not return the coefficients of the polynomial, they only provide interpolated data as far as I understand. I constructed this code before posting, but as said, the result is a more precised grid, not the polynomial coefficients

close all
clear
clc

[X,Y] = ndgrid(1:10,1:10);
V = X.^2 + 3*(Y).^2;
F = griddedInterpolant(X,Y,V,'cubic');
[Xq,Yq] = ndgrid(1:0.5:10,1:0.5:10);
Vq = F(Xq,Yq);
mesh(Xq,Yq,Vq)
figure
mesh(X, Y, V)

Answer 1

[David K.]

First off, I looked a little bit on wikis to find https://en.wikipedia.org/wiki/Multivariate_interpolation which might have what you want but I am not sure. However, the way I wanted to do it is with fmincon. I attached the objective function I created. Basically it creates a function of the form

a1x^2+a2y^2+a3xy+a4x+a5x+a6

then it loads in the X,Y, and wind data that you have and determines how closely it fits by a basic difference and summing.

This is utilized by fmincon as such:

a = fmincon(@objfun, rand(6,1), eye(6),ones(6,1).*20)

I decided to make the initial conditions random just because. Then, because this is a constrained minimum search, I just made it so that no value could go over 20, it would also likely work with large constraints as well but I doubt they will come into affect much. I tried using other optimization algorithms that did not require constraints but they did not work as well. For the example you gave the output should be

a = [0 1 0 1 0 0]

It generally gets pretty close, sometimes the constant value messes up though.

[Jaime De La Mota Sanchis]

Thanks again for your answer and dedication.

I am however lost, since there seems to be an error that I cannot solve.

I have looked inside of objfun and defined testData.mat as

close all
clear
clc
x = 0:440;
y = 0:680;
[X,Y] = meshgrid(x,y);
wind=X+Y.^2;
%surf(X,Y,wind,'EdgeColor','none','LineStyle','none','FaceLighting','phong')
save('testData', 'wind')

And now, I call a = fmincon(@objfun, rand(6,1), eye(6),ones(6,1).*20) from the comand window; however, there is an error which reads as:

Reference to non-existent field 'X'.

Error in objfun (line 8)

currGuess = currPoly(data.X,data.Y); % eval the coefficients

Error in fmincon (line 535)

initVals.f = feval(funfcn{3},X,varargin{:});

Caused by:

Failure in initial objective function evaluation. FMINCON cannot continue.

This is extrange since the variables X and Y are already defined.

I am also attaching a screenshot from my computer.

[David K.]

I am loading in testData.mat so if you change your save to be

save('testData.mat', 'wind','X','Y')

it should work fine.

Or, we can remove saving in general to make it a little nicer. Change the function as seen

function output = objfun(coeffs,wind,X,Y)

    % create the function ax^2+by^2+cxy+dx+ey+f
    currPoly = @(x,y) coeffs(1).*x.^2 + coeffs(2).*y.^2 + coeffs(3).*x.*y + coeffs(4).*x + coeffs(5).*y + coeffs(6);
    
%     data = load('testData.mat');
    
    currGuess = currPoly(X,Y); % eval the coefficients
    
    error = abs(wind-currGuess);
    output = sum(error,[1,2]); % sum the error

Then change the way you are calling it as such:

a = fmincon(@(x) objfun(x,wind,X,Y), rand(6,1), eye(6),ones(6,1).*20)

This will actually also be faster since calling a load function multiple times is pretty slow.

[Jaime De La Mota Sanchis]

I am very sorry, but the code still does not work. The error is

Error in objfun (line 10)

output = sum(error,[1,2]); % sum the error

Error in @(x)objfun(x,wind,X,Y)

Error in fmincon (line 535)

initVals.f = feval(funfcn{3},X,varargin{:});

Caused by:

Failure in initial objective function evaluation. FMINCON cannot continue.

I have tested in two diferent computers and I still get the same error.

I really don't understand what is wrong. I generate x, y, X, Y and wind as I have done before (see the attached file) and thencalling fmincon from the comand window.

The sequence is as follows:

first, I run grid_generator_2.m and then I call fmincon from the comand window as

>> a = fmincon(@(x) objfun(x,wind,X,Y), rand(6,1), eye(6),ones(6,1).*20)

I can only think of one reason for this not to work. In the function objfun, there is the following:

output = objfun(coeffs,wind,X,Y), x takes the place of the coeff and x is defined as x = 0:440; Perhaps should I use a vector of zeros of length 6 instead?

[David K.]

That is so weird. I added the fmincon call to your grid gen function and ran it on mine. It looks like we have the exact same thing. Perhaps the second argument for sum doesnt work in your version? try just replacing it with sum(sum(error)) I guess.

Answer 2

[John D'Errico]

I don't think you realize that an interpolation polynomial that passes exctly through 300321 points will be impossible to evaluate in double precision arithmetic. In fact, it may require a precision that is on the order of many thousand of decimal digits to get any thing out if it. Possibly hundreds of thousands of digits will be required.

And then, expect complete crapola anyway. Why? Because the polynomial might actually predict the points, but give you meaningless garbage between the points!

Answer 3

[Sharon Dolberg]

I wrote a multidimensional polyval function, a bit different from the original polyval. May it assist the users:

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function [a,ErrV]=polyvalN(MatIn,TargetIn)
[n,m]=size(MatIn);
MatOut=sum(pagemtimes(reshape(MatIn',m,1,n),reshape(MatIn',1,m,n)),3);
TargetOut=sum(MatIn.*TargetIn)';
a=MatOut\TargetOut;
ErrV=MatIn*a-TargetIn;
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

This function takes in a matrix MatIn which already includes the approximated structure for a linear approximated function (polynom) which the coefficient (a) of the approximated function needs to be found according to the target values TargetIn

for example:

------------

We have 8 points on a box

X=[-1;1;-1;1;-1;1;-1;1];

Y=[-1;-1;1;1;-1;-1;1;1]*2;

Z=[-1;-1;-1;-1;1;1;1;1]*0.5;

The approximated function which we intend to find the coefficients for is:

MatIn=[ones(8,1), X , Y , Z , X.*Y , X.*Z , Y.*Z ];

The value in each point is given. We can use rand or any known values in each point. Just in order to check the function polyvalN, we'll use known a's

intentionally:

TargetIn= 7 + 6*X + 5*Y + 4*Z + 3*X.*Y + 2*X.*Z + 1*Y.*Z;

Now, is we pose MatIn,TargetIn to polyvalN , we expect to get a's

[7;6;5;...1]

-------------

Notes:1. MatIn should be nXm wherein n>=m

TargetIn should be nx1

2. The point represented by X,Y,Z are not necessarily be in a box

structure as in the example. Actually, they can also represent any

single-dimension/multidimensional location, as long as note 1 is kept, and that MatIn includes

the vectors in the approximation.