# Avoid allocation of small arrays in median, extrema, etc

**URL:** <https://discourse.julialang.org/t/avoid-allocation-of-small-arrays-in-median-extrema-etc/58075>\
**Category:** Performance\
**Tags:** memory-allocation\
**Created:** [March 27, 2021, 12:22pm UTC](https://discourse.julialang.org/t/avoid-allocation-of-small-arrays-in-median-extrema-etc/58075 "2021-03-27T12:22:46Z")\
**Posts on this page:** 20\
**Page:** 1

<div class="post-metadata">

**Author:** ![weymouth](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/weymouth/32/15839_2.png) [@weymouth](https://discourse.julialang.org/u/weymouth)\
**Post date:** [March 27, 2021, 12:22pm UTC](https://discourse.julialang.org/t/avoid-allocation-of-small-arrays-in-median-extrema-etc/58075/1 "2021-03-27T12:22:46Z")

</div>

I have a performance critical call to set bounds on a number. This code works like magic, but I want to know why it doesn’t need to allocate

```julia
julia> @benchmark clamp(v[1],extrema((v[2],v[3],v[4]))...) setup=(v=rand(4))
BenchmarkTools.Trial: 
  memory estimate: 0 bytes
  allocs estimate: 0
  --------------
  minimum time: 8.571 ns (0.00% GC)
  median time: 9.320 ns (0.00% GC)
  mean time: 9.318 ns (0.00% GC)
  maximum time: 27.229 ns (0.00% GC)
  --------------
  samples: 10000
  evals/sample: 999

```

So again, I love this solution. It is extremely readable. But I need to create a tuple for `extrema` which seems like it should require an allocation. For example, this is mathematically equivalent

```julia
julia> @benchmark median((v[1],v[1],v[2],v[3],v[4])) setup=(v=rand(4))
BenchmarkTools.Trial: 
  memory estimate: 128 bytes
  allocs estimate: 1
  --------------
  minimum time: 94.199 ns (0.00% GC)
  median time: 103.288 ns (0.00% GC)
  mean time: 108.575 ns (1.64% GC)
  maximum time: 1.196 μs (81.55% GC)
  --------------
  samples: 10000
  evals/sample: 957

```

but requires an allocation and so is 10x slower. I tried `median(SArray(v[1],v[1],v[2],v[3],v[4]))` and `median((v[1],v[2],median((v[1],v[3],v[4]))))` but they all allocate.

Is this something special about `extrema` or something general about Julia I should understand? Or maybe @benchmarktools is lying to me?

---

<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:** [March 27, 2021, 12:59pm UTC](https://discourse.julialang.org/t/avoid-allocation-of-small-arrays-in-median-extrema-etc/58075/3 "2021-03-27T12:59:06Z")

</div>

`extrema` iterates over elements and keeps track of the min and max. So it does not allocate. See line 1711 of [https://github.com/JuliaLang/julia/blob/master/base/multidimensional.jl](https://github.com/JuliaLang/julia/blob/master/base/multidimensional.jl)

`median` needs to sort and therefore allocate a vector.

If you can pass `v` as a `Vector`, `median!` does not allocate:

```julia
julia> @benchmark median!(v) setup=(v=rand(4))
BenchmarkTools.Trial:
  memory estimate: 0 bytes
  allocs estimate: 0
  --------------
  minimum time: 15.952 ns (0.00% GC)
  median time: 16.424 ns (0.00% GC)
  mean time: 17.729 ns (0.00% GC)
  maximum time: 91.392 ns (0.00% GC)
  --------------
  samples: 10000
  evals/sample: 997

```

---

<div class="post-metadata">

**Author:** ![weymouth](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/weymouth/32/15839_2.png) [@weymouth](https://discourse.julialang.org/u/weymouth)\
**Post date:** [March 27, 2021, 2:25pm UTC](https://discourse.julialang.org/t/avoid-allocation-of-small-arrays-in-median-extrema-etc/58075/4 "2021-03-27T14:25:53Z")

</div>

Thanks, that makes sense. Unfortunately, the bounds aren’t constant, so I would need to refill the whole array every call.

In this case it looks like a custom function wins the race!

```julia
@fastmath function median(a,b,c)
    x = a-b
    if x*(b-c) ≥ 0
        return b
    elseif x*(a-c) > 0
        return c
    else
        return a
    end
end
julia> @benchmark median(v[2],v[1],median(v[3],v[1],v[4])) setup=(v=rand(4))
BenchmarkTools.Trial: 
  memory estimate: 0 bytes
  allocs estimate: 0
  --------------
  minimum time: 2.399 ns (0.00% GC)
  median time: 3.400 ns (0.00% GC)
  mean time: 3.264 ns (0.00% GC)
  maximum time: 26.301 ns (0.00% GC)
  --------------
  samples: 10000
  evals/sample: 1000

```

Plus this code is especially fast if the value **is** within the bounds (which is usually the case).

```julia
julia> @benchmark median(v[1],m,median(v[2],m,v[3])) setup=(v=rand(3);m=mean(v))
BenchmarkTools.Trial: 
  memory estimate: 0 bytes
  allocs estimate: 0
  --------------
  minimum time: 2.399 ns (0.00% GC)
  median time: 2.500 ns (0.00% GC)
  mean time: 2.769 ns (0.00% GC)
  maximum time: 75.501 ns (0.00% GC)
  --------------
  samples: 10000
  evals/sample: 1000

```

---

<div class="post-metadata">

**Author:** ![gustaphe](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/gustaphe/32/18174_2.png) [@gustaphe](https://discourse.julialang.org/u/gustaphe)\
**Post date:** [March 27, 2021, 3:16pm UTC](https://discourse.julialang.org/t/avoid-allocation-of-small-arrays-in-median-extrema-etc/58075/5 "2021-03-27T15:16:06Z")

</div>

I think you can make that function faster and more self-explanatory by replacing the multiplications by comparisons.

```julia
if a>b
    b>=c && return b
    a>c && return c
else
    b<=c && return b
    a<c && return c
end
return a

```

(Disclaimer: I’m on my phone, and I’m not entirely sure what your function is supposed to do so this might be a bit misspelled, but the idea holds)

---

<div class="post-metadata">

**Author:** ![weymouth](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/weymouth/32/15839_2.png) [@weymouth](https://discourse.julialang.org/u/weymouth)\
**Post date:** [March 27, 2021, 4:34pm UTC](https://discourse.julialang.org/t/avoid-allocation-of-small-arrays-in-median-extrema-etc/58075/6 "2021-03-27T16:34:39Z")

</div>

That’s the same speed on my system, but you’re right that it is more readable. Thanks!

---

<div class="post-metadata">

**Author:** ![gustaphe](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/gustaphe/32/18174_2.png) [@gustaphe](https://discourse.julialang.org/u/gustaphe)\
**Post date:** [March 28, 2021, 10:40am UTC](https://discourse.julialang.org/t/avoid-allocation-of-small-arrays-in-median-extrema-etc/58075/7 "2021-03-28T10:40:06Z")

</div>

If you rewrite it

```julia
function fun(a, b, c)
    (a>b)*(b>=c) && return b
    (a>b)*(a>c) && return c
    return a
end

```

it only has to perform boolean mutliplications;  
 ![medians](https://global.discourse-cdn.com/julialang/original/3X/8/4/84dace263af280ed1160796727aa79aa96645b0c.png)

“newcompmedian” is the new one. I’m guessing the three modes there correspond to shortcutting after 1, 2 or 3 comparisons.

---

<div class="post-metadata">

**Author:** ![weymouth](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/weymouth/32/15839_2.png) [@weymouth](https://discourse.julialang.org/u/weymouth)\
**Post date:** [March 28, 2021, 10:54am UTC](https://discourse.julialang.org/t/avoid-allocation-of-small-arrays-in-median-extrema-etc/58075/8 "2021-03-28T10:54:19Z")

</div>

Cool! Can you post the code for this? I would love to take a look.

---

<div class="post-metadata">

**Author:** ![gustaphe](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/gustaphe/32/18174_2.png) [@gustaphe](https://discourse.julialang.org/u/gustaphe)\
**Post date:** [March 28, 2021, 11:09am UTC](https://discourse.julialang.org/t/avoid-allocation-of-small-arrays-in-median-extrema-etc/58075/9 "2021-03-28T11:09:28Z")

</div>

Did it in the REPL (stupidly), but

```julia
using BenchmarkTools, StatsPlots, BenchmarkPlots # Requires this pull request: https://github.com/JuliaCI/BenchmarkTools.jl/pull/194/commits/230663872002eae1f9513590f68322636f2fa05d
default(fontfamily="Computer Modern")

function mulmedian(a, b, c)
   x = a-b
   if x*(b-c) >= 0
       return b
   elseif x*(a-c) > 0
       return c
   else
       return a
   end
end

function compmedian(a, b, c)
    if a>b
        b>=c && return b
        a>c && return c
    else
        b<=c && return b
        a<c && return c
    end
    return a
end
function newcompmedian(a, b, c)
    xor((a>b), (b<=c)) && return b
    xor((a>b), (a<c)) && return c
    return a
end

clampfun(v) = clamp(v[1], extrema((v[2],v[3],v[4]))...)
medianfun(v) = median((v[1],v[1],v[2],v[3],v[4]))
mulmedianfun(v) = mulmedian(v[2], v[1], mulmedian(v[3], v[1], v[4]))
compmedianfun(v) = compmedian(v[2], v[1], compmedian(v[3], v[1], v[4]))
newcompmedianfun(v) = newcompmedian(v[2], v[1], newcompmedian(v[3], v[1], v[4]))

suite = BenchmarkGroup()
suite["clamp"] = @benchmarkable clampfun(v) setup=(v=rand(4))
suite["median"] = @benchmarkable medianfun(v) setup=(v=rand(4))
suite["mulmedian"] = @benchmarkable mulmedianfun(v) setup=(v=rand(4))
suite["compmedian"] = @benchmarkable compmedianfun(v) setup=(v=rand(4))
suite["newcompmedian"] = @benchmarkable newcompmedianfun(v) setup=(v=rand(4))

results = run(suite)
plot(results; yscale=:log10)

```

---

<div class="post-metadata">

**Author:** ![Skoffer](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/skoffer/32/378_2.png) [@Skoffer](https://discourse.julialang.org/u/Skoffer)\
**Post date:** [March 28, 2021, 11:17am UTC](https://discourse.julialang.org/t/avoid-allocation-of-small-arrays-in-median-extrema-etc/58075/10 "2021-03-28T11:17:35Z")

</div>

Of note here, though, `(a>b)*(b >= c)` is not equal to `(a-b)*(b-c) >= 0`

```julia
julia> a, b, c = 1, 2, 3
(1, 2, 3)

julia> (a > b) *(b >= c)
false

julia> (a - b)*(b - c) >=0
true

```

---

<div class="post-metadata">

**Author:** ![gustaphe](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/gustaphe/32/18174_2.png) [@gustaphe](https://discourse.julialang.org/u/gustaphe)\
**Post date:** [March 28, 2021, 11:44am UTC](https://discourse.julialang.org/t/avoid-allocation-of-small-arrays-in-median-extrema-etc/58075/11 "2021-03-28T11:44:28Z")

</div>

Yeah, you’re right. I meant `xor`. Which actually makes it a bit faster. Edited.

---

<div class="post-metadata">

**Author:** ![Skoffer](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/skoffer/32/378_2.png) [@Skoffer](https://discourse.julialang.org/u/Skoffer)\
**Post date:** [March 28, 2021, 11:48am UTC](https://discourse.julialang.org/t/avoid-allocation-of-small-arrays-in-median-extrema-etc/58075/12 "2021-03-28T11:48:04Z")

</div>

I do not quite understand, why `xor` is better

```julia
julia> xor((a>b), (b>=c))
false

```

I think proper formula is

```julia
julia> (a == b) | (b == c) | ((a > b)&(b>c)) | ((a < b)&(b < c))
true

```

but it is rather slow.

---

<div class="post-metadata">

**Author:** ![gustaphe](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/gustaphe/32/18174_2.png) [@gustaphe](https://discourse.julialang.org/u/gustaphe)\
**Post date:** [March 28, 2021, 12:20pm UTC](https://discourse.julialang.org/t/avoid-allocation-of-small-arrays-in-median-extrema-etc/58075/13 "2021-03-28T12:20:34Z")

</div>

If `(a-b)*(b-c) >= 0`, then both parentheses are `>=0`, or they are both `<=0`. Depending on how sensitive the application is to the difference between `<=` and `<` (I still don’t quite understand what we are doing here, but it’s some kind of sorting/clamping thing so it shouldn’t matter), you can use `xor((a>=b), (b<c))`. Or some other combination of comparisons. Shouldn’t matter much for the timings.

---

<div class="post-metadata">

**Author:** ![weymouth](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/weymouth/32/15839_2.png) [@weymouth](https://discourse.julialang.org/u/weymouth)\
**Post date:** [March 28, 2021, 12:41pm UTC](https://discourse.julialang.org/t/avoid-allocation-of-small-arrays-in-median-extrema-etc/58075/15 "2021-03-28T12:41:02Z")

</div>

We’re ensuring that a number is within the extrema of three other numbers. The clamp function version is certainly the most direct and since it doesn’t use any custom code, it could be used to unit test any of the other versions.

I’m happy to do so but (dumb question) how do I get your PR above? Do I need to `git fetch https://github.com/JuliaCI/BenchmarkTools.jl/pull/194/commits/230663872002eae1f9513590f68322636f2fa05d`

---

<div class="post-metadata">

**Author:** ![gustaphe](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/gustaphe/32/18174_2.png) [@gustaphe](https://discourse.julialang.org/u/gustaphe)\
**Post date:** [March 28, 2021, 12:48pm UTC](https://discourse.julialang.org/t/avoid-allocation-of-small-arrays-in-median-extrema-etc/58075/16 "2021-03-28T12:48:59Z")

</div>

I’m not sure if there’s a native way to select a specific commit, but `]dev https://github.com/gustaphe/BenchmarkTools.jl/` should do the job for now because it’s on my master branch.

---

<div class="post-metadata">

**Author:** ![mcabbott](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/mcabbott/32/6603_2.png) [@mcabbott](https://discourse.julialang.org/u/mcabbott)\
**Post date:** [March 28, 2021, 1:13pm UTC](https://discourse.julialang.org/t/avoid-allocation-of-small-arrays-in-median-extrema-etc/58075/17 "2021-03-28T13:13:21Z")

</div>

Perhaps a bit of a tangent, but here’s an attempt at at a non-allocating `median` for tuples, by simply listing all the options. This will end up making the same comparison many times, and only works for short tuples, surely there is a better algorithm:

```julia
using Combinatorics, Statistics
function _medex(n)
    out = quote end
    for p in permutations(1:n)
        z = zip(Iterators.repeated(:(<=)), [:(x[$i]) for i in p])
        c = Expr(:comparison, Iterators.drop(Iterators.flatten(z),1)...)
        v = isodd(n) ? :(x[$(p[1+n÷2])]) : :((x[$(p[n÷2])] + x[$(p[1+n÷2])])/2)
        push!(out.args, :($c && return $v))
    end
    out
end
_medex(3) # to see what this generates
@generated unrollmedian(x::Tuple) = _medex(length(x.parameters))
unrollmedian(xs::Number...) = unrollmedian(xs)
for n in 2:6
    x = Tuple(rand(Int8,n))
    @show n
    @btime median($(Ref(x))[])
    @btime unrollmedian($(Ref(x))[]) # at n=7 I have to kill julia 
end

using Random, Statistics, BenchmarkTools
Statistics.median(xs::Number...) = median(xs) # to match @gustaphe's functions
Random.seed!(1);
xs = Tuple(rand(Int8, 3)) # (64, 82, -103)
for f in [median, unrollmedian, fastmath_median, fun, mulmedian, compmedian, newcompmedian]    
    @show f f(xs...)
    @btime $f($(Ref(xs))[]...)
end

```

This comparison prints the following, note that several functions above give the wrong answer, I think they assume that all numbers are positive. Not sure how meaningful these tiny times are:

```julia
f = Statistics.median
f(xs...) = -82.0
  26.439 ns (1 allocation: 96 bytes)
f = unrollmedian
f(xs...) = -82
  1.208 ns (0 allocations: 0 bytes)
f = fastmath_median
f(xs...) = 103
  0.916 ns (0 allocations: 0 bytes)
f = fun
f(xs...) = -82
  1.292 ns (0 allocations: 0 bytes)
f = mulmedian
f(xs...) = 103
  0.916 ns (0 allocations: 0 bytes)
f = compmedian
f(xs...) = -82
  1.166 ns (0 allocations: 0 bytes)
f = newcompmedian
f(xs...) = -82
  1.333 ns (0 allocations: 0 bytes)

```

---

<div class="post-metadata">

**Author:** ![weymouth](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/weymouth/32/15839_2.png) [@weymouth](https://discourse.julialang.org/u/weymouth)\
**Post date:** [March 28, 2021, 3:31pm UTC](https://discourse.julialang.org/t/avoid-allocation-of-small-arrays-in-median-extrema-etc/58075/18 "2021-03-28T15:31:54Z")

</div>

I can confirm that `newcompmedian` doesn’t work. If `a<b<c` then `xor( a>b, b>=c ) = xor(false,false)=false` which means it misses `b` as the median.

This one seems to work

```julia
function eqcompmedian(a, b, c)
    (a>b) == (b>=c) && return b
    (a>b) == (a>c) && return c
    return a
end

```

I think this is exactly the same logical chain as `mulmedian`. When doing the timing, I get

```julia
julia> run(suite)
5-element BenchmarkTools.BenchmarkGroup:
  tags: []
  "median" => Trial(144.000 ns)
  "compmedian" => Trial(43.000 ns)
  "eqcompmedian" => Trial(45.000 ns)
  "clamp" => Trial(63.000 ns) # this one seems quite variable??
  "mulmedian" => Trial(47.000 ns)

```

So `compmedian` is the best by a hair.

---

<div class="post-metadata">

**Author:** ![gustaphe](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/gustaphe/32/18174_2.png) [@gustaphe](https://discourse.julialang.org/u/gustaphe)\
**Post date:** [March 28, 2021, 3:49pm UTC](https://discourse.julialang.org/t/avoid-allocation-of-small-arrays-in-median-extrema-etc/58075/19 "2021-03-28T15:49:27Z")

</div>

Yeah, it would have to be opposite comparisons. `(a>b) == (b>=c) <==> ~xor((a>b), (b<c))`. Don’t know if there’s any performance implication.

---

<div class="post-metadata">

**Author:** ![weymouth](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/weymouth/32/15839_2.png) [@weymouth](https://discourse.julialang.org/u/weymouth)\
**Post date:** [March 28, 2021, 3:56pm UTC](https://discourse.julialang.org/t/avoid-allocation-of-small-arrays-in-median-extrema-etc/58075/20 "2021-03-28T15:56:37Z")

</div>

That works! On my machine all of the double median version are essentially the same. Clamp is a little slower (actually, still pretty much the same) and the single vectorized median is much worse due to the allocation.

---

<div class="post-metadata">

**Author:** ![rafael.guerra](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/rafael.guerra/32/216610_2.png) [@rafael.guerra](https://discourse.julialang.org/u/rafael.guerra)\
**Post date:** [March 28, 2021, 4:38pm UTC](https://discourse.julialang.org/t/avoid-allocation-of-small-arrays-in-median-extrema-etc/58075/21 "2021-03-28T16:38:54Z")

</div>

@mcabbot, in case it might interest, for larger lists of numbers, a median finding algorithm, seemingly faster than sort-based methods, is presented in a simple way [here](https://rcoh.me/posts/linear-time-median-finding/) and the author points also to a more advanced [reference](http://erdani.com/research/sea2017.pdf) on this topic by Andrei Alexandrescu.

---

<div class="post-metadata">

**Author:** ![Skoffer](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/skoffer/32/378_2.png) [@Skoffer](https://discourse.julialang.org/u/Skoffer)\
**Post date:** [March 28, 2021, 4:44pm UTC](https://discourse.julialang.org/t/avoid-allocation-of-small-arrays-in-median-extrema-etc/58075/22 "2021-03-28T16:44:46Z")

</div>

There is also approximate median algorithm, which is promised to be effective: [An Efficient Algorithm for the Approximate Median Selection Problem | SpringerLink](https://link.springer.com/chapter/10.1007%2F3-540-46521-9_19)

[Next page](https://discourse.julialang.org/t/avoid-allocation-of-small-arrays-in-median-extrema-etc/58075.md?page=2)