# Speed up column looping in matrix

**URL:** <https://discourse.julialang.org/t/speed-up-column-looping-in-matrix/65295>\
**Category:** Performance\
**Created:** [July 26, 2021, 8:34am UTC](https://discourse.julialang.org/t/speed-up-column-looping-in-matrix/65295 "2021-07-26T08:34:11Z")\
**Posts on this page:** 7\
**Page:** 1

<div class="post-metadata">

**Author:** ![HJW019](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/hjw019/32/22102_2.png) [@HJW019](https://discourse.julialang.org/u/HJW019)\
**Post date:** [July 26, 2021, 8:34am UTC](https://discourse.julialang.org/t/speed-up-column-looping-in-matrix/65295/1 "2021-07-26T08:34:11Z")

</div>

I am working on a project where one of the functions is the bottleneck of speed. The problem seems to do with accessing a large matrix. I’d appreciate any advice to improve the code.

The problem is represented by the following MWE. `test3()` is what I need to do in the project and I found it slow. `test1()` and `test2()` are for comparisons.

```julia
using Statistics, BenchmarkTools

function test1()
    c = rand(1000)
    M = rand(256, 1000)
    a = zero(eltype(c))
    @inbounds for i in 1:1000
       @views a += mean(c[i] .+ M[:,1]) # a fixed column of M
    end
    return a
end     

function test2()
    c = rand(1000)
    M = rand(256, 1000)
    a = zero(eltype(c))
    @inbounds for i in 1:1000
        @views a += mean(c[i] .+ M[:, 1 + i % 10]) # cycling through the first 10 columns of M
    end
    return a
end

function test3()
    c = rand(1000)
    M = rand(256, 1000)
    a = zero(eltype(c))
    @inbounds for i in 1:1000
       @views a += mean(c[i] .+ M[:, i]) # using each column of M
    end
    return a
end

```

The execution time and memory allocation is:

```julia
julia> @btime test1();
  789.300 μs (1003 allocations: 4.04 MiB)

julia> @btime test2();
  801.600 μs (1003 allocations: 4.04 MiB)

julia> @btime test3();
  985.400 μs (1003 allocations: 4.04 MiB)

```

`test3()` is more than 20% slower than `test1()` and `test2()`, though memory allocation is the same. If there a way to improve `test3()`'s performance?

Many thanks.

---

<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:** [July 26, 2021, 10:08am UTC](https://discourse.julialang.org/t/speed-up-column-looping-in-matrix/65295/2 "2021-07-26T10:08:41Z")

</div>

> [@HJW019](#):
>
> `c[i] .+ M[:, i]`

This is allocating a new array at every iteration. Preallocating the array of results, evaluating it and then computing the mean will be faster. (I’m assuming here that this is not exactly the computation you need, the other options posted below can be faster for _this_ specific calculation).

```julia
julia> @btime test3() #original
  661.438 μs (1003 allocations: 4.04 MiB)
1007.7474721530921

julia> function test3()
           c = rand(1000)
           M = rand(256, 1000)
           a = zero(eltype(c))
           d = zeros(256)
           @inbounds for i in 1:1000
              @views d .= c[i] .+ M[:,i]
              a += mean(d) # using each column of M
           end
           return a
       end
test3 (generic function with 1 method)

julia> @btime test3()
  310.763 μs (4 allocations: 1.96 MiB)
991.5350998679983

```

---

<div class="post-metadata">

**Author:** ![xiaodai](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/xiaodai/32/15937_2.png) [@xiaodai](https://discourse.julialang.org/u/xiaodai)\
**Post date:** [July 26, 2021, 10:10am UTC](https://discourse.julialang.org/t/speed-up-column-looping-in-matrix/65295/3 "2021-07-26T10:10:12Z")

</div>

I just used algebra. But the slowest part is actually about generating the random numbers

```julia
using Random
function test3()
    Random.seed!(1)
    c = rand(1000)
    M = rand(256, 1000)
    a = zero(eltype(c))
    @inbounds for i in 1:1000
       @views a += mean(c[i] .+ M[:, i]) # using each column of M
    end
    return a
end

using BenchmarkTools

test3()
@benchmark test3()

function test4()
    Random.seed!(1)
    c = rand(1000)
    M = rand(256, 1000)
    a = sum(c) + mean(M)*1000
end

test4()
@benchmark test4()

test3() ≈ test4() #true

```

---

<div class="post-metadata">

**Author:** ![DataFrames](https://avatars.discourse-cdn.com/v4/letter/d/e19b73/32.png) [@DataFrames](https://discourse.julialang.org/u/DataFrames)\
**Post date:** [July 26, 2021, 10:45am UTC](https://discourse.julialang.org/t/speed-up-column-looping-in-matrix/65295/4 "2021-07-26T10:45:48Z")

</div>

```julia
function test3()
           c = rand(1000)
           M = rand(256, 1000)
           a = zero(eltype(c))
           @inbounds for i in 1:1000
               s1 = 0.0
               for j in 1:size(M, 1)
                   s1 += c[i] + M[j,i]
               end
               a += s1/size(M,1)
           end
           return a
       end

```

---

<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:** [July 26, 2021, 11:55am UTC](https://discourse.julialang.org/t/speed-up-column-looping-in-matrix/65295/5 "2021-07-26T11:55:50Z")

</div>

If that is the actual calculation you are performing @DataFrames approach is probably the best. But let us move the random number generation out of the function, and checking the loop speed:

```julia
julia> function test3(M,c)
                  #c = rand(1000)
                  #M = rand(256, 1000)
                  a = zero(eltype(c))
                  @inbounds for i in 1:1000
                      s1 = 0.0
                      for j in 1:size(M, 1)
                          s1 += c[i] + M[j,i]
                      end
                      a += s1/size(M,1)
                  end
                  return a
              end
test3 (generic function with 3 methods)

julia> @btime test3($M,$c) # DataFrames version above
  230.734 μs (0 allocations: 0 bytes)
1020.9477509086946

```

This is probably close to the fastest you can get without multi-threading:

```julia
julia> using LoopVectorization 

julia> function test3(M,c)
         a = zero(eltype(c))
         for i in axes(M,2)
           s1 = 0.0
           @turbo for j in axes(M,1)
             s1 += c[i] + M[j,i]
           end
           a += s1
         end
         return a/size(M,1)
       end
test3 (generic function with 2 methods)

julia> @btime test3($M,$c)
  36.215 μs (0 allocations: 0 bytes)
1020.9477509086943

```

You can get some speedup with muliti-threading:

```julia
julia> using Base.Threads # 4 cores here

julia> function test3(M,c)
         a = zeros(nthreads())
         @threads for i in axes(M,2)
           s = 0.
           @turbo for j in axes(M,1)
             s += c[i] + M[j,i]
           end
           a[threadid()] += s
         end
         return sum(a)/size(M,1)
       end
test3 (generic function with 2 methods)

julia> @btime test3($M,$c)
  18.697 μs (42 allocations: 3.80 KiB)
1020.9477509086946

```

Using the cheap threads from Polyester is even better:

```julia
julia> using Polyester

julia> function test3(M,c)
         a = zeros(nthreads())
         @batch for i in axes(M,2)
           s = 0.
           @turbo for j in axes(M,1)
             s += c[i] + M[j,i]
           end
           a[threadid()] += s
         end
         return sum(a)/size(M,1)
       end
test3 (generic function with 3 methods)

julia> @btime test3($M,$c)
  12.903 μs (1 allocation: 144 bytes)
1020.947750908695

```

---

<div class="post-metadata">

**Author:** ![cgeoga](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/cgeoga/32/216186_2.png) [@cgeoga](https://discourse.julialang.org/u/cgeoga)\
**Post date:** [July 26, 2021, 1:07pm UTC](https://discourse.julialang.org/t/speed-up-column-looping-in-matrix/65295/6 "2021-07-26T13:07:14Z")

</div>

As @lmiq also prefaces, this advice may not be relevant unless this is literally your function, but I think re-arranging some of your operations can really benefit a lot. Here is your function as a one-liner that is pretty close to the speed of the `@turbo` version Leandro gives you:

```julia
test4(M,c) = sum(z->z[1]+mean(z[2]), zip(c, eachcol(M)))

```

Still pretty readable, in my opinion, and I think for operations like this the Julia compiler is pretty good at applying SIMD-like optimizations. This of course depends a lot on the fact that you’re using `+` everywhere, though, so that we can rearrange things so aggressively. This can be reduced even more to

```julia
test5(M,c) = sum(c) + sum(M)/size(M,1)

```

which is even a little bit faster, kind of like @xiaodai is saying.

In general, I think the point that people are making here is that it might make sense to try and optimize the actual math operations a bit more before trying to optimize the code, because there’s a lot of performance on the table just by doing that, even if you don’t want to think about SIMD/threads/eliding bounds checking/etc at all.

---

<div class="post-metadata">

**Author:** ![HJW019](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/hjw019/32/22102_2.png) [@HJW019](https://discourse.julialang.org/u/HJW019)\
**Post date:** [July 26, 2021, 1:54pm UTC](https://discourse.julialang.org/t/speed-up-column-looping-in-matrix/65295/7 "2021-07-26T13:54:23Z")

</div>

I thank all of you for the information and suggestions; they are all _very_ helpful to me! Though my real task is not the literal sum() job, the suggestions pointed me to promising directions.

@lmiq showed the importance of preallocating array and the power of multi-threading; @xiaodai proved that a little math goes a long way; @DataFrames showed the strength of the for-loop in Julia; and @cgeoga summarized it all and the one-liner is amazing!

A good lesson to me about Julia as well as programming in general.