### Author:  [@Wooheon](https://discourse.julialang.org/u/Wooheon)
#### Post date: [December 17, 2022, 9:01pm UTC](https://discourse.julialang.org/t/my-julia-code-is-somehow-much-slower-than-the-matlab-code/91792/1 "2022-12-17T21:01:15Z")
Hi, I’m new to Julia. So I used my existing MatLab code to make a counterpart for Julia.
I heard Julia prefers the performance-demanding part inside a function, so I did it.
However, still, the performance is far lagging behind that of MatLab. I want to know what makes my Julia code inefficient.
The followings are MatLab code and Julia code each.
```julia
%First, set up given data%
alpha=0.261;beta=0.99;delta=0.0176;n_star=1;
z_grid=[0.9526;0.9760;1;1.0246;1.0497];
Pi_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 k_star and k_L and k_H%
k_star=(1/alpha*(1/beta-1+delta))^(1/(alpha-1));
k_L=k_star*0.85;k_H=1.15*k_star;
disp(round(k_star,4));disp(round(k_L,4));disp(round(k_H,4));
%Setting up an 499 points of kgrid%
k_grid=(k_L:((k_H-k_L)/498):k_H).';
disp(median(k_grid));
%%(b)%%
%First, set up the length variables%
n_k=length(k_grid);n_z=length(z_grid);
%Now initialize a V array%
V=zeros(n_k,n_z);
%Make settings for the iteration%
%First, make an array for each z_i%
k_prime_grid_mat_proto=repmat(k_grid,[1 n_z n_k]);%build prototype k' array with 3-d block array%
k_grid_mat=permute(k_prime_grid_mat_proto,[3 2 1]);%build k array with 3-d block array%;
k_prime_grid_mat=k_prime_grid_mat_proto.*(k_grid_mat.^alpha.*z_grid.'+k_grid_mat.*(1-delta)-k_prime_grid_mat_proto>=0);
k_prime_grid_mat(k_prime_grid_mat==0)=nan;
logC=log(k_grid_mat.^alpha.*z_grid.'+k_grid_mat.*(1-delta)-k_prime_grid_mat);
%build k' array that assigns nan to k' that is out of bound%
precision=0.01;distance=100;sz=size(k_prime_grid_mat);%precision and initial distance%
i=1;%iteration number%
%Start iteration%
while (distance > precision)
[TV,Tk_index]=max(logC+repmat(V*Pi_z.',[1 1 n_k]).*beta,[],'omitnan');%max within column omitting nan%
TV=reshape(permute(TV,[1,3,2]),[],size(TV,2));%change 1x5x499 array into 499x5 matrix%
distance=max(max(TV-V));
i=i+1;disp(i);disp(distance);
if distance0);
k_prime_grid_mat[k_prime_grid_mat.==0].=NaN;#build k' array that is out of bound#
logC=log.(k_grid_mat.^alpha.*transpose(z_grid)+k_grid_mat.*(1-delta)-k_prime_grid_mat);
logC[isnan.(logC)].=-Inf;
const precision=0.01;
const startdistance=100;#precision and initial distance#
#iteration number#
#start iteration#
function loop_over_global(distance, precision)
#Now initialize a V array#
global V=zeros(n_k,n_z);global i=1;
while distance>precision
global (TV,Tk_index)=findmax(logC::Array{Float64, 3}+repeat(V*transpose(Pi_z),1, 1, n_k).*beta,dims=1);
global TV=transpose(reshape(TV,5,499));
distance=maximum(TV-V);
println(i);println(distance);
if distance≤precision
break
else
global V=TV; global k_index=Tk_index; i=i+1;
end
end
return (TV,V,Tk_index,k_index,distance);
end
global (TV,V,Tk_index,k_index,distance)=loop_over_global(startdistance,precision);
#Reviving k-rules from the k_index array#
k_prime_sol_mat=transpose(reshape(k_prime_grid_mat_proto[Tk_index],5,499));
k_prime_sol_mat_previous=transpose(reshape(k_prime_grid_mat_proto[k_index],5,499));#matrix of previous capital to compare with solutions#
distance_k_prime_sol=maximum(abs.(k_prime_sol_mat-k_prime_sol_mat_previous));
##for information; how to convert a cartesian index array into matrix with only necessary info(first entry of the cartesian index)##
Tk_index=first.(Tuple.(Tk_index));Tk_index=transpose(reshape(Tk_index,5,499));
#answers#
println(i);println(distance);println(distance_k_prime_sol);
```
Thanks for reading this cumbersome codes. I wish I could get fluent in this fantastic language.
---
### Author:  [@tfiers](https://discourse.julialang.org/u/tfiers)
#### Post date: [December 17, 2022, 9:05pm UTC](https://discourse.julialang.org/t/my-julia-code-is-somehow-much-slower-than-the-matlab-code/91792/2 "2022-12-17T21:05:39Z")
I haven’t grokked the code, but as a first idea here: it’s the inner loop _body_ that should mostly be in a function. Can you try that (extract the body of the `while` to a function) and see if it runs faster?
---
### Author:  [@Wooheon](https://discourse.julialang.org/u/Wooheon)
#### Post date: [December 17, 2022, 9:41pm UTC](https://discourse.julialang.org/t/my-julia-code-is-somehow-much-slower-than-the-matlab-code/91792/4 "2022-12-17T21:41:19Z")
> [@Wooheon](#):
>
> I didn’t quite get the idea of extracting the body of while loop into a function. I thought this was what I’d done.
> Are you talking about the second code? First one is MatLabs.
I didn’t quite get the idea of extracting the body of while loop into a function. I thought this was what I’d done.
Are you talking about the second code? First one is MatLabs.
---
### Author:  [@tfiers](https://discourse.julialang.org/u/tfiers)
#### Post date: [December 18, 2022, 7:36am UTC](https://discourse.julialang.org/t/my-julia-code-is-somehow-much-slower-than-the-matlab-code/91792/5 "2022-12-18T07:36:05Z")
The idea is:
```julia
function main(…)
while {condition}
{result} = loop_body(…)
end
end
function loop_body(…)
…
end
```
---
### Author:  [@Krastanov](https://discourse.julialang.org/u/Krastanov)
#### Post date: [December 18, 2022, 7:44am UTC](https://discourse.julialang.org/t/my-julia-code-is-somehow-much-slower-than-the-matlab-code/91792/6 "2022-12-18T07:44:51Z")
On first sight: You are creating global variables in the inner loop. That would be bad in any language on purely stylistic grounds, but in Julia working with non-const untyped global variables is particularly bad for performance.
Have you had a chance to read [Performance Tips · The Julia Language](https://docs.julialang.org/en/v1/manual/performance-tips/) ? It is a particularly important part of the documentation for people coming from other languages, as Julia slightly different habits than Python and Matlab.
Side note: you are printing from inside your inner loop. That would throw off the measurement results both for the Matlab and Julia examples.
---
### Author:  [@Krastanov](https://discourse.julialang.org/u/Krastanov)
#### Post date: [December 18, 2022, 7:51am UTC](https://discourse.julialang.org/t/my-julia-code-is-somehow-much-slower-than-the-matlab-code/91792/7 "2022-12-18T07:51:26Z")
A quick guess: This line `distance=maximum(TV-V);` probably causes significant allocations that will be slowing down your code. `TV-V` would be allocating a new array each time, which would be much slower than the rest of the code combined. I.e. reserving the memory would be slower than actually running the computation of the values to put in that memory.
Also, could you share how you were benchmarking? I would strongly suggest using `@btime` from `BenchmarkTools.jl` when measuring the speed of julia code, in order to make sure you are not measuring compilation time.
---
### Author:  [@mbauman](https://discourse.julialang.org/u/mbauman)
#### Post date: [December 19, 2022, 3:33pm UTC](https://discourse.julialang.org/t/my-julia-code-is-somehow-much-slower-than-the-matlab-code/91792/8 "2022-12-19T15:33:29Z")
10 posts were split to a new topic: [Digression on short variable names](https://discourse.julialang.org/t/digression-on-short-variable-names/91859)
---
### Author:  [@ufechner7](https://discourse.julialang.org/u/ufechner7)
#### Post date: [December 18, 2022, 8:58am UTC](https://discourse.julialang.org/t/my-julia-code-is-somehow-much-slower-than-the-matlab-code/91792/11 "2022-12-18T08:58:52Z")
In Julia there is no need for a semicolon at the end of each line.
---
### Author:  [@Wooheon](https://discourse.julialang.org/u/Wooheon)
#### Post date: [December 19, 2022, 5:36am UTC](https://discourse.julialang.org/t/my-julia-code-is-somehow-much-slower-than-the-matlab-code/91792/19 "2022-12-19T05:36:39Z")
Is there a way to calculate the supnorm of two functions avoiding such a way like `distance=maximum(TV-V)`?
I just measured with my phone, but the difference was so significant that the measurement didn’t matter.
But I would like to use the package from now on, Thanks!
---
### Author:  [@Wooheon](https://discourse.julialang.org/u/Wooheon)
#### Post date: [December 19, 2022, 6:21am UTC](https://discourse.julialang.org/t/my-julia-code-is-somehow-much-slower-than-the-matlab-code/91792/20 "2022-12-19T06:21:48Z")
Thank you for the suggestion. I tried to change the code as below.
```plaintext
#Now initialize a V array#
V=zeros(n_k,n_z);i=1;k_index=Array{CartesianIndex{3}, 3}
while distance>precision
k_index=k_index;
(TV,Tk_index,distance)=loop_body(V);
println(i);println(distance);
if distance≤precision
return (TV,V,Tk_index,distance,k_index,i);
break
else
V=TV; k_index=Tk_index::Array{CartesianIndex{3}, 3}; i=i+1;
end
end
end
function loop_body(V)
(A,B)=findmax(logC::Array{Float64,3}+repeat(V*transpose(Pi_z),1,1,n_k).*beta,dims=1);
A1=transpose(reshape(A,size(V)[2],size(V)[1]));
C=maximum(A1-V);
return (A1,B,C)
(TV,V,Tk_index,distance,k_index,i)=loop_over_global(startdistance,precision);
end
```
However, I think there was not much difference from the previous one. Is there something I missed?
---
### Author:  [@DNF](https://discourse.julialang.org/u/DNF)
#### Post date: [December 19, 2022, 10:04am UTC](https://discourse.julialang.org/t/my-julia-code-is-somehow-much-slower-than-the-matlab-code/91792/21 "2022-12-19T10:04:37Z")
> [@Wooheon](#):
>
> `(A,B)=findmax(logC::Array{Float64,3}+repeat(V*transpose(Pi_z),1,1,n_k).*beta,dims=1)`
I’m on holiday without my computer, so can’t be very helpful. But I can tell that this is inefficient. You never need to use `repeat`, and should not. This is one of the good things about broadcasting. You should not create arrays with many repeated entries. Also, get rid of the `Array{Float64,3}` annotation, it can only do harm.
Try something like
```julia
(A, B) = findmax(logC .+ V .* transpose(Pi_z) .* beta; dims=1)
```
I cannot test it, so cannot guarantee it works, maybe you need a `reshape`, for example.
Also, you don’t need to use `repeat` in Matlab either. They also support broadcasting, though implicitly, without dots.
Generally, just remove `repeat` from your programming vocabulary (in any programming language). Imho, you will never need it for anything, it’s just a waste of time and memory.
(Btw, your code would be more readable if you put spaces in your code, around operators, etc.)
---
### Author:  [@ufechner7](https://discourse.julialang.org/u/ufechner7)
#### Post date: [December 19, 2022, 1:02pm UTC](https://discourse.julialang.org/t/my-julia-code-is-somehow-much-slower-than-the-matlab-code/91792/22 "2022-12-19T13:02:07Z")
And please do not put ; at the end of the lines…
---
### Author:  [@giordano](https://discourse.julialang.org/u/giordano)
#### Post date: [December 19, 2022, 1:09pm UTC](https://discourse.julialang.org/t/my-julia-code-is-somehow-much-slower-than-the-matlab-code/91792/23 "2022-12-19T13:09:41Z")
> [@Wooheon](#):
>
> ` V=zeros(n_k,n_z);i=1;k_index=Array{CartesianIndex{3}, 3}`
If I’m reading this line correctly, you’re setting `k_index` equal to a datatype and then in the hot loop you’re changing its value to an _instance_ of that type.
---
### Author:  [@lmiq](https://discourse.julialang.org/u/lmiq)
#### Post date: [December 19, 2022, 1:42pm UTC](https://discourse.julialang.org/t/my-julia-code-is-somehow-much-slower-than-the-matlab-code/91792/24 "2022-12-19T13:42:35Z")
Can you post a complete running example of your code, and the time you are getting?
> [@Wooheon](#):
>
> ` C=maximum(A1-V);`
In terms of this, here are some options:
```julia
maximum(a1[i] - v[i] for i in eachindex(a1,v))
```
This is not faster for a single call, but does not allocate the intermediate array `A1-V`, and since this is being called from within the loop, that can be faster for the whole timing.
The fastest you can get for that operation is, probably something like this:
```julia
julia> using LoopVectorization
julia> function f2(a1,v)
m = zero(eltype(a1))
@turbo for i in eachindex(a1,v)
val = a1[i] - v[i]
m = ifelse(val > m, val, m)
end
return m
end
f2 (generic function with 1 method)
julia> @btime f2($a1, $v)
86.956 ns (0 allocations: 0 bytes)
0.9590868980321948
```
This is about 15 times faster than `maximum(a1-v)`.
(I don’t know how much time is being spent on this opeartion, though)
---
### Author:  [@DNF](https://discourse.julialang.org/u/DNF)
#### Post date: [December 19, 2022, 2:06pm UTC](https://discourse.julialang.org/t/my-julia-code-is-somehow-much-slower-than-the-matlab-code/91792/25 "2022-12-19T14:06:27Z")
> [@lmiq](#):
>
> `m = ifelse(val > m, val, m)`
Is this better than `m = max(m, val)`?
---
### Author:  [@TheLateKronos](https://discourse.julialang.org/u/TheLateKronos)
#### Post date: [December 19, 2022, 2:28pm UTC](https://discourse.julialang.org/t/my-julia-code-is-somehow-much-slower-than-the-matlab-code/91792/26 "2022-12-19T14:28:35Z")
A more readable version of the code, with more line breaks and no useless semicolons:
```julia
using Statistics
const alpha=0.261
const beta=0.99
const delta=0.0176
const n_star=1
const z_grid=[0.9526;0.9760;1;1.0246;1.0497]
const Pi_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 k_star and k_L and k_H#
const k_star=(1/alpha*(1/beta-1+delta))^(1/(alpha-1))
const k_L=k_star*0.85
const k_H=1.15*k_star
println(round(k_star,digits=4));println(round(k_L,digits=4));println(round(k_H,digits=4))
#Setting up an 499 points of kgrid#
const k_grid=[k_L:((k_H-k_L)/498):k_H;]
println(median(k_grid))
##(b)##
#First, set up the length variables#
const n_k=length(k_grid)
const n_z=length(z_grid)
#Make settings for the iteration#
#First, make an array for each z_i#
const k_prime_grid_mat_proto=repeat(k_grid,1, n_z, n_k)
const k_grid_mat=permutedims(k_prime_grid_mat_proto,(3,2,1))
k_prime_grid_mat=k_prime_grid_mat_proto.*(k_grid_mat.^alpha.*transpose(z_grid)+k_grid_mat.*(1-delta)-k_prime_grid_mat_proto.>0)
k_prime_grid_mat[k_prime_grid_mat.==0].=NaN#build k' array that is out of bound#
logC=log.(k_grid_mat.^alpha.*transpose(z_grid)+k_grid_mat.*(1-delta)-k_prime_grid_mat)
logC[isnan.(logC)].=-Inf
const precision=0.01
const startdistance=100#precision and initial distance#
#iteration number#
#start iteration#
function loop_over_global(distance, precision)
#Now initialize a V array#
global V=zeros(n_k,n_z)
global i=1
while distance>precision
global (TV,Tk_index)=findmax(logC::Array{Float64, 3}+repeat(V*transpose(Pi_z),1, 1, n_k).*beta,dims=1)
global TV=transpose(reshape(TV,5,499))
distance=maximum(TV-V)
println(i)
println(distance)
if distance≤precision
break
else
global V=TV
global k_index=Tk_index
i=i+1
end
end
return (TV,V,Tk_index,k_index,distance)
end
global (TV,V,Tk_index,k_index,distance)=loop_over_global(startdistance,precision)
#Reviving k-rules from the k_index array#
k_prime_sol_mat=transpose(reshape(k_prime_grid_mat_proto[Tk_index],5,499))
k_prime_sol_mat_previous=transpose(reshape(k_prime_grid_mat_proto[k_index],5,499))#matrix of previous capital to compare with solutions#
distance_k_prime_sol=maximum(abs.(k_prime_sol_mat-k_prime_sol_mat_previous))
##for information; how to convert a cartesian index array into matrix with only necessary info(first entry of the cartesian index)##
Tk_index=first.(Tuple.(Tk_index))
Tk_index=transpose(reshape(Tk_index,5,499))
#answers#
println(i)
println(distance)
println(distance_k_prime_sol)
```
Timing all of it takes around 11.43 seconds for me (on occation as low as 8)
Putting all that (except `using Statistics`) in a `begin`-`end` block, and annotating with `@profview`, lets me know that 88% of the time is spent in `poptask`:
By the way, removing all `const` annotations changes the runtime from around 11.43 to around 11.6 for me, so I will do that, as it allows `BenchmarkTools`’s `@btime`.
Removing the hot loop printing takes the time down to 10.3 s the first time, and 8.648 s the next time. There is clearly some variation in the runs, but it seems to cut about a second.
Simply removing the global annotations (which had no functionality!) takes it down to 5.8 seconds. A new `@profview` without constant and global annotations, and without inner loop printing, shows that _still_ 88% of the time is spendt in poptask
I am not sure what that is, nor how to fix it.
Zooming in on the 12% of the time spent on `loop_over_global` shows the following “flamegraph”. It can be a bit cryptic at first, but tells you everything about where the time is spent. It can probably help others spot the performance problems in the inner loop.
---
### Author:  [@lmiq](https://discourse.julialang.org/u/lmiq)
#### Post date: [December 19, 2022, 3:38pm UTC](https://discourse.julialang.org/t/my-julia-code-is-somehow-much-slower-than-the-matlab-code/91792/27 "2022-12-19T15:38:08Z")
In my tests that allowed vectorization, and `max` didn’t.
---
### Author:  [@dlakelan](https://discourse.julialang.org/u/dlakelan)
#### Post date: [December 19, 2022, 6:01pm UTC](https://discourse.julialang.org/t/my-julia-code-is-somehow-much-slower-than-the-matlab-code/91792/28 "2022-12-19T18:01:18Z")
This “pop\_task” stuff is because of unused threads.
---
### Author:  [@TheLateKronos](https://discourse.julialang.org/u/TheLateKronos)
#### Post date: [December 20, 2022, 9:45am UTC](https://discourse.julialang.org/t/my-julia-code-is-somehow-much-slower-than-the-matlab-code/91792/29 "2022-12-20T09:45:09Z")
I would also say that your code is clearly written in a matlab kind of way, and that a Julian formulation would probably look very different. But it takes some effort to see what your code is actually trying to do, which is a key part of making an improved rewrite. A general tip to make it easy to help you is to boil down the code as much as possible in a clear example, which is then easier to optimize.
One of the big problems is that the timing varies significantly - you need to get rid of this randomness to make it easy to optimize.
One smaller thing is e.g. the following line:
```julia
k_grid=[k_L:((k_H-k_L)/498):k_H;]
```
You are in fact (WRONG, I did not know the syntax well enough in writing this) making a 1-element vector containing a range as its only element, which does not seem to be what you want. Espesially concidering your calls to `length(k_grid)`. But if I remove the semicolon, your code breaks. It is just a messy code to understand and optimize, which is unrelated to Julia and Matlab.
---
### Author:  [@GunnarFarneback](https://discourse.julialang.org/u/GunnarFarneback)
#### Post date: [December 20, 2022, 10:00am UTC](https://discourse.julialang.org/t/my-julia-code-is-somehow-much-slower-than-the-matlab-code/91792/30 "2022-12-20T10:00:49Z")
> [@TheLateKronos](#):
>
> But if I remove the semicolon, your code breaks.
That’s because the semicolon is syntactically significant. Cf. the construction with two elements:
```julia
julia> [1:2, 3:5]
2-element Vector{UnitRange{Int64}}:
1:2
3:5
julia> [1:2; 3:5]
5-element Vector{Int64}:
1
2
3
4
5
```
Effectively `[a:b;]` is equivalent to `collect(a:b)` but most likely directly using the range `a:b` would work equally well and be more efficient.
[Next page](https://discourse.julialang.org/t/my-julia-code-is-somehow-much-slower-than-the-matlab-code/91792.md?page=2)