# Weighted histogram ~2x as slow in Julia vs. Python

**URL:** <https://discourse.julialang.org/t/weighted-histogram-2x-as-slow-in-julia-vs-python/74953>\
**Category:** Performance\
**Tags:** question, statistics, python\
**Created:** [January 20, 2022, 6:55pm UTC](https://discourse.julialang.org/t/weighted-histogram-2x-as-slow-in-julia-vs-python/74953 "2022-01-20T18:55:17Z")\
**Posts on this page:** 14\
**Page:** 1

<div class="post-metadata">

**Author:** ![kirklong](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/kirklong/32/19433_2.png) [@kirklong](https://discourse.julialang.org/u/kirklong)\
**Post date:** [January 20, 2022, 6:55pm UTC](https://discourse.julialang.org/t/weighted-histogram-2x-as-slow-in-julia-vs-python/74953/1 "2022-01-20T18:55:17Z")

</div>

I’m working on a project where I need to compute the weighted histogram of something many times. Basically I need the sum of each of the bins of the histogram, and I’ve been using PyCall to call the SciPy “binned statistic” [method](https://docs.scipy.org/doc/scipy/reference/generated/scipy.stats.binned_statistic.html) until this point as there wasn’t an obvious equivalent I could find in Julia.

While this is fine, it’s always annoyed me, and I recently decided I could just write it myself in native Julia and stop being annoyed by it. So I took this naive stab at it, which works, but is twice as slow as calling the Python version (~30 vs. ~60 ms):

```julia
function histSum(x::Array{Float64,},y::Array{Float64,},bins::Int=100) #61 ms
    x = vec(x); y = vec(y)
    binMax = maximum(x); binMin = minimum(x)
    binEdges = range(binMin,stop=binMax,length=bins+1)
    result = zeros(length(binEdges)-1)
    for i=1:bins-1
        result[i] = sum(y[(x .>= binEdges[i]) .& (x .< binEdges[i+1])])
    end
    result[end] = sum(y[(x .>= binEdges[end-1]) .& (x .<= binEdges[end])])
    return binEdges, result
end

```

I posted the problem on the Julia slack and there were some great suggestions from Dave MacMahon, including this much more elegant solution that uses the built-in StatsBase histogram function:

```julia
using StatsBase
function histSum(x::Array,y::Array,bins::Int=100)
    h = fit(Histogram, vec(x), aweights(y), nbins=bins)
    h.edges, h.weights
end

```

But this way is nearly identical in speed to my first naive way, so still twice as slow as using PyCall. The relevant Python source code for binned\_statistic is [here](https://github.com/scipy/scipy/blob/47bb6febaa10658c72962b9615d5d5aa2513fa3a/scipy/stats/_binned_statistic.py#L361), which in the summation case relies on the [numpy bincount function](https://github.com/numpy/numpy/blob/8dbd507fb6c854b362c26a0dd056cd04c9c10f25/numpy/core/multiarray.py#L883) which is calling some C code I can’t find.

One thought is that maybe the numpy version is automatically multi-threaded behind the scenes and this is responsible for the discrepancy? If I place Threads.@threads in front of the for loop in my first version of the code I can beat the Python version by a few ms. Any suggestions much appreciated!

---

<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, 2022, 7:10pm UTC](https://discourse.julialang.org/t/weighted-histogram-2x-as-slow-in-julia-vs-python/74953/2 "2022-01-20T19:10:36Z")

</div>

> [@kirklong](#):
>
> ```julia
> result[i] = sum(y[(x .>= binEdges[i]) .& (x .< binEdges[i+1])])
> 
> ```

When you have performance issues, looking for unnecessary allocations is often a good idea. In this line you are creating one BitArray and one regular array for each bin.

If you can write this out using loop, using only scalar operations, there should be a decent speedup, I think. Maybe for each element search for the right bin using one of the sorted search functions.

---

<div class="post-metadata">

**Author:** ![kirklong](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/kirklong/32/19433_2.png) [@kirklong](https://discourse.julialang.org/u/kirklong)\
**Post date:** [January 20, 2022, 7:18pm UTC](https://discourse.julialang.org/t/weighted-histogram-2x-as-slow-in-julia-vs-python/74953/3 "2022-01-20T19:18:14Z")

</div>

This was a suggestion I got on the Slack as well (to try sorting instead) and we came up with this implementation:

```julia
function histSumSort(x::Array{Float64,},y::Array{Float64,},bins::Int=100) #87 ms
    x = vec(x); y = vec(y)
    p = sortperm(x)
    #sortedX = x[p]; sortedY = y[p]
    bins = range(minimum(x),stop=maximum(x),length=bins+1)
    I = 1
    result = zeros(length(bins)-1)
    for i in p
        while x[i] > bins[I+1]
            I += 1
        end
        result[I] += y[i]
    end
    return bins,result
end

```

Unfortunately this way is actually _slower_ (by about 1/3) than my initial naive guess…but maybe you can think of a better way to write this?

---

<div class="post-metadata">

**Author:** ![GunnarFarneback](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/gunnarfarneback/32/1827_2.png) [@GunnarFarneback](https://discourse.julialang.org/u/GunnarFarneback)\
**Post date:** [January 20, 2022, 7:19pm UTC](https://discourse.julialang.org/t/weighted-histogram-2x-as-slow-in-julia-vs-python/74953/4 "2022-01-20T19:19:59Z")

</div>

Try this:

```julia
function histSum(x::Array{Float64,}, y::Array{Float64,}, bins::Int=100)
    binMin, binMax = extrema(x)
    result = zeros(bins)
    α = bins / (binMax - binMin)
    for (x, y) in zip(x, y)
        i = min(bins, 1 + floor(Int, α * (x - binMin)))
        result[i] += y
    end
    return range(binMin, stop=binMax, length=bins+1), result
end

```

---

<div class="post-metadata">

**Author:** ![kirklong](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/kirklong/32/19433_2.png) [@kirklong](https://discourse.julialang.org/u/kirklong)\
**Post date:** [January 20, 2022, 7:24pm UTC](https://discourse.julialang.org/t/weighted-histogram-2x-as-slow-in-julia-vs-python/74953/5 "2022-01-20T19:24:14Z")

</div>

This does it! What about this makes it so much faster? I notice if I place the @threads macro in front of this for loop I don’t get any speedup, and the time it takes is basically identical to my first way with the @threads macro, so I’m wondering if somehow this formulation is by default multi-threaded or something?

---

<div class="post-metadata">

**Author:** ![GunnarFarneback](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/gunnarfarneback/32/1827_2.png) [@GunnarFarneback](https://discourse.julialang.org/u/GunnarFarneback)\
**Post date:** [January 20, 2022, 7:28pm UTC](https://discourse.julialang.org/t/weighted-histogram-2x-as-slow-in-julia-vs-python/74953/6 "2022-01-20T19:28:14Z")

</div>

Since the bins are uniform you can directly compute into which bin each element go, which is much faster than searching for the bin.

---

<div class="post-metadata">

**Author:** ![kirklong](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/kirklong/32/19433_2.png) [@kirklong](https://discourse.julialang.org/u/kirklong)\
**Post date:** [January 20, 2022, 7:29pm UTC](https://discourse.julialang.org/t/weighted-histogram-2x-as-slow-in-julia-vs-python/74953/7 "2022-01-20T19:29:07Z")

</div>

Thanks so much!

---

<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, 2022, 7:35pm UTC](https://discourse.julialang.org/t/weighted-histogram-2x-as-slow-in-julia-vs-python/74953/8 "2022-01-20T19:35:02Z")

</div>

> [@kirklong](#):
>
> This was a suggestion I got on the Slack as well (to try sorting instead)

I didn’t actually mean that you should sort, but exploit that the _bins_ are sorted. @GunnarFarneback’s solution is an even faster version of that.

---

<div class="post-metadata">

**Author:** ![goerch](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/goerch/32/29122_2.png) [@goerch](https://discourse.julialang.org/u/goerch)\
**Post date:** [January 20, 2022, 7:47pm UTC](https://discourse.julialang.org/t/weighted-histogram-2x-as-slow-in-julia-vs-python/74953/9 "2022-01-20T19:47:00Z")

</div>

> [@GunnarFarneback](#):
>
> ```julia
> function histSum(x::Array{Float64,}, y::Array{Float64,}, bins::Int=100)
> binMin, binMax = extrema(x)
> result = zeros(bins)
> α = bins / (binMax - binMin)
> for (x, y) in zip(x, y)
> i = min(bins, 1 + floor(Int, α * (x - binMin)))
> result[i] += y
> end
> return range(binMin, stop=binMax, length=bins+1), result
> end
> 
> ```

Interestingly, changing `extrema` back to `minimum` and `maximum` seems quite a bit faster on my system (Julia 1.7.1)? Otherwise very nice!

---

<div class="post-metadata">

**Author:** ![jling](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/jling/32/212909_2.png) [@jling](https://discourse.julialang.org/u/jling)\
**Post date:** [January 20, 2022, 8:13pm UTC](https://discourse.julialang.org/t/weighted-histogram-2x-as-slow-in-julia-vs-python/74953/10 "2022-01-20T20:13:49Z")

</div>

[https://github.com/Moelf/FHist.jl](https://github.com/Moelf/FHist.jl)

give this a try if you want

---

<div class="post-metadata">

**Author:** ![albheim](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/albheim/32/34660_2.png) [@albheim](https://discourse.julialang.org/u/albheim)\
**Post date:** [January 20, 2022, 9:19pm UTC](https://discourse.julialang.org/t/weighted-histogram-2x-as-slow-in-julia-vs-python/74953/11 "2022-01-20T21:19:29Z")

</div>

I can see a similar difference here, also 1.7.1. Seems like something must be wrong with the `extrema` implementation?

```julia
julia> a = randn(1000);

julia> @btime extrema($a)
  5.849 μs (0 allocations: 0 bytes)
(-3.155119094596879, 3.0526046427540967)

julia> @btime minimum($a), maximum($a)
  1.744 μs (0 allocations: 0 bytes)
(-3.155119094596879, 3.0526046427540967)

```

---

<div class="post-metadata">

**Author:** ![kirklong](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/kirklong/32/19433_2.png) [@kirklong](https://discourse.julialang.org/u/kirklong)\
**Post date:** [January 20, 2022, 9:22pm UTC](https://discourse.julialang.org/t/weighted-histogram-2x-as-slow-in-julia-vs-python/74953/12 "2022-01-20T21:22:12Z")

</div>

I’m only using 1.6.1, and there’s less of a difference but indeed the minimum / maximum way is faster by ~40% which does seem strange.

---

<div class="post-metadata">

**Author:** ![lmiq](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/lmiq/32/18314_2.png) [@lmiq](https://discourse.julialang.org/u/lmiq)\
**Post date:** [January 20, 2022, 9:36pm UTC](https://discourse.julialang.org/t/weighted-histogram-2x-as-slow-in-julia-vs-python/74953/13 "2022-01-20T21:36:02Z")

</div>

[https://github.com/JuliaLang/julia/issues/31442](https://github.com/JuliaLang/julia/issues/31442)

---

<div class="post-metadata">

**Author:** ![tbeason](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/tbeason/32/15898_2.png) [@tbeason](https://discourse.julialang.org/u/tbeason)\
**Post date:** [January 20, 2022, 9:36pm UTC](https://discourse.julialang.org/t/weighted-histogram-2x-as-slow-in-julia-vs-python/74953/14 "2022-01-20T21:36:35Z")

</div>

The fact that `extrema` is slow compared to separate calls of `minimum` and `maximum` has had multiple github issues raised I believe.