# Why Julia and linear algebra module are incredibly slow, compare to C++ and Eigen

**URL:** <https://discourse.julialang.org/t/why-julia-and-linear-algebra-module-are-incredibly-slow-compare-to-c-and-eigen/21304>\
**Category:** Performance\
**Tags:** question, package, performance, linearalgebra\
**Created:** [February 28, 2019, 2:32pm UTC](https://discourse.julialang.org/t/why-julia-and-linear-algebra-module-are-incredibly-slow-compare-to-c-and-eigen/21304 "2019-02-28T14:32:20Z")\
**Posts on this page:** 15\
**Page:** 1

<div class="post-metadata">

**Author:** ![Amin\_Abouee](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/amin_abouee/32/7275_2.png) [@Amin\_Abouee](https://discourse.julialang.org/u/Amin_Abouee)\
**Post date:** [February 28, 2019, 2:32pm UTC](https://discourse.julialang.org/t/why-julia-and-linear-algebra-module-are-incredibly-slow-compare-to-c-and-eigen/21304/1 "2019-02-28T14:32:20Z")

</div>

Hi,  
I developed a julia file to compare its performance with C++. I implemented the [homography estimation based on DLT algortihm](https://en.wikipedia.org/wiki/Direct_linear_transformation). In brief, I created a 8\*9 matrix A from two bundle of four 3D points and solve Ah=0. So we just need to compute nullspace vectors of matrix A and finish.  
For this one, I implemented the following julia code and one equivalence C++ code that you can find [here](https://github.com/amin-abouee/test-cplusplus-and-julia/blob/master/main.cpp). I tried to have same number of variables and function’s calls in both languages. the linear algebra part of C++ was implemented by eigen library.

```julia
using LinearAlgebra

function position_in_world(x::Float64, y::Float64)
      target_width_px = 1123
      target_height_px = 791
      target_width_mm = 297
      target_height_mm = 210.025
      X = ((x + 0.5) / target_width_px) * target_width_mm - 0.5 * target_width_mm
      Y = ((y + 0.5) / target_height_px) * target_height_mm - 0.5 * target_height_mm
      Z = 0.0
      return [X, Y, Z]
end

# function position_in_world2!(x::Float64, y::Float64, vec::Array{Float64,1})
# target_width_px = 1123
# target_height_px = 791
# target_width_mm = 297
# target_height_mm = 210.025
# vec[1] = ((x + 0.5) / target_width_px) * target_width_mm - 0.5 * target_width_mm
# vec[2] = ((y + 0.5) / target_height_px) * target_height_mm - 0.5 * target_height_mm
# vec[3] = 0.0
# end

function π_projection_function(K::Array{Float64,2}, R::Array{Float64,2}, t::Array{Float64,1}, point::Array{Float64,1})
      point_in_camera = R * point + t
      point_in_image = K * point_in_camera
      return point_in_image / point_in_camera[end]
end

function compute_homograpy_DLT(M, P)
      A = zeros(Float64, 8, 9)
      for i in 1:size(M)[1]
            A[2 * i - 1, 1:3] = M[i,:]
            A[2 * i - 1, 7:9] = - P[i,1] * M[i,:]
            A[2 * i, 4:6] = M[i,:]
            A[2 * i, 7:9] = - P[i,2] * M[i,:]
      end
      return nullspace(A)
end

@timev begin
K = [ 1169.19630 0.0 652.98743;
      0.0 1169.61014 528.83429;
      0.0 0.0 1.0]

T = [0.961255 -0.275448 0.0108487 112.79;
 0.171961 0.629936 0.75737 -217.627;
 -0.21545 -0.72616 0.652895 1385.13]

R = T[1:3, 1:3]
t = T[:,end]

tl = position_in_world(0.0, 0.0)
tr = position_in_world(1123.0, 0.0)
br = position_in_world(1123.0, 791.0)
bl = position_in_world(0.0, 791.0)

# println("pos in world 2")
# @timev begin
# tl = Float64[0.0, 0.0, 0.0]
# tr = Float64[0.0, 0.0, 0.0]
# br = Float64[0.0, 0.0, 0.0]
# bl = Float64[0.0, 0.0, 0.0]
# position_in_world2!(0.0, 0.0,tl)
# position_in_world2!(1123.0, 0.0,tr)
# position_in_world2!(1123.0, 791.0,br)
# position_in_world2!(0.0, 791.0,bl)
# end

p1 = π_projection_function(K, R, t, tl)
p2 = π_projection_function(K, R, t, tr)
p3 = π_projection_function(K, R, t, br)
p4 = π_projection_function(K, R, t, bl)
P = [p1 p2 p3 p4]'

m1 = [0.0, 0.0, 1.0]
m2 = [1123.0, 0.0, 1.0]
m3 = [1123.0, 791.0, 1.0]
m4 = [0.0, 791.0, 1.0]
M = [m1 m2 m3 m4]'

res = compute_homograpy_DLT(M , P)
homography = reshape(res, (3,3))'
println("Homography: ", homography)
end

```

And the result:

```julia
Julia (version 1.1.0)
Homography: [0.000244416 -0.000198718 0.915494; 2.16745e-5 8.80403e-5 0.40233; -5.35586e-8 -1.81231e-7 0.00140359]
  0.073007 seconds (218.81 k allocations: 11.512 MiB)
elapsed time (ns): 73007284
bytes allocated: 12071266
pool allocs: 218781
non-pool GC allocs:33

C++ ( clang++ version: 7.0.1, release-flags: "-march=native -O3 -DNDEBUG", Eigen version: 3.3.4) 
Homography: 
 0.000244416 -0.000198718 0.915494
 2.16745e-05 8.80403e-05 0.40233
-5.35586e-08 -1.81231e-07 0.00140359
elapsed time (ns): 275887

```

As can be seen, the homography results are same but c++ and eigen are approximately 250! times faster.  
It is noticeable that I am beginner in julia and there is a possibility that I did not implement my code in efficient way. So the main question is, how can I improve the speed of my julia code?

---

<div class="post-metadata">

**Author:** ![carstenbauer](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/carstenbauer/32/4981_2.png) [@carstenbauer](https://discourse.julialang.org/u/carstenbauer)\
**Post date:** [February 28, 2019, 2:47pm UTC](https://discourse.julialang.org/t/why-julia-and-linear-algebra-module-are-incredibly-slow-compare-to-c-and-eigen/21304/2 "2019-02-28T14:47:38Z")

</div>

Before anything you should take a look at the [Performance tips](https://docs.julialang.org/en/latest/manual/performance-tips/).

To get you started, you should really try to avoid all those temporary allocations. For example, a slice, like `M[i,:]`, on the rhs of an assignment allocates a whole new array. Here you might want to consider `@views`.

---

<div class="post-metadata">

**Author:** ![Yifan\_Liu](https://avatars.discourse-cdn.com/v4/letter/y/4da419/32.png) [@Yifan\_Liu](https://discourse.julialang.org/u/Yifan_Liu)\
**Post date:** [February 28, 2019, 3:27pm UTC](https://discourse.julialang.org/t/why-julia-and-linear-algebra-module-are-incredibly-slow-compare-to-c-and-eigen/21304/3 "2019-02-28T15:27:12Z")

</div>

I remember reading somewhere that Eigen is even faster than MKL in small matrix computations, however such advantage will disappear when the dimension of matrix increases. Could you try a larger matrix?

---

<div class="post-metadata">

**Author:** ![Egwene\_al\_Vere](https://avatars.discourse-cdn.com/v4/letter/e/df788c/32.png) [@Egwene\_al\_Vere](https://discourse.julialang.org/u/Egwene_al_Vere)\
**Post date:** [February 28, 2019, 4:02pm UTC](https://discourse.julialang.org/t/why-julia-and-linear-algebra-module-are-incredibly-slow-compare-to-c-and-eigen/21304/4 "2019-02-28T16:02:40Z")

</div>

Apart from pre-allocating, consider using `@btime` from `BenchmarkTools`, also looks like your test may have included in the compile time of the jit. You may try running `@btime` or just run your @timev block a second time to exclude the compile time. Refer to e.g [here](https://github.com/JuliaCI/BenchmarkTools.jl) for tips on how to use `@btime`. Additionally, looks like your matrix size is small, you may try [https://github.com/JuliaArrays/StaticArrays.jl](https://github.com/JuliaArrays/StaticArrays.jl)

---

<div class="post-metadata">

**Author:** ![tkoolen](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/tkoolen/32/1603_2.png) [@tkoolen](https://discourse.julialang.org/u/tkoolen)\
**Post date:** [February 28, 2019, 4:27pm UTC](https://discourse.julialang.org/t/why-julia-and-linear-algebra-module-are-incredibly-slow-compare-to-c-and-eigen/21304/5 "2019-02-28T16:27:49Z")

</div>

You’re using `Eigen::Vector3d` for points in C++, yet you’re using the equivalent of `Eigen::VectorXd` in Julia (i.e., a `Vector{Float64}`). You’ll definitely want to use StaticArrays.

---

<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:** [February 28, 2019, 7:59pm UTC](https://discourse.julialang.org/t/why-julia-and-linear-algebra-module-are-incredibly-slow-compare-to-c-and-eigen/21304/6 "2019-02-28T19:59:07Z")

</div>

Just wrapping your main code in a function and timing it with `@btime` from BenchmarkTools I got a runtime of 21us, which is several thousand times faster than what you’re seeing.

I also tried using StaticArrays, in which case runtime is entirely dominated by by the call to `nullspace` (btw. `nullspace` doesn’t work for static arrays, unfortunately): all the code, except `nullspace`, took 150ns to run, but then `nullspace(A)` takes 18us.

Here’s a simple benchmark without static arrays, but using BenchmarkTools:

```julia
function foo()
    K = [ 1169.19630 0.0 652.98743;
             0.0 1169.61014 528.83429;
             0.0 0.0 1.0]

    T = [0.961255 -0.275448 0.0108487 112.79;
         0.171961 0.629936 0.75737 -217.627;
        -0.21545 -0.72616 0.652895 1385.13]

    R = T[:, SVector(1,2,3)]
    t = T[:,end]

    tl = position_in_world(0.0, 0.0)
    tr = position_in_world(1123.0, 0.0)
    br = position_in_world(1123.0, 791.0)
    bl = position_in_world(0.0, 791.0)

    p1 = π_projection_function(K, R, t, tl)
    p2 = π_projection_function(K, R, t, tr)
    p3 = π_projection_function(K, R, t, br)
    p4 = π_projection_function(K, R, t, bl)
    P = [p1 p2 p3 p4]'

    M = [0.0 0.0 1.0;
       1123.0 0.0 1.0;
       1123.0 791.0 1.0;
             0.0 791.0 1.0]

    res = compute_homograpy_DLT(M, P)
    homography = reshape(res, (3,3))'
    return homography
end

```

Then benchmark it:

```julia
julia> @btime bar()
  20.643 μs (65 allocations: 14.48 KiB)
3×3 Adjoint{Float64,Array{Float64,2}}:
  0.000244416 -0.000198718 0.915494  
  2.16745e-5 8.80403e-5 0.40233   
 -5.35586e-8 -1.81231e-7 0.00140359

```

**Edit:** To clarify, the problem isn’t with your code, but with your benchmarking.

---

<div class="post-metadata">

**Author:** ![Zach\_Christensen](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/zach_christensen/32/7220_2.png) [@Zach\_Christensen](https://discourse.julialang.org/u/Zach_Christensen)\
**Post date:** [March 1, 2019, 1:16am UTC](https://discourse.julialang.org/t/why-julia-and-linear-algebra-module-are-incredibly-slow-compare-to-c-and-eigen/21304/7 "2019-03-01T01:16:33Z")

</div>

I think the comment about StaticArrays may still stand for a thorough comparison. Eigen is pretty explicit about this sort of thing and I believe it is where libraries like tensorflow get a lot of their speed because I read earlier today about someone having trouble getting past the 32 dimension limit that’s hard coded into tensorflow.

---

<div class="post-metadata">

**Author:** ![c42f](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/c42f/32/52842_2.png) [@c42f](https://discourse.julialang.org/u/c42f)\
**Post date:** [March 1, 2019, 2:56am UTC](https://discourse.julialang.org/t/why-julia-and-linear-algebra-module-are-incredibly-slow-compare-to-c-and-eigen/21304/8 "2019-03-01T02:56:43Z")

</div>

> [@DNF](#):
>
> I also tried using StaticArrays, in which case runtime is entirely dominated by by the call to `nullspace` (btw. `nullspace` doesn’t work for static arrays, unfortunately): all the code, except `nullspace` , took 150ns to run, but then `nullspace(A)` takes 18us.

Yes it it’s not possible to make a super fast type stable `nullspace` because the output size depends on the elements rather than the dimension of the input. Eigen would have the same problem of course. Even so, it’s a pity that the `Base` version of `nullspace` doesn’t work out of the box with `SMatrix`. Here, I fixed this:

> <https://github.com/JuliaArrays/StaticArrays.jl/pull/596>
>
> A small update to the \`svd\` wrapper so that \`full=true\` is supported. Code occas…ionally uses this, for example \`nullspace\`... in fact this is the only place \`full=true\` occurs in \`LinearAlgebra\`.
> 
> XRef: https://discourse.julialang.org/t/why-julia-and-linear-algebra-module-are-incredibly-slow-compare-to-c-and-eigen/21304/6

(Hmm. It would be possible to return an SMatrix wrapped in a view from `nullspace` if somebody cared enough to try it out. I’m not sure that would be useful though.)

---

<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:** [March 1, 2019, 8:30am UTC](https://discourse.julialang.org/t/why-julia-and-linear-algebra-module-are-incredibly-slow-compare-to-c-and-eigen/21304/9 "2019-03-01T08:30:11Z")

</div>

> [@c42f](#):
>
> Yes it it’s not possible to make a super fast type stable `nullspace` because the output size depends on the elements rather than the dimension of the input.

Would it help if the output were a `Vector` of `SVector{N}`s, where the length `N` is known?

---

<div class="post-metadata">

**Author:** ![c42f](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/c42f/32/52842_2.png) [@c42f](https://discourse.julialang.org/u/c42f)\
**Post date:** [March 1, 2019, 8:52am UTC](https://discourse.julialang.org/t/why-julia-and-linear-algebra-module-are-incredibly-slow-compare-to-c-and-eigen/21304/10 "2019-03-01T08:52:45Z")

</div>

> [@DNF](#):
>
> Would it help if the output were a `Vector` of `SVector{N}` s, where the length `N` is known?

I don’t think so. But regardless of that, we don’t have a statically sized svd implementation anyway so there’s much more room for improvement there. Currently we fall back to `svd(::AbstractMatrix)` and simply wrap the output to make sure we preserve the size information. We could probably do much better by using `MMatrix` and passing that into `LinearAlgebra.LAPACK.gesdd!` directly for `BlasFloat` types, but nobody has bothered to dig in that deep yet.

---

<div class="post-metadata">

**Author:** ![Tamas\_Papp](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/tamas_papp/32/25949_2.png) [@Tamas\_Papp](https://discourse.julialang.org/u/Tamas_Papp)\
**Post date:** [March 1, 2019, 9:09am UTC](https://discourse.julialang.org/t/why-julia-and-linear-algebra-module-are-incredibly-slow-compare-to-c-and-eigen/21304/11 "2019-03-01T09:09:17Z")

</div>

> [@c42f](#):
>
> a statically sized svd implementation

I think that 3x3 has a more or less manageable closed form, at least I have seen it hardcoded.

---

<div class="post-metadata">

**Author:** ![Amin\_Abouee](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/amin_abouee/32/7275_2.png) [@Amin\_Abouee](https://discourse.julialang.org/u/Amin_Abouee)\
**Post date:** [March 1, 2019, 1:03pm UTC](https://discourse.julialang.org/t/why-julia-and-linear-algebra-module-are-incredibly-slow-compare-to-c-and-eigen/21304/12 "2019-03-01T13:03:24Z")

</div>

To be honest, before this post, I didn’t know about the overhead of JIT in julia. I read all of your comments and fixed my code. First of all, the BenchmarkTools is applied to evaluate the execution speed and I had significant improvement after JIT warm-up. Secondly, all possible variables of type `AbstractArray` (dynamic memory allocation) are replaced with `SVector` and `SMatrix`. Finnaly, I used `@view` macro for `SubArray` slicing.

```julia
function compute_homograpy_DLT(M, P)
      A = zeros(8,9)
      for i in 1:size(M)[1]
            A[2 * i - 1, 1:3] = @view M[i,:]
            A[2 * i - 1, 7:9] = - P[i,1] * @view M[i,:]
            A[2 * i, 4:6] = @view M[i,:]
            A[2 * i, 7:9] = - P[i,2] * @view M[i,:]
      end
      return nullspace(A)
end

```

and

```julia
K = @SMatrix [ 1169.19630 0.0 652.98743;
            0.0 1169.61014 528.83429;
            0.0 0.0 1.0]
T = @SMatrix [0.961255 -0.275448 0.0108487 112.79;
       0.171961 0.629936 0.75737 -217.627;
       -0.21545 -0.72616 0.652895 1385.13]

R = @view T[1:3, 1:3]
t = @view T[:,end]

tl = position_in_world(0.0, 0.0)
tr = position_in_world(1123.0, 0.0)
br = position_in_world(1123.0, 791.0)
bl = position_in_world(0.0, 791.0)

p1 = π_projection_function(K, R, t, tl)
p2 = π_projection_function(K, R, t, tr)
p3 = π_projection_function(K, R, t, br)
p4 = π_projection_function(K, R, t, bl)
P = [p1 p2 p3 p4]'

M = @SMatrix [0.0 0.0 1.0;
                    1123.0 0.0 1.0;
                    1123.0 791.0 1.0;
                    0.0 791.0 1.0]

```

and re-execute the code:

```julia
151.783 μs (133 allocations: 13.27 KiB)

```

BUT, for a fair comparison between julia and eigen, I also modified my c++ code. For this one, I ran the code for 100,000 and the output average time is:

```julia
elapsed time (microseconds): 34.2004

```

So the C++ code is still 5 times faster and I think the reason is coming from null space computation.

---

<div class="post-metadata">

**Author:** ![Amin\_Abouee](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/amin_abouee/32/7275_2.png) [@Amin\_Abouee](https://discourse.julialang.org/u/Amin_Abouee)\
**Post date:** [March 1, 2019, 2:09pm UTC](https://discourse.julialang.org/t/why-julia-and-linear-algebra-module-are-incredibly-slow-compare-to-c-and-eigen/21304/13 "2019-03-01T14:09:06Z")

</div>

I have never installed MKL and consequently, I don’t have any estimation about the performance of it. But for Eigen, there is one [benchmark](https://eigen.tuxfamily.org/dox/group__DenseDecompositionBenchmark.html), which is done by eigen team. They evaluates the performance of Eigen for different dense decomposition methods, based on the size of input matrix.

---

<div class="post-metadata">

**Author:** ![tkoolen](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/tkoolen/32/1603_2.png) [@tkoolen](https://discourse.julialang.org/u/tkoolen)\
**Post date:** [March 1, 2019, 2:51pm UTC](https://discourse.julialang.org/t/why-julia-and-linear-algebra-module-are-incredibly-slow-compare-to-c-and-eigen/21304/14 "2019-03-01T14:51:25Z")

</div>

The following runs in ~11 microseconds on my machine:

```julia
using LinearAlgebra
using StaticArrays

function position_in_world(x::Float64, y::Float64)
    target_width_px = 1123
    target_height_px = 791
    target_width_mm = 297
    target_height_mm = 210.025
    X = ((x + 0.5) / target_width_px) * target_width_mm - 0.5 * target_width_mm
    Y = ((y + 0.5) / target_height_px) * target_height_mm - 0.5 * target_height_mm
    Z = 0.0
    return SVector(X, Y, Z)
end

function compute_homograpy_DLT(M, P)
    T = promote_type(eltype(M), eltype(P))
    A = zeros(T, 8, 9)
    for i in 1:size(M, 1)
        Mi = M[i,:]
        A[2 * i - 1, 1:3] = Mi
        A[2 * i - 1, 7:9] = - P[i,1] * Mi
        A[2 * i, 4:6] = Mi
        A[2 * i, 7:9] = - P[i,2] * Mi
    end
    return nullspace(A)
end

function π_projection_function(K::AbstractMatrix, R::AbstractMatrix, t::AbstractVector, point::AbstractVector)
    point_in_camera = R * point + t
    point_in_image = K * point_in_camera
    return point_in_image / point_in_camera[end]
end
                  
function foo()
    K = @SMatrix [ 1169.19630 0.0 652.98743;
                0.0 1169.61014 528.83429;
                0.0 0.0 1.0]

    T = @SMatrix [0.961255 -0.275448 0.0108487 112.79;
            0.171961 0.629936 0.75737 -217.627;
        -0.21545 -0.72616 0.652895 1385.13]

    R = T[:, SOneTo(3)]
    t = T[:,end]

    tl = position_in_world(0.0, 0.0)
    tr = position_in_world(1123.0, 0.0)
    br = position_in_world(1123.0, 791.0)
    bl = position_in_world(0.0, 791.0)

    p1 = π_projection_function(K, R, t, tl)
    p2 = π_projection_function(K, R, t, tr)
    p3 = π_projection_function(K, R, t, br)
    p4 = π_projection_function(K, R, t, bl)
    P = [p1 p2 p3 p4]'

    M = @SMatrix [0.0 0.0 1.0;
        1123.0 0.0 1.0;
        1123.0 791.0 1.0;
                0.0 791.0 1.0]

    res = compute_homograpy_DLT(M, P)
    homography = reshape(res, (3,3))'
    return homography
end

using BenchmarkTools
@btime foo()

```

Note that everything up to `compute_homography_DLT` is completely constant-folded (computed once, at compile time), which is fair game for the compiler since `foo` takes no arguments and all of the functions called on the constants defined in `foo` are pure. In fact, ideally the whole thing would get constant-folded, but the creation of dynamically-sized `Array`s in `compute_homograpy_DLT` is preventing this from happening. It depends on your actual use case whether this is a good benchmark or not.

---

<div class="post-metadata">

**Author:** ![Amin\_Abouee](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/amin_abouee/32/7275_2.png) [@Amin\_Abouee](https://discourse.julialang.org/u/Amin_Abouee)\
**Post date:** [March 1, 2019, 3:37pm UTC](https://discourse.julialang.org/t/why-julia-and-linear-algebra-module-are-incredibly-slow-compare-to-c-and-eigen/21304/15 "2019-03-01T15:37:29Z")

</div>

In the previous evaluation, I did not remove the print line which you did. By the way, for a practical and more fair test, I also removed the print line from both julia and C++ codes.  
Julia : ~ 17 micro  
C++ : ~ 13 micro  
and we can ignore this 4 microseconds. Acceptable. Thanks.