Interpolate a ND vector at once

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
	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])

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]]

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?


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]]
method1 (generic function with 1 method)

julia> function method2(FF, xy_data)
         diag(FF[xy_data[:, 1], xy_data[:, 2]])
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”.


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]]
function method2(FF, xy_data)  #ENTIRE MATRIX
         output_matrix = FF[xy_data[:, 1], xy_data[:, 2]]

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,

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.

1 Like