Why is dot so slow for small vectors?

General Usage Performance
1

I used LinearAlgebra.dot for small to medium sized vectors but discovered recently that it is much slower for short vectors than an simple for loop. The following code

using Printf
using LinearAlgebra
using BenchmarkTools

function mydot(x, y)
    s = zero(eltype(x))
    axes(x) === axes(y) || error("axes must match")
    @inbounds @simd for i in eachindex(x, y)
        s += x[i]*y[i]
    end
    return s
end

lens = [3, 11, 21, 51, 501, 2501]
mtdot = similar(lens, Float64)
mtmydot = similar(lens, Float64)

for i in eachindex(lens)
    global mt
    l = lens[i]
    x = randn(l)
    y = randn(l)
    #x = @view(randn(l+2)[begin+1:end-1])
    #y = @view(randn(l+2)[begin+1:end-1])
    #x = SVector{length(x)}(x)
    #y = SVector{length(x)}(y)
    bm = @benchmark dot($x, $y)
    mtdot[i] = minimum(bm.times)
    bm = @benchmark mydot($x, $y)
    mtmydot[i] = minimum(bm.times)
end

println("length  dot/ns  mydot/ns")
for i in eachindex(lens)
    @printf("%4d   %1.2f ns      %1.2f ns\n", lens[i], mtdot[i]/lens[i], mtmydot[i]/lens[i])
end

gives this output for the times needed per bin:

length dot/ns mydot/ns
3 2.00 ns 0.70 ns
11 0.76 ns 0.32 ns
21 0.34 ns 0.27 ns
51 0.17 ns 0.15 ns
501 0.06 ns 0.06 ns
2501 0.05 ns 0.05 ns

So for small, β€œgeometrical” vectors dot is nearly three times slower than the loop. It needs more than 50 bins until they run equally fast. The picture is the same if I use views (see code in comments). If I use static vectors both are equally fast. So is it recommended to use dot only for long vectors? Actually it would be nice to have a function that is equally fast for all sizes.

My versionfo() output is:

Julia Version 1.12.7
Commit 6d172b025e (2026-08-15 08:05 UTC)
Build Info:
  Official https://julialang.org release
Platform Info:
  OS: Windows (x86_64-w64-mingw32)
  CPU: 16 Γ— AMD Ryzen AI 7 PRO 350 w/ Radeon 860M
  WORD_SIZE: 64
  LLVM: libLLVM-18.1.7 (ORCJIT, znver5)
  GC: Built with stock GC
Threads: 7 default, 1 interactive, 7 GC (on 16 virtual cores)
Environment:
  JULIA_NUM_THREADS = 7,1
  JULIA_EDITOR = code
  JULIA_VSCODE_REPL = 1
2

LinearAlgebra.dot uses BLAS. This is a library of foreign functions. This means that there is more call overhead than something native. The overhead is small in absolute terms but can be significant in relative terms.

I prefer to look at the absolute timings:

length  dot  mydot
   3   4.70 ns      2.42 ns
  11   7.46 ns      4.77 ns
  21   6.83 ns      4.05 ns
  51   7.47 ns      4.68 ns
 501   27.33 ns      25.65 ns
2501   132.83 ns      129.90 ns

So the clearer picture is that dot is 2-3ns slower than mydot at all sizes.

3

Looking at this as the time/length is masking what is actually going on:

julia> for i in eachindex(lens)
           @printf("%4d   %1.2f ns      %1.2f ns\n", lens[i], mtdot[i], mtmydot[i])
       end
   3   4.69 ns      2.34 ns
  11   6.46 ns      4.19 ns
  21   6.27 ns      7.13 ns
  51   7.19 ns      8.82 ns
 501   34.87 ns      38.04 ns
2501   262.67 ns      261.39 ns

There is a small, constant overhead in the regular dot function. If you look at how the function is implemented, it’s because it dispatches to a ccall, which basically means there’s a layer of indirection here.

This could be fixed by making a julia native version of this function. I’m actually rather surprised we use BLAS for this, I thought usually we only used BLAS for things like matmul.

Edit: lol @mikmoore beat me to it by a few seconds. Sorry for the double post

4

Small followup, but perhaps this makes more sense @mwlidar when you account for the fact that BLAS will start multithreading these calls beyond a certain length cutoff, so a performance gap opens up again:

julia> let len = 100_000
           u, v = randn(len), randn(len)
           @info "len = $len" dot=@benchmark(dot($u, $v)) mydot=@benchmark(mydot($u, $v))
       end
β”Œ Info: len = 100000
β”‚   dot =
β”‚    BenchmarkTools.Trial: 10000 samples with 8 evaluations per sample.
β”‚     Range (min … max):  3.446 ΞΌs … 54.877 ΞΌs  β”Š GC (min … max): 0.00% … 0.00%
β”‚     Time  (median):     3.632 ΞΌs              β”Š GC (median):    0.00%
β”‚     Time  (mean Β± Οƒ):   3.732 ΞΌs Β±  1.265 ΞΌs  β”Š GC (mean Β± Οƒ):  0.00% Β± 0.00%
β”‚    
β”‚      β–ƒβ–„β–…β–ˆβ–‡β–‚β–‚                                                    ▁
β”‚      β–ˆβ–ˆβ–ˆβ–ˆβ–ˆβ–ˆβ–ˆβ–‡β–‡β–†β–…β–‡β–†β–ˆβ–†β–†β–…β–†β–…β–…β–…β–…β–†β–†β–…β–…β–…β–…β–ƒβ–„β–ƒβ–…β–„β–„β–„β–„β–„β–ƒβ–„β–…β–β–„β–β–„β–β–ƒβ–ƒβ–„β–β–ƒβ–β–ƒβ–β–ƒβ–„β–„β–„β–… β–ˆ
β”‚      3.45 ΞΌs      Histogram: log(frequency) by time     6.33 ΞΌs <
β”‚    
β”‚     Memory estimate: 0 bytes, allocs estimate: 0.
β”‚   mydot =
β”‚    BenchmarkTools.Trial: 10000 samples with 1 evaluation per sample.
β”‚     Range (min … max):  12.904 ΞΌs … 68.269 ΞΌs  β”Š GC (min … max): 0.00% … 0.00%
β”‚     Time  (median):     14.698 ΞΌs              β”Š GC (median):    0.00%
β”‚     Time  (mean Β± Οƒ):   15.240 ΞΌs Β±  2.256 ΞΌs  β”Š GC (mean Β± Οƒ):  0.00% Β± 0.00%
β”‚    
β”‚      ▂▁   β–„β–ˆβ–ˆβ–‡β–†β–„β–ƒβ–‚β–‚β–β–                                            β–‚
β”‚      β–ˆβ–ˆβ–‡β–‡β–†β–ˆβ–ˆβ–ˆβ–ˆβ–ˆβ–ˆβ–ˆβ–ˆβ–ˆβ–ˆβ–ˆβ–ˆβ–ˆβ–ˆβ–‡β–†β–‡β–‡β–†β–ˆβ–ˆβ–‡β–‡β–†β–†β–†β–†β–‡β–‡β–†β–‡β–†β–†β–†β–†β–†β–†β–…β–…β–†β–†β–…β–…β–†β–…β–…β–†β–†β–„β–„β–…β–…β–„β–„ β–ˆ
β”‚      12.9 ΞΌs      Histogram: log(frequency) by time      26.5 ΞΌs <
β”‚    
β””     Memory estimate: 0 bytes, allocs estimate: 0.

I suspect we can handle this alright though with native julia multithreading if someone was interested in taking a stab at it.

5

If you are dealing with β€œgeometrical” vectors, typically 2- and 3-component vectors, where the dimensionality is effectively known ahead of time (statically), then you are usually much better off using StaticArrays.jl.

Not only will operations like dot be much faster than loop-based code (the loop is completely unrolled to 3 multiplications and 2 additions, and no branches or indexing), but many other operations will be faster, and you don’t need to worry about allocating temporary vectors so the code becomes much cleaner (e.g. you can just do things like vec1 = 2vec2 + vec3 * (vec4'vec5) in the most natural way, without messing around with pre-allocation and in-place operations like you would for larger arrays).

See also When Static Arrays may be useful.