# Why matrix multiplication is much slower than PyTorch

**URL:** <https://discourse.julialang.org/t/why-matrix-multiplication-is-much-slower-than-pytorch/63661>\
**Category:** Performance\
**Tags:** linearalgebra, python, pytorch\
**Created:** [June 27, 2021, 7:43pm UTC](https://discourse.julialang.org/t/why-matrix-multiplication-is-much-slower-than-pytorch/63661 "2021-06-27T19:43:16Z")\
**Posts on this page:** 5\
**Page:** 1

<div class="post-metadata">

**Author:** ![yingqiuz](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/yingqiuz/32/25165_2.png) [@yingqiuz](https://discourse.julialang.org/u/yingqiuz)\
**Post date:** [June 27, 2021, 7:43pm UTC](https://discourse.julialang.org/t/why-matrix-multiplication-is-much-slower-than-pytorch/63661/1 "2021-06-27T19:43:16Z")

</div>

For example, matrix multiplication of 10,000 x 10,100 matrices, single threaded

In julia:

```julia
BLAS.set_num_threads(1)
A = randn(10000, 10000)
B = randn(10000, 10000)
C = Matrix{Float64}(undef, 10000, 10000)

@benchmark mul!($C, $A, $B)

```

The result shows

```julia
BenchmarkTools.Trial: 
  memory estimate: 0 bytes
  allocs estimate: 0
  --------------
  minimum time: 18.192 s (0.00% GC)
  median time: 18.192 s (0.00% GC)
  mean time: 18.192 s (0.00% GC)
  maximum time: 18.192 s (0.00% GC)
  --------------
  samples: 1
  evals/sample: 1

```

I also tried LoopVectorization:

```julia
function mygemmavx!(C, A, B)
    @turbo for m ∈ axes(A,1), n ∈ axes(B,2)
        Cmn = zero(eltype(C))
        for k ∈ axes(A,2)
            Cmn += A[m,k] * B[k,n]
        end
        C[m,n] = Cmn
    end
end
@benchmark mygemmavx!($C, $A, $B)

```

Results:

```julia
BenchmarkTools.Trial: 
  memory estimate: 0 bytes
  allocs estimate: 0
  --------------
  minimum time: 149.960 s (0.00% GC)
  median time: 149.960 s (0.00% GC)
  mean time: 149.960 s (0.00% GC)
  maximum time: 149.960 s (0.00% GC)
  --------------
  samples: 1
  evals/sample: 1

```

However, using Pytorch:

```julia
import torch

torch.set_num_threads(1)
A = torch.randn(10000, 10000)
B = torch.randn(10000, 10000)
%timeit -n 3 torch.matmul(A, B)
torch.matmul(A, B)

```

It took only ~9s in average (7 runs, 3 loops each)

I’m so confused. Could anyone shed some light on this? Thank you!

---

<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:** [June 27, 2021, 7:46pm UTC](https://discourse.julialang.org/t/why-matrix-multiplication-is-much-slower-than-pytorch/63661/2 "2021-06-27T19:46:14Z")

</div>

> [@yingqiuz](#):
>
> `set_num_threads`

because of this maybe: [pytorch not respecting torch.set\_num\_threads · Issue #975 · pytorch/pytorch · GitHub](https://github.com/pytorch/pytorch/issues/975)

---

<div class="post-metadata">

**Author:** ![yingqiuz](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/yingqiuz/32/25165_2.png) [@yingqiuz](https://discourse.julialang.org/u/yingqiuz)\
**Post date:** [June 27, 2021, 7:50pm UTC](https://discourse.julialang.org/t/why-matrix-multiplication-is-much-slower-than-pytorch/63661/3 "2021-06-27T19:50:59Z")

</div>

Thanks - but I monitored the usage and pytorch indeed took up to 100% cpu. Does it suggest it was single-threaded?

---

<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:** [June 27, 2021, 7:51pm UTC](https://discourse.julialang.org/t/why-matrix-multiplication-is-much-slower-than-pytorch/63661/4 "2021-06-27T19:51:54Z")

</div>

also, eh, pytorch default type is float32…

```julia
In [8]: torch.set_num_threads(1)
   ...: A = torch.randn(1000, 1000, dtype=torch.float64)
   ...: B = torch.randn(1000, 1000, dtype=torch.float64)
   ...: %timeit -n 5 torch.matmul(A, B)
43.1 ms ± 810 µs per loop (mean ± std. dev. of 7 runs, 5 loops each)

In [9]: torch.set_num_threads(1)
   ...: A = torch.randn(1000, 1000)
   ...: B = torch.randn(1000, 1000)
   ...: %timeit -n 5 torch.matmul(A, B)
22.1 ms ± 245 µs per loop (mean ± std. dev. of 7 runs, 5 loops each)

```

```julia
julia> A = randn(Float32, 1000, 1000);

julia> B = randn(Float32, 1000, 1000);

julia> C = Matrix{Float32}(undef, 1000, 1000);

julia> @benchmark mul!($C, $A, $B)
BenchmarkTools.Trial: 
  memory estimate: 0 bytes
  allocs estimate: 0
  --------------
  minimum time: 20.951 ms (0.00% GC)
  median time: 22.395 ms (0.00% GC)
  mean time: 22.528 ms (0.00% GC)
  maximum time: 25.443 ms (0.00% GC)
  --------------
  samples: 222
  evals/sample: 1

```

* * *

we really need some FAQ saying please don’t benchmark a single BLAS instruction because they literally should be the same. In the case they are not, it’s either OpenBLAS vs. MKL, or one of them is not respecting # of threads setting

---

<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:** [June 27, 2021, 8:49pm UTC](https://discourse.julialang.org/t/why-matrix-multiplication-is-much-slower-than-pytorch/63661/5 "2021-06-27T20:49:27Z")

</div>

> [@yingqiuz](#):
>
> I also tried LoopVectorization:

It’s missing optimizations for handling such large sizes. Single threaded, those optimizations are needed above 300x300, or so.  
`Octavian.jl` will do much better, but will probably be a little behind OpenBLAS.

```julia
julia> M = K = N = 10_000; T = Float32; A = rand(T,M,K); B = rand(T,K,N); C0 = Array{T}(undef, M, N);

julia> @time mul!(C0, A, B);
  9.055349 seconds (2.27 M allocations: 120.548 MiB, 0.92% gc time, 5.70% compilation time)

julia> @time matmul_serial!(C0, A, B);
 21.109921 seconds (23.04 M allocations: 1.302 GiB, 0.51% gc time, 52.73% compilation time)

julia> @time mul!(C0, A, B);
  8.565136 seconds

julia> @time matmul_serial!(C0, A, B);
  9.987270 seconds

```

Octavian’s single threaded performance could use some work; the multithreaded performance is closer.  
I think one reason why the performance is worse is because Octavian will only use a fraction of the L3 cache equal to the number of threads it’s using.  
It (incorrectly, here) assumes that you aren’t running single threaded matmul benchmarks, and that if you’re running it single threaded, that’s because the other cores are busy doing something else that they’re going to need L3 cache for. Like maybe they’re also calling single threaded matmuls.