# Optimizing Direct 2D Convolution Code

**URL:** https://discourse.julialang.org/t/optimizing-direct-2d-convolution-code/106598
**Category:** Performance
**Tags:** performance, simd, loopvectorization, image-processing, convolution
**Created:** [November 22, 2023, 9:25pm UTC](https://discourse.julialang.org/t/optimizing-direct-2d-convolution-code/106598 "2023-11-22T21:25:10Z")
**Posts on this page:** 15
**Page:** 1

<div class="post-metadata">

### Author: ![RoyiAvital](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/royiavital/32/571_2.png) [@RoyiAvital](https://discourse.julialang.org/u/RoyiAvital)
#### Post date: [November 22, 2023, 9:25pm UTC](https://discourse.julialang.org/t/optimizing-direct-2d-convolution-code/106598/1 "2023-11-22T21:25:11Z")

</div>

I have implemented 2D convolution using direct calculation:

```julia
using BenchmarkTools;

function _Conv2D!( mO :: Matrix{T}, mI :: Matrix{T}, mK :: Matrix{T} ) where {T <: AbstractFloat}

    numRowsI, numColsI = size(mI);
    numRowsK, numColsK = size(mK);

    for jj ∈ 1:(numColsK - 1), ii ∈ 1:(numRowsK - 1) #<! Top Left
        sumVal = zero(T);
        for nn ∈ 1:numColsK, mm ∈ 1:numRowsK
            ib0 = (jj >= nn) && (ii >= mm);
            @inbounds oa = ib0 ? mI[ii - mm + 1, jj - nn + 1] : zero(T);
            @inbounds sumVal += mK[mm, nn] * oa;
        end
        mO[ii, jj] = sumVal;
    end

    for jj ∈ 1:(numColsK - 1), ii ∈ numRowsK:(numRowsI - 1) #<! Middle Left
        sumVal = zero(T);
        for nn ∈ 1:numColsK, mm ∈ 1:numRowsK
            ib0 = (jj >= nn);
            @inbounds oa = ib0 ? mI[ii - mm + 1, jj - nn + 1] : zero(T);
            @inbounds sumVal += mK[mm, nn] * oa;
        end
        mO[ii, jj] = sumVal;
    end

    for jj ∈ 1:(numColsK - 1), ii ∈ numRowsI:(numRowsI + numRowsK - 1) #<! Bottom Left
        sumVal = zero(T);
        for nn ∈ 1:numColsK, mm ∈ 1:numRowsK
            ib0 = (jj >= nn) && (ii < numRowsI + mm);;
            @inbounds oa = ib0 ? mI[ii - mm + 1, jj - nn + 1] : zero(T);
            @inbounds sumVal += mK[mm, nn] * oa;
        end
        mO[ii, jj] = sumVal;
    end

    for jj ∈ numColsK:(numColsI - 1), ii ∈ 1:(numRowsK - 1) #<! Top Middle
        sumVal = zero(T);
        for nn ∈ 1:numColsK, mm ∈ 1:numRowsK
            ib0 = (ii >= mm);
            @inbounds oa = ib0 ? mI[ii - mm + 1, jj - nn + 1] : zero(T);
            @inbounds sumVal += mK[mm, nn] * oa;
        end
        mO[ii, jj] = sumVal;
    end

    for jj ∈ numColsK:(numColsI - 1), ii ∈ numRowsK:(numRowsI - 1) #<! Middle Middle
        sumVal = zero(T);
        for nn ∈ 1:numColsK, mm ∈ 1:numRowsK
            @inbounds sumVal += mK[mm, nn] * mI[ii - mm + 1, jj - nn + 1];
        end
        mO[ii, jj] = sumVal;
    end

    for jj ∈ numColsK:(numColsI - 1), ii ∈ numRowsI:(numRowsI + numRowsK - 1) #<! Bottom Middle
        sumVal = zero(T);
        for nn ∈ 1:numColsK, mm ∈ 1:numRowsK
            ib0 = (ii < numRowsI + mm);;
            @inbounds oa = ib0 ? mI[ii - mm + 1, jj - nn + 1] : zero(T);
            @inbounds sumVal += mK[mm, nn] * oa;
        end
        mO[ii, jj] = sumVal;
    end

    for jj ∈ numColsI:(numColsI + numColsK - 1), ii ∈ 1:(numRowsK - 1) #<! Top Right
        sumVal = zero(T);
        for nn ∈ 1:numColsK, mm ∈ 1:numRowsK
            ib0 = (jj < numColsI + nn) && (ii >= mm);
            @inbounds oa = ib0 ? mI[ii - mm + 1, jj - nn + 1] : zero(T);
            @inbounds sumVal += mK[mm, nn] * oa;
        end
        mO[ii, jj] = sumVal;
    end

    for jj ∈ numColsI:(numColsI + numColsK - 1), ii ∈ numRowsK:(numRowsI - 1) #<! Middle Right
        sumVal = zero(T);
        for nn ∈ 1:numColsK, mm ∈ 1:numRowsK
            ib0 = (jj < numColsI + nn);
            @inbounds oa = ib0 ? mI[ii - mm + 1, jj - nn + 1] : zero(T);
            @inbounds sumVal += mK[mm, nn] * oa;
        end
        mO[ii, jj] = sumVal;
    end

    for jj ∈ numColsI:(numColsI + numColsK - 1), ii ∈ numRowsI:(numRowsI + numRowsK - 1) #<! Bottom Right
        sumVal = zero(T);
        for nn ∈ 1:numColsK, mm ∈ 1:numRowsK
            ib0 = (jj < numColsI + nn) && (ii < numRowsI + mm);;
            @inbounds oa = ib0 ? mI[ii - mm + 1, jj - nn + 1] : zero(T);
            @inbounds sumVal += mK[mm, nn] * oa;
        end
        mO[ii, jj] = sumVal;
    end

end

mA = rand(1000, 1000);
mK = rand(9, 9);
mO = zeros(size(mA) .+ size(mK) .- (1, 1));

@btime _Conv2D!($mO, $mA, $mK);

```

I wonder how can I optimize the code farther?  
On MATLAB, the direct **built in** implementation `conv2()` executes the above in `6 [Milli Sec]`.  
The Julia implementation above takes `60 [Milli Sec]`.  
I’d be happy for tips to narrow down the difference (MATLAB’s code is probably _SIMD_ vectroized).

~~Moreover, by `@btime` output it seems to allocate while there is no allocation in the code. How is that?~~

~~ 

```julia
julia> @benchmark _Conv2D!($mO, $mA, $mK)
BenchmarkTools.Trial: 1 sample with 1 evaluation.
 Single result which took 6.086 s (1.46% GC) to evaluate,
 with a memory estimate of 6.07 GiB, over 407121552 allocations.

```

 ~~

I’m using `Julia 1.9.4` on Windows.

**Remark** : The concept of the loop is heavily influenced by a past code of @Elrod .  
**Remark** : The allocation issue was a typo of mine as noted by @Dan and @Elrod . Fixed now.

---

<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: [November 22, 2023, 10:07pm UTC](https://discourse.julialang.org/t/optimizing-direct-2d-convolution-code/106598/2 "2023-11-22T22:07:48Z")

</div>

the easy speedup is to move your boundary handling outside the main loop. also if you are convolving with a small array, you might want to use a StaticArray so the compiler knows the size.

---

<div class="post-metadata">

### Author: ![RoyiAvital](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/royiavital/32/571_2.png) [@RoyiAvital](https://discourse.julialang.org/u/RoyiAvital)
#### Post date: [November 22, 2023, 10:17pm UTC](https://discourse.julialang.org/t/optimizing-direct-2d-convolution-code/106598/3 "2023-11-22T22:17:36Z")

</div>

> [@Oscar\_Smith](#):
>
> the easy speedup is to move your boundary handling outside the main loop.

Could you elaborate?

---

<div class="post-metadata">

### Author: ![Dan](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/dan/32/42581_2.png) [@Dan](https://discourse.julialang.org/u/Dan)
#### Post date: [November 22, 2023, 11:04pm UTC](https://discourse.julialang.org/t/optimizing-direct-2d-convolution-code/106598/4 "2023-11-22T23:04:29Z")

</div>

> [@RoyiAvital](#):
>
> Could you elaborate?

I think Oscar means it is better to keep the loops going over the bulk of the 2D arrays as simple as possible (and as fast as possible) without any conditionals checking for edge cases.  
The edges can be handled seperately outside the ‘hot’ loops.

---

<div class="post-metadata">

### Author: ![RoyiAvital](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/royiavital/32/571_2.png) [@RoyiAvital](https://discourse.julialang.org/u/RoyiAvital)
#### Post date: [November 22, 2023, 11:28pm UTC](https://discourse.julialang.org/t/optimizing-direct-2d-convolution-code/106598/5 "2023-11-22T23:28:17Z")

</div>

Could you give a concrete example?  
I can see how to make the loop on `mm` and `nn` depends on `ii` and `jj`.

Yet I refrain from this as it was a guideline by @Elrod that doing so will hurt the vectorization.

---

<div class="post-metadata">

### Author: ![Dan](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/dan/32/42581_2.png) [@Dan](https://discourse.julialang.org/u/Dan)
#### Post date: [November 23, 2023, 12:09am UTC](https://discourse.julialang.org/t/optimizing-direct-2d-convolution-code/106598/6 "2023-11-23T00:09:15Z")

</div>

There is a reference to `mA` in the function which is a global. It’s a typo.

---

<div class="post-metadata">

### Author: ![Raf](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/raf/32/3383_2.png) [@Raf](https://discourse.julialang.org/u/Raf)
#### Post date: [November 23, 2023, 12:13am UTC](https://discourse.julialang.org/t/optimizing-direct-2d-convolution-code/106598/7 "2023-11-23T00:13:54Z")

</div>

1. The compiler knows nothing about your arrays, so the loops will be slow. Use StaticArrays for an easy win.

2. Use conditionals only for the corners of each matrix you are looping over. You don’t need to check every position. Knowing their size at compile time helps calculations here too.

Stencils.jl does these things, and your code would be something like:

```julia
mapstensil(kernelproduct, Kernel(Window(4), mK), mO)

```

(I’m not sure that’s exactly what your code does but it looks like it - kernelproduct just does non-recursive dot product on the kernel and each window)

Its threaded too, so on a laptop I get:

```julia
julia> @time mapstencil(kernelproduct, Kernel(Window(4), mK), mO)
  0.034262 seconds (241 allocations: 7.765 MiB, 0.00% compilation time)

```

and that is also allocating the new array.

It seems there is something else wrong with your code (like type instability) to get this timing, it shouldn’t be that much slower.

---

<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: [November 23, 2023, 3:29am UTC](https://discourse.julialang.org/t/optimizing-direct-2d-convolution-code/106598/8 "2023-11-23T03:29:40Z")

</div>

> [@Dan](#):
>
> The edges can be handled seperately outside the ‘hot’ loops.

Just handle the edges by padding the arrays with the desired boundary conditions.

---

<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: [November 23, 2023, 3:37am UTC](https://discourse.julialang.org/t/optimizing-direct-2d-convolution-code/106598/9 "2023-11-23T03:37:03Z")

</div>

It allocates because `mA` is a global variable.  
You should probably rename the second input argument, `mI`, to `mA`.

---

<div class="post-metadata">

### Author: ![RoyiAvital](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/royiavital/32/571_2.png) [@RoyiAvital](https://discourse.julialang.org/u/RoyiAvital)
#### Post date: [November 23, 2023, 6:27am UTC](https://discourse.julialang.org/t/optimizing-direct-2d-convolution-code/106598/10 "2023-11-23T06:27:31Z")

</div>

> [@Dan](#):
>
> There is a reference to `mA` in the function which is a global. It’s a typo.

I missed that! Now the allocations are solved.  
Thanks for catching that.

---

<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: [November 23, 2023, 6:45am UTC](https://discourse.julialang.org/t/optimizing-direct-2d-convolution-code/106598/11 "2023-11-23T06:45:16Z")

</div>

> [@RoyiAvital](#):
>
> Could you give a concrete example?

Most of your hot loops look like this, with a branch in the central part, preventing vectorization.

> [@RoyiAvital](#):
>
> ```julia
> for nn ∈ 1:numColsK, mm ∈ 1:numRowsK
> ib0 = (jj >= nn) && (ii >= mm);
> @inbounds oa = ib0 ? mI[ii - mm + 1, jj - nn + 1] : zero(T);
> @inbounds sumVal += mK[mm, nn] * oa;
> end
> 
> ```

---

<div class="post-metadata">

### Author: ![RoyiAvital](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/royiavital/32/571_2.png) [@RoyiAvital](https://discourse.julialang.org/u/RoyiAvital)
#### Post date: [November 23, 2023, 6:46am UTC](https://discourse.julialang.org/t/optimizing-direct-2d-convolution-code/106598/12 "2023-11-23T06:46:39Z")

</div>

> [@stevengj](#):
>
> > [@Dan](#):
> >
> > The edges can be handled seperately outside the ‘hot’ loops.
> 
> Just handle the edges by padding the arrays with the desired boundary conditions.

Indeed.  
For that I have the `valid` variant of the code:

```julia
function _Conv2DValid!( mO :: Matrix{T}, mI :: Matrix{T}, mK :: Matrix{T} ) where {T <: AbstractFloat}

    numRowsI, numColsI = size(mI);
    numRowsK, numColsK = size(mK);

    for jj ∈ 1:(numColsI - numColsK + 1)
        @turbo for ii in 1:(numRowsI - numRowsK + 1)
            sumVal = zero(T);
            for nn ∈ 1:numColsK, mm ∈ 1:numRowsK
                @inbounds sumVal += mK[mm, nn] * mI[ii - mm + numRowsK, jj - nn + numColsK];
            end
            mO[ii, jj] = sumVal;
        end
    end

end

```

Yet I wonder if there is something to optimize for the `full` case.  
Specifically beyond adding `@turbo` on the `ii` loop.

---

<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: [November 23, 2023, 7:18am UTC](https://discourse.julialang.org/t/optimizing-direct-2d-convolution-code/106598/13 "2023-11-23T07:18:47Z")

</div>

maybe using linear indices can speed this `_Conv2DValid!` up. that is, if the input are vectors

---

<div class="post-metadata">

### Author: ![Raf](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/raf/32/3383_2.png) [@Raf](https://discourse.julialang.org/u/Raf)
#### Post date: [November 23, 2023, 8:20pm UTC](https://discourse.julialang.org/t/optimizing-direct-2d-convolution-code/106598/14 "2023-11-23T20:20:20Z")

</div>

> [@stevengj](#):
>
> Just handle the edges by padding the arrays with the desired boundary conditions.

Its not always clear that this is better after the allocation. In Stencils.jl `Conditional` (runtime checks without padding) can sometimes beat or be very close to `Halo` (physical padding)

---

<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: [November 23, 2023, 8:48pm UTC](https://discourse.julialang.org/t/optimizing-direct-2d-convolution-code/106598/15 "2023-11-23T20:48:40Z")

</div>

> [@Raf](#):
>
> Its not always clear that this is better after the allocation.

I’m assuming you are doing this over and over, if you care about performance. Whereas you only need to allocate padding once. Ideally you adjust your code so that the padding is built-in from the beginning.