# slow performance of eig
**URL:** https://discourse.julialang.org/t/slow-performance-of-eig/13074
**Category:** Performance
**Tags:** question
**Created:** [August 9, 2018, 9:48am UTC](https://discourse.julialang.org/t/slow-performance-of-eig/13074 "2018-08-09T09:48:32Z")
**Posts on this page:** 15
**Page:** 1
### Author:  [@itzui](https://discourse.julialang.org/u/itzui)
#### Post date: [August 9, 2018, 9:48am UTC](https://discourse.julialang.org/t/slow-performance-of-eig/13074/1 "2018-08-09T09:48:32Z")
Hello,
the performance of the determination of the eigenvalues and the eigenvector for small complex matrices in Julia seems for me quite slow. If I use the following code
> N = 100000
> dim = 3
> function eigen()
> for i = 1:N
> A = rand(dim, dim) + im \* rand(dim, dim)
> eig(A)
> end
> end
> @time eigen()
it takes about 3sec.
> 3.226431 seconds (5.76 M allocations: 436.489 MiB, 2.34% gc time)
If I use the same code in Python (with numpy) the code is slighty slower (4.4s), however by using numba I can speed up the loop and it takes then only 0.7s which is much faster than Julia.
Is there a workaround or package to enhance the performance for these matrices?
PS: versioninfo in Julia gives the following
Julia Version 0.6.2
Commit d386e40c17\* (2017-12-13 18:08 UTC)
Platform Info:
OS: Windows (x86\_64-w64-mingw32)
CPU: Intel(R) Core™ i5-6500 CPU @ 3.20GHz
WORD\_SIZE: 64
BLAS: mkl\_rt
LAPACK: mkl\_rt
LIBM: libopenlibm
LLVM: libLLVM-3.9.1 (ORCJIT, skylake)
---
### Author:  [@ufechner7](https://discourse.julialang.org/u/ufechner7)
#### Post date: [August 9, 2018, 3:06pm UTC](https://discourse.julialang.org/t/slow-performance-of-eig/13074/2 "2018-08-09T15:06:10Z")
Most likely you are benchmarking the performance of rand() and not of eig(). Did you try to pre-calculate the random arrays and @time only the eig() calculations?
---
### Author:  [@Tero\_Frondelius](https://discourse.julialang.org/u/Tero_Frondelius)
#### Post date: [August 9, 2018, 5:04pm UTC](https://discourse.julialang.org/t/slow-performance-of-eig/13074/4 "2018-08-09T17:04:37Z")
Can you try BenchmarkTools? For me it looks like you are measuring compile time.
---
### Author:  [@DNF](https://discourse.julialang.org/u/DNF)
#### Post date: [August 9, 2018, 5:19pm UTC](https://discourse.julialang.org/t/slow-performance-of-eig/13074/5 "2018-08-09T17:19:00Z")
If you want to benchmark `eig`, you should do that without mixing in the performance of `rand`, and also avoid using global variables (`N` and `dim` are globals), which violates the “first guideline of Julia benchmarking”.
Here’s the “right” way:
```julia
using BenchmarkTools
A = rand(Complex{Float64}, dim, dim)
@btime eig($A)
```
---
### Author:  [@ufechner7](https://discourse.julialang.org/u/ufechner7)
#### Post date: [August 9, 2018, 5:21pm UTC](https://discourse.julialang.org/t/slow-performance-of-eig/13074/6 "2018-08-09T17:21:31Z")
> [@DNF](#):
>
> and also avoid using global variables ( `N` and `dim` are globals)
Well, if you define them as const this should also be fine:
const N=100000
const dim=3
---
### Author:  [@itzui](https://discourse.julialang.org/u/itzui)
#### Post date: [August 10, 2018, 6:56am UTC](https://discourse.julialang.org/t/slow-performance-of-eig/13074/7 "2018-08-10T06:56:33Z")
Dear All,
thanks a lot for the comments and help. If I use BenchmarkTools and define the matrix before the loop I can reduce the computation time by 1s. But its still about 3 times slower than using numpy + numba.
The code I used now is the following:
> using BenchmarkTools
> const mat = rand(3, 3) + im \* rand(3, 3)
> @btime for i = 1:100000 eig(mat) end;
The output is then:
2.132 s (3800000 allocations: 289.92 MiB)
---
### Author:  [@mohamed82008](https://discourse.julialang.org/u/mohamed82008)
#### Post date: [August 10, 2018, 7:21am UTC](https://discourse.julialang.org/t/slow-performance-of-eig/13074/8 "2018-08-10T07:21:00Z")
You don’t need the loop since `btime` will take care of this for you and you should interpolate the variables in a `btime` expression as such `@btime eig($mat)`.
---
### Author:  [@itzui](https://discourse.julialang.org/u/itzui)
#### Post date: [August 10, 2018, 7:28am UTC](https://discourse.julialang.org/t/slow-performance-of-eig/13074/9 "2018-08-10T07:28:10Z")
> [@mohamed82008](#):
>
> You don’t need the loop since `btime` will take care of this for you and you should interpolate the variables in a `btime` expression as such `@btime eig($mat)`
Thanks. Sorry, I’m still a beginner in Julia. If I remove the loop I’ll end up with 20µs for one call. But still, for 10^5 calls I need 2s. So there is only a small enhancement of the time. The main result, that its much slower than numpy+numba, is not changing.
---
### Author:  [@DNF](https://discourse.julialang.org/u/DNF)
#### Post date: [August 10, 2018, 8:48am UTC](https://discourse.julialang.org/t/slow-performance-of-eig/13074/10 "2018-08-10T08:48:02Z")
> [@itzui](#):
>
> `mat = rand(3, 3) + im * rand(3, 3)`
Note that this creates several unnecessary temporary arrays. You should do this instead:
```julia
rand(Complex{Float64}, 3, 3)
```
which creates just one array.
> [@itzui](#):
>
> The main result, that its much slower than numpy+numba, is not changing.
Yes, the speedups from better benchmarking was not as dramatic as expected. Can you provide the python code you use for benchmarking? I am seeing a factor of 2x difference. Perhaps an issue should be opened.
Python:
```python
In [12]: A = np.random.rand(3,3) + 1j*np.random.rand(3,3) # I don't know how to do this correctly in Python! :/
In [13]: %timeit np.linalg.eig(A)
30.4 µs ± 539 ns per loop (mean ± std. dev. of 7 runs, 10000 loops each)
@jit
...: def myeig(A):
...: return np.linalg.eig(A)
In [27]: %timeit myeig(A)
10.7 µs ± 25.2 ns per loop (mean ± std. dev. of 7 runs, 100000 loops each)
```
Julia
```julia
julia> A = rand(Complex{Float64}, 3, 3);
julia> @btime eigen($A);
20.381 μs (35 allocations: 8.80 KiB)
```
---
### Author:  [@ufechner7](https://discourse.julialang.org/u/ufechner7)
#### Post date: [August 10, 2018, 8:51am UTC](https://discourse.julialang.org/t/slow-performance-of-eig/13074/11 "2018-08-10T08:51:11Z")
Are you using Python and Numba with MKL or with OpenBLAS?
---
### Author:  [@DNF](https://discourse.julialang.org/u/DNF)
#### Post date: [August 10, 2018, 9:10am UTC](https://discourse.julialang.org/t/slow-performance-of-eig/13074/12 "2018-08-10T09:10:05Z")
Are you asking me or @itzui? I’m not sure, my `np. __config__ ` seems to indicate MKL. It’s whatever ships with Anaconda/Python 3.6.
I guess this is all up to the BLAS library then. But I am surprised how slow this all is. Even for 2x2 matrices it runs in 16us, while I would have expected 16 _nano_ seconds, or something. There must be a tremendous overhead due to BLAS.
So I tried StaticArrays, and to my surprise it was no faster. (In fact, StaticArrays gives an error for non-hermitian matrices! But even for hermitian ones it was very slow.)
---
### Author:  [@itzui](https://discourse.julialang.org/u/itzui)
#### Post date: [August 10, 2018, 10:11am UTC](https://discourse.julialang.org/t/slow-performance-of-eig/13074/13 "2018-08-10T10:11:47Z")
My code for testing in python was the following (i know that it is optimized):
> from numba import njit, prange
> import numpy as np
> import time
> @njit
> def test():
> for i in prange(0,100000):
> A = np.random.rand(3, 3) + 1j \* np.random.rand(3,3)
> np.linalg.eig(A)
> return
As result I got
773 ms ± 8.37 ms per loop (mean ± std. dev. of 7 runs, 1 loop each)
`np. __config__ ` tells me, that MKL is used. But this should also be true for Julia if I interpreted the versioninfo() (BLAS: mkl\_rt) correctly
---
### Author:  [@antoine-levitt](https://discourse.julialang.org/u/antoine-levitt)
#### Post date: [August 10, 2018, 10:29am UTC](https://discourse.julialang.org/t/slow-performance-of-eig/13074/14 "2018-08-10T10:29:47Z")
Static arrays don’t have optimized implementations of complex hermitian eigenproblems, and fall back to regular eig : [https://github.com/JuliaArrays/StaticArrays.jl/blob/master/src/eigen.jl](https://github.com/JuliaArrays/StaticArrays.jl/blob/master/src/eigen.jl), although it possibly shouldn’t be too hard to adapt the existing code? Lapack and Blas libraries are pretty slow for small matrices (newer mkl have a special “fast” mode for such operations, but I’m not sure how fast it is)
---
### Author:  [@stabbles](https://discourse.julialang.org/u/stabbles)
#### Post date: [August 10, 2018, 1:16pm UTC](https://discourse.julialang.org/t/slow-performance-of-eig/13074/15 "2018-08-10T13:16:13Z")
One of the issues with `eigen` is that it is type unstable, which is causing overhead for tiny matrices.
```julia
julia> @benchmark LinearAlgebra.LAPACK.geevx!('B', 'N', 'V', 'N', copy($A))
BenchmarkTools.Trial:
memory estimate: 7.59 KiB
allocs estimate: 13
--------------
minimum time: 6.627 μs (0.00% GC)
median time: 6.833 μs (0.00% GC)
julia> @benchmark eigen($A)
BenchmarkTools.Trial:
memory estimate: 8.98 KiB
allocs estimate: 38
--------------
minimum time: 16.252 μs (0.00% GC)
median time: 17.885 μs (0.00% GC)
```
Type instability stems from the fact that it sometimes uses real arithmetic when the eigenvalues and eigenvectors are real, and complex arithmetic otherwise.
Edit: I just realized that the above might look unintuitive when dealing with complex-valued matrices, so here’s an example. When the matrix is Hermitian (`A' == A`) the eigenvalues are real, and so is the vector that Julia returns. For non-Hermitian matrices vice-versa:
```julia
julia> eigen(begin
A = rand(ComplexF64, 3, 3)
A + A'
end).values |> typeof
Array{Float64,1}
julia> eigen(rand(ComplexF64, 3, 3)).values |> typeof
Array{Complex{Float64},1}
```
---
### Author:  [@itzui](https://discourse.julialang.org/u/itzui)
#### Post date: [August 14, 2018, 5:39am UTC](https://discourse.julialang.org/t/slow-performance-of-eig/13074/16 "2018-08-14T05:39:40Z")
Thanks for the comment. Does it mean, that for small matrices it is better to test the type of the matrix and call the LAPACK routine directly?