# Acceleration of a function

**URL:** <https://discourse.julialang.org/t/acceleration-of-a-function/50181>\
**Category:** Performance\
**Tags:** question\
**Created:** [November 15, 2020, 11:35am UTC](https://discourse.julialang.org/t/acceleration-of-a-function/50181 "2020-11-15T11:35:22Z")\
**Posts on this page:** 17\
**Page:** 1

<div class="post-metadata">

**Author:** ![Jonas](https://avatars.discourse-cdn.com/v4/letter/j/aeb1de/32.png) [@Jonas](https://discourse.julialang.org/u/Jonas)\
**Post date:** [November 15, 2020, 11:35am UTC](https://discourse.julialang.org/t/acceleration-of-a-function/50181/1 "2020-11-15T11:35:22Z")

</div>

Hello everyone,

I am currently transforming some Matlab code into Julia. Everything works fine, but I recognized a big difference in runtime. I don’t have much experience in that regard, so I started reading the guidelines of Julia. Further, I know that Matlab is highly parallelized, which I plan for my Julia code to some extent, too. But first I want my cde to run as fast as possible without parallelization.  
In my code, the following function is used repeatedly and it takes rather long time to calculate (around 36 ms, evaluated with @btime)

```julia
gsm_cascade(Sc::Array{Array{Complex{Float64},2},2}, S0::Array{Array{Complex{Float64},2},2}, S1::Array{Array{Complex{Float64},2},2}, D::Diagonal{Complex{Float64},Array{Complex{Float64},1}})

     S111 = S1[1,1]
     S121 = S1[1,2]
     S211 = S1[2,1]
     S221 = S1[2,2]

     S110 = S0[1,1]
     S120 = S0[1,2]
     S210 = S0[2,1]
     S220 = S0[2,2]

     I = sparse(1.0LinearAlgebra.I, size(D,1), size(D,1))
     U1 = inv(I - S221*D*S110*D))
     U2 = inv(I - S110*D*S221*D))

     Sc[1,1] = S111 + S121*D*U2*S110*D*S211
     Sc[1,2] = S121*D*U2*S120
     Sc[2,1] = S210*D*U1*S211
     Sc[2,2] = S220 + S210*D*U1*S221*D*S120

     return Sc

end

```

For explanation, I have two Matrices S0 and S1, which are used to get matrix Sc. Since the matrices are constructed by four submatrices, I decided to use an Array of Arrays to represent each of S0, S1 and Sc, in order to keep track of the submatrices.  
At first, I thought that taking the inverses is the most time consuming part, but actually the matrix multiplications and additions of

```julia
     Sc[1,1] = S111 + S121*D*U2*S110*D*S211
     Sc[1,2] = S121*D*U2*S120
     Sc[2,1] = S210*D*U1*S211
     Sc[2,2] = S220 + S210*D*U1*S221*D*S120

```

are the slow part of the function. Does anybody know how I can accelerate this function? The biggest matrix is S111, which is currently 32 by 32. Or would be using parallelized code the only option here?

Further, I recognized that

```julia
S111 = S1[1,1]

```

does not allocate anything on the heap . Why is that the case? Is it because the matrix S1 is aleady allocated and fed into the function?

Thank you in advance and have a nice day 😃  
Jonas

---

<div class="post-metadata">

**Author:** ![antoine-levitt](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/antoine-levitt/32/4008_2.png) [@antoine-levitt](https://discourse.julialang.org/u/antoine-levitt)\
**Post date:** [November 15, 2020, 12:51pm UTC](https://discourse.julialang.org/t/acceleration-of-a-function/50181/2 "2020-11-15T12:51:59Z")

</div>

> [@Jonas](#):
>
> Further, I recognized that
> 
> ```julia
> S111 = S1[1,1]
> 
> ```
> 
> does not allocate anything on the heap

This is just a pointer assignment: S1 is a matrix of matrices, stored in practice as four pointers.

- that line `I = sparse(1.0LinearAlgebra.I, size(D,1), size(D,1))` is unnecessary I think; just using `I` is fine
- you probably shouldn’t use `inv` but rather factor the matrix; probably not going to do much for speed if that’s not the bottleneck though.
- Lines like `S111 + S121*D*U2*S110*D*S211` allocate a lot of temporary arrays. That might be the problem here. You can use in-place operations (`mul!`, dot assignment and fusion)
- If your arrays are small, look up StaticArrays.jl.
- If you want to compare against matlab be careful that matlab uses MKL, while julia uses OpenBLAS by default. You can make julia use MKL with the MKL.jl package.
- You don’t need to type your input arguments
- In all cases, profile to see what takes time and what impact your modifications have

---

<div class="post-metadata">

**Author:** ![antoine-levitt](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/antoine-levitt/32/4008_2.png) [@antoine-levitt](https://discourse.julialang.org/u/antoine-levitt)\
**Post date:** [November 15, 2020, 1:03pm UTC](https://discourse.julialang.org/t/acceleration-of-a-function/50181/3 "2020-11-15T13:03:50Z")

</div>

> [@antoine-levitt](#):
>
> - Lines like `Sc[1,1] = S111 + S121*D*U2*S110*D*S211` allocate a lot of temporary arrays. That might be the problem here. You can use in-place operations ( `mul!` , dot assignment and fusion)

If that is indeed the case the manual version (5 mul! and one add) is pretty annoying to write down. It would be cool if there was a package that could do figure it out by itself in the general case where matrices might not be of the same size.

---

<div class="post-metadata">

**Author:** ![Jonas](https://avatars.discourse-cdn.com/v4/letter/j/aeb1de/32.png) [@Jonas](https://discourse.julialang.org/u/Jonas)\
**Post date:** [November 15, 2020, 1:58pm UTC](https://discourse.julialang.org/t/acceleration-of-a-function/50181/4 "2020-11-15T13:58:26Z")

</div>

Hello Mr. Levitt,

thank you very much for your answer. I will look into everything you said. Since my code is dominated by matrix multiplications, maybe MKL.jl is the better choice for me.

What do you mean by " You don’t need to type your input arguments"? 😃

Best regards,  
Jonas

---

<div class="post-metadata">

**Author:** ![antoine-levitt](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/antoine-levitt/32/4008_2.png) [@antoine-levitt](https://discourse.julialang.org/u/antoine-levitt)\
**Post date:** [November 15, 2020, 2:05pm UTC](https://discourse.julialang.org/t/acceleration-of-a-function/50181/5 "2020-11-15T14:05:27Z")

</div>

This

```julia
gsm_cascade(Sc::Array{Array{Complex{Float64},2},2}, S0::Array{Array{Complex{Float64},2},2}, S1::Array{Array{Complex{Float64},2},2}, D::Diagonal{Complex{Float64},Array{Complex{Float64},1}})

```

can just be replaced by `gsm_cascade(Sc, S0, S1, D)`. In Julia type annotations are used for dispatch and/or clarity, there’s no performance impact.

---

<div class="post-metadata">

**Author:** ![Jonas](https://avatars.discourse-cdn.com/v4/letter/j/aeb1de/32.png) [@Jonas](https://discourse.julialang.org/u/Jonas)\
**Post date:** [November 15, 2020, 4:13pm UTC](https://discourse.julialang.org/t/acceleration-of-a-function/50181/6 "2020-11-15T16:13:46Z")

</div>

If I call the functions I don’t use the type annotations. But if I define a function for example in a Module, aren’t the type annotations somehow important?

By the way, I installed MKL.jl but unfortunalely the package Arpack doen’t work anymore for solving needed generalized eigenvalue problems. So I have to find a solution for this now haha.

---

<div class="post-metadata">

**Author:** ![DNF](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/dnf/32/10191_2.png) [@DNF](https://discourse.julialang.org/u/DNF)\
**Post date:** [November 15, 2020, 4:20pm UTC](https://discourse.julialang.org/t/acceleration-of-a-function/50181/7 "2020-11-15T16:20:59Z")

</div>

> [@Jonas](#):
>
> aren’t the type annotations somehow important?

Nope. They can be useful, but you mostly need them if you want to make several versions of the same function (also called ‘methods’).

---

<div class="post-metadata">

**Author:** ![Jonas](https://avatars.discourse-cdn.com/v4/letter/j/aeb1de/32.png) [@Jonas](https://discourse.julialang.org/u/Jonas)\
**Post date:** [November 15, 2020, 4:27pm UTC](https://discourse.julialang.org/t/acceleration-of-a-function/50181/8 "2020-11-15T16:27:24Z")

</div>

Ok thank you for your reply. But if the annotations aren’t destructive in any way, I think I will leave them for clarity 😃

---

<div class="post-metadata">

**Author:** ![DNF](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/dnf/32/10191_2.png) [@DNF](https://discourse.julialang.org/u/DNF)\
**Post date:** [November 15, 2020, 4:56pm UTC](https://discourse.julialang.org/t/acceleration-of-a-function/50181/9 "2020-11-15T16:56:19Z")

</div>

Well, it depends what you mean by destructive. Those annotations mean that your code will fail for real-valued matrices, or if they are Diagonal, Hermitian, SharedArrays, StaticArrays, etc etc. Every array that is not strictly of type Array, and with precisely 64 bit floating point complex elements, which is _very_ restrictive.

It also makes the function signature very verbose, which does not necessarily help clarity.

---

<div class="post-metadata">

**Author:** ![antoine-levitt](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/antoine-levitt/32/4008_2.png) [@antoine-levitt](https://discourse.julialang.org/u/antoine-levitt)\
**Post date:** [November 15, 2020, 6:04pm UTC](https://discourse.julialang.org/t/acceleration-of-a-function/50181/10 "2020-11-15T18:04:32Z")

</div>

> [@Jonas](#):
>
> By the way, I installed MKL.jl but unfortunalely the package Arpack doen’t work anymore for solving needed generalized eigenvalue problems. So I have to find a solution for this now haha.

Really? There were issues like that that I thought were fixed recently by switching to 32bits MKL. Also worth noting that there are several good pure-julia alternatives to Arpack, eg IterativeSolvers.jl and KrylovKit.jl

> Well, it depends what you mean by destructive. Those annotations mean that your code will fail for real-valued matrices, or if they are Diagonal, Hermitian, SharedArrrays, StaticArrays, etc etc. Every array that is not strictly of type Array, and with precisely 64 bit floating point complex elements, which is _very_ restrictive.

To expand on that, the standard way of doing things in julia is 1) do not type input arguments unless you use it for dispatch or it helps understand what the function does 2) if you have to type, type it in the least possible restrictive way. This makes code less verbose, and allows things like GPU, autodiff, multiple floating point types, etc. to work without modifications to the code.

---

<div class="post-metadata">

**Author:** ![Jonas](https://avatars.discourse-cdn.com/v4/letter/j/aeb1de/32.png) [@Jonas](https://discourse.julialang.org/u/Jonas)\
**Post date:** [November 15, 2020, 7:15pm UTC](https://discourse.julialang.org/t/acceleration-of-a-function/50181/11 "2020-11-15T19:15:37Z")

</div>

> [@antoine-levitt](#):
>
> Really? There were issues like that that I thought were fixed recently by switching to 32bits MKL. Also worth noting that there are several good pure-julia alternatives to Arpack, eg IterativeSolvers.jl and KrylovKit.jl

Yes, I get this error:  
ERROR: LoadError: ARPACKException: unspecified ARPACK error: -10  
Stacktrace:  
[1] aupd\_wrapper(::Type{T} where T, ::Arpack.var"#matvecA!#24"{SparseMatrixCSC{Float64,Int64}}, ::Arpack.var"#21#28"{SparseMatrixCSC{Float64,Int64}}, ::Arpack.var"#23#30", ::Int64, ::Bool, ::Bool, ::String,  
::Int64, ::Int32, ::String, ::Float64, ::Int64, ::Int64, ::Array{Float64,1}) at path.julia\packages\Arpack\o35I5\src\libarpack.jl:76  
[2] \_eigs(::SparseMatrixCSC{Float64,Int64}, ::SparseMatrixCSC{Float64,Int64}; nev::Int64, ncv::Int64, which::Symbol, tol::Float64, maxiter::Int64, sigma::Nothing, v0::Array{Float64,1}, ritzvec::Bool) at path.julia\packages\Arpack\o35I5\src\Arpack.jl:181  
[3] #eigs#11 at path.julia\packages\Arpack\o35I5\src\Arpack.jl:48 [inlined]  
[4] top-level scope at …path…\Filter Tool\myfile.jl:19  
[5] include\_string(::Function, ::Module, ::String, ::String) at .\loading.jl:1088  
in expression starting at …path…\myfile.jl:19

Before switching to MKL.jl, using the “eigs” function worked like a charm. I tried KrylovKit.jl and it works, even with Sparse Arrays. I am not sure if it is as fast as eigs() before though. Thank you for the suggestion.

---

<div class="post-metadata">

**Author:** ![antoine-levitt](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/antoine-levitt/32/4008_2.png) [@antoine-levitt](https://discourse.julialang.org/u/antoine-levitt)\
**Post date:** [November 15, 2020, 7:28pm UTC](https://discourse.julialang.org/t/acceleration-of-a-function/50181/12 "2020-11-15T19:28:29Z")

</div>

Possibly you need to rebuild ARPACK? You can try opening an issue at ARPACK.jl.

> I am not sure if it is as fast as eigs() before though

They both implement variants of the same broad class of algorithms, so I would imagine their performance is similar. Comparisons would be interesting though. @juthohaegeman may have some thoughts on this.

---

<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:** [November 15, 2020, 8:53pm UTC](https://discourse.julialang.org/t/acceleration-of-a-function/50181/13 "2020-11-15T20:53:19Z")

</div>

I often find that you should supply the keyword argument `sigma` for `eigs` to work well. Of course, it depends on the matrix but…

---

<div class="post-metadata">

**Author:** ![Jonas](https://avatars.discourse-cdn.com/v4/letter/j/aeb1de/32.png) [@Jonas](https://discourse.julialang.org/u/Jonas)\
**Post date:** [November 16, 2020, 10:01am UTC](https://discourse.julialang.org/t/acceleration-of-a-function/50181/14 "2020-11-16T10:01:14Z")

</div>

Ok thank you for your reply, I will try it!  
best regards.

---

<div class="post-metadata">

**Author:** ![Jonas](https://avatars.discourse-cdn.com/v4/letter/j/aeb1de/32.png) [@Jonas](https://discourse.julialang.org/u/Jonas)\
**Post date:** [November 16, 2020, 10:20am UTC](https://discourse.julialang.org/t/acceleration-of-a-function/50181/15 "2020-11-16T10:20:26Z")

</div>

I rebuilt Arpack with no success. But thank you very much for your help. I learned a lot 😃

---

<div class="post-metadata">

**Author:** ![colintbowers](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/colintbowers/32/8033_2.png) [@colintbowers](https://discourse.julialang.org/u/colintbowers)\
**Post date:** [November 16, 2020, 10:28am UTC](https://discourse.julialang.org/t/acceleration-of-a-function/50181/16 "2020-11-16T10:28:12Z")

</div>

You do you! Personally I like to leave them in too - when I’m looking back at the code in 3 months time the type annotations often get me back in the right headspace quickly.

---

<div class="post-metadata">

**Author:** ![Jonas](https://avatars.discourse-cdn.com/v4/letter/j/aeb1de/32.png) [@Jonas](https://discourse.julialang.org/u/Jonas)\
**Post date:** [November 16, 2020, 2:09pm UTC](https://discourse.julialang.org/t/acceleration-of-a-function/50181/17 "2020-11-16T14:09:49Z")

</div>

For anyone interested, I got it to work now. Instead of Arpack I used the IterativeSolvers package. And everything works with MKL now. I got a **massiv** boost in execution time. My function before took about 800ms and now with MKS it only takes about 2.5ms. Back to work, thank you to everyone for helping me!

Best regards,  
Jonas