# Why are my linear interpolations 10x faster in MATLAB?

**URL:** <https://discourse.julialang.org/t/why-are-my-linear-interpolations-10x-faster-in-matlab/41232>\
**Category:** Performance\
**Tags:** interpolations\
**Created:** [June 11, 2020, 10:27pm UTC](https://discourse.julialang.org/t/why-are-my-linear-interpolations-10x-faster-in-matlab/41232 "2020-06-11T22:27:20Z")\
**Posts on this page:** 20\
**Page:** 1

<div class="post-metadata">

**Author:** ![elenev](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/elenev/32/18440_2.png) [@elenev](https://discourse.julialang.org/u/elenev)\
**Post date:** [June 11, 2020, 10:27pm UTC](https://discourse.julialang.org/t/why-are-my-linear-interpolations-10x-faster-in-matlab/41232/1 "2020-06-11T22:27:20Z")

</div>

The consensus on this forum seems to be that Interpolations.jl is extremely fast, so I must not be using it in a smart way.

Here’s an MWE of the Matlab code I’m trying to adapt to Julia. If some stuff looks weird (e.g. generating random numbers inside a loop), it’s because I’m using that stuff as a runnable placeholder for more complex operations in my application that would just distract.

```julia
f=@(a,b) a.^2 + b.^2;

xs = linspace(0,5,100);
ys = linspace(-3,3,200);

[XX,YY] = meshgrid(xs,ys);
vals = f(XX,YY);

fct = griddedInterpolant( {xs, ys} ,vals', 'linear','none');

Npt = 1e4;
pts = [rand(Npt,1)*5, zeros(Npt,1)];

Ntrials = 10;
tic;
for i=1:Ntrials
	N = 1e3;
	results = zeros(Npt,N);
	for j=1:1e3
		pts(:,2) = rand(Npt,1) * 6 - 3;
		results(:,j) = fct(pts);
	end
end
elapsedtime=toc/Ntrials;
fprintf('mean time: %0.2f ms\n', elapsedtime*1000 );

```

When I run this, I get:

```julia
mean time: 284.85 ms

```

Now, here’s how I implemented this in Julia:

```julia
using Interpolations, BenchmarkTools
f(a,b) = a^2 + b^2
const xs = range(0,5,length=100)
const ys = range(-3,3,length=200)
const vals = [f(x,y) for x in xs, y in ys]
const fct = LinearInterpolation( (xs,ys), vals)
const Npt = 10_000;
pts = [rand(1,Npt); zeros(1,Npt)]

function evaluate!(pts,fct)
    N = 1_000
    result = Matrix{Float64}(undef,size(pts,2),N)
    for j in 1:N
        pts[2,:] = rand(1,Npt)*6 .- 3
        # Need this loop b/c unlike Matlab fct here requires scalar args
        for p in axes(pts,2)
            result[p,j] = fct(pts[:,p]...)
        end
    end
    return result
end

benchres = @benchmark evaluate!(pts,fct)

```

The output of `@show benchres` is

```julia
benchres = Trial(2.393 s)
BenchmarkTools.Trial:
  memory estimate: 1.64 GiB
  allocs estimate: 40006002
  --------------
  minimum time: 2.393 s (27.67% GC)
  median time: 2.415 s (27.72% GC)
  mean time: 2.408 s (27.91% GC)
  maximum time: 2.416 s (28.33% GC)
  --------------
  samples: 3
  evals/sample: 1

```

I can reduce this ~3x down to 840 ms and only 6000 allocations by replacing `pts[:,p]...` with `pts[1,p],pts[2,p]` but (1) that’s still \>3x slower than Matlab, and (2) the evaluate! function now only works for 2-D functions (yes, I know I can replace splats with dispatching on the type parameter of the interpolant object, but that’s pretty ugly and the resulting time is still bounded below by 840 ms)

Where did I mess up?

---

<div class="post-metadata">

**Author:** ![sgaure](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/sgaure/32/14779_2.png) [@sgaure](https://discourse.julialang.org/u/sgaure)\
**Post date:** [June 11, 2020, 10:49pm UTC](https://discourse.julialang.org/t/why-are-my-linear-interpolations-10x-faster-in-matlab/41232/2 "2020-06-11T22:49:00Z")

</div>

You’re splatting an array, i.e. the call `fct(pts[:,p]...)` allocates a copy, and array splatting is generally slow. Do `fct(pts[1,p], pts[2,p])` instead. An throw in an `@inbounds` to avoid array bounds checking. This speeds it up by a factor of 10 on my machine.

```julia

function evaluate!(pts,fct)
    N = 1_000
    result = Matrix{Float64}(undef,size(pts,2),N)
    @inbounds for j in 1:N
        pts[2,:] = rand(1,Npt)*6 .- 3
        # Need this loop b/c unlike Matlab fct here requires scalar args
        for p in axes(pts,2)
            result[p,j] = fct(pts[1,p],pts[2,p])
        end
    end
    return result
end

```

---

<div class="post-metadata">

**Author:** ![elenev](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/elenev/32/18440_2.png) [@elenev](https://discourse.julialang.org/u/elenev)\
**Post date:** [June 11, 2020, 10:53pm UTC](https://discourse.julialang.org/t/why-are-my-linear-interpolations-10x-faster-in-matlab/41232/3 "2020-06-11T22:53:10Z")

</div>

Thanks. As I mentioned in my post, I tried that:

> I can reduce this ~3x down to 840 ms and only 6000 allocations by replacing pts[:,p]… with pts[1,p],pts[2,p] but (1) that’s still \>3x slower than Matlab, and (2) the evaluate! function now only works for 2-D functions (yes, I know I can replace splats with dispatching on the type parameter of the interpolant object, but that’s pretty ugly and the resulting time is still bounded below by 840 ms)

In other words, it does help considerably but it still leaves me with 3x the time of Matlab, and comes at the expense of generic applicability (i.e. I won’t always have 2 dimensions).

Adding `@inbounds` in addition lowers it to 694 ms, still about 2.5x slower.

---

<div class="post-metadata">

**Author:** ![Elrod](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/elrod/32/22461_2.png) [@Elrod](https://discourse.julialang.org/u/Elrod)\
**Post date:** [June 11, 2020, 10:53pm UTC](https://discourse.julialang.org/t/why-are-my-linear-interpolations-10x-faster-in-matlab/41232/4 "2020-06-11T22:53:55Z")

</div>

MATLAB:

> `pts(:,2) = rand(Npt,1) * 6 - 3;`

Julia:

> `pts[2,:] = rand(1,Npt)*6 .- 3`

Given that `pts` is `2` x `10_000`, assigning to 10\_000 elements on each information is going to hurt performance (and give you the wrong answer).

---

<div class="post-metadata">

**Author:** ![elenev](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/elenev/32/18440_2.png) [@elenev](https://discourse.julialang.org/u/elenev)\
**Post date:** [June 11, 2020, 10:59pm UTC](https://discourse.julialang.org/t/why-are-my-linear-interpolations-10x-faster-in-matlab/41232/5 "2020-06-11T22:59:52Z")

</div>

Ugh duh, sorry. The vector vs. one-row matrix always confuses me when I’m going back and forth between languages. Fixed:

```julia
pts[2,:] = rand(Npt)*6 .- 3

```

With no splatting and with `@inbounds`, it’s now 650 ms. If I further replace `=` with `.=`, I can shave off another 2000 allocations and 50 ms. But we’re still 2x MATLAB time, and the code is still less flexible b/c it’s hard-coding two dimensions.

---

<div class="post-metadata">

**Author:** ![sgaure](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/sgaure/32/14779_2.png) [@sgaure](https://discourse.julialang.org/u/sgaure)\
**Post date:** [June 11, 2020, 11:06pm UTC](https://discourse.julialang.org/t/why-are-my-linear-interpolations-10x-faster-in-matlab/41232/6 "2020-06-11T23:06:14Z")

</div>

Ok, try replacing the inner loop with

`result[:,j] .= fct.(eachrow(pts)...)`

---

<div class="post-metadata">

**Author:** ![Elrod](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/elrod/32/22461_2.png) [@Elrod](https://discourse.julialang.org/u/Elrod)\
**Post date:** [June 11, 2020, 11:07pm UTC](https://discourse.julialang.org/t/why-are-my-linear-interpolations-10x-faster-in-matlab/41232/7 "2020-06-11T23:07:18Z")

</div>

What are the dimensions of `pts` in Julia and MATLAB?  
If in MATLAB it is 10\_000 x 2, and in Julia it is 2 x 10\_000, I’d like to point out that:

```julia
julia> x1 = rand(2,10_000);

julia> x2 = rand(10_000,2);

julia> @benchmark $x1[1,:] .= 0.5
BenchmarkTools.Trial:
  memory estimate: 0 bytes
  allocs estimate: 0
  --------------
  minimum time: 3.142 μs (0.00% GC)
  median time: 3.431 μs (0.00% GC)
  mean time: 3.454 μs (0.00% GC)
  maximum time: 6.034 μs (0.00% GC)
  --------------
  samples: 10000
  evals/sample: 8

julia> @benchmark $x2[:,1] .= 0.5
BenchmarkTools.Trial:
  memory estimate: 0 bytes
  allocs estimate: 0
  --------------
  minimum time: 1.132 μs (0.00% GC)
  median time: 1.175 μs (0.00% GC)
  mean time: 1.179 μs (0.00% GC)
  maximum time: 2.933 μs (0.00% GC)
  --------------
  samples: 10000
  evals/sample: 10

```

Of course, doing it a thousand times will only make up another couple of milliseconds.

---

<div class="post-metadata">

**Author:** ![chakravala](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/chakravala/32/6832_2.png) [@chakravala](https://discourse.julialang.org/u/chakravala)\
**Post date:** [June 11, 2020, 11:27pm UTC](https://discourse.julialang.org/t/why-are-my-linear-interpolations-10x-faster-in-matlab/41232/8 "2020-06-11T23:27:59Z")

</div>

Why would x1 and x2 have different performance? Seems like that shouldnt be a thing, and if it is a thing then it should be documented somewhere.

---

<div class="post-metadata">

**Author:** ![tbeason](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/tbeason/32/15898_2.png) [@tbeason](https://discourse.julialang.org/u/tbeason)\
**Post date:** [June 11, 2020, 11:43pm UTC](https://discourse.julialang.org/t/why-are-my-linear-interpolations-10x-faster-in-matlab/41232/9 "2020-06-11T23:43:43Z")

</div>

Row-major vs column-major storage. It is faster to access items in the array in the sequence that they are stored.

---

<div class="post-metadata">

**Author:** ![Elrod](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/elrod/32/22461_2.png) [@Elrod](https://discourse.julialang.org/u/Elrod)\
**Post date:** [June 11, 2020, 11:46pm UTC](https://discourse.julialang.org/t/why-are-my-linear-interpolations-10x-faster-in-matlab/41232/10 "2020-06-11T23:46:02Z")

</div>

This is general (other than requiring at least 2 dimensions, which the MATLAB code also does), and fairly fast:

```julia
using Interpolations, Random, BenchmarkTools
f(a,b) = a^2 + b^2
const xs = range(0,5,length=100);
const ys = range(-3,3,length=200);
const vals = [f(x,y) for x in xs, y in ys];
const fct = LinearInterpolation( (xs,ys), vals);
const Npt = 10_000;
pts1 = rand(Npt);
pts2 = zeros(Npt);
function evaluate!(fct, pts::Vararg{<:AbstractVector,N}) where {N}
    Niter = 1_000
    pts2 = pts[2]
    result = Matrix{Float64}(undef,length(pts2),Niter)
    @inbounds for j in 1:Niter
        rand!(pts2)
        # Need this loop b/c unlike Matlab fct here requires scalar args
        for p in axes(pts2,1)
            pts2[p] = 6pts2[p] - 3
            result[p,j] = fct(ntuple(n -> pts[n][p], Val(N))...)
        end
    end
    return result
end

```

This yields:

```julia
julia> @benchmark evaluate!(fct, $pts1, $pts2)
BenchmarkTools.Trial:
  memory estimate: 76.29 MiB
  allocs estimate: 2
  --------------
  minimum time: 138.638 ms (0.00% GC)
  median time: 140.861 ms (0.66% GC)
  mean time: 144.204 ms (3.04% GC)
  maximum time: 208.384 ms (32.85% GC)
  --------------
  samples: 35
  evals/sample: 1

```

@chakravala  
[Here](https://docs.julialang.org/en/v1/manual/performance-tips/index.html#man-performance-column-major-1) is the documentation.  
Also worth pointing out that how matrices/arrays are stored varies by language. Python is row-major, for example, and C++ libraries vary or let you choose.

---

<div class="post-metadata">

**Author:** ![sgaure](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/sgaure/32/14779_2.png) [@sgaure](https://discourse.julialang.org/u/sgaure)\
**Post date:** [June 11, 2020, 11:52pm UTC](https://discourse.julialang.org/t/why-are-my-linear-interpolations-10x-faster-in-matlab/41232/11 "2020-06-11T23:52:09Z")

</div>

> [@chakravala](#):
>
> Why would x1 and x2 have different performance? Seems like that shouldnt be a thing, and if it is a thing then it should be documented somewhere.

It’s a cpu-cache issue. Memory is transferred between the cpu and memory in chunks of 64 consecutive bytes or so (a cache line). Each such transfer is a “cache fill”. If possible you should use all the bytes in a cache line before you advance to the next. That is what you do in x2, but in x1 only half of them are used. So you need twice as many cache fills (which may take a significant number of cpu clock cycles, like 50 or 100 or 300 depending on your cpu-architecture. Then “prefetching” is important, but may not alleviate the problem entirely).

---

<div class="post-metadata">

**Author:** ![Elrod](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/elrod/32/22461_2.png) [@Elrod](https://discourse.julialang.org/u/Elrod)\
**Post date:** [June 11, 2020, 11:55pm UTC](https://discourse.julialang.org/t/why-are-my-linear-interpolations-10x-faster-in-matlab/41232/12 "2020-06-11T23:55:57Z")

</div>

It’s also a vectorization issue.

If you have AVX2, you can do:

```julia
@fastmath begin
x[1,1] = 6x[1,1] - 3
x[2,1] = 6x[2,1] - 3
x[3,1] = 6x[3,1] - 3
x[4,1] = 6x[4,1] - 3
end

```

with 1 load instruction, 1 fused multiply add, and 1 store (you’ll also have to load the constants, but that’s only necessary once if this happens in a loop).  
Meanwhile, this:

```julia
@fastmath begin
x[1,1] = 6x[1,1] - 3
x[1,2] = 6x[1,2] - 3
x[1,3] = 6x[1,3] - 3
x[1,4] = 6x[1,4] - 3
end

```

will probably use 4 load instructions, 4 fused multiply add instructions, and 4 stores.

With big matrices, the first one is going to be memory bound rather than compute bound, so the cache issue you describe will be more important.

---

<div class="post-metadata">

**Author:** ![elenev](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/elenev/32/18440_2.png) [@elenev](https://discourse.julialang.org/u/elenev)\
**Post date:** [June 12, 2020, 3:56am UTC](https://discourse.julialang.org/t/why-are-my-linear-interpolations-10x-faster-in-matlab/41232/13 "2020-06-12T03:56:29Z")

</div>

Man, this is a powerful demonstration of how important it is to store the dimensions in a type rather than have it be determined at runtime. Since your code has a couple of different optimizations (including the mutating `rand`, which I hadn’t thought of), I decided to go from my code to yours step by step. I also added a line inside each function to fix the seed so that I can make sure all of them produce the same results (they do).

```julia
using Interpolations, BenchmarkTools, Random
f(a,b) = a^2 + b^2
const xs = range(0,5,length=100)
const ys = range(-3,3,length=200)
const vals = [f(x,y) for x in xs, y in ys]
const fct = LinearInterpolation( (xs,ys), vals)
const Npt = 10_000;
pts = [rand(1,Npt); zeros(1,Npt)]

function evaluate!(pts,fct)
    Random.seed!(100)
    Niter = 1_000
    nd = size(pts,1)
    result = Matrix{Float64}(undef,size(pts,2),Niter)
    @inbounds for j in 1:Niter
        pts[2,:] .= rand(Npt)*6 .- 3
        for p in axes(pts,2)
            result[p,j] = fct(pts[:,p]...)
            #result[p,j] = fct(pts[1,p],pts[2,p])
        end
    end
    return result
end

function evaluate2!(pts::NTuple{N,<:AbstractVector},fct) where {N}
    Random.seed!(100)
    Niter = 1_000
    Npt=length(pts[2])
    result = Matrix{Float64}(undef,Npt,Niter)
    @inbounds for j in 1:Niter
        pts[2] .= 6rand(Npt) .- 3
        for p in 1:Npt
            result[p,j] = fct(ntuple(n -> pts[n][p], Val(N))...)
        end
    end
    return result
end

function evaluate3!(pts::NTuple{N,<:AbstractVector},fct) where {N}
    Random.seed!(100)
    Niter = 1_000
    pts2=pts[2]
    result = Matrix{Float64}(undef,length(pts2),Niter)
    @inbounds for j in 1:Niter
        rand!(pts2)
        for p in axes(pts2,1)
            pts2[p] = 6pts2[p] - 3
            result[p,j] = fct(ntuple(n -> pts[n][p], Val(N))...)
        end
    end
    return result
end

function evaluateElrod!(fct, pts::Vararg{<:AbstractVector,N}) where {N}
    Random.seed!(100)
    Niter = 1_000
    pts2 = pts[2]
    result = Matrix{Float64}(undef,length(pts2),Niter)
    @inbounds for j in 1:Niter
        rand!(pts2)
        # Need this loop b/c unlike Matlab fct here requires scalar args
        for p in axes(pts2,1)
            pts2[p] = 6pts2[p] - 3
            result[p,j] = fct(ntuple(n -> pts[n][p], Val(N))...)
        end
    end
    return result
end

```

1. `evaluate!` is my original function.
2. `evaluate2!` replaces the matrix with a tuple of vectors and splats a tuple instead of of a SubArray (btw, I also tried a vector of vectors instead of a tuple of vectors, and that was 10x slower than my original code, so I’m not even showing that result)
3. `evaluate3!` is `evaluate2!` but with the mutating `rand()`
4. `evaluateElrod!` is @Elrod 's function, which is like `evaluate3!` except it collects dimensions using `Vararg` instead of just taking a tuple as an argument.

Here are the results:

```julia
julia> @benchmark evaluate!($pts,fct)
BenchmarkTools.Trial:
  memory estimate: 1.56 GiB
  allocs estimate: 40004004
  --------------
  minimum time: 1.961 s (6.56% GC)
  median time: 2.023 s (15.61% GC)
  mean time: 2.095 s (13.33% GC)
  maximum time: 2.303 s (17.07% GC)
  --------------
  samples: 3
  evals/sample: 1

```

```julia
julia> @benchmark evaluate2!(ntuple(n->$pts[n,:],2),fct)
@show benchres3
benchres3 = Trial(391.896 ms)
BenchmarkTools.Trial:
  memory estimate: 229.19 MiB
  allocs estimate: 4010
  --------------
  minimum time: 391.896 ms (3.32% GC)
  median time: 396.804 ms (3.35% GC)
  mean time: 414.187 ms (7.33% GC)
  maximum time: 586.854 ms (36.01% GC)
  --------------
  samples: 13
  evals/sample: 1

```

```julia
julia> @benchmark evaluate3!(ntuple(n->$pts[n,:],2),fct)
BenchmarkTools.Trial:
  memory estimate: 76.45 MiB
  allocs estimate: 10
  --------------
  minimum time: 337.228 ms (0.00% GC)
  median time: 350.747 ms (1.46% GC)
  mean time: 361.871 ms (5.70% GC)
  maximum time: 460.437 ms (26.81% GC)
  --------------
  samples: 14
  evals/sample: 1

```

```julia
julia> @benchmark evaluateElrod!(fct,$pts[1,:],$pts[2,:])
BenchmarkTools.Trial:
  memory estimate: 76.45 MiB
  allocs estimate: 8
  --------------
  minimum time: 342.385 ms (0.18% GC)
  median time: 352.621 ms (1.47% GC)
  mean time: 366.353 ms (5.62% GC)
  maximum time: 468.002 ms (26.52% GC)
  --------------
  samples: 14
  evals/sample: 1

```

The extra allocations in my benchmark of @Elrod 's code relative to his are due to (1) slicing `pts` inside the benchmark as opposed to outside (4 allocs), and (2) fixing the seed (2 allocs). These extra allocations do not affect the speed by more than a few milliseconds i.e. within a 95% CI of the mean.

Upshot: by far the most important improvement is from `evaluate!` to `evaluate2!`, where the matrix is replaced with the tuple of vectors. This is all about `N` being in the type declaration. As I mentioned, if I replace the tuple with a vector of vectors and `Val(N)` with `Val(ND)` where `ND = length(pts)` somewhere outside the loops, that takes 20 seconds! But given `N` in the type, it doesn’t matter much whether there are variable number of arguments or a tuple passed in. I prefer the tuple because it allows the code calling the function to also be agnostic about the number of dimensions, without splatting.

What would be really nice is if there was some cousin of `StaticArray` that would allow you to “type” all but one dimension. E.g. for an `MMatrix`, you could say that it will have 5 columns and an unspecified number of rows. Otherwise, I have to keep doing stuff like `ntuple(n->$pts[n,:],2)`, which is kinda boilerplate.

Btw, all of these are still marginally slower than MATLAB. I guess the `griddedInterpolant` object is just crazy fast. It’s definitely not written in MATLAB itself.

---

<div class="post-metadata">

**Author:** ![jling](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/jling/32/212909_2.png) [@jling](https://discourse.julialang.org/u/jling)\
**Post date:** [June 12, 2020, 4:12am UTC](https://discourse.julialang.org/t/why-are-my-linear-interpolations-10x-faster-in-matlab/41232/14 "2020-06-12T04:12:58Z")

</div>

I suspect someone can make this faster but interestingly, nobody knows what it’s written in / what algorithm is used:

> **[what is the source code of gridinterpolant function? what is the al...](https://www.mathworks.com/matlabcentral/answers/485336-what-is-the-source-code-of-gridinterpolant-function-what-is-the-algorithm)**
>
> I want to know what is behind griddedinterpolant function?
> what is the original code? can I know? what is the algorithm?
> best regards
> Behzad- Nashville, Tennessee

---

<div class="post-metadata">

**Author:** ![baggepinnen](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/baggepinnen/32/693_2.png) [@baggepinnen](https://discourse.julialang.org/u/baggepinnen)\
**Post date:** [June 12, 2020, 5:50am UTC](https://discourse.julialang.org/t/why-are-my-linear-interpolations-10x-faster-in-matlab/41232/15 "2020-06-12T05:50:45Z")

</div>

> [@elenev](#):
>
> What would be really nice is if there was some cousin of `StaticArray` that would allow you to “type” all but one dimension.

> **[GitHub - JuliaArrays/HybridArrays.jl: Arrays with both statically and...](https://github.com/JuliaArrays/HybridArrays.jl)**
>
> Arrays with both statically and dynamically sized axes in Julia - GitHub - JuliaArrays/HybridArrays.jl: Arrays with both statically and dynamically sized axes in Julia

---

<div class="post-metadata">

**Author:** ![giordano](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/giordano/32/2166_2.png) [@giordano](https://discourse.julialang.org/u/giordano)\
**Post date:** [June 12, 2020, 9:57am UTC](https://discourse.julialang.org/t/why-are-my-linear-interpolations-10x-faster-in-matlab/41232/16 "2020-06-12T09:57:16Z")

</div>

> [@elenev](#):
>
> ```julia
> pts[2,:] .= rand(Npt)*6 .- 3
> 
> ```

Small nitpicking: this line in `evaluate!` can be replaced by

```julia
pts[2,:] .= rand.() .* 6 .- 3

```

which is equivalent (there is a new call to `rand()` in every iteration) and doesn’t allocate an `Npt`-element array in each iteration. This doesn’t speedup much your code, just saves a few allocations, but it’s a good trick to keep in mind in general

---

<div class="post-metadata">

**Author:** ![Elrod](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/elrod/32/22461_2.png) [@Elrod](https://discourse.julialang.org/u/Elrod)\
**Post date:** [June 12, 2020, 2:05pm UTC](https://discourse.julialang.org/t/why-are-my-linear-interpolations-10x-faster-in-matlab/41232/17 "2020-06-12T14:05:23Z")

</div>

> [@elenev](#):
>
> Btw, all of these are still marginally slower than MATLAB. I guess the `griddedInterpolant` object is just crazy fast. It’s definitely not written in MATLAB itself.

I’ve heard that broadcasting is multithreaded in MATLAB by default. Any chance this is?  
Either way, we could multithread it in Julia:

```julia
julia> using Interpolations, Random, BenchmarkTools

julia> f(a,b) = a^2 + b^2
f (generic function with 1 method)

julia> const xs = range(0,5,length=100);

julia> const ys = range(-3,3,length=200);

julia> const vals = [f(x,y) for x in xs, y in ys];

julia> const fct = LinearInterpolation( (xs,ys), vals);

julia> const Npt = 10_000;

julia> pts1 = rand(Npt);

julia> pts2 = zeros(Npt);

julia> function evaluateElrod!(fct, pts::Vararg{<:AbstractVector,N}) where {N}
           Random.seed!(100)
           Niter = 1_000
           pts2 = pts[2]
           result = Matrix{Float64}(undef,length(pts2),Niter)
           @inbounds for j in 1:Niter
               rand!(pts2)
               # Need this loop b/c unlike Matlab fct here requires scalar args
               Threads.@threads for p in axes(pts2,1)
                   pts2[p] = 6pts2[p] - 3
                   result[p,j] = fct(ntuple(n -> pts[n][p], Val(N))...)
               end
           end
           return result
       end
evaluateElrod! (generic function with 1 method)

julia> @benchmark evaluateElrod!(fct,$pts1,$pts2)
BenchmarkTools.Trial:
  memory estimate: 83.62 MiB
  allocs estimate: 51517
  --------------
  minimum time: 41.556 ms (0.00% GC)
  median time: 43.659 ms (1.21% GC)
  mean time: 44.964 ms (4.54% GC)
  maximum time: 116.069 ms (60.29% GC)
  --------------
  samples: 112
  evals/sample: 1

julia> versioninfo()
Julia Version 1.6.0-DEV.205
Commit a018dc0ad5* (2020-06-10 04:25 UTC)
Platform Info:
  OS: Linux (x86_64-generic-linux)
  CPU: Intel(R) Core(TM) i9-7900X CPU @ 3.30GHz
  WORD_SIZE: 64
  LIBM: libopenlibm
  LLVM: libLLVM-10.0.0 (ORCJIT, skylake)
Environment:
  JULIA_NUM_THREADS = 10

```

I saw much better performance by setting the number of threads to equal the number of physical cores than the number of logical cores, but you could benchmark both as results may vary.

Also, the original `evaluate!` was actually about the same fast on my computer as yours, or maybe a little slower:

```julia
julia> @benchmark evaluate!($pts,fct)
BenchmarkTools.Trial:
  memory estimate: 2.76 GiB
  allocs estimate: 50004004
  --------------
  minimum time: 1.875 s (9.27% GC)
  median time: 2.223 s (13.94% GC)
  mean time: 2.149 s (12.63% GC)
  maximum time: 2.349 s (14.07% GC)
  --------------
  samples: 3
  evals/sample: 1

```

So I’d hoped you’d see a similar speed up:

```julia
julia> @benchmark evaluate3!(ntuple(n->$pts[n,:],2),fct)
BenchmarkTools.Trial:
  memory estimate: 76.45 MiB
  allocs estimate: 9
  --------------
  minimum time: 136.669 ms (0.00% GC)
  median time: 139.478 ms (0.00% GC)
  mean time: 143.595 ms (3.01% GC)
  maximum time: 208.128 ms (33.02% GC)
  --------------
  samples: 35
  evals/sample: 1

```

If you had, then that’d already have been faster than MATLAB’s time, without the need for threading.  
I do not have a MATLAB license, so I have no idea how fast MATLAB’s `griddedInterpolant` would be on this computer.

---

<div class="post-metadata">

**Author:** ![tbeason](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/tbeason/32/15898_2.png) [@tbeason](https://discourse.julialang.org/u/tbeason)\
**Post date:** [June 12, 2020, 2:33pm UTC](https://discourse.julialang.org/t/why-are-my-linear-interpolations-10x-faster-in-matlab/41232/18 "2020-06-12T14:33:25Z")

</div>

This has been a minor but still weird “pain point” of mine while moving form MATLAB to Julia. I kind of liked `griddedInterpolant` but there isn’t any direct equivalent in Julia and there isn’t enough information that exists about its source to remake it! There is a longstanding open issue in Interpolations.jl regarding cubic splines on gridded data but nothing concrete has emerged from it.

---

<div class="post-metadata">

**Author:** ![elenev](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/elenev/32/18440_2.png) [@elenev](https://discourse.julialang.org/u/elenev)\
**Post date:** [June 12, 2020, 3:31pm UTC](https://discourse.julialang.org/t/why-are-my-linear-interpolations-10x-faster-in-matlab/41232/19 "2020-06-12T15:31:54Z")

</div>

Thanks for your response.

Here’s my benchmark for the multi-threaded version of your function. Is it correct to compare Julia’s `Threads.@threads for` to MATLAB’s `parfor`?

```julia
BenchmarkTools.Trial: 
  memory estimate: 79.75 MiB
  allocs estimate: 23664
  --------------
  minimum time: 137.605 ms (0.00% GC)
  median time: 166.774 ms (0.35% GC)
  mean time: 175.249 ms (9.40% GC)
  maximum time: 301.986 ms (50.10% GC)
  --------------
  samples: 29
  evals/sample: 1

```

But our setups are quite different, so these numbers should be compared with my earlier ones, not with yours.

```julia
Julia Version 1.4.2
Commit 44fa15b150* (2020-05-23 18:35 UTC)
Platform Info:
  OS: Windows (x86_64-w64-mingw32)
  CPU: Intel(R) Xeon(R) CPU E3-1240 v3 @ 3.40GHz
  WORD_SIZE: 64
  LIBM: libopenlibm
  LLVM: libLLVM-8.0.1 (ORCJIT, haswell)
Environment:
  JULIA_NUM_THREADS = 4 #physical cores

```

If MATLAB default is indeed to multi-thread vector query (“evaluate”) calls, how would that work if the queries themselves happen in a multi-threaded loop e.g. in `parfor`? I guess I don’t actually care that much about MATLAB’s internals – I care about what would happen if I did this in Julia i.e. if I nested multithreaded loops (across a function barrier, presumably).

---

<div class="post-metadata">

**Author:** ![elenev](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/elenev/32/18440_2.png) [@elenev](https://discourse.julialang.org/u/elenev)\
**Post date:** [June 12, 2020, 3:39pm UTC](https://discourse.julialang.org/t/why-are-my-linear-interpolations-10x-faster-in-matlab/41232/20 "2020-06-12T15:39:07Z")

</div>

Re irregular grid splines: We came across this a few years ago to speed up our gridded spline interpolations in MATLAB.

> **[Extrimely fast general n-dimensional interpolators](https://www.mathworks.com/matlabcentral/fileexchange/30836-extrimely-fast-general-n-dimensional-interpolators)**
>
> Mex C++ interpolator routines for general pp-forms in any dimension. Multithreaded.

I have basically no knowledge of C or C++ so I don’t think I’d be able to adapt this to Julia myself, but if lack of a good algorithm is what’s holding folks back, maybe this helps. This is an order of magnitude faster than the built-in MATLAB functions it’s replacing.

[Next page](https://discourse.julialang.org/t/why-are-my-linear-interpolations-10x-faster-in-matlab/41232.md?page=2)