interpn vs griddata, not producing the same results when interpolation is only in 2 dimensions

Hi, I am trying to understand why these two codes are not producing the same results.

in one I am doing interpolation in the third and fourth dimension and going over the arrays manually using for loops.

in the other I am doing the entire interpolation without for loops using interpn.

Because the default interpolation method is linear, I would expect that if the interpolated points fall exactly on the original grid it wouldnt use the information from neighbouring points and if so the two methods should produce the exact same thing, but interpolatedLF != interpolatedLFn.

%% First Code
[u_index, v_index, s_index, t_index] = ndgrid(1:size(LF,1), 1:size(LF,2), 1:size(LF,3), 1:size(LF,4));
 
shifted_s_index = min(s_max, s_index + d * u_index);
shifted_t_index = min(t_max, t_index + d * v_index);

for uu = 1:u_max
     for vv = 1:v_max
         for cc = 1:colors
            interpolatedLF(uu, vv, :, :, cc) = griddata(squeeze(s_index(uu,vv,:,:)), squeeze(t_index(uu,vv,:,:)), squeeze(LF(uu,vv,:,:,cc)), squeeze(shifted_s_index(uu,vv,:,:)), squeeze(shifted_t_index(uu,vv,:,:)));
        end
     end
end

%% second code
[u_index, v_index, s_index, t_index, color_index] = ndgrid(1:size(LF,1), 1:size(LF,2), 1:size(LF,3), 1:size(LF,4), 1:size(LF,5));
shifted_s_index = min(s_max, s_index + d * u_index);
shifted_t_index = min(t_max, t_index + d * v_index);
interpolatedLFn = interpn(u_index, v_index, s_index, t_index, color_index, LF, u_index, v_index, shifted_s_index, shifted_t_index, color_index);

Answer 1

[John D'Errico]

Are you testing for EXACT equality, which is something you should absolutely never do in any floating point environment? Instead, you should test the difference. Remember that any pair of computations, done in a a different sequence need not produce the same result, when that computation is done in floating point arithmetic.

For example, compare these results:

0.3 -0.1 -0.2

ans = -2.7756e-17

-0.1 -0.2 + 0.3

ans = -5.5511e-17

which surely must produce the same (zero) results. Yet they are not so, because the two computations, done in a different sequence, experience subtly different errors in the least significant bits. (If you look carefully at the binary expansions of those numbers and how they are stored, you can see where the subtle errors arise.)

Instead, you should test the DIFFERENCE between the two computations. This is the first thing you should look at, instead of asking if the two computations are identical.

[Roee Mazor]

Thanks, but i checked the actual difference, at times it was 0.01 and sometimes as big as 0.134. For values between 0 to 1 that is unacceptable as equal :)

[Bruno Luong]

Until you post an example so we can see, there is no point to discuss.

Actually I beilieve you make a mistake somewhere.

[John D'Errico]

Agreed - that difference is significant for numbers on the order of 1 or less. Therefore it is pretty certain there is a mistake in what you did. Those squeezes make me fear for what you have done though. And they make the code amazingly difficult to read and understand. If I had to make my next guess at the most probable source of this problem, it would be in the squeezes, and how the arrays got implicitly reshaped.

[Roee Mazor]

do you have a suggestion on how to do his the right way? permute to reorder the dimensions?

The code really doesnt do much so i assume the idea of what i was trying to do is clear, can you suggest the right way of doing that?

[Stephen23]

Use permute, it will reliably permute only the specified dimensions.

Avoid squeeze for exactly the same reason for avoiding awful length:

https://www.mathworks.com/matlabcentral/discussions/highlights/133276-new-in-r2020b-new-way-to-return-height-or-width-of-arrays/36288#reply_36288

[Roee Mazor]

Thank you, I will rewrite and update ASAP.

[Roee Mazor]

Btw, this is exactly what i am trying yo understand, why am i not supposed to be getting bit equivalence?

Answer 2

[Stephen23]

Your code seems to have some bugs:

• the first and second dimensions should NOT be included in the indices for griddata (i.e. they should NOT be inputs/outputs of ndgrid).
• the inputs to griddata must be vectors.

A simple example using random data shows that the outputs are exactly identical:

inp = rand(5,4,3,2); % fake input data
v3o = [1,3,2]; % output sample locations
v4o = [2,1,2]; % output sample locations
% GRIDDATA
[m3o,m4o] = ndgrid(v3o,v4o); % Exclude 1st and 2nd dimensions!
[m3i,m4i] = ndgrid(1:3,1:2); % Exclude 1st and 2nd dimensions!
one = nan(5,4,numel(v3o),numel(v4o));
for ii = 1:5 % d1
    for jj = 1:4 % d2
        tmp = inp(ii,jj,:); % all inputs are vectors:
        mat = griddata(m3i(:),m4i(:),tmp(:),m3o(:),m4o(:));
        one(ii,jj,:) = mat(:);
    end
end
% INTERPN:
[m1o,m2o,m3o,m4o] = ndgrid(1:5,1:4,v3o,v4o);
two = interpn(inp,m1o,m2o,m3o,m4o);
% Checking:
isequal(one,two)

ans = logical
   1

[Roee Mazor]

Hi, thanks for your answer!

I made small modifications to your code to better represent my original problem

inp = rand(5,4,3,2); % fake input data
% v3o = [1.1,2.5,1.9]; % output sample locations
% v4o = [1.5,1,1.9]; % output sample locations
% GRIDDATA
[m3i,m4i] = ndgrid(1:3,1:2); % Exclude 1st and 2nd dimensions!
m3o = m3i + 0.12;
m4o = m4i + 0.56;
%[m3o,m4o] = ndgrid(v3o,v4o); % Exclude 1st and 2nd dimensions!
one = nan(5,4,3,2);
for ii = 1:5 % d1
    for jj = 1:4 % d2
        tmp = inp(ii,jj,:); % all inputs are vectors:
        mat = griddata(m3i(:),m4i(:),tmp(:),m3o(:),m4o(:));
        one(ii,jj,:) = mat(:);
    end
end
% INTERPN:
[m1o,m2o,m3i,m4i] = ndgrid(1:5,1:4,1:3,1:2);
m3o = m3i + 0.12;
m4o = m4i + 0.56;
two = interpn(inp,m1o,m2o,m3o,m4o);
% Checking:
isequal(one,two)

the main difference is to use a simple alteration of the original grid as the output grid and to use non-integers as the indices. I am no longer getting equal results.

[Roee Mazor]

an even smaller modification (non-integers) also doesnt produce equality:

inp = rand(5,4,3,2); % fake input data
v3o = [1.1,3,2]; % output sample locations
v4o = [2.1,1,2]; % output sample locations
% GRIDDATA
[m3o,m4o] = ndgrid(v3o,v4o); % Exclude 1st and 2nd dimensions!
[m3i,m4i] = ndgrid(1:3,1:2); % Exclude 1st and 2nd dimensions!
one = nan(5,4,numel(v3o),numel(v4o));
for ii = 1:5 % d1
    for jj = 1:4 % d2
        tmp = inp(ii,jj,:); % all inputs are vectors:
        mat = griddata(m3i(:),m4i(:),tmp(:),m3o(:),m4o(:));
        one(ii,jj,:) = mat(:);
    end
end
% INTERPN:
[m1o,m2o,m3o,m4o] = ndgrid(1:5,1:4,v3o,v4o);
two = interpn(inp,m1o,m2o,m3o,m4o);
% Checking:
isequal(one,two)

[Stephen23]

"to use non-integers as the indices. I am no longer getting equal results."

In general I would not expect them to be equal: you are performing completely different numeric operations on different binary floating point data, so getting exactly the same output values would be nothing short of magic.

griddata performs delauney triangulation between the input sample points and interpolates along those edges. interpn presumably applies a repeated linear interpolation along each non-scalar dimension. Given non-integer input data supplied to two totally different algorithms which therefore accumulate totally different floating point error, it would be nothing short of a miracle to get exact numeric equivalence of the output data.**

In your original question you wrote "... if the interpolated points fall exactly on the original grid it wouldnt use the information from neighbouring points and if so the two methods should produce the exact same thing", which is exactly what my answer shows: the outputs are exactly the same when the sample point fall exactly on the original grid. But for the examples in your comments above the sample points do NOT fall on the original grid, that is quite a different situation.

Your expectation that "...the two methods should produce the exact same thing..." for any aribitrary non-integer values does not reflect the reality of how binary floating point numbers actually behave when numeric operations are performed:

https://docs.oracle.com/cd/E19957-01/806-3568/ncg_goldberg.html

** excepting some special cases: low-value integers and perhaps powers of two.

[Roee Mazor]

The different methods should explain the difference. The data type is not to be blamed here in my opinion though.

I will try to use interp2 or interpn instead of griddata and see if that produces identical results.

[Stephen23]

"I will try to use interp2 or interpn instead of griddata and see if that produces identical results."

That the results might or might not be identical is meaningless, because that floating point error is effectively just noise within the accuracy we would expect from that data type for those (unknown but rather irrelevant) numeric operations.

"The data type is not to be blamed here in my opinion though."

No one is "blaming" anything, just several people have explained that it is entirely expected that numeric operations on binary floating point numbers will accumulate floating point error. This is completely dependent on the data type.

[Roee Mazor]

Ok, thank you very much!