# Help with poor LinearInterpolation performance

**URL:** <https://discourse.julialang.org/t/help-with-poor-linearinterpolation-performance/22325>\
**Category:** General Usage\
**Created:** [March 25, 2019, 11:01pm UTC](https://discourse.julialang.org/t/help-with-poor-linearinterpolation-performance/22325 "2019-03-25T23:01:14Z")\
**Posts on this page:** 11\
**Page:** 1

<div class="post-metadata">

**Author:** ![angelv](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/angelv/32/133_2.png) [@angelv](https://discourse.julialang.org/u/angelv)\
**Post date:** [March 25, 2019, 11:01pm UTC](https://discourse.julialang.org/t/help-with-poor-linearinterpolation-performance/22325/1 "2019-03-25T23:01:14Z")

</div>

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:

```julia
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:

```julia
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
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

---

<div class="post-metadata">

**Author:** ![JeffreySarnoff](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/jeffreysarnoff/32/1980_2.png) [@JeffreySarnoff](https://discourse.julialang.org/u/JeffreySarnoff)\
**Post date:** [March 26, 2019, 12:18am UTC](https://discourse.julialang.org/t/help-with-poor-linearinterpolation-performance/22325/2 "2019-03-26T00:18:29Z")

</div>

Does this help?

```julia
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);

```

---

<div class="post-metadata">

**Author:** ![NiclasMattsson](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/niclasmattsson/32/21988_2.png) [@NiclasMattsson](https://discourse.julialang.org/u/NiclasMattsson)\
**Post date:** [March 26, 2019, 12:48am UTC](https://discourse.julialang.org/t/help-with-poor-linearinterpolation-performance/22325/3 "2019-03-26T00:48:20Z")

</div>

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.

---

<div class="post-metadata">

**Author:** ![aharoun](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/aharoun/32/6887_2.png) [@aharoun](https://discourse.julialang.org/u/aharoun)\
**Post date:** [March 26, 2019, 12:57am UTC](https://discourse.julialang.org/t/help-with-poor-linearinterpolation-performance/22325/4 "2019-03-26T00:57:50Z")

</div>

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.

```julia
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

```

---

<div class="post-metadata">

**Author:** ![NiclasMattsson](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/niclasmattsson/32/21988_2.png) [@NiclasMattsson](https://discourse.julialang.org/u/NiclasMattsson)\
**Post date:** [March 26, 2019, 1:12am UTC](https://discourse.julialang.org/t/help-with-poor-linearinterpolation-performance/22325/5 "2019-03-26T01:12:10Z")

</div>

True, if all he wants is the two endpoints. But judging from the [IDL interpolate manpage](https://www.harrisgeospatial.com/docs/INTERPOLATE.html), 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.

---

<div class="post-metadata">

**Author:** ![JeffreySarnoff](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/jeffreysarnoff/32/1980_2.png) [@JeffreySarnoff](https://discourse.julialang.org/u/JeffreySarnoff)\
**Post date:** [March 26, 2019, 1:26am UTC](https://discourse.julialang.org/t/help-with-poor-linearinterpolation-performance/22325/6 "2019-03-26T01:26:56Z")

</div>

```julia
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
```

---

<div class="post-metadata">

**Author:** ![angelv](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/angelv/32/133_2.png) [@angelv](https://discourse.julialang.org/u/angelv)\
**Post date:** [March 26, 2019, 8:49am UTC](https://discourse.julialang.org/t/help-with-poor-linearinterpolation-performance/22325/7 "2019-03-26T08:49:35Z")

</div>

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

---

<div class="post-metadata">

**Author:** ![angelv](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/angelv/32/133_2.png) [@angelv](https://discourse.julialang.org/u/angelv)\
**Post date:** [March 26, 2019, 8:53am UTC](https://discourse.julialang.org/t/help-with-poor-linearinterpolation-performance/22325/8 "2019-03-26T08:53:43Z")

</div>

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.

---

<div class="post-metadata">

**Author:** ![angelv](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/angelv/32/133_2.png) [@angelv](https://discourse.julialang.org/u/angelv)\
**Post date:** [March 26, 2019, 10:44am UTC](https://discourse.julialang.org/t/help-with-poor-linearinterpolation-performance/22325/9 "2019-03-26T10:44:32Z")

</div>

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:

```julia
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:

```julia
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 ☹

AdV

---

<div class="post-metadata">

**Author:** ![tim.holy](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/tim.holy/32/52_2.png) [@tim.holy](https://discourse.julialang.org/u/tim.holy)\
**Post date:** [March 26, 2019, 11:12am UTC](https://discourse.julialang.org/t/help-with-poor-linearinterpolation-performance/22325/10 "2019-03-26T11:12:52Z")

</div>

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.

---

<div class="post-metadata">

**Author:** ![angelv](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/angelv/32/133_2.png) [@angelv](https://discourse.julialang.org/u/angelv)\
**Post date:** [March 26, 2019, 2:10pm UTC](https://discourse.julialang.org/t/help-with-poor-linearinterpolation-performance/22325/11 "2019-03-26T14:10:38Z")

</div>

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

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

Many thanks,