Bilinear interpolation for large grids

Question

I need to make a bilinear interpolation of a function $Y(i,j)$ tabulated on a $100\times 100$ grid. I tried to do it with the Fortran polin2.f and polint.f subroutines of Numerical Recipes. These routines seem to be written for small grids (maximum expected values are $20\times 20$). They work fine in the code when I consider a small grid.

Can anybody suggest me a subroutine that can be used for larger grids?

Answer 1

You could try the GNU Scientific Library: https://www.gnu.org/software/gsl/doc/html/interp.html#id5

Answer 2

I finally found a couple of function in f90 that can perform a bilinear interpolation for a mesh of any size quite nicely. These routines can be found here.

Answer 3

Here's an idea that works if you have 1D interpolation routines and the query points are on a regular grid.

You can compute the SVD of the data matrix, $Y = U\Sigma V^T$, interpolate the columns of $U$ and $V$, to obtain, $\hat{U}$ and $\hat{V}$, and then compute the interpolated values as $\hat{U} \Sigma \hat{V}^T$. If $Y$ is approximately low-rank, you can discard singular vectors corresponding to small singular values.

Note that when interpolating the columns of $U$ and $V$, determining which pair of nodes is closest to the query point only needs to be done once for all columns.

If the number of query points is large, forming the output via this matrix-matrix product can be faster than direct bilinear interpolation.

Here's some matlab code that demonstrates this idea:

x = linspace(-1, 1, 100);
[X,Y] = meshgrid(x);
A = exp(-(X.^2+Y.^2)).*sin(2*pi*X.*Y);
xq = linspace(-1, 1, 1000);
[Xq,Yq] = meshgrid(xq);
tic
Aq = interp2(X, Y, A, Xq, Yq);
toc

tic
[U,S,V] = svd(A);
sv = diag(S);
Uq = interp1(x, U, xq);
Vq = interp1(x, V, xq);
AqSVD = Uq*(sv.*Vq');
toc
norm(AqSVD - Aq, inf)/norm(Aq, inf)

Elapsed time is 0.033942 seconds.
Elapsed time is 0.012933 seconds.

ans =

      1.29298230462322e-15