# My julia code is somehow much slower than the matlab code

**URL:** <https://discourse.julialang.org/t/my-julia-code-is-somehow-much-slower-than-the-matlab-code/91792>\
**Category:** New to Julia\
**Tags:** question, performance, matlab\
**Created:** [December 17, 2022, 9:01pm UTC](https://discourse.julialang.org/t/my-julia-code-is-somehow-much-slower-than-the-matlab-code/91792 "2022-12-17T21:01:15Z")\
**Posts on this page:** 16\
**Page:** 3

<div class="post-metadata">

**Author:** ![Wooheon](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/wooheon/32/206569_2.png) [@Wooheon](https://discourse.julialang.org/u/Wooheon)\
**Post date:** [December 23, 2022, 9:58pm UTC](https://discourse.julialang.org/t/my-julia-code-is-somehow-much-slower-than-the-matlab-code/91792/52 "2022-12-23T21:58:48Z")

</div>

Thank you for a lot of relevant information about my major! A lot of things to study, phew.  
I also found out that your code with a non-generic version is as fast as @Seif_Shebl 's code above. It seemed like the generic version needed some loading time to operate.  
Honestly, I have no idea how the codes work differently from mine in inner operation, but I’ll try to understand it!🙂  
Thanks a lot for your support.

---

<div class="post-metadata">

**Author:** ![Wooheon](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/wooheon/32/206569_2.png) [@Wooheon](https://discourse.julialang.org/u/Wooheon)\
**Post date:** [December 23, 2022, 10:35pm UTC](https://discourse.julialang.org/t/my-julia-code-is-somehow-much-slower-than-the-matlab-code/91792/53 "2022-12-23T22:35:12Z")

</div>

I think my computer has some differences from yours. On my computer, this code was even a little bit slower than your previous one with 1.870 s while the previous one showed 1.530 s with `@btime`. But if the performance difference comes from the part about multithreading, I don’t know why is matlab code so consistently faster on my computer.  
I’m learning a lot, thx!!

---

<div class="post-metadata">

**Author:** ![Seif\_Shebl](https://avatars.discourse-cdn.com/v4/letter/s/eada6e/32.png) [@Seif\_Shebl](https://discourse.julialang.org/u/Seif_Shebl)\
**Post date:** [December 24, 2022, 12:04am UTC](https://discourse.julialang.org/t/my-julia-code-is-somehow-much-slower-than-the-matlab-code/91792/54 "2022-12-24T00:04:18Z")

</div>

This is weird, make sure you set the `--check-bounds=no` option when you call your Julia script. It’s well-known in the Julia community that any code relying on vectorized loops will be quite fast. Make sure the loops get vectorized by inspecting the emitted LLVM code via `@code_llvm main()`, for instance. You should see llvm code like `4×double` which means simd vectorization works.

Something strange is happening with your Julia code. Notice that the MATLAB code runs at the same speed on our computers, so the Julia code should also run at similar speeds. Maybe some people here can confirm before we look further into what is happening.

Also, what version of Julia you are using? CPU thread count? CPU type? I’m using latest Julia 1.9 alpha, and my computer has 4-core Intel i7-4790k CPU.

---

<div class="post-metadata">

**Author:** ![Wooheon](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/wooheon/32/206569_2.png) [@Wooheon](https://discourse.julialang.org/u/Wooheon)\
**Post date:** [December 24, 2022, 1:52am UTC](https://discourse.julialang.org/t/my-julia-code-is-somehow-much-slower-than-the-matlab-code/91792/55 "2022-12-24T01:52:51Z")

</div>

Can you give me an example of how to use `--check-bounds=no`? I don’t know ways to use it.  
Yes, I think it is weird since mine is a laptop that will perform worse than yours but has a similar run time for Matlab and not for Julia. My laptop spec is as below.

 ![image](https://global.discourse-cdn.com/julialang/original/3X/c/0/c0c4299f4f6b8203852ec135c8039f5a45f0ae8c.jpeg)  
Also, my Julia version is 1.73. Can this be a source of the difference? I am looking for ways to update Julia, but it needs a bit of time to understand this, too… 😅  
I also don’t see `4×double` from `@code_llvm main()`. But I did `using LoopVectorization`. What may be the problem?

---

<div class="post-metadata">

**Author:** ![Wooheon](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/wooheon/32/206569_2.png) [@Wooheon](https://discourse.julialang.org/u/Wooheon)\
**Post date:** [December 24, 2022, 1:59am UTC](https://discourse.julialang.org/t/my-julia-code-is-somehow-much-slower-than-the-matlab-code/91792/56 "2022-12-24T01:59:07Z")

</div>

My current code:  
I added the `println(i)` inside the loop to make the comparison fair.

```julia
using LinearAlgebra
using Statistics
using BenchmarkTools;
using LoopVectorization

function main()
    ## (a) ##
    # First, set up the given data#
    α = 0.261
    β = 0.99
    δ = 0.0176
    zGrid = [0.9526, 0.9760, 1, 1.0246, 1.0497]
    Π_z = [0.9149 0.0823 0.0028 0 0
           0.0206 0.9163 0.0618 0.0014 0
           0.0005 0.0412 0.9167 0.0412 0.0005
           0 0.0014 0.0618 0.9163 0.0206
           0 0 0.0028 0.0823 0.9149]
    # Next, get the kStar and kL and kH #
    kStar = (1 / α * (1 / β - 1 + δ))^(1 / (α - 1))
    kL = kStar * 0.85
    kH = 1.15 * kStar
    println(kStar)
    println(kL)
    println(kH)

    # Setting up an 499 points of kgrid #
    kGrid = range(kL, kH, 499)
    println(median(kGrid))

    ## (b) ##
    # First, set up the length variables #
    nk = length(kGrid)
    nz = length(zGrid)

    # Make settings for the iteration #
    # First, make an array for each z_i #
    k′Grid = repeat(kGrid, 1, nz, nk)
    kGridMat = permutedims(k′Grid, (3, 2, 1))
    k′GridMat = k′Grid .* (kGridMat .^ α .* zGrid' .+ kGridMat .* (1 - δ) .- k′Grid .> 0)
    replace!(k′GridMat, 0 => NaN) 
    logC = log.(kGridMat .^ α .* zGrid' .+ kGridMat .* (1 - δ) .- k′GridMat)
    replace!(logC, NaN => -Inf)

    # precision and initial distance #
    precision = 0.01
    distance = 100    

    # start iteration #
    i = 1
    V = zeros(nk, nz)
    VΠ_z = similar(V)
    logCC = similar(logC)    
    while distance > precision
        mul!(VΠ_z, V, Π_z')
        d = -Inf
        for k = 1:nk, j = 1:nz
            m = -Inf
            for i = 1:nk
                logCC[i,j,k] = logC[i,j,k] + VΠ_z[i,j] * β
                m = ifelse(m > logCC[i,j,k], m, logCC[i,j,k])
            end 
            distance = max(d, m-V[k,j])
            d = m - V[k,j]
            V[k,j] = m
        end
        i += 1
        println(i);
    end
    # Add one more iteration to find Tk_index
    _, k_index = findmax(logCC, dims=1)
    mul!(VΠ_z, V, Π_z')
    logCC .= logC .+ VΠ_z .* β
    _, Tk_index = findmax(logCC, dims=1)

    # Reviving k-rules from the k_index array #
    k′sol = transpose(reshape(k′Grid[Tk_index], 5, 499))
    k′sol_prev = transpose(reshape(k′Grid[k_index], 5, 499)) # matrix of previous capital to compare with solutions #
    dist_k′sol = maximum(abs(i-j) for (i,j) in zip(k′sol,k′sol_prev))

    # answers #
    println(i)
    println(distance)
    println(dist_k′sol)
end
@btime main()
@code_llvm main()

```

Thanks a lot for the whole thing, everyone.

---

<div class="post-metadata">

**Author:** ![Wooheon](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/wooheon/32/206569_2.png) [@Wooheon](https://discourse.julialang.org/u/Wooheon)\
**Post date:** [December 24, 2022, 2:27am UTC](https://discourse.julialang.org/t/my-julia-code-is-somehow-much-slower-than-the-matlab-code/91792/57 "2022-12-24T02:27:26Z")

</div>

I was not able to run your code but found out minor typo(changed `k_index = CartesianIndex{3}[;;;]` into `k_index = CartesianIndex{3}[]`) and finally found out that your `loop_over_withcache!` is also very fast (around 1.5 s in my laptop).  
I honestly don’t know what part made your code very different from mine. Can you explain what role `similar` has or anything that is essential in the difference of allocation?  
Thanks a lot!

---

<div class="post-metadata">

**Author:** ![Seif\_Shebl](https://avatars.discourse-cdn.com/v4/letter/s/eada6e/32.png) [@Seif\_Shebl](https://discourse.julialang.org/u/Seif_Shebl)\
**Post date:** [December 24, 2022, 2:48am UTC](https://discourse.julialang.org/t/my-julia-code-is-somehow-much-slower-than-the-matlab-code/91792/58 "2022-12-24T02:48:15Z")

</div>

> [@Wooheon](#):
>
> Can you give me an example of how to use `--check-bounds=no`?

I mean to start Julia like this:

```julia
julia.exe --optimize=3 -t auto --check-bounds=no 

```

Then, `include` your file as: `include("\path\to\file.jl")`

This will turn off bounds checking in the whole code. So, in debugging phase, it’s not recommended to call Julia like this. The other easier option for you now is just to add `@inbounds` in front of `while` and keep your workflow as is.

```julia
@inbounds while ..
    ...
end

```

The version of Julia is not so important, it makes a negligible difference.

This is how it looks when you issue `@code_llvm main()` in the REPL. Notice these `<4 x double>` in the middle of a lot of emitted code.

 ![image](https://global.discourse-cdn.com/julialang/original/3X/a/8/a809a944ac45bb40c8db1258830c05d769e50d21.png)

Finally my cpu is 4th gen, that’s 10 years old, but yours is 11th gen, so it should be powerful enough to compete with a desktop cpu that age.

---

<div class="post-metadata">

**Author:** ![Wooheon](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/wooheon/32/206569_2.png) [@Wooheon](https://discourse.julialang.org/u/Wooheon)\
**Post date:** [December 24, 2022, 6:58am UTC](https://discourse.julialang.org/t/my-julia-code-is-somehow-much-slower-than-the-matlab-code/91792/59 "2022-12-24T06:58:35Z")

</div>

Wow, I don’t know why but as I put `@inbounds while ` as you said, the `@btime` showed me 887.031 ms, which is definitely faster than my Matlab code.  
Thanks a lot!  
As a novice in Julia, I am very confused since, in my eyes, your code contains a lot of **for loop** (of course, it is a vectorized loop but still it is a loop), which until now, with my knowledge in R and Matlab, was something that should be avoided. Could you explain what is happening? Also, why did my first code be so heavy, but your code is so light in memory?

---

<div class="post-metadata">

**Author:** ![nilshg](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/nilshg/32/2283_2.png) [@nilshg](https://discourse.julialang.org/u/nilshg)\
**Post date:** [December 24, 2022, 8:27am UTC](https://discourse.julialang.org/t/my-julia-code-is-somehow-much-slower-than-the-matlab-code/91792/60 "2022-12-24T08:27:58Z")

</div>

For loops are fast in Julia, like they are in compiled languages like C or C++. There’s almost never a reason to avoid them, and they can often be the simplest way of writing highly performant code.

---

<div class="post-metadata">

**Author:** ![jmair](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/jmair/32/35117_2.png) [@jmair](https://discourse.julialang.org/u/jmair)\
**Post date:** [December 24, 2022, 9:56am UTC](https://discourse.julialang.org/t/my-julia-code-is-somehow-much-slower-than-the-matlab-code/91792/61 "2022-12-24T09:56:18Z")

</div>

> [@Seif\_Shebl](#):
>
> The version of Julia is not so important, it makes a negligible difference.

Sometimes, this is worth checking. Just yesterday, updating to the latest Julia version gave me a 70% performance increase on one of my tasks. This is obviously not always guaranteed but sometimes it can make a difference.

> [@Wooheon](#):
>
> I am looking for ways to update Julia, but it needs a bit of time to understand this, too… 😅

There’s a tool called `juliaup` which lets you update julia seamlessly. You can read about it on [the GitHub page](https://github.com/JuliaLang/juliaup). `juliaup` lets you have multiple versions of Julia installed at the same time, so you can always try the latest version with the option of reverting back if something in your code breaks/is slower (unlikely but it could happen).

---

<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:** [December 24, 2022, 11:17am UTC](https://discourse.julialang.org/t/my-julia-code-is-somehow-much-slower-than-the-matlab-code/91792/62 "2022-12-24T11:17:54Z")

</div>

> [@Wooheon](#):
>
> **for loop** (of course, it is a vectorized loop but still it is a loop), which until now, with my knowledge in R and Matlab, was something that should be avoided.

Those ‘vectorized’ calls in Matlab or R just send the inputs to a C/Fortran library that uses a loop. It’s all loops, Matlab/R/numpy just hide them from you.

If you use a _fast_ language, like C/Fortran/Julia, you don’t need to send your data somewhere else to do the loops for you, you just write the loops yourself, directly.

---

<div class="post-metadata">

**Author:** ![fatteneder](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/fatteneder/32/33991_2.png) [@fatteneder](https://discourse.julialang.org/u/fatteneder)\
**Post date:** [December 24, 2022, 2:30pm UTC](https://discourse.julialang.org/t/my-julia-code-is-somehow-much-slower-than-the-matlab-code/91792/63 "2022-12-24T14:30:41Z")

</div>

> [@Wooheon](#):
>
> changed `k_index = CartesianIndex{3}[;;;]` into `k_index = CartesianIndex{3}[]`

Not really a typo, but this is a feature that is only available with Julia `v1.8` or newer: [add syntax for empty n-dimensional arrays by simeonschaub · Pull Request #41618 · JuliaLang/julia · GitHub](https://github.com/JuliaLang/julia/issues/41618)

> Can you explain what role `similar` has or anything that is essential in the difference of allocation?

`similar(V)` alone is not responsible for the performance gain. I could have also used `zeros(n_k,n_z)` and its purpose is just to create another array that has the same size and element type as `V`. Try `julia> ?similar` to learn more about it.

What I tried to do in my version is to avoid computations like `logC::Array{Float64,3}.+V*transpose(Pi_z).*beta` from the `loop_body` functions to create a new array, hence, an allocation, whenever you call these functions, which you do in every step of the `while` loop. So my version sets up an array called `tmp_product` with the right size before the `while` loop and I then forward it to `loop_body_withcache!`. In there I use the `.=` assignment (notice the dot!!!) to tell julia that it should store the result of the rhs computation into `tmp_product`. Using only a `=` would have discarded the buffer I set up before.

The effect of this is that the program needs around a factor of x100 less memory to do the same work. When I wrote the code I knew one could shave off some more allocations this way by avoiding `findmax` to create new arrays all the time, but I was too lazy on my train 😄. Luckily @Seif_Shebl showed how it can be done.

---

<div class="post-metadata">

**Author:** ![Wooheon](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/wooheon/32/206569_2.png) [@Wooheon](https://discourse.julialang.org/u/Wooheon)\
**Post date:** [December 25, 2022, 11:59pm UTC](https://discourse.julialang.org/t/my-julia-code-is-somehow-much-slower-than-the-matlab-code/91792/64 "2022-12-25T23:59:52Z")

</div>

That’s awesome, considering that I struggling to turn ‘for loop easy’ problems into arrays. Thank you for explaining it.

---

<div class="post-metadata">

**Author:** ![Wooheon](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/wooheon/32/206569_2.png) [@Wooheon](https://discourse.julialang.org/u/Wooheon)\
**Post date:** [December 27, 2022, 8:07am UTC](https://discourse.julialang.org/t/my-julia-code-is-somehow-much-slower-than-the-matlab-code/91792/65 "2022-12-27T08:07:18Z")

</div>

Thank you for the explanation. I feel shame about my past time spent on changing loops into vectors. 😂

---

<div class="post-metadata">

**Author:** ![Wooheon](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/wooheon/32/206569_2.png) [@Wooheon](https://discourse.julialang.org/u/Wooheon)\
**Post date:** [December 30, 2022, 2:38pm UTC](https://discourse.julialang.org/t/my-julia-code-is-somehow-much-slower-than-the-matlab-code/91792/66 "2022-12-30T14:38:36Z")

</div>

Thanks a lot for your explanation. I am getting to understand some stuff (though I’m still working on this understanding 😂).

I feel desperate to understand the whole thing and reproduce these myself since, if not, this means I cannot benefit from learning one of the fastest languages. In this regard, I thank you for giving me valuable materials to learn.

I have some questions, by the way. In your code, I found out that you set `make_cache` function for the `V, TV, tmp_product ` and then again changed them with all zeros in `loop_withcache!` function. Is there a particular reason for this? (The reason why you made a separate function for generating the right sized `V, TV, tmp_product`.)

Also, I am confused about the relationship between allocation and arrays. First, I thought avoiding `findmax` with an array inside was the key. But next, I found out that @Seif_Shebl 's code also happened to have two of the `findmax` that contains an array in lines 74 and 75. Why are these not damaging the performance? Are these arrays not consuming data allocation as my original arrays?

Lastly, what actually makes the difference between `.=` and `=` between arrays? Is the former consuming less data? or is it just related to performance? (or both?)

---

<div class="post-metadata">

**Author:** ![fatteneder](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/fatteneder/32/33991_2.png) [@fatteneder](https://discourse.julialang.org/u/fatteneder)\
**Post date:** [December 30, 2022, 4:26pm UTC](https://discourse.julialang.org/t/my-julia-code-is-somehow-much-slower-than-the-matlab-code/91792/67 "2022-12-30T16:26:35Z")

</div>

> [@Wooheon](#):
>
> I have some questions, by the way. In your code, I found out that you set `make_cache` function for the `V, TV, tmp_product ` and then again changed them with all zeros in `loop_withcache!` function. Is there a particular reason for this? (The reason why you made a separate function for generating the right sized `V, TV, tmp_product`.)

To be fair, my first version did not have the `fill!` in there. If you leave that out you will see that running

```julia
cache = make_cache(n_k, n_z)
(Tk_index,distance,k_index,i) = loop_withcache!(startdistance,precision,cache)

```

will only be slow the first time you run it, but all subsequent runs will be much faster. The reason is that `V, TV, tmp_product` will contain the results from the previous run and, hence, there will be no need to iterate it again. If you look at my EDIT of that post (click the little pencile at the beginning of my answer) you see that initially I reported a x1000 speedup, because of that.  
To make the benchmark fair again, I then added the `fill!`s so that the buffers are reset every time when the function is run and, thus, every call does ‘real’ work.

> (The reason why you made a separate function for generating the right sized `V, TV, tmp_product`.)

The reason is that I can create the cache before I call `loop_withcache!`. Otherwise `@benchmark` will count those allocations too.  
Let me say that you could also have created the cache inside `loop_withcache!` and the benchmark would have still displayed less allocations overall compared to your versions, because the really saving comes from using `.=` in the inner loop. However, creating the cache outside of the function allows you to reuse it multiple times. E.g. imagine calling `loop_withcache!` inside another loop many times, each iteration with with slightly modified parameters. Preparing the cache separately would allow you to reuse it here too.

> Also, I am confused about the relationship between allocation and arrays. First, I thought avoiding `findmax` with an array inside was the key. But next, I found out that @Seif_Shebl 's code also happened to have two of the `findmax` that contains an array in lines 74 and 75. Why are these not damaging the performance? Are these arrays not consuming data allocation as my original arrays?

I guess you are referring to this one: [My julia code is somehow much slower than the matlab code - #49 by Seif\_Shebl](https://discourse.julialang.org/t/my-julia-code-is-somehow-much-slower-than-the-matlab-code/91792/49)

In my version `findmax` is called in every iteration of the `while` loop, of which we need 412 (see the benchmark results). On the other hand, @Seif_Shebl s version calls `findmax` function exactly once after the _expensive_ `while` loop is done. It is not that `findmax` needs to be avoided at all costs in Julia, but it should be avoided running it in a _hot_ loop (a loop that needs to do lots of work) where you cannot cache the result. As his example demonstrates, the computational costs of the `while` loop clearly exceeds that of running `findmax` once, when compared to my version.

For the next points I wanna say that I am not 100% an expert with Julia’s memory model and my vocabulary might be a bit off. Others should please correct me if I am using wrong language here.

> Lastly, what actually makes the difference between `.=` and `=` between arrays? Is the former consuming less data? or is it just related to performance? (or both?)

Consider

```julia
julia> N = 2;

julia> a = zeros(N)
2-element Vector{Float64}:
 0.0
 0.0

julia> b = ones(N)
2-element Vector{Float64}:
 1.0
 1.0

julia> a = b
2-element Vector{Float64}:
 1.0
 1.0

julia> a[1] = 2;

julia> a
2-element Vector{Float64}:
 2.0
 1.0

julia> b
2-element Vector{Float64}:
 2.0
 1.0

julia> a = zeros(N)
2-element Vector{Float64}:
 0.0
 0.0

julia> b = ones(N)
2-element Vector{Float64}:
 1.0
 1.0

julia> a .= b
2-element Vector{Float64}:
 1.0
 1.0

julia> a[1] = 2;

julia> a
2-element Vector{Float64}:
 2.0
 1.0

julia> b
2-element Vector{Float64}:
 1.0
 1.0

```

As you can see, using `=` will make `a, b` refer to the same array: If you change an element of `a` this will be reflected in `b` and vice versa. Note that in this example the array to which `a` referred to _before_ the assignment (the zeros) are now lost (they might be still around if another variable refers to them, but this not the case in our example). In fact, you would not have needed to setup `a` with `zeros(N)`.  
However, if you use `.=` instead, then the contents of `b` will be copied over to `a`: changing elements of `a` will not be reflected in `b` and vice versa. For this to work, it is important that before running `a .= b` there exists an `a` such that `size(a) == size(b)`. The benefit of this is that Julia does not have to dispose of the storage to which `a` referred to before, but instead just overwrites its contents.

If you understood this concept than it should be easy to see why using `.=` above caused less allocations.

[Previous page](https://discourse.julialang.org/t/my-julia-code-is-somehow-much-slower-than-the-matlab-code/91792.md?page=2)