# Julian way to write this code
**URL:** https://discourse.julialang.org/t/julian-way-to-write-this-code/119348
**Category:** Performance
**Tags:** question
**Created:** [September 12, 2024, 7:09pm UTC](https://discourse.julialang.org/t/julian-way-to-write-this-code/119348 "2024-09-12T19:09:48Z")
**Posts on this page:** 19
**Page:** 1
### Author:  [@gdkrmr](https://discourse.julialang.org/u/gdkrmr)
#### Post date: [September 12, 2024, 7:09pm UTC](https://discourse.julialang.org/t/julian-way-to-write-this-code/119348/1 "2024-09-12T19:09:48Z")
Is there a pretty way to write this _and_ have good performance? `g1` produces runtime dispatch, `g2` produces an allocation and `g3` doesn’t look very “Julian” and is also unsafe.
```julia
function g1(x, y, w)
sum(i -> x[i] * y[i] * w[i], 1:length(x))
end
function g2(x, y, w)
sum(x .* y .* w)
end
function g3(x::AbstractArray{T}, y::AbstractArray{T}, w::AbstractArray{T}) where T
l = length(x)
r = T(0)
@inbounds for i in 1:l
r += x[i] * y[i] * w[i]
end
return r
end
function g4(x, y, w)
x .*= y
x .*= w
sum(x)
end
x = rand(3)
y = rand(3)
w = rand(3)
```
```julia-repl
julia> @benchmark g1($x, $y, $w)
BenchmarkTools.Trial: 10000 samples with 999 evaluations.
Range (min … max): 9.969 ns … 46.624 ns ┊ GC (min … max): 0.00% … 0.00%
Time (median): 12.405 ns ┊ GC (median): 0.00%
Time (mean ± σ): 12.423 ns ± 0.565 ns ┊ GC (mean ± σ): 0.00% ± 0.00%
▃▁ ▂ ▁ ▄▇█▇ ▂
▃▃▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁██▁█▁█▄▆████▅▆▇▅▆▄▄▁▃▄▃▄▆▅▅▅▅▅▆▆▆▇ █
9.97 ns Histogram: log(frequency) by time 14.1 ns <
Memory estimate: 0 bytes, allocs estimate: 0.
julia> @benchmark g2($x, $y, $w)
BenchmarkTools.Trial: 10000 samples with 991 evaluations.
Range (min … max): 38.842 ns … 115.352 ns ┊ GC (min … max): 0.00% … 0.00%
Time (median): 49.113 ns ┊ GC (median): 0.00%
Time (mean ± σ): 48.779 ns ± 4.394 ns ┊ GC (mean ± σ): 0.00% ± 0.00%
▁▅█▇▃
▁▁▁▂▂▂▂▂▃▃▂▂▂▂▂▃▃▃▃▄▆█████▇▃▂▂▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁ ▂
38.8 ns Histogram: frequency by time 65.6 ns <
Memory estimate: 80 bytes, allocs estimate: 1.
julia> @benchmark g3($x, $y, $w)
BenchmarkTools.Trial: 10000 samples with 1000 evaluations.
Range (min … max): 3.135 ns … 22.953 ns ┊ GC (min … max): 0.00% … 0.00%
Time (median): 3.367 ns ┊ GC (median): 0.00%
Time (mean ± σ): 3.433 ns ± 0.426 ns ┊ GC (mean ± σ): 0.00% ± 0.00%
▂█
▅▃▁▃▁▁▃▅▄▅▅▆███▃▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▁▁▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▃▂▂▂▂▃▃▂▃ ▃
3.14 ns Histogram: frequency by time 4.25 ns <
Memory estimate: 0 bytes, allocs estimate: 0.
julia> @benchmark g4($x, $y, $w)
BenchmarkTools.Trial: 10000 samples with 999 evaluations.
Range (min … max): 8.153 ns … 35.191 ns ┊ GC (min … max): 0.00% … 0.00%
Time (median): 8.876 ns ┊ GC (median): 0.00%
Time (mean ± σ): 8.961 ns ± 0.570 ns ┊ GC (mean ± σ): 0.00% ± 0.00%
▄ ▁ ▅▆▇▇█▅▄▃▁ ▁▁ ▁▂▃▃▂▂▃ ▁▁▁ ▁ ▂
█▄▁▁▇█▄▆▁▄█▇▄▁███████████████████████████▆▆▅▅▅▅▅▄▅▅▄▄▅▄▅▆▆ █
8.15 ns Histogram: log(frequency) by time 10.5 ns <
Memory estimate: 0 bytes, allocs estimate: 0.
```
---
### Author:  [@gdalle](https://discourse.julialang.org/u/gdalle)
#### Post date: [September 12, 2024, 7:17pm UTC](https://discourse.julialang.org/t/julian-way-to-write-this-code/119348/2 "2024-09-12T19:17:32Z")
I think `g3` comes closest but you can make it safer by removing `@inbounds` and using `eachindex` instead:
```julia
function g5(x, y, w)
T = promote_type(eltype(x), eltype(y), eltype(w))
r = zero(T)
for i in eachindex(x, y, w)
r += x[i] * y[i] * w[i]
end
return r
end
```
EDIT: weirdly enough this one is rather slow too, unsure why
---
### Author:  [@gdkrmr](https://discourse.julialang.org/u/gdkrmr)
#### Post date: [September 12, 2024, 7:26pm UTC](https://discourse.julialang.org/t/julian-way-to-write-this-code/119348/3 "2024-09-12T19:26:35Z")
didn’t know that `eachindex` could take more than one array.
it isn’t too bad:
```julia-rep
julia> @benchmark g5($x1, $y1, $w1)
BenchmarkTools.Trial: 10000 samples with 1000 evaluations.
Range (min … max): 5.240 ns … 39.784 ns ┊ GC (min … max): 0.00% … 0.00%
Time (median): 5.290 ns ┊ GC (median): 0.00%
Time (mean ± σ): 5.314 ns ± 0.444 ns ┊ GC (mean ± σ): 0.00% ± 0.00%
▅█▂
▂▄▃▆▃███▅▂▄▂▃▂▁▂▂▂▂▁▂▂▂▁▁▂▂▂▁▁▂▂▂▂▂▂▁▂▁▁▂▂▂▁▂▂▂▂▁▂▂▁▂▁▁▁▂▂ ▂
5.24 ns Histogram: frequency by time 5.7 ns <
Memory estimate: 0 bytes, allocs estimate: 0.
```
---
### Author:  [@gdalle](https://discourse.julialang.org/u/gdalle)
#### Post date: [September 12, 2024, 7:30pm UTC](https://discourse.julialang.org/t/julian-way-to-write-this-code/119348/4 "2024-09-12T19:30:36Z")
Also be careful not to benchmark `g4` that way, because it will modify your `x` to store very small numbers, and then subsequent benchmarks will be slower (presumably because of `NaN`s appearing?). I recommend you use a setup phase instead, see [Manual · BenchmarkTools.jl](https://juliaci.github.io/BenchmarkTools.jl/stable/manual/#Setup-and-teardown-phases)
---
### Author:  [@DanielVandH](https://discourse.julialang.org/u/DanielVandH)
#### Post date: [September 12, 2024, 7:35pm UTC](https://discourse.julialang.org/t/julian-way-to-write-this-code/119348/5 "2024-09-12T19:35:16Z")
You can also a `@simd` to your `g5` to improve it further.
```julia
julia> function g5(x, y, w)
T = promote_type(eltype(x), eltype(y), eltype(w))
r = zero(T)
for i in eachindex(x, y, w)
r += x[i] * y[i] * w[i]
end
return r
end
g5 (generic function with 1 method)
julia> function g6(x, y, w)
T = promote_type(eltype(x), eltype(y), eltype(w))
r = zero(T)
@simd for i in eachindex(x, y, w)
r += x[i] * y[i] * w[i]
end
return r
end
g6 (generic function with 1 method)
julia> @benchmark g5(x, y, z) (setup=(x = rand(100); y = rand(100); z = rand(100)))
BenchmarkTools.Trial: 10000 samples with 993 evaluations.
Range (min … max): 37.865 ns … 268.983 ns ┊ GC (min … max): 0.00% … 0.00%
Time (median): 38.469 ns ┊ GC (median): 0.00%
Time (mean ± σ): 40.368 ns ± 5.999 ns ┊ GC (mean ± σ): 0.00% ± 0.00%
▃█▄ ▂ ▁▅▅ ▁
██████▇██▇▇▇██████▇█▇▇▇▆▇▇▇▆▆▆▆▅▆▆▅▆▅▄▅▅▃▄▆▅▄▄▅▆▄▃▄▅▅▄▅▄▅▄▆▅ █
37.9 ns Histogram: log(frequency) by time 56.7 ns <
Memory estimate: 0 bytes, allocs estimate: 0.
julia> @benchmark g6(x, y, z) (setup=(x = rand(100); y = rand(100); z = rand(100)))
BenchmarkTools.Trial: 10000 samples with 999 evaluations.
Range (min … max): 12.112 ns … 137.337 ns ┊ GC (min … max): 0.00% … 0.00%
Time (median): 12.412 ns ┊ GC (median): 0.00%
Time (mean ± σ): 13.260 ns ± 2.603 ns ┊ GC (mean ± σ): 0.00% ± 0.00%
█▃ ▄ ▂▇▁ ▁▃▄▂▁ ▂
██▇▃▅█▇█████▇▅█████▆▇██▆▆▇▆▆▇▇▁▇▅▅▆▅▅▅▇▅▄▄▄▁▃▆▄▅▄▅▁▃▃▁▄▃▃▃▁▄ █
12.1 ns Histogram: log(frequency) by time 22.3 ns <
Memory estimate: 0 bytes, allocs estimate: 0.
```
---
### Author:  [@Satvik](https://discourse.julialang.org/u/Satvik)
#### Post date: [September 12, 2024, 8:21pm UTC](https://discourse.julialang.org/t/julian-way-to-write-this-code/119348/6 "2024-09-12T20:21:42Z")
Note the caveat that @simd can produce slightly different results because floating point addition is not associative, and it’s allowed to reorder the inputs:
```julia
x = rand(1000)
y = rand(1000)
w = rand(1000)
g5(x, y, w)
> 136.6347749999703
g6(x, y, w)
> 136.6347749999705
```
---
### Author:  [@Sevi](https://discourse.julialang.org/u/Sevi)
#### Post date: [September 12, 2024, 8:28pm UTC](https://discourse.julialang.org/t/julian-way-to-write-this-code/119348/7 "2024-09-12T20:28:07Z")
> [@gdkrmr](#):
>
> `g1` produces runtime dispatch
Is that so? Neither `@code_warntype` nor `@descent` (Cthulhu.jl) seem to indicate that for me. Can you share how to test that?
The performance difference might simply be due to the fact that `sum` does not exactly do the same as `g3`, but as far as I know uses a more accurate algorithm for summing a (large) sequence of numbers (please correct me if I’m wrong).
For small vectors like in your benchmarks, there might be no difference, but for longer ones you can check that the results are slightly different. Also `@simd` can change the output (ah, @Satvik beat me to it 😁 ).
Whether you need the extra accuracy or speed depends on your use case of course. If the vectors stay small, maybe static vectors are another option? Not sure if they will improve the summation performance though, just throwing it out there.
---
### Author:  [@bertschi](https://discourse.julialang.org/u/bertschi)
#### Post date: [September 12, 2024, 8:53pm UTC](https://discourse.julialang.org/t/julian-way-to-write-this-code/119348/8 "2024-09-12T20:53:39Z")
Imho, this is nice and quite fast:
```julia
g7(x, y, z) = sum(splat(*), zip(x, y, z))
```
```julia-repl
julia> @btime g5(x, y, z) (setup=(x = rand(100); y = rand(100); z = rand(100)));
210.963 ns (0 allocations: 0 bytes)
julia> @btime g7(x, y, z) (setup=(x = rand(100); y = rand(100); z = rand(100)));
223.047 ns (0 allocations: 0 bytes)
```
---
### Author:  [@alfaromartino](https://discourse.julialang.org/u/alfaromartino)
#### Post date: [September 13, 2024, 4:10am UTC](https://discourse.julialang.org/t/julian-way-to-write-this-code/119348/9 "2024-09-13T04:10:07Z")
One more option with `LazyArrays` (with the caveat of `@simd`)
```julia
using LazyArrays
g8(x, y, w) = sum(@~ x .* y .* w)
@btime g5(x, y, z) (setup=(x = rand(100); y = rand(100); z = rand(100)))
@btime g8(x, y, z) (setup=(x = rand(100); y = rand(100); z = rand(100)))
```
with times
```julia
#g5 32.425 ns (0 allocations: 0 bytes)
#g8 17.512 ns (0 allocations: 0 bytes)
```
---
### Author:  [@GunnarFarneback](https://discourse.julialang.org/u/GunnarFarneback)
#### Post date: [September 13, 2024, 5:43am UTC](https://discourse.julialang.org/t/julian-way-to-write-this-code/119348/10 "2024-09-13T05:43:51Z")
> [@Satvik](#):
>
> Note the caveat that @simd can produce slightly different results because floating point addition is not associative, and it’s allowed to reorder the inputs:
This is good to keep in mind but at a point of investigation where the alternatives involve calling `sum`, the `@simd` reordering is a minor effect and more likely than not to be numerically beneficial compared to the straight accumulation.
---
### Author:  [@gdkrmr](https://discourse.julialang.org/u/gdkrmr)
#### Post date: [September 13, 2024, 6:19am UTC](https://discourse.julialang.org/t/julian-way-to-write-this-code/119348/11 "2024-09-13T06:19:51Z")
This is what I get for `g1`. `@code_warntype` says it is fine, haven’t tested `@descent`
```julia-repl
julia> @report_opt g1(x1, y1, w1)
═════ 1 possible error found ═════
┌ g1(x::Vector{Float64}, y::Vector{Float64}, w::Vector{Float64}) @ Main ./REPL[32]:2
│┌ sum(f::var"#5#6"{Vector{Float64}, Vector{Float64}, Vector{Float64}}, a::UnitRange{Int64}) @ Base ./reducedim.jl:1011
││┌ sum(f::var"#5#6"{Vector{Float64}, Vector{Float64}, Vector{Float64}}, a::UnitRange{Int64}; dims::Colon, kw::@Kwargs
{}) @ Base ./reducedim.jl:1011
│││┌ _sum(f::var"#5#6"{Vector{Float64}, Vector{Float64}, Vector{Float64}}, a::UnitRange{Int64}, ::Colon) @ Base ./reducedim.jl:1015
││││┌ _sum(f::var"#5#6"{Vector{Float64}, Vector{Float64}, Vector{Float64}}, a::UnitRange{Int64}, ::Colon; kw::@Kwargs{}) @ Base ./reducedim.jl:1015
│││││┌ mapreduce(f::var"#5#6"{Vector{…}, Vector{…}, Vector{…}}, op::typeof(Base.add_sum), A::UnitRange{Int64}) @ Base ./reducedim.jl:357
││││││┌ mapreduce(f::var"#5#6"{…}, op::typeof(Base.add_sum), A::UnitRange{…}; dims::Colon, init::Base._InitialValue) @ Base ./reducedim.jl:357
│││││││┌ _mapreduce_dim(f::var"#5#6"{…}, op::typeof(Base.add_sum), ::Base._InitialValue, A::UnitRange{…}, ::Colon) @ Base ./reducedim.jl:365
││││││││┌ _mapreduce(f::var"#5#6"{Vector{…}, Vector{…}, Vector{…}}, op::typeof(Base.add_sum), ::IndexLinear, A::UnitRange{Int64}) @ Base ./reduce.jl:432
│││││││││┌ mapreduce_empty_iter(f::var"#5#6"{…}, op::typeof(Base.add_sum), itr::UnitRange{…}, ItrEltype::Base.HasEltype) @ Base ./reduce.jl:380
││││││││││┌ reduce_empty_iter(op::Base.MappingRF{var"#5#6"{…}, typeof(Base.add_sum)}, itr::UnitRange{Int64}, ::Base.HasEltype) @ Base ./reduce.jl:384
│││││││││││┌ reduce_empty(op::Base.MappingRF{var"#5#6"{Vector{…}, Vector{…}, Vector{…}}, typeof(Base.add_sum)}, ::Type
{Int64}) @ Base ./reduce.jl:361
││││││││││││ runtime dispatch detected: Base.mapreduce_empty(%1::var"#5#6"{Vector{Float64}, Vector{Float64}, Vector{Float64}}, add_sum, ::Int64)
│││││││││││└────────────────────
```
The vectors are quite short, so I don’t think I will benefit from `@simd`.
---
### Author:  [@LaurentPlagne](https://discourse.julialang.org/u/LaurentPlagne)
#### Post date: [September 13, 2024, 6:29am UTC](https://discourse.julialang.org/t/julian-way-to-write-this-code/119348/12 "2024-09-13T06:29:25Z")
> [@gdkrmr](#):
>
> The vectors are quite short, so I don’t think I will benefit from `@simd`
simd instructions operate on very small pack of data (2,4,8) and can be beneficial for small vectors (no overhead).
For this operation, if the data is not already in cache the major price will be to read the data from the RAM. It may be non trivial to evaluate because the benchmark tends to create an artificial temporal locality that may or may not represents your actual computation.
---
### Author:  [@lmiq](https://discourse.julialang.org/u/lmiq)
#### Post date: [September 13, 2024, 7:01am UTC](https://discourse.julialang.org/t/julian-way-to-write-this-code/119348/13 "2024-09-13T07:01:58Z")
> [@gdalle](#):
>
> EDIT: weirdly enough this one is rather slow too, unsure why
afaik, you should keep the `@inbounds` there, which is then safe because it is guaranteed. As of today the each index used that way is not used by the compiler to avoid bounds checking.
---
### Author:  [@abraemer](https://discourse.julialang.org/u/abraemer)
#### Post date: [September 13, 2024, 7:10am UTC](https://discourse.julialang.org/t/julian-way-to-write-this-code/119348/14 "2024-09-13T07:10:20Z")
> [@lmiq](#):
>
> afaik, you should keep the `@inbounds` there, which is then safe because it is guaranteed. As of today the each index used that way is not used by the compiler to avoid bounds checking.
That’s incorrect (at least for my Julia 1.10.5). Generally using `eachindex` signals the compiler that the index access are in fact `inbounds`. This is the great advantage of using it in the first place 🙂 You get **safe** `@inbounds` from it!
I checked this example and adding `@inbounds` gives the exact same LLVM-code.
---
### Author:  [@ffevotte](https://discourse.julialang.org/u/ffevotte)
#### Post date: [September 13, 2024, 9:48am UTC](https://discourse.julialang.org/t/julian-way-to-write-this-code/119348/15 "2024-09-13T09:48:27Z")
Here is yet another suggestion along the same lines as `g1` and `g7`:
```julia-repl
julia> g9(x,y,z) = sum(x[i] * y[i] * z[i] for i in eachindex(x, y, z))
g9 (generic function with 1 method)
julia> x, y, z = rand(3), rand(3), rand(3);
julia> @btime g9($x, $y, $z)
4.472 ns (0 allocations: 0 bytes)
0.7977237806641253
```
---
### Author:  [@lmiq](https://discourse.julialang.org/u/lmiq)
#### Post date: [September 13, 2024, 11:53am UTC](https://discourse.julialang.org/t/julian-way-to-write-this-code/119348/16 "2024-09-13T11:53:53Z")
> [@abraemer](#):
>
> I checked this example and adding `@inbounds` gives the exact same LLVM-code.
Good to know, the last time I checked that was not the case (it isn´t, for instance, in Julia 1.6).
> **code**
>
> ```julia-repl
> julia> using BenchmarkTools
>
> julia> function g6(x, y, w)
> T = promote_type(eltype(x), eltype(y), eltype(w))
> r = zero(T)
> for i in eachindex(x, y, w)
> r += x[i] * y[i] * w[i]
> end
> return r
> end
> g6 (generic function with 1 method)
>
> julia> @btime g6($(rand(1000)),$(rand(1000)),$(rand(1000)))
> 518.948 ns (0 allocations: 0 bytes)
> 116.36617338389264
>
> julia> function g6(x, y, w)
> T = promote_type(eltype(x), eltype(y), eltype(w))
> r = zero(T)
> for i in eachindex(x, y, w)
> @inbounds r += x[i] * y[i] * w[i]
> end
> return r
> end
> g6 (generic function with 1 method)
>
> julia> @btime g6($(rand(1000)),$(rand(1000)),$(rand(1000)))
> 416.673 ns (0 allocations: 0 bytes)
> 120.4890370653176
>
> ```
ps: nobody mentioned, I think, that `LoopVectorization.@turbo` is still around:
```julia-repl
julia> x, y, z = (rand(10^3), rand(10^3), rand(10^3));
julia> function g6(x, y, w)
T = promote_type(eltype(x), eltype(y), eltype(w))
r = zero(T)
@turbo for i in eachindex(x, y, w)
r += x[i] * y[i] * w[i]
end
return r
end
g6 (generic function with 1 method)
julia> @btime g6($x,$y,$z)
65.375 ns (0 allocations: 0 bytes)
117.85528649756924
julia> function g6(x, y, w)
T = promote_type(eltype(x), eltype(y), eltype(w))
r = zero(T)
@simd for i in eachindex(x, y, w)
r += x[i] * y[i] * w[i]
end
return r
end
g6 (generic function with 1 method)
julia> @btime g6($x,$y,$z)
79.421 ns (0 allocations: 0 bytes)
117.85528649756922
```
but also the OP should probably clarify if the goal is to _actually_ perform these operations in those short (3-element) arrays.
---
### Author:  [@sbuercklin](https://discourse.julialang.org/u/sbuercklin)
#### Post date: [September 13, 2024, 1:31pm UTC](https://discourse.julialang.org/t/julian-way-to-write-this-code/119348/17 "2024-09-13T13:31:25Z")
The length 3 seems odd here. If you don’t need anything longer than 3, `StaticArrays.jl` yields a comparable (if not marginally faster) benefit to SIMD for me, and is generally preferred for vectors of this size. Is 3 a representative length for this problem, or should we benchmark with a larger vector?
```julia
julia> using Chairmarks, StaticArrays
julia> function g6(x, y, w)
T = promote_type(eltype(x), eltype(y), eltype(w))
r = zero(T)
@simd for i in eachindex(x, y, w)
r += x[i] * y[i] * w[i]
end
return r
end
g6 (generic function with 1 method)
julia> x = rand(3); y = rand(3); w = rand(3);
julia> xstatic = SVector{3}(x); ystatic = SVector{3}(y); wstatic = SVector{3}(w);
julia> g_mapreduce(x, y, w) = mapreduce(*, +, x, y, w)
g_mapreduce (generic function with 1 method)
julia> @b g6(x, y, w)
20.227 ns (1.01 allocs: 16.175 bytes)
julia> @b g6(xstatic, ystatic, wstatic)
18.295 ns (1.01 allocs: 16.135 bytes)
julia> @b g_mapreduce(x, y, w)
754.138 ns (6.28 allocs: 214.621 bytes)
julia> @b g_mapreduce(xstatic, ystatic, wstatic)
17.724 ns (1.01 allocs: 16.127 bytes)
```
---
### Author:  [@ArthurW](https://discourse.julialang.org/u/ArthurW)
#### Post date: [September 18, 2024, 12:10pm UTC](https://discourse.julialang.org/t/julian-way-to-write-this-code/119348/19 "2024-09-18T12:10:41Z")
Do you know why `splat(*)` behaves differently from `prod` here?
I tried it and while `prod` triggers some allocations, it’s still faster for larger array sizes.
---
### Author:  [@bertschi](https://discourse.julialang.org/u/bertschi)
#### Post date: [September 22, 2024, 6:32pm UTC](https://discourse.julialang.org/t/julian-way-to-write-this-code/119348/20 "2024-09-22T18:32:22Z")
Sorry, don’t know where the allocations are coming from … it seems that there are several special paths/fallbacks in `prod` though:
```julia-repl
julia> g8(x, y, z) = sum(prod, zip(x, y, z))
g8 (generic function with 1 method)
julia> @btime g8(x, y, z) (setup=(x = rand(100); y = rand(100); z = rand(100)));
201.941 ns (3 allocations: 80 bytes)
julia> @btime g8(x, y, z) (setup=(x = rand(Int, 100); y = rand(Int, 100); z = rand(Int, 100)));
124.463 ns (0 allocations: 0 bytes)
```