Interpolate a ND vector at once

General Usage
1

Hi all,

I am trying to perform linear interpolation on some gridded data using the Interpolations.jl package. Suppose that I have 2 variables (x,y) and a function z=f(x,y). I want to interpolate a vector of values V={ (x1,y1), (x2,y2), … (xm,ym) } at once, without using a loop across the m points.

Suppose that I have the following:

using Interpolations
x        = collect(linspace(0,1,10))  #x grid
y        = collect(linspace(0,1,10))  #y grid
function fxy(x,y)
	z = x.^2+y.^2
	return z
end
	z = zeros(length(x), length(y))
	for i_x=1:length(x)
	for i_y=1:length(y)
	z[i_x,i_y] = fxy(x[i_x],y[i_y])
	end
	end

I want to interpolate a vector of data points at once, without the use of a loop across the different points. To be clear:

#Create interpolation object
FF  = interpolate((x,y),z,Gridded(Linear()))

#I want to interpolate fxy on the following points:
xy_data = [collect(linspace(1,5,5)) collect(linspace(5,10,5))] #Object that I want to interpolate

#This is the output that I want (a 5x1 vector with the values of fxy at each of the data points)
output_desired = zeros(size(xy_data,1))
for i=1:size(xy_data,1)
output_desired[i] = FF[xy_data[i,1], xy_data[i,2]]
end

I cannot find a way of obtaining the same result without the use of a loop.

For my application, I need to repeat these interpolations many times as the interpolation object will be inside a non-linear solver. So performance is key and I need to avoid looping across the M points of the xy_data vector.

For example, the following will interpolate all possible combinations of (x,y) in the xy_data matrix (it returns a 5x5 matrix)

output_matrix = FF[xy_data[:,1], xy_data[:,2]]

The above line is in fact faster than using the loop across points. But I just want the diagonal of output_matrix.

Is there any way to do this?

Thanks,

2

Can I ask why? There is no performance reason not to use a loop.

Are you sure? Let’s try it:

julia> function method1(FF, xy_data)
         output_desired = zeros(size(xy_data,1))
         for i=1:size(xy_data,1)
           output_desired[i] = FF[xy_data[i,1], xy_data[i,2]]
         end
         output_desired
       end
method1 (generic function with 1 method)

julia> function method2(FF, xy_data)
         diag(FF[xy_data[:, 1], xy_data[:, 2]])
       end
method2 (generic function with 1 method)

julia> using BenchmarkTools

julia> @btime method1($FF, $xy_data);
  165.027 ns (1 allocation: 128 bytes)

julia> @btime method2($FF, $xy_data);
  341.631 ns (10 allocations: 1.05 KiB)

The loop is significantly faster and will be vastly faster than the matrix version as your data gets larger.

Note that if you are operating with (non-constant) global variables, as your original code does, you may get timing results that are both slower and somewhat uninformative. That’s why the very first performance tip is “Avoid Global Variables”.

3

Thanks for you reply, rdeits.

I’m not using global variables in my code. When I mentioned the performance, I was comparing these 2 functions actually:

function method1(FF, xy_data)  #ENTIRE MATRIX
         output_matrix = zeros(size(xy_data,1),size(xy_data,1))
         for i=1:size(xy_data,1)
         for j=1:size(xy_data,2)
         output_matrix[i,j] = FF[xy_data[i,1], xy_data[j,2]]
         end
         end
         output_matrix
       end
         
function method2(FF, xy_data)  #ENTIRE MATRIX
         output_matrix = FF[xy_data[:, 1], xy_data[:, 2]]
       end

This is the time it takes to compile and run for the first time:

Method 1:  0.024658 seconds (5.05 k allocations: 224.457 KiB)
Method 2:  0.007001 seconds (1.95 k allocations: 95.635 KiB)

In subsequent runs, as you said, it is true that the method2 fx is slower once the function has been compiled. In my code, however, I need to re-define the “methodX fx” along the main loop of the code and therefore I need to compile this function in each loop. At least this is how I was doing it so far. I will make sure to change this.

Thanks again,

4

Can you explain more about why this is necessary? That’s a very unusual design pattern, and there’s likely to be a better way.