Help with poor LinearInterpolation performance

General Usage
1

Hi,

I’m trying to use the function LinearInterpolation from Interpolations.jl for the first time, but I guess I’m doing many things wrong, since I get very bad performance when compared with a similar IDL routine.

The julia function I have now is:

function interp_test(size)
    xs = range(0,10,length=size); ys = range(0,10,length=size); zs = range(0,10,length=size)
    xsi = xs .+ 0.5 ; ysi = ys .+ 0.5 ; zsi = zs .+ 0.5

    dat = randn(size,size,size)
    dat_out = zeros(size,size,size)

    dat_int = LinearInterpolation((xs,ys,zs),dat,extrapolation_bc = Line());

    for C in CartesianIndices((length(xsi),length(ysi),length(zsi)))
        dat_out[C] = dat_int(xsi[C[1]],ysi[C[2]],zsi[C[3]])
    end

    println(dat_out[1,1,1],dat_out[end,end,end])
end

The similar IDL that I have is:

PRO interp_test, size
  
  xs = findgen(size)
  ys = findgen(size)
  zs = findgen(size)
  
  xsi = xs + 0.5
  ysi = ys + 0.5
  zsi = zs + 0.5
  
  dat = randomu(1,size,size,size)
   
  dat_out = interpolate(dat,xsi,ysi,zsi,/grid)
  
  print,dat_out[0,0,0],dat_out[-1,-1,-1]
end

But the times in Julia are much worse:

julia> @time interp_test(500)                                                                                                                                                                                        
1.171280353464721146.0781291105493                                                                                                                                                                                   
 27.314111 seconds (125.00 M allocations: 4.657 GiB, 9.00% gc time)

IDL> interp_test, 500
interp_test, 500
     0.457461    0.0822094
% Time elapsed: 1.5881102 seconds.
IDL>

Any help on getting this to work better?

Many thanks,
AdV

2

Does this help?

using Interpolations
n = 20;
xs = range(0,10,length=n); ys = range(0,10,length=n); zs = range(0,10,length=n);
dat = randn(n,n,n);
result = LinearInterpolation((xs,ys,zs),dat);
3

I think that only sets up the trilinear interpolation without evaluating it at any new grid points. To perform the evaluation, I think you have to add result(xs.+0.5,ys.+0.5,zs.+0.5) which shouldn’t be faster than the original example (it throws an error unless you add the extrapolation option as in the OP).

A (small) part of the speed difference lies in the fact that IDL is “cheating” by doing the calculations in single precision.

4

Your Julia function spends most of its time at the for loop, which looks unnecessary. You can simply evaluate the interpolation result at the desired point you want to return.

function interp_test(size)
    xs = range(0,10,length=size); ys = range(0,10,length=size); zs = range(0,10,length=size)
    xsi = xs .+ 0.5 ; ysi = ys .+ 0.5 ; zsi = zs .+ 0.5

    dat = randn(size,size,size)
    dat_out = zeros(size,size,size)

    dat_int = LinearInterpolation((xs,ys,zs),dat,extrapolation_bc = Line());

    return dat_int(xsi[1],ysi[1],zsi[1]), dat_int(xsi[end],ysi[end],zsi[end])
end
5

True, if all he wants is the two endpoints. But judging from the IDL interpolate manpage, the IDL function evaluates the interpolation at all 500^3 grid points, as does @angelv’s Julia code. So the large performance difference remains. Or am I wrong? Is IDL doing lazy evaluation?

EDIT: Another contribution to the difference: the IDL function is multithreaded.

6 
it seems to slow down with n>70
 
function interp_test(n)
           xs = range(0,10,length=n); ys = range(0,10,length=n); zs = range(0,10,length=n)
           xsi = xs .+ 0.5 ; ysi = ys .+ 0.5 ; zsi = zs .+ 0.5

           dat = randn(n, n, n)
           dat_out = zeros(n, n, n)

           dat_int = LinearInterpolation((xs,ys,zs),dat,extrapolation_bc = Line());

           for xi=1:n
             xv = xsi[xi]
             for yi=1:n
              yv = ysi[yi] 
               for zi=1:n
                  dat_out[xi, yi, zi] = dat_int(xv, yv, zsi[zi])
               end
             end  
           end
           
           return dat_out
       end
7

Hi,
as NiclasMattsson points, this only sets the interpolation, but does not evaluate it at any points, and I need to evaluate the interpolated values at the whole 500^3 grid.
Thanks

8

AFAIK, IDL doesn’t do lazy evaluation. Just in case, I changed the last line to be
return sum(dat_out) in my Julia code and print, total(sum) in the IDL code to make sure that all the values have to be evaluated before returning, and the time difference between Julia and IDL stays basically the same.

9

In my original post I was using different computers (I didn’t have IDL in my machine), etc. so here it goes a fairer comparison: modified IDL to use doubles everywhere, and to use only 1 CPU, plus now both Julia and IDL running in the same box.

Julia time has not changed: ~25 s

IDL code now is:

PRO interp_test, size
  TIC
  
  xs = DOUBLE(findgen(size))
  ys = DOUBLE(findgen(size))
  zs = DOUBLE(findgen(size))
  
  xsi = xs + 0.5
  ysi = ys + 0.5
  zsi = zs + 0.5
  
  dat = DOUBLE(randomu(1,size,size,size))
  
  dat_out = interpolate(dat,xsi,ysi,zsi,/double,/grid)
  
  help,xs
  help,xsi
  help,dat
  help,dat_out
  
  print,total(dat_out)

  TOC
end

Run like:

IDL> CPU, TPOOL_NTHREADS = 1
CPU, TPOOL_NTHREADS = 1
IDL> interp_test, 500
interp_test, 500
XS              DOUBLE    = Array[500]
XSI             DOUBLE    = Array[500]
DAT             DOUBLE    = Array[500, 500, 500]
DAT_OUT         DOUBLE    = Array[500, 500, 500]
       62495677.
% Time elapsed: 6.0175951 seconds.
IDL>

Still, IDL beats Julia badly :frowning:

AdV

10

The fact that so much memory is allocated (one allocation per iteration) is a clear sign of some problem in the Julia implementation—it shouldn’t allocate at all. I don’t have time to debug this now, but perhaps start playing with dropping certain elements and seeing if the allocation goes away. The two obvious candidates are the ranges (just use, e.g., itp = extrapolate(interpolate(randn(sz, sz, sz), BSpline(Linear())), Line())) and the extrapolation. (Does the IDL implementation set the ranges to 0…10 or are you only doing that on the Julia side?) Try it all with ProfileView and then file an issue with your findings at Interpolations.jl.

1 Like
11

OK, I filed an issue at Interpolations.jl (Performance issues with linear interpolations · Issue #311 · JuliaMath/Interpolations.jl · GitHub).

Let’s see if there is some quick way to fix it.

Many thanks,

2 Likes