5

I need to perform many (thousands) of look-up operations, where the break-points in the look-up do not change. A simple example would be,

% create some dummy data
% In practice
% - the grid will be denser, and not as regular
% - each page of v will be different
% - v will have thousands of pages
[x,y] = ndgrid(-5:0.1:5);
n_surfaces = 10;
v = repmat(sin(x.^2 + y.^2) ./ (x.^2 + y.^2),1,1,n_surfaces);

[xq,yq] = ndgrid(-5:0.2:5);
vq = nan([size(xq),n_surfaces]);
for idx = 1:n_surfaces
   F = griddedInterpolant(x,y,v(:,:,idx));
   vq(:,:,idx) = F(xq,yq);
end

Note that the above code can be sped up slightly by doing,

F = griddedInterpolant(x,y,v(:,:,1));
for idx = 1:n_surfaces
   F.Values = v(:,:,idx);
   vq(:,:,idx) = F(xq,yq);
end

However, in general interpolation is a two step process,

  1. Determining the index and interval fraction for each new point
  2. Performing the interpolation to obtain the new value

and in the above code both of these steps are being performed during every loop. However Step 1 will be identical in every loop and hence performing it thousands of times is inefficient. I'm wondering if anyone has a workaround to split the two steps and only perform the first step once, and just perform the second step in the loop?

(For those familiar with Simulink, this is equivalent to using the Prelookup in conjunction with multiple Interpolation Using Prelookup blocks.)

Edit:

The question linked in the comment by @rahnema1 (Precompute weights for multidimensional linear interpolation) is pretty much what I am looking for. However, on converting that code to run on the CPU (rather than a GPU), and using double arithmetic, it is about 3 times slower than using the m-code at the start of my question. That timing seems to hold irrespective of the number of surfaces being interpolated (I have tried values from 10 through to 1000.)

The problem is in performing the indexing operation V(I) used in the linked code. Even when the complete operation sum(W.*V(I),2) is implemented in a mex file, the execution times are slower than the above m-code.

Phil Goddard
  • 10,571
  • 1
  • 16
  • 28
  • I guess you need to implement your own interpolation algorithm, so you can split the components up. By default `griddedInterpolant` uses linear interpolation, so this should be quite simple to implement. -- Actually, are the `x` and `xq` vectors you show here the ones you actually intend to use? For the example code you could just do `vq=v(1:2:end,1:2:end,:)`. – Cris Luengo Jul 24 '19 at 15:24
  • I have started implementing my own algorithm, I'm just wondering if anyone has already solved the problem as it is a pretty standard thing to want to do. (Given that the Simulink development team are obviously aware that this can provide a significant efficiency, I'm a little surprised that the functionality is missing in MATLAB). Unfortunately in practice the `v` and `vq` points aren't as well behaved as in the example so indexing won't work. – Phil Goddard Jul 24 '19 at 15:49
  • 1
    Related question: https://stackoverflow.com/q/46126019/6579744 – rahnema1 Jul 24 '19 at 15:57
  • Note also that, going from a grid to a grid, you could skip the `ndgrid` calls and compute weights per row and column. This saves computations as well as memory. – Cris Luengo Jul 24 '19 at 16:01

0 Answers0