# Matrix product \`M \* transpose(M)\`

**URL:** <https://discourse.julialang.org/t/matrix-product-m-transpose-m/61212>\
**Category:** General Usage\
**Tags:** performance, linearalgebra, tullio, loopvectorization\
**Created:** [May 15, 2021, 7:11pm UTC](https://discourse.julialang.org/t/matrix-product-m-transpose-m/61212 "2021-05-15T19:11:00Z")\
**Posts on this page:** 11\
**Page:** 1

<div class="post-metadata">

**Author:** ![stla](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/stla/32/20382_2.png) [@stla](https://discourse.julialang.org/u/stla)\
**Post date:** [May 15, 2021, 7:11pm UTC](https://discourse.julialang.org/t/matrix-product-m-transpose-m/61212/1 "2021-05-15T19:11:01Z")

</div>

Hello,

Is there a function in Julia to efficiently compute the matrix product `M * transpose(M)`?

This matrix is symmetric, so there should be a more efficient way than doing `M * transpose(M)`. In R one can use the `tcrossprod` function to perform this product.

---

<div class="post-metadata">

**Author:** ![Mason](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/mason/32/2423_2.png) [@Mason](https://discourse.julialang.org/u/Mason)\
**Post date:** [May 15, 2021, 7:27pm UTC](https://discourse.julialang.org/t/matrix-product-m-transpose-m/61212/2 "2021-05-15T19:27:15Z")

</div>

One thing you could do is

```julia
using Tullio, LoopVectorization
@tullio out[i, j] := M[i, k] * M[j, k]

```

This may not take advantage of symmetry, but it’ll take advantage of the fact that it’s just one matrix involved. Here’s how it compares to regular matrix multiplication:

```julia
using Tullio, LoopVectorization

f1(M) = M * transpose(M)
f2(M) = @tullio out[i, j] := M[i, k] * M[j, k]

foreach((10, 100, 1000)) do N
    M = rand(N, N)
    @show N
    @btime f1($M)
    @btime f2($M)
    println()
end

#+RESULTS:
N = 10
  535.363 ns (1 allocation: 896 bytes)
  197.315 ns (1 allocation: 896 bytes)

N = 100
  58.979 μs (2 allocations: 78.20 KiB)
  32.970 μs (50 allocations: 80.52 KiB)

N = 1000
  10.812 ms (2 allocations: 7.63 MiB)
  16.134 ms (76 allocations: 7.63 MiB)

```

As you can see, Tullio eventually starts losing to OpenBLAS for very large matrices, mostly because it doesn’t do array packing.

---

<div class="post-metadata">

**Author:** ![carstenbauer](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/carstenbauer/32/4981_2.png) [@carstenbauer](https://discourse.julialang.org/u/carstenbauer)\
**Post date:** [May 15, 2021, 8:25pm UTC](https://discourse.julialang.org/t/matrix-product-m-transpose-m/61212/3 "2021-05-15T20:25:59Z")

</div>

At the risk of mentioning the obvious: note that in Julia multiple dispatch automatically selects efficient methods of generic functions, like matrix multiply, based on the input types. Hence, by indicating that `M` is symmetric at the type level, i.e. making it a `Symmetric`, Julia will be smart about choosing an efficient multiply method:

```julia
julia> M = rand(100,100); M = (M+M')/2; # M is symmetric but of type Matrix{Float64}

julia> @which M * transpose(M) # dispatches to a generic implementation in matmul.jl
*(A::AbstractMatrix{T} where T, B::AbstractMatrix{T} where T) in LinearAlgebra at /Applications/Julia-1.6.app/Contents/Resources/julia/share/julia/stdlib/v1.6/LinearAlgebra/src/matmul.jl:151

julia> M = Symmetric(M); # M is of type Symmetric{Float64, Matrix{Float64}}

julia> @which M * transpose(M) # dispatches to a specialised method in symmetric.jl
*(A::Union{Hermitian{T, S}, Symmetric{T, S}} where {T, S}, B::Union{Hermitian{T, S}, Symmetric{T, S}} where {T, S}) in LinearAlgebra at /Applications/Julia-1.6.app/Contents/Resources/julia/share/julia/stdlib/v1.6/LinearAlgebra/src/symmetric.jl:572

```

---

<div class="post-metadata">

**Author:** ![stla](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/stla/32/20382_2.png) [@stla](https://discourse.julialang.org/u/stla)\
**Post date:** [May 15, 2021, 8:29pm UTC](https://discourse.julialang.org/t/matrix-product-m-transpose-m/61212/4 "2021-05-15T20:29:47Z")

</div>

Thanks everyone.

@carstenbauer I don’t mean that `M` is symmetric. I mean that the product `M * transpose(M)` is symmetric (whatever `M` is). So in order to get this product it is enough to calculate the entries of the upper triangular part.

---

<div class="post-metadata">

**Author:** ![carstenbauer](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/carstenbauer/32/4981_2.png) [@carstenbauer](https://discourse.julialang.org/u/carstenbauer)\
**Post date:** [May 15, 2021, 8:32pm UTC](https://discourse.julialang.org/t/matrix-product-m-transpose-m/61212/5 "2021-05-15T20:32:03Z")

</div>

Oh, I misunderstood. Sorry for the noise 🙂

---

<div class="post-metadata">

**Author:** ![fredrikekre](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/fredrikekre/32/1688_2.png) [@fredrikekre](https://discourse.julialang.org/u/fredrikekre)\
**Post date:** [May 15, 2021, 8:56pm UTC](https://discourse.julialang.org/t/matrix-product-m-transpose-m/61212/6 "2021-05-15T20:56:51Z")

</div>

You are still correct, although the example you used where wrong. `A * transpose(A)` dispatches to [`syrk`](http://www.netlib.org/lapack/explore-html/d1/d54/group __double__ blas__level3_gae0ba56279ae3fa27c75fefbc4cc73ddf.html) instead of [`gemm`](http://www.netlib.org/lapack/explore-html/d1/d54/group__double __blas__ level3_gaeda3cbd99c8fb834a60a6412878226e1.html#gaeda3cbd99c8fb834a60a6412878226e1).

---

<div class="post-metadata">

**Author:** ![lhnguyen-vn](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/lhnguyen-vn/32/15449_2.png) [@lhnguyen-vn](https://discourse.julialang.org/u/lhnguyen-vn)\
**Post date:** [May 15, 2021, 9:01pm UTC](https://discourse.julialang.org/t/matrix-product-m-transpose-m/61212/7 "2021-05-15T21:01:22Z")

</div>

Julia does this automatically under the hood by calling BLAS’s SYRK routine. You can examine the relevant code snippet in the standard library `LinearAlgebra` [here](https://github.com/JuliaLang/julia/blob/5aca7a37be174ad9f4b28d3586ba9cc112c12d75/stdlib/LinearAlgebra/src/matmul.jl#L371-L379).

---

<div class="post-metadata">

**Author:** ![rafael.guerra](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/rafael.guerra/32/216610_2.png) [@rafael.guerra](https://discourse.julialang.org/u/rafael.guerra)\
**Post date:** [May 16, 2021, 5:32pm UTC](https://discourse.julialang.org/t/matrix-product-m-transpose-m/61212/8 "2021-05-16T17:32:12Z")

</div>

> [@Mason](#):
>
> ```julia
> using Tullio, LoopVectorization
> @tullio out[i, j] := M[i, k] * M[j, k]
> 
> ```

Sorry if it is a dumb question: _why doesn’t Tullio load LoopVectorization as default?_

---

<div class="post-metadata">

**Author:** ![mcabbott](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/mcabbott/32/6603_2.png) [@mcabbott](https://discourse.julialang.org/u/mcabbott)\
**Post date:** [May 16, 2021, 5:56pm UTC](https://discourse.julialang.org/t/matrix-product-m-transpose-m/61212/9 "2021-05-16T17:56:57Z")

</div>

> [@rafael.guerra](#):
>
> why doesn’t Tullio load LoopVectorization as default?

Mostly for startup & compilation times. Startup keeps getting better, but the time to generate & compile these fast versions is still significant:

```julia
julia> @time using Tullio, LoopVectorization
  1.325795 seconds (3.77 M allocations: 230.068 MiB, 2.90% gc time)

julia> M = rand(100,100); @time @tullio out[i, j] := M[i, k] * M[j, k];
  6.619863 seconds (19.40 M allocations: 1.078 GiB, 4.57% gc time, 0.00% compilation time)

```

vs

```julia
julia> @time using Tullio
  0.277620 seconds (766.49 k allocations: 46.810 MiB, 3.51% gc time)

julia> M = rand(100,100); @time @tullio out[i, j] := M[i, k] * M[j, k];
  0.133953 seconds (790.37 k allocations: 43.063 MiB, 3.37% gc time, 97.81% compilation time)

```

---

<div class="post-metadata">

**Author:** ![Euhan](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/euhan/32/36548_2.png) [@Euhan](https://discourse.julialang.org/u/Euhan)\
**Post date:** [December 8, 2022, 4:11pm UTC](https://discourse.julialang.org/t/matrix-product-m-transpose-m/61212/10 "2022-12-08T16:11:27Z")

</div>

Hello [fredrikekre](https://discourse.julialang.org/u/fredrikekre), [lhnguyen-vn](https://discourse.julialang.org/u/lhnguyen-vn)

While I think it is great that this test is in place it doesn’t replace a presumptive dedicated tcrossprod (or whatever name the principle of least surprise dictates) or perhaps tcrossprod! for a few reasons:

1. You don’t get a `Symmetric` in return from `A * transpose(A)` and this can’t be handled automatically.

2. You can’t use a preallocated `Symmetric` with `LinearAlgebra.mul!`, thereby having to allocate about twice as much memory for the result. I suppose that this could potentially change because it would only require dispatching to another method, rather than causing type instability.

3. Such a function could make for clearer, shorter (depending on how long the name will be, of course) and easier-to-read code because one would write `f(A)` rather than `f(A, g(A))`.

4. The presence of this function would allow people to profit from time and memory optimisations you get with symmetric matrices without even having to think about it.

I think this functionality should be available or if it is in some way already I would be curious to know how to access it.

---

<div class="post-metadata">

**Author:** ![jishnub](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/jishnub/32/33620_2.png) [@jishnub](https://discourse.julialang.org/u/jishnub)\
**Post date:** [December 8, 2022, 6:02pm UTC](https://discourse.julialang.org/t/matrix-product-m-transpose-m/61212/11 "2022-12-08T18:02:45Z")

</div>

> You can’t use a preallocated `Symmetric` with `LinearAlgebra.mul!`, thereby having to allocate about twice as much memory for the result

Perhaps `mul!(S::Symmetric, A, transpose(A))` should be special-cased, instead of hitting a generic fallback?