# High performance vector/matrix/tensor linear algebra operations

**URL:** <https://discourse.julialang.org/t/high-performance-vector-matrix-tensor-linear-algebra-operations/88756>\
**Category:** Performance\
**Tags:** question, performance, linearalgebra\
**Created:** [October 14, 2022, 9:30pm UTC](https://discourse.julialang.org/t/high-performance-vector-matrix-tensor-linear-algebra-operations/88756 "2022-10-14T21:30:01Z")\
**Posts on this page:** 10\
**Page:** 1

<div class="post-metadata">

**Author:** ![FireCrumb](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/firecrumb/32/35370_2.png) [@FireCrumb](https://discourse.julialang.org/u/FireCrumb)\
**Post date:** [October 14, 2022, 9:30pm UTC](https://discourse.julialang.org/t/high-performance-vector-matrix-tensor-linear-algebra-operations/88756/1 "2022-10-14T21:30:01Z")

</div>

Hi, I have vectors (say 20 of them) of the same size (say 256 Float64).  
Now I want to add two of them and put the result into the third one, with the index of which one given as input.

_TL;DR:_ what is the fastest way to perform linear algebra operations on such vectors on CPU?

Sample code for adding two vectors:

```julia
vectors1 = Vector{Vector{Float64}}([zeros(Float64, 256) for _ = 1:20]);
@inline function add_vec(mem, out, in1, in2)
    @inbounds mem[out] .= mem[in1] .+ mem[in2]
end

```

When I check the performance, it’s unsatisfactory – because I’m doing this operation _a lot_ so it adds up.

```julia
using BenchmarkTools
@benchmark add_vec(vectors1, 3, 1, 2)

```

```julia
BenchmarkTools.Trial: 10000 samples with 994 evaluations.
 Range (min … max): 30.430 ns … 95.037 ns ┊ GC (min … max): 0.00% … 0.00%
 Time (median): 31.739 ns ┊ GC (median): 0.00%
 Time (mean ± σ): 31.967 ns ± 1.484 ns ┊ GC (mean ± σ): 0.00% ± 0.00%

        ▁▂▄▅▇███▇▅▄▄▃▂ ▂
  ▂▂▂▇▇▇███████████████▇▆▅▄▂▄▅▄▆▄▅▄▄▅▄▂▂▄▅▆▆▇▇█▇▇▇▇▇▆▆▆▅▄▅▅▄▅ █
  30.4 ns Histogram: log(frequency) by time 36.6 ns <

 Memory estimate: 0 bytes, allocs estimate: 0.

```

Although obviously can always be worse, this is what I started with:

```julia
function add_vec_unoptimized(mem, out, in1, in2)
    mem[out] = mem[in1] + mem[in2]
end

@benchmark add_vec_unoptimized(vectors1, 3, 1, 2)

```

```julia
BenchmarkTools.Trial: 10000 samples with 681 evaluations.
 Range (min … max): 164.420 ns … 4.327 μs ┊ GC (min … max): 0.00% … 86.74%
 Time (median): 214.507 ns ┊ GC (median): 0.00%
 Time (mean ± σ): 304.479 ns ± 341.656 ns ┊ GC (mean ± σ): 17.15% ± 13.98%

  ▄█▅▃▂▁ ▃▂ ▁
  ███████▇▇▆▆▆▇████▇▆▆▅▁▃▃▁▁▁▁▁▁▁▁▁▄▃▁▃▃▄▅▄▅▆▆▆▆▆▆▆▆▆▄▅▅▅▆▆▅▇▆▇ █
  164 ns Histogram: log(frequency) by time 2.15 μs <

 Memory estimate: 2.12 KiB, allocs estimate: 1.

```

When looking at the underlying code using `code_warntype` and `code_native`, I see it’s doing lots of things internally – checking dimensions, `materialize`, `broadcast` (in the `.+` version), `axes`, `shape`, `convert`, `alias` etc.

Tried also using [StaticArrays.jl](https://github.com/JuliaArrays/StaticArrays.jl) to hint the compiler about the vector size, but it was slower.

What is the correct way to do this?  
I’m asking for this one operation for sake of simplicity and analysis, but I’m looking for a solution for all linear algebra operations.  
Preferably single-threaded CPU.

Thanks!

---

<div class="post-metadata">

**Author:** ![Elrod](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/elrod/32/22461_2.png) [@Elrod](https://discourse.julialang.org/u/Elrod)\
**Post date:** [October 14, 2022, 10:21pm UTC](https://discourse.julialang.org/t/high-performance-vector-matrix-tensor-linear-algebra-operations/88756/2 "2022-10-14T22:21:11Z")

</div>

Loops are going to be reliable. `Fastbroadcast.jl` should also shave off a little overhead and match the loop.

The loop would also make it easy to pass in static size information.

---

<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:** [October 14, 2022, 10:54pm UTC](https://discourse.julialang.org/t/high-performance-vector-matrix-tensor-linear-algebra-operations/88756/3 "2022-10-14T22:54:23Z")

</div>

> [@FireCrumb](#):
>
> When I check the performance, it’s unsatisfactory

What sort of performance are you expecting? 256 additions in 30ns seems really fast to me. Maybe there’s still more to be had, but probably not a lot.

Multi-threading a single vector addition won’t help for such short time spans, but you could parallelize different vector operations.

---

<div class="post-metadata">

**Author:** ![matthias314](https://avatars.discourse-cdn.com/v4/letter/m/a88e4f/32.png) [@matthias314](https://discourse.julialang.org/u/matthias314)\
**Post date:** [October 15, 2022, 1:31am UTC](https://discourse.julialang.org/t/high-performance-vector-matrix-tensor-linear-algebra-operations/88756/4 "2022-10-15T01:31:36Z")

</div>

Just for the record: If you use `$` in your benchmarks, then the time is reduced by 1/3:

```julia
julia> @btime add_vec(vectors1, 3, 1, 2);
  46.584 ns (0 allocations: 0 bytes)

julia> @btime add_vec($vectors1, 3, 1, 2);
  31.934 ns (0 allocations: 0 bytes)

```

---

<div class="post-metadata">

**Author:** ![Elrod](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/elrod/32/22461_2.png) [@Elrod](https://discourse.julialang.org/u/Elrod)\
**Post date:** [October 15, 2022, 2:23am UTC](https://discourse.julialang.org/t/high-performance-vector-matrix-tensor-linear-algebra-operations/88756/5 "2022-10-15T02:23:50Z")

</div>

FWIW

```julia
julia> vectors1 = Vector{Vector{Float64}}([zeros(Float64, 256) for _ = 1:20]);

julia> @inline function add_vec_loop(mem, out, in1, in2)
           @inbounds begin
             mo = mem[out]; mi1 = mem[in1]; mi2 = mem[in2];
           end
           for i = eachindex(mo)
               @inbounds mo[i] = mi1[i] + mi2[i]
           end
       end
add_vec_loop (generic function with 1 method)

julia> @benchmark add_vec_loop($vectors1, 3, 1, 2)
BenchmarkTools.Trial: 10000 samples with 999 evaluations.
 Range (min … max): 12.189 ns … 32.116 ns ┊ GC (min … max): 0.00% … 0.00%
 Time (median): 13.191 ns ┊ GC (median): 0.00%
 Time (mean ± σ): 13.088 ns ± 0.500 ns ┊ GC (mean ± σ): 0.00% ± 0.00%

                            ▇ █▁
  ▂▂▅▆▃▃▂▂▂▂▂▂▂▂▁▂▃▃▃▂▂▃▄▂▂▄█▃▄██▂▂▁▂▂▂▂▂▂▂▂▂▂▁▂▂▂▂▂▂▂▂▂▂▂▂▂▂ ▃
  12.2 ns Histogram: frequency by time 14.4 ns <

 Memory estimate: 0 bytes, allocs estimate: 0.

julia> @inline function add_vec_fb(mem, out, in1, in2)
           @inbounds begin
             mo = mem[out]; mi1 = mem[in1]; mi2 = mem[in2];
           end
           @.. mo = mi1 + mi2
       end
add_vec_fb (generic function with 1 method)

julia> @benchmark add_vec_fb($vectors1, 3, 1, 2)
BenchmarkTools.Trial: 10000 samples with 999 evaluations.
 Range (min … max): 12.191 ns … 39.547 ns ┊ GC (min … max): 0.00% … 0.00%
 Time (median): 12.462 ns ┊ GC (median): 0.00%
 Time (mean ± σ): 12.492 ns ± 0.364 ns ┊ GC (mean ± σ): 0.00% ± 0.00%

         ▂▂ ▄█ ▁
  ▂▃▃▃▄▅▅██▆██▇▃▅█▆▂▂▁▂▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂ ▃
  12.2 ns Histogram: frequency by time 13.7 ns <

 Memory estimate: 0 bytes, allocs estimate: 0.

julia> @inline function add_vec2(mem, out, in1, in2)
           @inbounds begin
             mo = mem[out]; mi1 = mem[in1]; mi2 = mem[in2];
           end
           @. mo = mi1 + mi2
       end
add_vec2 (generic function with 1 method)

julia> @benchmark add_vec2($vectors1, 3, 1, 2)
BenchmarkTools.Trial: 10000 samples with 998 evaluations.
 Range (min … max): 14.466 ns … 33.422 ns ┊ GC (min … max): 0.00% … 0.00%
 Time (median): 15.247 ns ┊ GC (median): 0.00%
 Time (mean ± σ): 15.318 ns ± 0.419 ns ┊ GC (mean ± σ): 0.00% ± 0.00%

   ▁ ▁▃▆██▇▇▅▂ ▁ ▂
  ██▆▁▃▁▁▁▁▇██████████▅▁▁▁▁▃▃▁▃▄▃▆▇███▆▅▄▁▄▆▆██▅▆▅▃▃▃▅▆▅▄▅▆▅▅ █
  14.5 ns Histogram: log(frequency) by time 17.4 ns <

 Memory estimate: 0 bytes, allocs estimate: 0.

```

I get a decent speedup from 15 → 13ns when switching from `@.` to `FastBroadcast.@..` or a loop.

---

<div class="post-metadata">

**Author:** ![FireCrumb](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/firecrumb/32/35370_2.png) [@FireCrumb](https://discourse.julialang.org/u/FireCrumb)\
**Post date:** [January 20, 2023, 8:29am UTC](https://discourse.julialang.org/t/high-performance-vector-matrix-tensor-linear-algebra-operations/88756/6 "2023-01-20T08:29:26Z")

</div>

Just back to update on the solution

Using [LoopVectorization.jl](https://docs.juliahub.com/LoopVectorization/)

Example:  
[https://docs.juliahub.com/LoopVectorization/4TogI/0.12.147/examples/matrix\_multiplication/](https://docs.juliahub.com/LoopVectorization/4TogI/0.12.147/examples/matrix_multiplication/)

---

<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:** [January 20, 2023, 8:58am UTC](https://discourse.julialang.org/t/high-performance-vector-matrix-tensor-linear-algebra-operations/88756/7 "2023-01-20T08:58:02Z")

</div>

> [@FireCrumb](#):
>
> `Vector{Vector{Float64}}([zeros(Float64, 256) for _ = 1:20]);`

One thing I forgot to mention. In the above code you are taking a `Vector{Vector{Float64}}` and passing it into a `Vector{Vector{Float64}}` constructor. That is redundant, since `[zeros(Float64, 256) for _ = 1:20]` is _already_ a `Vector{Vector{Float64}}`.

So you can just write

```julia
vectors1 = [zeros(Float64, 256) for _ = 1:20]

```

This is considered idiomatic.

---

<div class="post-metadata">

**Author:** ![FireCrumb](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/firecrumb/32/35370_2.png) [@FireCrumb](https://discourse.julialang.org/u/FireCrumb)\
**Post date:** [January 20, 2023, 9:01am UTC](https://discourse.julialang.org/t/high-performance-vector-matrix-tensor-linear-algebra-operations/88756/8 "2023-01-20T09:01:40Z")

</div>

Thanks for the reminder about the data structure.

So also update about that – after investigation, much better performance is by using one big two dimensional [MArray](https://juliaarrays.github.io/StaticArrays.jl/stable/pages/api/#Mutable-arrays:-MVector,-MMatrix-and-MArray) and subsetting it, rather than using a _vector of vectors_.

---

<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:** [January 20, 2023, 9:06am UTC](https://discourse.julialang.org/t/high-performance-vector-matrix-tensor-linear-algebra-operations/88756/9 "2023-01-20T09:06:15Z")

</div>

That is a bit surprising. `MArray`s are specifically intended for small arrays, and their performance degrades as their sizes grow – normally much beyond 100 elements is not worthwhile. Also, compilation time gets gradually worse.

---

<div class="post-metadata">

**Author:** ![FireCrumb](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/firecrumb/32/35370_2.png) [@FireCrumb](https://discourse.julialang.org/u/FireCrumb)\
**Post date:** [January 20, 2023, 9:12am UTC](https://discourse.julialang.org/t/high-performance-vector-matrix-tensor-linear-algebra-operations/88756/10 "2023-01-20T09:12:59Z")

</div>

It’s true that after some size the performance degrades, but it seems to be way beyond 100 elements.

Probably worth checking on case by case basis.