# Efficient finite difference operators

**URL:** <https://discourse.julialang.org/t/efficient-finite-difference-operators/12439>\
**Category:** Performance\
**Tags:** diffeq\
**Created:** [July 17, 2018, 3:47pm UTC](https://discourse.julialang.org/t/efficient-finite-difference-operators/12439 "2018-07-17T15:47:56Z")\
**Posts on this page:** 20\
**Page:** 1

<div class="post-metadata">

**Author:** ![milankl](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/milankl/32/4198_2.png) [@milankl](https://discourse.julialang.org/u/milankl)\
**Post date:** [July 17, 2018, 3:47pm UTC](https://discourse.julialang.org/t/efficient-finite-difference-operators/12439/1 "2018-07-17T15:47:56Z")

</div>

Hey there,

I have a model that spends most time in functions like this

```julia
function Gux!(dudx::Matrix{Numtype},u::Matrix{Numtype})

    dudx[1:end-1,:] = u[2:end,:]-u[1:end-1,:]
    dudx[end,:] = u[1,:]-u[end,:]`

end

```

that will be evaluated over and over again with changing input matrices u. Preallocating the result dudx and reusing it already gives a speed advantage. I am currently testing a couple of other ways to write essentially the same operation, but I was wondering whether there are any performance tips that are obvious for someone with more insight into Julia. For example, does the order of the lines inside the function matter (regarding reading and writing from memory)? Is Matrix{Numtype} where Numtype could be Float32, Float64 or BigFloat etc. the way to give the compiler enough information what input arguments to expect? The size of u will not change throughtout a computation, hence would it be advantageous to also pass on that information? Furthermore, is the matrix slicing that I use to write, as I find it more convenient than writing loops, a bottleneck? I am aware of the column-major order for instance but not sure whether writing with matrix slices always respects that.

Thanks for any hints!

---

<div class="post-metadata">

**Author:** ![Evizero](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/evizero/32/10118_2.png) [@Evizero](https://discourse.julialang.org/u/Evizero)\
**Post date:** [July 17, 2018, 4:03pm UTC](https://discourse.julialang.org/t/efficient-finite-difference-operators/12439/2 "2018-07-17T16:03:12Z")

</div>

you could try the following

```julia
julia> function Gux2!(dudx::Matrix{T},u::Matrix{T}) where T
           @views dudx[1:end-1,:] .= u[2:end,:] .- u[1:end-1,:]
           @views dudx[end,:] .= u[1,:] .- u[end,:]
           dudx
       end

```

```julia
julia> u = rand(100,100);

julia> du = zeros(u);

julia> @btime Gux!($du, $u);
  26.539 μs (18 allocations: 235.17 KiB)

julia> @btime Gux2!($du, $u);
  3.274 μs (2 allocations: 112 bytes)

```

---

<div class="post-metadata">

**Author:** ![Evizero](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/evizero/32/10118_2.png) [@Evizero](https://discourse.julialang.org/u/Evizero)\
**Post date:** [July 17, 2018, 4:21pm UTC](https://discourse.julialang.org/t/efficient-finite-difference-operators/12439/3 "2018-07-17T16:21:32Z")

</div>

> [@milankl](#):
>
> Is Matrix{Numtype} where Numtype could be Float32, Float64 or BigFloat etc. the way to give the compiler enough information what input arguments to expect?

Don’t think of the types you specify in the method signature as “helping the compiler” in terms of performance - which AFAIK it doesn’t. Think of it as “restricting this specific method to just the specified subset of types”.

For example the solution I posted above is probably overly restrictive in its signature. As it is written, the method will only actually be called when

1. both `dudx` and `u` are actually a `Matrix` (so no `SubArray` etc), and
2. both `dudx` and `u` have the same `eltype` (which in this case is probably expected, but not actually necesarry for the code to work)

everything else will simply throw a `MethodError`.

So what would one typically do? Since the code itself assumes that both parameters are matrices, I’d probably just write the signature as `Gux!(dudx::AbstractMatrix, u::AbstractMatrix)`.

---

<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:** [July 17, 2018, 4:32pm UTC](https://discourse.julialang.org/t/efficient-finite-difference-operators/12439/4 "2018-07-17T16:32:19Z")

</div>

> Don’t think of the types you specify in the method signature as “helping the compiler” in terms of performance - which AFAIK it doesn’t.

That’s correct.  
If you call `Gux!` with Numtype = Float32, again with Numtype = Float64, and then again with Numtype = BigFloat, the compiler will have created three separate versions of the function – one for when Numtype == Float32, another for when Numtype == Float64, another…

Whenever you call a function with a different combination of input types, it under the hood creates a version specialized on that combination of input types.  
Applying restrictions just forbids it from accepting certain things.

Of course, if the algorithm or implementation changes based on input types, then annotations are useful.

EDIT:  
Also, unfortunately this is still a case where for loops are fastest:

```julia
julia> @btime Gux2!($du, $u);
  3.014 μs (6 allocations: 336 bytes)

julia> function Gux3!(dudx::AbstractMatrix,u::AbstractMatrix)
           m, n = size(dudx)
           @boundscheck (m,n) == size(u) || throw(BoundsError())
       
           @inbounds for i ∈ 1:n
               for j ∈ 1:m-1
                   dudx[j,i] = u[j+1,i]-u[j,i]
               end
               dudx[m,i] = u[1,i] - u[m,i]
           end
       
       
       end
Gux3! (generic function with 1 method)

julia> @btime Gux3!($du, $u);
  2.177 μs (0 allocations: 0 bytes)

```

The difference was smaller on Julia 0.6:

```julia
julia> @btime Gux2!($du, $u);
  2.957 μs (2 allocations: 112 bytes)

julia> @btime Gux3!($du, $u);
  2.472 μs (0 allocations: 0 bytes)

```

(I think the broadcasting regression we see here will be fixed before 0.7 is released.)

---

<div class="post-metadata">

**Author:** ![milankl](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/milankl/32/4198_2.png) [@milankl](https://discourse.julialang.org/u/milankl)\
**Post date:** [July 17, 2018, 4:44pm UTC](https://discourse.julialang.org/t/efficient-finite-difference-operators/12439/5 "2018-07-17T16:44:39Z")

</div>

Great. Thanks for the explanation! Is there any general advice whether writing the same operation as a loop is beneficial? I know that is says somewhere in the manual that you should only use matrix operations when it feels natural … whatever that means.

---

<div class="post-metadata">

**Author:** ![milankl](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/milankl/32/4198_2.png) [@milankl](https://discourse.julialang.org/u/milankl)\
**Post date:** [July 17, 2018, 4:55pm UTC](https://discourse.julialang.org/t/efficient-finite-difference-operators/12439/6 "2018-07-17T16:55:59Z")

</div>

Thanks so much guys! That saves me a lot of computing time 😉

---

<div class="post-metadata">

**Author:** ![ChrisRackauckas](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/chrisrackauckas/32/77_2.png) [@ChrisRackauckas](https://discourse.julialang.org/u/ChrisRackauckas)\
**Post date:** [July 17, 2018, 5:00pm UTC](https://discourse.julialang.org/t/efficient-finite-difference-operators/12439/7 "2018-07-17T17:00:20Z")

</div>

I just want to point out that DiffEqOperators.jl has matrix-free operators that make `*` or `mul!` write this loop for you.

[https://github.com/JuliaDiffEq/DiffEqOperators.jl/blob/master/docs/HeatEquation.md](https://github.com/JuliaDiffEq/DiffEqOperators.jl/blob/master/docs/HeatEquation.md)

I have found the full loop to be better, as seen in this tutorial notebook:

[http://nbviewer.jupyter.org/github/JuliaDiffEq/DiffEqTutorials.jl/blob/master/Introduction/OptimizingDiffEqCode.ipynb#Optimizing-Large-Systems](http://nbviewer.jupyter.org/github/JuliaDiffEq/DiffEqTutorials.jl/blob/master/Introduction/OptimizingDiffEqCode.ipynb#Optimizing-Large-Systems)

so DiffEqOperators basically tries to make those.

---

<div class="post-metadata">

**Author:** ![milankl](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/milankl/32/4198_2.png) [@milankl](https://discourse.julialang.org/u/milankl)\
**Post date:** [July 17, 2018, 5:11pm UTC](https://discourse.julialang.org/t/efficient-finite-difference-operators/12439/8 "2018-07-17T17:11:55Z")

</div>

Yeah, I heard of this package, however, I will need a couple of other unusual stencils with unusual datatypes without promotion from one type to the other. I.e. something like `0.5*u[i,j]` needs to be written as `one_half*u[i,j]` where the type of one\_half corresponds to the type in `u` etc… Thanks for mentioning it though!

---

<div class="post-metadata">

**Author:** ![Evizero](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/evizero/32/10118_2.png) [@Evizero](https://discourse.julialang.org/u/Evizero)\
**Post date:** [July 17, 2018, 5:14pm UTC](https://discourse.julialang.org/t/efficient-finite-difference-operators/12439/9 "2018-07-17T17:14:10Z")

</div>

> [@milankl](#):
>
> Is there any general advice whether writing the same operation as a loop is beneficial? I know that is says somewhere in the manual that you should only use matrix operations when it feels natural … whatever that means.

Personally I think the best advice here is to get into to the habit of benchmarking your code to build up your own intuition first hand. Using [GitHub - JuliaCI/BenchmarkTools.jl: A benchmarking framework for the Julia language](https://github.com/JuliaCI/BenchmarkTools.jl) has helped me countless times and continues to do so even after years of coding in julia

---

<div class="post-metadata">

**Author:** ![ChrisRackauckas](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/chrisrackauckas/32/77_2.png) [@ChrisRackauckas](https://discourse.julialang.org/u/ChrisRackauckas)\
**Post date:** [July 17, 2018, 5:18pm UTC](https://discourse.julialang.org/t/efficient-finite-difference-operators/12439/10 "2018-07-17T17:18:13Z")

</div>

> [@milankl](#):
>
> Yeah, I heard of this package, however, I will need a couple of other unusual stencils with unusual datatypes without promotion from one type to the other. I.e. something like `0.5*u[i,j]` needs to be written as `one_half*u[i,j]` where the type of one\_half corresponds to the type in `u` etc… Thanks for mentioning it though!

Feel free to open an issue on this. We are still actively developing it so it would be great to know the use cases.

---

<div class="post-metadata">

**Author:** ![dpsanders](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/dpsanders/32/3573_2.png) [@dpsanders](https://discourse.julialang.org/u/dpsanders)\
**Post date:** [July 17, 2018, 9:19pm UTC](https://discourse.julialang.org/t/efficient-finite-difference-operators/12439/11 "2018-07-17T21:19:12Z")

</div>

The explicit loop is much more readable IMO.

---

<div class="post-metadata">

**Author:** ![milankl](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/milankl/32/4198_2.png) [@milankl](https://discourse.julialang.org/u/milankl)\
**Post date:** [August 7, 2018, 11:49am UTC](https://discourse.julialang.org/t/efficient-finite-difference-operators/12439/12 "2018-08-07T11:49:08Z")

</div>

As a follow up on that topic: Once there is a multiplication with a constant involved, the matrix version is 3x faster. Or am I missing here some obvious optimization?

```julia
function Ix!(ux::AbstractMatrix,u::AbstractMatrix)
    
    m, n = size(ux)
    @boundscheck (m+1,n) == size(u) || throw(BoundsError())

    @inbounds for i ∈ 1:n
        for j ∈ 1:m
            ux[j,i] .= 0.5*(u[j+1,i] .+ u[j,i])
        end
    end
end

function Ix2!(ux::AbstractMatrix,u::AbstractMatrix)
    
    m, n = size(ux)
    @boundscheck (m+1,n) == size(u) || throw(BoundsError())

    @inbounds @views ux[:,:] .= 0.5*(u[2:end,:] .+ u[1:end-1,:])
end

```

```julia
julia> u = rand(500,500);

julia> ux = zeros(499,500);

julia> @btime Ix!(ux,u);
  2.064 ms (249500 allocations: 11.42 MiB)

julia> @btime Ix2!(ux,u);
  594.146 μs (5 allocations: 3.81 MiB)

```

---

<div class="post-metadata">

**Author:** ![rveltz](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/rveltz/32/2707_2.png) [@rveltz](https://discourse.julialang.org/u/rveltz)\
**Post date:** [August 7, 2018, 1:13pm UTC](https://discourse.julialang.org/t/efficient-finite-difference-operators/12439/13 "2018-08-07T13:13:12Z")

</div>

you dont need the dot here

```julia
ux[j,i] .=

```

I then get

```julia
julia> @btime Ix!(ux,u);
  80.177 μs (0 allocations: 0 bytes)

```

---

<div class="post-metadata">

**Author:** ![kristoffer.carlsson](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/kristoffer.carlsson/32/22_2.png) [@kristoffer.carlsson](https://discourse.julialang.org/u/kristoffer.carlsson)\
**Post date:** [August 7, 2018, 1:16pm UTC](https://discourse.julialang.org/t/efficient-finite-difference-operators/12439/14 "2018-08-07T13:16:52Z")

</div>

Should also use `0.5 .* ... `

---

<div class="post-metadata">

**Author:** ![RoyiAvital](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/royiavital/32/571_2.png) [@RoyiAvital](https://discourse.julialang.org/u/RoyiAvital)\
**Post date:** [August 7, 2018, 1:18pm UTC](https://discourse.julialang.org/t/efficient-finite-difference-operators/12439/15 "2018-08-07T13:18:05Z")

</div>

If the `.=` and `.+` creates such an overhead it is really starnge.  
Could you show the timing results of all 3 (The 2 codes written above by @milankl and yours without the dots)?

---

<div class="post-metadata">

**Author:** ![RoyiAvital](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/royiavital/32/571_2.png) [@RoyiAvital](https://discourse.julialang.org/u/RoyiAvital)\
**Post date:** [August 7, 2018, 1:28pm UTC](https://discourse.julialang.org/t/efficient-finite-difference-operators/12439/17 "2018-08-07T13:28:23Z")

</div>

Could you show the timing for with and without the dots?

---

<div class="post-metadata">

**Author:** ![kristoffer.carlsson](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/kristoffer.carlsson/32/22_2.png) [@kristoffer.carlsson](https://discourse.julialang.org/u/kristoffer.carlsson)\
**Post date:** [August 7, 2018, 1:28pm UTC](https://discourse.julialang.org/t/efficient-finite-difference-operators/12439/18 "2018-08-07T13:28:25Z")

</div>

I meant on the second version to have it fuse.

---

<div class="post-metadata">

**Author:** ![milankl](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/milankl/32/4198_2.png) [@milankl](https://discourse.julialang.org/u/milankl)\
**Post date:** [August 7, 2018, 4:39pm UTC](https://discourse.julialang.org/t/efficient-finite-difference-operators/12439/19 "2018-08-07T16:39:36Z")

</div>

A loop version, a matrix version, and a loop version with .=

```julia
function Ix_loop!(ux::AbstractMatrix,u::AbstractMatrix)

    m, n = size(ux)
    @boundscheck (m+1,n) == size(u) || throw(BoundsError())

    @inbounds for i ∈ 1:n
        for j ∈ 1:m
            ux[j,i] = 0.5*(u[j+1,i] + u[j,i])
        end
    end
end

function Ix_mat!(ux::AbstractMatrix,u::AbstractMatrix)

    m, n = size(ux)
    @boundscheck (m+1,n) == size(u) || throw(BoundsError())

    @inbounds @views ux[:,:] .= 0.5*(u[2:end,:] .+ u[1:end-1,:])
end

function Ix_loop_dot!(ux::AbstractMatrix,u::AbstractMatrix)

    m, n = size(ux)
    @boundscheck (m+1,n) == size(u) || throw(BoundsError())

    @inbounds for i ∈ 1:n
        for j ∈ 1:m
            ux[j,i] .= 0.5*(u[j+1,i] + u[j,i])
        end
    end
end

```

and the timings are

```julia
julia> @btime Ix_loop!($ux,$u);
  85.346 μs (0 allocations: 0 bytes)

julia> @btime Ix_mat!($ux,$u);
  565.237 μs (5 allocations: 3.81 MiB)

julia> @btime Ix_loop_dot!($ux,$u);
  2.108 ms (249500 allocations: 11.42 MiB)

```

u and ux are as above.

---

<div class="post-metadata">

**Author:** ![milankl](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/milankl/32/4198_2.png) [@milankl](https://discourse.julialang.org/u/milankl)\
**Post date:** [August 7, 2018, 4:45pm UTC](https://discourse.julialang.org/t/efficient-finite-difference-operators/12439/20 "2018-08-07T16:45:29Z")

</div>

And probably the worst is if you put dots everywhere

```julia
function Ix_loop_dotdotdot!(ux::AbstractMatrix,u::AbstractMatrix)

    m, n = size(ux)
    @boundscheck (m+1,n) == size(u) || throw(BoundsError())

    @inbounds for i ∈ 1:n
        for j ∈ 1:m
            ux[j,i] .= 0.5.*(u[j+1,i] .+ u[j,i])
        end
    end
end

```

```julia
julia> @btime Ix_loop_dotdotdot!($ux,$u);
  55.061 ms (1247500 allocations: 26.65 MiB)

```

---

<div class="post-metadata">

**Author:** ![RoyiAvital](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/royiavital/32/571_2.png) [@RoyiAvital](https://discourse.julialang.org/u/RoyiAvital)\
**Post date:** [August 8, 2018, 4:17am UTC](https://discourse.julialang.org/t/efficient-finite-difference-operators/12439/21 "2018-08-08T04:17:50Z")

</div>

Could anyone say why the dots add such an overhead?