# What is idiomatic way of vectorizing longer operations in Julia?

**URL:** <https://discourse.julialang.org/t/what-is-idiomatic-way-of-vectorizing-longer-operations-in-julia/104117>\
**Category:** Performance\
**Tags:** performance, loopvectorization\
**Created:** [September 21, 2023, 5:43pm UTC](https://discourse.julialang.org/t/what-is-idiomatic-way-of-vectorizing-longer-operations-in-julia/104117 "2023-09-21T17:43:25Z")\
**Posts on this page:** 6\
**Page:** 1

<div class="post-metadata">

**Author:** ![Devetak](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/devetak/32/50611_2.png) [@Devetak](https://discourse.julialang.org/u/Devetak)\
**Post date:** [September 21, 2023, 5:43pm UTC](https://discourse.julialang.org/t/what-is-idiomatic-way-of-vectorizing-longer-operations-in-julia/104117/1 "2023-09-21T17:43:25Z")

</div>

Hello.

I have recently tried writing some code and need help with vectorization. I like that Julia’s loops are efficient, but sometimes I want to write less verbose code. How to do it in this case? What are the best practices? Below is an example of what I mean:

Imagine you have a vector `a` and you are interested in knowing the sum of `a[i]/(1 + a[i])`. If you want to get it done, you might be tempted to write:

```julia
function f1(a)
    return sum( a ./ (1. .+ a))
end

```

But you might also write

```julia
function f2(a)
    Val = 0.0
    for i in each index(a)
        Val += a[i] / (1. + a[i])
    end
    return value
end

```

Ignoring any typing issues, that is, `a` is an array of Float64.

Then you have

```julia
using BenchmarkTools

a = rand(9000000)

@btime f1(a)
# 22.670 ms (3 allocations: 68.66 MiB)

@btime f2(a)
# 15.136 ms (1 allocation: 16 bytes)

```

Which is a slight difference, but it quickly adds up.

How to make the `f1` as efficient memory and time-wise?

---

<div class="post-metadata">

**Author:** ![adienes](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/adienes/32/37459_2.png) [@adienes](https://discourse.julialang.org/u/adienes)\
**Post date:** [September 21, 2023, 5:54pm UTC](https://discourse.julialang.org/t/what-is-idiomatic-way-of-vectorizing-longer-operations-in-julia/104117/2 "2023-09-21T17:54:43Z")

</div>

```julia
julia> @btime mapreduce(x -> x/(1+x), +, a)
  2.042 ms (1 allocation: 16 bytes)
2.7611581397048105e6

```

`mapreduce` seems pretty good. probably there is more wizardry to be done for another factor of 2 or so if you’re willing to try a lot harder

---

<div class="post-metadata">

**Author:** ![hendri54](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/hendri54/32/9621_2.png) [@hendri54](https://discourse.julialang.org/u/hendri54)\
**Post date:** [September 21, 2023, 5:57pm UTC](https://discourse.julialang.org/t/what-is-idiomatic-way-of-vectorizing-longer-operations-in-julia/104117/3 "2023-09-21T17:57:37Z")

</div>

For your particular example:

```julia
function f3(a)
  return sum(x -> (x / (1+x)), a)
end

@btime f3($a)
2.043 ms (0 allocations: 0 bytes)

```

But your question is probably more general. `Tullio.jl` is one way of writing complex loops in a very compact format.

---

<div class="post-metadata">

**Author:** ![Devetak](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/devetak/32/50611_2.png) [@Devetak](https://discourse.julialang.org/u/Devetak)\
**Post date:** [September 21, 2023, 6:00pm UTC](https://discourse.julialang.org/t/what-is-idiomatic-way-of-vectorizing-longer-operations-in-julia/104117/4 "2023-09-21T18:00:20Z")

</div>

What??? Very cool, thank you!

Why is `mapreduce` faster than the for loop?

Are there other functions like that I should know about?

Also I feel it loses a bit of the readability that the vectorization has

---

<div class="post-metadata">

**Author:** ![hendri54](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/hendri54/32/9621_2.png) [@hendri54](https://discourse.julialang.org/u/hendri54)\
**Post date:** [September 21, 2023, 6:10pm UTC](https://discourse.julialang.org/t/what-is-idiomatic-way-of-vectorizing-longer-operations-in-julia/104117/5 "2023-09-21T18:10:34Z")

</div>

`mapreduce` likely uses optimizations along the lines of

```julia
julia> function f2(a)
           vOut = 0.0
           @simd for a1 in a
               vOut += a1 / (1.0 + a1)
           end
           return vOut
       end
f2 (generic function with 1 method)

julia> @btime f2($a)
  2.114 ms (0 allocations: 0 bytes)
2.76175963830546e6

```

Now it’s the same speed at f3

---

<div class="post-metadata">

**Author:** ![alfaromartino](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/alfaromartino/32/52986_2.png) [@alfaromartino](https://discourse.julialang.org/u/alfaromartino)\
**Post date:** [September 21, 2023, 6:22pm UTC](https://discourse.julialang.org/t/what-is-idiomatic-way-of-vectorizing-longer-operations-in-julia/104117/6 "2023-09-21T18:22:06Z")

</div>

If you want to stick to one liners and vectorized operations, I suggest using `LoopVectorization` and `LazyArrays`. When you have to sum a broadcasted operation, LazyArrays is usually the best option, it achieves almost the same speed as `mapreduce`.

Specifically:

```julia
using BenchmarkTools
a = rand(1_000)

#baseline
f1(a) = sum( a ./ (1. .+ a))
@btime f1($a)

#with a generator
f2(a) = sum(x / (1. + x) for x in a)
@btime f2($a)

# with lazyarray (notice @~ in the broadcasted part)
using LazyArrays
f3(a) = sum(@~ a ./ (1. .+ a))
@btime f3($a)

#with @turbo
using LoopVectorization
f4(a) = sum(@turbo a ./ (1. .+ a))
@btime f4($a)

# with mapreduce
f5(a) = mapreduce(x -> x/(1+x), +, a)
@btime f5($a)

```

with time

```julia
  693.377 ns (1 allocation: 7.94 KiB) # baseline
  889.583 ns (0 allocations: 0 bytes) # generator
  517.188 ns (0 allocations: 0 bytes) # Lazy Array option
  551.875 ns (1 allocation: 7.94 KiB) # LoopVectorization
  454.315 ns (0 allocations: 0 bytes) # mapreduce

```

I wrote a week ago how to keep a balance between speed and clear code, with emphasis on vectorized operations. See [here](https://alfaromartino.github.io/blog/PAGES/08_structural/)