# Performance of naive convolution against Python Numpy

**URL:** <https://discourse.julialang.org/t/performance-of-naive-convolution-against-python-numpy/75603>\
**Category:** New to Julia\
**Tags:** performance, loopvectorization\
**Created:** [February 1, 2022, 8:33pm UTC](https://discourse.julialang.org/t/performance-of-naive-convolution-against-python-numpy/75603 "2022-02-01T20:33:59Z")\
**Posts on this page:** 20\
**Page:** 1

<div class="post-metadata">

**Author:** ![romainvieme](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/romainvieme/32/33349_2.png) [@romainvieme](https://discourse.julialang.org/u/romainvieme)\
**Post date:** [February 1, 2022, 8:34pm UTC](https://discourse.julialang.org/t/performance-of-naive-convolution-against-python-numpy/75603/1 "2022-02-01T20:34:00Z")

</div>

I wrote a naive convolution function for some personal purpose and tested it against Python numpy’s function (which does not use FFT).  
On my computer the Julia script is about 7x slower than Python numpy’s version. Could anyone help me understand why it’s slow or how to accelerate the code?  
Here is my tentative:

```julia
using PyCall
np = pyimport("numpy")

function naive_convol_full!(w, u, v)
	if length(u)<length(v)
		naive_convol_full!(w, v, u)
	else
		n=length(u)
		m=length(v)
		for i=0:n+m-2
			k=max(i+1-m, 0)
			l=min(n-1,i) 

			w[i+1]=zero(eltype(u))
			for j=k:l
				w[i+1]+=u[j+1]*v[i+1-j]
			end
		end
	end
end
A = rand(10000); B=rand(10000); D = zeros(Float64, (length(A)+ length(B)-1));

```

and the REPL output:

```julia
julia> @time naive_convol_full!(D, A, B)
  0.126681 seconds
julia> @time C=np.convolve(A, B, "full");
  0.018948 seconds (41 allocations: 158.094 KiB)
julia> C≈D
true

```

---

<div class="post-metadata">

**Author:** ![juliohm](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/juliohm/32/215266_2.png) [@juliohm](https://discourse.julialang.org/u/juliohm)\
**Post date:** [February 1, 2022, 8:40pm UTC](https://discourse.julialang.org/t/performance-of-naive-convolution-against-python-numpy/75603/2 "2022-02-01T20:40:55Z")

</div>

Try to add @inbounds or @simd calls to the hot loops. Currently your implementation is checking if the indices are valid and that slow things down for example.

---

<div class="post-metadata">

**Author:** ![Oscar\_Smith](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/oscar_smith/32/25343_2.png) [@Oscar\_Smith](https://discourse.julialang.org/u/Oscar_Smith)\
**Post date:** [February 1, 2022, 8:45pm UTC](https://discourse.julialang.org/t/performance-of-naive-convolution-against-python-numpy/75603/3 "2022-02-01T20:45:24Z")

</div>

Also, it looks like your current implementation was copied from a language with 0 indexing.

```julia
function naive_convol_full!(w, u, v)
	if length(u)<length(v)
		return naive_convol_full!(w, v, u)
	end
        n=length(u)
	m=length(v)
	@inbounds @fastmath for i=1:n+m-1
		w[i]=zero(eltype(u))
		for j=max(i-m, 0):(min(n,i) - 1)
			w[i]+=u[j+1]*v[i-j]
		end
	end
end

```

This takes @juliohm’s suggestion and also changes `i+1` to `i` which makes things a little clearer in my opinion.

---

<div class="post-metadata">

**Author:** ![rveltz](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/rveltz/32/2707_2.png) [@rveltz](https://discourse.julialang.org/u/rveltz)\
**Post date:** [February 1, 2022, 8:46pm UTC](https://discourse.julialang.org/t/performance-of-naive-convolution-against-python-numpy/75603/4 "2022-02-01T20:46:51Z")

</div>

why would numpy version be this double for loop and not fft based (for example) ?

---

<div class="post-metadata">

**Author:** ![stevengj](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/stevengj/32/71_2.png) [@stevengj](https://discourse.julialang.org/u/stevengj)\
**Post date:** [February 1, 2022, 8:53pm UTC](https://discourse.julialang.org/t/performance-of-naive-convolution-against-python-numpy/75603/5 "2022-02-01T20:53:01Z")

</div>

> [@romainvieme](#):
>
> Python numpy’s function

If you drill down in the Numpy source code, I think its convolution ends up at [`multiarraymodule.c:_pyarray_correlate`](https://github.com/numpy/numpy/blob/ff4b33edcf6819f24f9a1b4a483d64873cea2d5c/numpy/core/src/multiarray/multiarraymodule.c#L1174-L1290). It looks like it just does a sequence of `dot` calls, which probably is some optimized SIMD-based function, and possibly uses threads.

To get similar SIMD utilization, you might want to try [GitHub - JuliaSIMD/LoopVectorization.jl: Macro(s) for vectorizing loops.](https://github.com/JuliaSIMD/LoopVectorization.jl)

---

<div class="post-metadata">

**Author:** ![romainvieme](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/romainvieme/32/33349_2.png) [@romainvieme](https://discourse.julialang.org/u/romainvieme)\
**Post date:** [February 1, 2022, 8:53pm UTC](https://discourse.julialang.org/t/performance-of-naive-convolution-against-python-numpy/75603/6 "2022-02-01T20:53:59Z")

</div>

Thanks, didn’t know this trick ! Doesn’t make much of a difference though…

```julia
function naive_convol_full!(w, u, v)
	if length(u)<length(v)
		naive_convol_full!(v, u, w)
	else
		n=length(u)
		m=length(v)
		@simd for i=0:n+m-2
			k=max(i+1-m, 0)
			l=min(n-1,i) 

			@inbounds w[i+1]=zero(eltype(u))
			@simd for j=k:l
				@inbounds w[i+1]+=u[j+1]*v[i+1-j]
			end
		end
	end
end

```

REPL

```julia
julia> @time naive_convol_full!(D, A,B)
  0.129330 seconds

```

---

<div class="post-metadata">

**Author:** ![romainvieme](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/romainvieme/32/33349_2.png) [@romainvieme](https://discourse.julialang.org/u/romainvieme)\
**Post date:** [February 1, 2022, 8:55pm UTC](https://discourse.julialang.org/t/performance-of-naive-convolution-against-python-numpy/75603/7 "2022-02-01T20:55:39Z")

</div>

Sorry, not copied, I wrote it myself 😅 Thanks for the code though !  
Just ran it, I gained a few milliseconds but not enough to make the difference

```julia
julia> function naive_convol_full!(w, u, v)
               if length(u)<length(v)
                       return naive_convol_full!(w, v, u)
               end
               n=length(u)
               m=length(v)
               @inbounds @fastmath for i=1:n+m-1
                       w[i]=zero(eltype(u))
                       for j=max(i-m, 0):(min(n,i) - 1)
                               w[i]+=u[j+1]*v[i-j]
                       end
               end
       end
naive_convol_full! (generic function with 1 method)

julia> @time naive_convol_full!(D, A,B)
  0.142548 seconds (35.08 k allocations: 1.943 MiB, 11.77% compilation time)

julia> @time naive_convol_full!(D, A,B)
  0.126336 seconds

```

---

<div class="post-metadata">

**Author:** ![romainvieme](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/romainvieme/32/33349_2.png) [@romainvieme](https://discourse.julialang.org/u/romainvieme)\
**Post date:** [February 1, 2022, 8:56pm UTC](https://discourse.julialang.org/t/performance-of-naive-convolution-against-python-numpy/75603/8 "2022-02-01T20:56:40Z")

</div>

I don’t actually know, but that’s what’s said in the docs !

---

<div class="post-metadata">

**Author:** ![rveltz](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/rveltz/32/2707_2.png) [@rveltz](https://discourse.julialang.org/u/rveltz)\
**Post date:** [February 1, 2022, 9:02pm UTC](https://discourse.julialang.org/t/performance-of-naive-convolution-against-python-numpy/75603/9 "2022-02-01T21:02:53Z")

</div>

I think @stevengj answer goes into that direction

---

<div class="post-metadata">

**Author:** ![romainvieme](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/romainvieme/32/33349_2.png) [@romainvieme](https://discourse.julialang.org/u/romainvieme)\
**Post date:** [February 1, 2022, 9:45pm UTC](https://discourse.julialang.org/t/performance-of-naive-convolution-against-python-numpy/75603/10 "2022-02-01T21:45:29Z")

</div>

Thanks, I couldn’t figure out how to use `@turbo` from `LoopVectorization` but I got a 3x speedup by replacing the innermost loop by a dot product (following your suggestion). That’ll have to do for now ! The numpy call is still about 2.2x faster than the Julia one though.  
Thank you all for your help — I’m still curious if anyone can figure out a better way!

```julia
using LinearAlgebra, PyCall
np = pyimport("numpy")

function naive_convol_full!(w, u, v)
	if length(u)<length(v)
		return naive_convol_full!(w, v, u)
	end
        n=length(u)
	m=length(v)
	@inbounds for i=1:n+m-1
		w[i]=zero(eltype(u))
		@simd for j=max(i-m, 0):(min(n,i) - 1)
			w[i]+=u[j+1]*v[i-j]
		end
	end
end

function matmul_convol_full!(w, u, v)
	if length(u)<length(v)
		matmul_convol_full!(w, v, u)
	else
		n=length(u)
		m=length(v)
		@inbounds for i=1:n+m-1
			k=max(i-m, 0)
			l=(min(n,i) - 1) 

			@views w[i]= dot(u[k+1:l+1], v[i-k:-1:i-l])
		end
	end
end
A = rand(10000); B=rand(10000); D = zeros(Float64, (length(A)+ length(B)-1));

```

```julia
julia> @time naive_convol_full!(D, A, B)
  0.128569 seconds
julia> @time matmul_convol_full!(D, A, B)
  0.041894 seconds
julia> @time C=np.convolve(A, B, "full");
  0.018964 seconds (41 allocations: 158.094 KiB)

```

---

<div class="post-metadata">

**Author:** ![Eben60](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/eben60/32/13475_2.png) [@Eben60](https://discourse.julialang.org/u/Eben60)\
**Post date:** [February 1, 2022, 10:20pm UTC](https://discourse.julialang.org/t/performance-of-naive-convolution-against-python-numpy/75603/11 "2022-02-01T22:20:18Z")

</div>

[romainvieme](https://discourse.julialang.org/u/romainvieme), without actually checking your code, may I point you to a discussion which followed this post:

> [@General questions from Python user](https://discourse.julialang.org/t/general-questions-from-python-user/55475/44):
>
> In the course of my weekend procrastinations I contrieved an example to check my statement: some home-made implementation of convolution in both languages. Make up some signal and kernel data in Python (unimportant part which is not benchmarked): import numpy as np def lorentz(x, gamma, mu=0): return (gamma\*\*2/(((x-mu)\*\*2)+gamma\*\*2))/(np.pi\*gamma) def signal(x): return lorentz(x, 0.1, -0.3) + lorentz(x, 0.03, 0.1) + lorentz(x, 0.2, 0.25) def gauss(x, sigma, mu=0.0): return np.e…

---

<div class="post-metadata">

**Author:** ![romainvieme](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/romainvieme/32/33349_2.png) [@romainvieme](https://discourse.julialang.org/u/romainvieme)\
**Post date:** [February 1, 2022, 10:28pm UTC](https://discourse.julialang.org/t/performance-of-naive-convolution-against-python-numpy/75603/12 "2022-02-01T22:28:15Z")

</div>

Looks relevant indeed 😄 I didn’t bump into that in my preliminary research ! I’ll be sure to have a look. Thanks !

---

<div class="post-metadata">

**Author:** ![Seif\_Shebl](https://avatars.discourse-cdn.com/v4/letter/s/eada6e/32.png) [@Seif\_Shebl](https://discourse.julialang.org/u/Seif_Shebl)\
**Post date:** [February 2, 2022, 12:58am UTC](https://discourse.julialang.org/t/performance-of-naive-convolution-against-python-numpy/75603/13 "2022-02-02T00:58:14Z")

</div>

The naive loop is faster than numpy, just make sure it gets SIMD’ed:

```julia
function naive_convol_full!(w, u, v)
    if length(u) < length(v)
        return naive_convol_full!(w, v, u)
    end
    n = length(u)
    m = length(v)
    for i = 1:n+m-1
        s = zero(eltype(u))
        @simd for j = max(i-m,0):min(n,i)-1
            s += u[j+1] * v[i-j]
        end
        w[i] = s
    end
end

@btime naive_convol_full!(D,A,B) setup=(A=rand(10000); B=rand(10000); D=zeros(length(A)+length(B)-1)) evals=1
  15.383 ms (0 allocations: 0 bytes)

```

And if you use `@tturbo` from `LoopVectorizations` , you get it even faster:

```julia
    ...
        @tturbo for j = max(i-m,0):min(n,i)-1
            s += u[j+1] * v[i-j]
        end
    ...

@btime naive_convol_full!(D,A,B) setup=(A=rand(10000); B=rand(10000); D=zeros(length(A)+length(B)-1)) evals=1
  9.398 ms (0 allocations: 0 bytes)

```

---

<div class="post-metadata">

**Author:** ![romainvieme](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/romainvieme/32/33349_2.png) [@romainvieme](https://discourse.julialang.org/u/romainvieme)\
**Post date:** [February 2, 2022, 11:33am UTC](https://discourse.julialang.org/t/performance-of-naive-convolution-against-python-numpy/75603/14 "2022-02-02T11:33:08Z")

</div>

Great ! Thanks ! That seems to work on my machine

```julia
function naive_convol_full!(w, u, v)
    if length(u) < length(v)
        return naive_convol_full!(w, v, u)
    end
    n = length(u)
    m = length(v)
    @inbounds for i = 1:n+m-1
        s = zero(eltype(u))
        @tturbo for j = max(i-m,0):min(n,i)-1
            s += u[j+1] * v[i-j]
        end
        w[i] = s
    end
end

```

```julia
julia> @btime naive_convol_full!(D,A,B) setup=(A=rand(10000); B=rand(10000); D=zeros(length(A)+length(B)-1)) evals=100;
  19.402 ms (0 allocations: 0 bytes)

julia> @btime C=np.convolve(A, B, "full") setup=(A=rand(10000); B=rand(10000)) evals=100;
  18.031 ms (41 allocations: 158.09 KiB)

```

So the difference is 1.5ms. Works great !

---

<div class="post-metadata">

**Author:** ![romainvieme](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/romainvieme/32/33349_2.png) [@romainvieme](https://discourse.julialang.org/u/romainvieme)\
**Post date:** [February 2, 2022, 11:48am UTC](https://discourse.julialang.org/t/performance-of-naive-convolution-against-python-numpy/75603/15 "2022-02-02T11:48:25Z")

</div>

Doing further homework it does not seem to scale nicely

```julia
julia> @btime C=np.convolve(A, B, "full") setup=(A=rand(100000); B=rand(10000)) evals=100;
  182.196 ms (41 allocations: 861.22 KiB)
julia> @btime naive_convol_full!(D,A,B) setup=(A=rand(100000); B=rand(10000); D=zeros(length(A)+length(B)-1)) evals=100;
  202.405 ms (0 allocations: 0 bytes)
julia> @btime naive_convol_full!(D,A,B) setup=(A=rand(100000); B=rand(100000); D=zeros(length(A)+length(B)-1)) evals=100;
  3.355 s (0 allocations: 0 bytes)
julia> @btime C=np.convolve(A, B, "full") setup=(A=rand(100000); B=rand(100000)) evals=100;
  724.136 ms (41 allocations: 1.53 MiB)

```

So I still don’t understand everything ! On another laptop however, the Julia code is clearly faster. I guess that really depends on hardware then…

Edit: The other laptop

```julia
julia> @btime naive_convol_full!(D,A,B) setup=(A=rand(10000); B=rand(10000); D=zeros(length(A)+length(B)-1)) evals=100;
  15.199 ms (0 allocations: 0 bytes)

julia> @btime C=np.convolve(A, B, "full") setup=(A=rand(10000); B=rand(10000)) evals=100;
  111.462 ms (41 allocations: 158.09 KiB)

```

---

<div class="post-metadata">

**Author:** ![romainvieme](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/romainvieme/32/33349_2.png) [@romainvieme](https://discourse.julialang.org/u/romainvieme)\
**Post date:** [February 2, 2022, 12:01pm UTC](https://discourse.julialang.org/t/performance-of-naive-convolution-against-python-numpy/75603/16 "2022-02-02T12:01:03Z")

</div>

By any chance did you benchmark the Python version on your laptop?

---

<div class="post-metadata">

**Author:** ![Seif\_Shebl](https://avatars.discourse-cdn.com/v4/letter/s/eada6e/32.png) [@Seif\_Shebl](https://discourse.julialang.org/u/Seif_Shebl)\
**Post date:** [February 2, 2022, 12:25pm UTC](https://discourse.julialang.org/t/performance-of-naive-convolution-against-python-numpy/75603/17 "2022-02-02T12:25:28Z")

</div>

> [@romainvieme](#):
>
> `100000`

Notice that Numpy uses multithreading, so you should be comparing with the multithreaded version in Julia too:

```julia
  9.295 ms (0 allocations: 0 bytes) # Julia (10000)
  17.367 ms (41 allocations: 158.09 KiB) # numpy (10000)

  873.935 ms (1 allocation: 32 bytes) # Julia (100000)
  1.213 s (41 allocations: 1.53 MiB) # numpy (100000)

```

---

<div class="post-metadata">

**Author:** ![romainvieme](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/romainvieme/32/33349_2.png) [@romainvieme](https://discourse.julialang.org/u/romainvieme)\
**Post date:** [February 2, 2022, 12:38pm UTC](https://discourse.julialang.org/t/performance-of-naive-convolution-against-python-numpy/75603/18 "2022-02-02T12:38:34Z")

</div>

> [@Seif\_Shebl](#):
>
> Notice that Numpy uses multithreading

I suspected as well, thanks for confirming 😀 I’ve tried to shut down this feature of numpy but didn’t succeed. Thanks again, I don’t have any further questions !

---

<div class="post-metadata">

**Author:** ![stevengj](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/stevengj/32/71_2.png) [@stevengj](https://discourse.julialang.org/u/stevengj)\
**Post date:** [February 2, 2022, 12:42pm UTC](https://discourse.julialang.org/t/performance-of-naive-convolution-against-python-numpy/75603/19 "2022-02-02T12:42:20Z")

</div>

> [@romainvieme](#):
>
> I’ve tried to shut down this feature of numpy but didn’t succeed.

> <https://stackoverflow.com/questions/17053671/how-do-you-stop-numpy-from-multithreading>

---

<div class="post-metadata">

**Author:** ![romainvieme](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/romainvieme/32/33349_2.png) [@romainvieme](https://discourse.julialang.org/u/romainvieme)\
**Post date:** [February 2, 2022, 12:45pm UTC](https://discourse.julialang.org/t/performance-of-naive-convolution-against-python-numpy/75603/20 "2022-02-02T12:45:42Z")

</div>

Yep, found that one but couldn’t figure out which BLAS library it’s using on my system (up-to-date Debian stable with Python 3.9.2). I’ll give it more thought.

```julia
In [3]: numpy. __config__.show()
blas_mkl_info:
  NOT AVAILABLE
blis_info:
  NOT AVAILABLE
openblas_info:
  NOT AVAILABLE
atlas_3_10_blas_threads_info:
  NOT AVAILABLE
atlas_3_10_blas_info:
  NOT AVAILABLE
atlas_blas_threads_info:
  NOT AVAILABLE
atlas_blas_info:
  NOT AVAILABLE
accelerate_info:
  NOT AVAILABLE
blas_info:
    libraries = ['blas', 'blas']
    library_dirs = ['/usr/lib/x86_64-linux-gnu']
    include_dirs = ['/usr/local/include', '/usr/include']
    language = c
    define_macros = [('HAVE_CBLAS', None)]
blas_opt_info:
    define_macros = [('NO_ATLAS_INFO', 1), ('HAVE_CBLAS', None)]
    libraries = ['blas', 'blas']
    library_dirs = ['/usr/lib/x86_64-linux-gnu']
    include_dirs = ['/usr/local/include', '/usr/include']
    language = c
lapack_mkl_info:
  NOT AVAILABLE
openblas_lapack_info:
  NOT AVAILABLE
openblas_clapack_info:
  NOT AVAILABLE
flame_info:
  NOT AVAILABLE
atlas_3_10_threads_info:
  NOT AVAILABLE
atlas_3_10_info:
  NOT AVAILABLE
atlas_threads_info:
  NOT AVAILABLE
atlas_info:
  NOT AVAILABLE
lapack_info:
    libraries = ['lapack', 'lapack']
    library_dirs = ['/usr/lib/x86_64-linux-gnu']
    language = f77
lapack_opt_info:
    libraries = ['lapack', 'lapack', 'blas', 'blas']
    library_dirs = ['/usr/lib/x86_64-linux-gnu']
    language = c
    define_macros = [('NO_ATLAS_INFO', 1), ('HAVE_CBLAS', None)]
    include_dirs = ['/usr/local/include', '/usr/include']

```

[Next page](https://discourse.julialang.org/t/performance-of-naive-convolution-against-python-numpy/75603.md?page=2)