# Improve the performance of multiplication of an arbitrary number of matrices

**URL:** <https://discourse.julialang.org/t/improve-the-performance-of-multiplication-of-an-arbitrary-number-of-matrices/10835>\
**Category:** Performance\
**Created:** [May 11, 2018, 5:01am UTC](https://discourse.julialang.org/t/improve-the-performance-of-multiplication-of-an-arbitrary-number-of-matrices/10835 "2018-05-11T05:01:49Z")\
**Posts on this page:** 20\
**Page:** 1

<div class="post-metadata">

**Author:** ![Ronis\_BR](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/ronis_br/32/50999_2.png) [@Ronis\_BR](https://discourse.julialang.org/u/Ronis_BR)\
**Post date:** [May 11, 2018, 5:01am UTC](https://discourse.julialang.org/t/improve-the-performance-of-multiplication-of-an-arbitrary-number-of-matrices/10835/1 "2018-05-11T05:01:49Z")

</div>

Hi guys!

I need to make a function that multiplies an arbitrary number of matrices. My current attempt is the following:

```julia
@inline function compose_rotation(D1::Matrix, Ds::Matrix...)
    for Di in Ds
        D1 = Di*D1
    end

    D1
end

```

The problem is that the performance is much worse than only the multiplication:

```julia
D1 = rand(3,3)
D2 = rand(3,3)

@time for i=1:10000000
    compose_rotation(D1,D2,D1,D2)
end

2.406881 seconds (30.00 M allocations: 4.470 GiB, 1.90% gc time)

@time for i=1:10000000
    D1*D2*D1*D2
end

2.004193 seconds (30.00 M allocations: 4.470 GiB, 2.44% gc time)

```

I am wondering: since I will always know how many matrices it is necessary to multiply at compile time, is there any way to “unroll” to loop so that I can improve the performance?

P.S.: This seems a little silly, but this function will be used to compute a composed rotation regardless the choice of the rotation description (matrix, quaternion, Euler axis/angle, etc.).

---

<div class="post-metadata">

**Author:** ![rdeits](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/rdeits/32/286_2.png) [@rdeits](https://discourse.julialang.org/u/rdeits)\
**Post date:** [May 11, 2018, 5:08am UTC](https://discourse.julialang.org/t/improve-the-performance-of-multiplication-of-an-arbitrary-number-of-matrices/10835/2 "2018-05-11T05:08:04Z")

</div>

1. You’re timing in global scope, so your results are not going to be informative. For a discussion of the same issue from earlier today, see: [How to assign real/imaginary part of complex array - #11 by rdeits](https://discourse.julialang.org/t/how-to-assign-real-imaginary-part-of-complex-array/10818/11)

2. You’ll get substantially better performance if you use [https://github.com/JuliaArrays/StaticArrays.jl](https://github.com/JuliaArrays/StaticArrays.jl) for your small fixed-size attrays.

---

<div class="post-metadata">

**Author:** ![rdeits](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/rdeits/32/286_2.png) [@rdeits](https://discourse.julialang.org/u/rdeits)\
**Post date:** [May 11, 2018, 5:12am UTC](https://discourse.julialang.org/t/improve-the-performance-of-multiplication-of-an-arbitrary-number-of-matrices/10835/3 "2018-05-11T05:12:28Z")

</div>

You also might be interested in [GitHub - JuliaGeometry/CoordinateTransformations.jl: A fresh approach to coordinate transformations...](https://github.com/FugroRoames/CoordinateTransformations.jl) which already supports compositions of transformations and has been stable and performant for quite a while.

---

<div class="post-metadata">

**Author:** ![Ronis\_BR](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/ronis_br/32/50999_2.png) [@Ronis\_BR](https://discourse.julialang.org/u/Ronis_BR)\
**Post date:** [May 11, 2018, 5:25am UTC](https://discourse.julialang.org/t/improve-the-performance-of-multiplication-of-an-arbitrary-number-of-matrices/10835/4 "2018-05-11T05:25:56Z")

</div>

Thanks!

But this is actually for my package:

[https://github.com/ronisbr/ReferenceFrameRotations.jl](https://github.com/ronisbr/ReferenceFrameRotations.jl)

I am almost finishing a toolbox with functions related to Satellite simulations. In this toolbox, there are a lot of functions to create rotations between reference frames (J2000, GCRF, PEF, TOD, MOD, TEME, etc.). I want that each function can compute the rotation using Quaternions ou DCMs. If I do not have this `compose_rotation` function, then I will have to create two functions for each rotation (which is a lot, really…).

I will try to benchmark inside a function and will post the results here!

---

<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:** [May 11, 2018, 6:40am UTC](https://discourse.julialang.org/t/improve-the-performance-of-multiplication-of-an-arbitrary-number-of-matrices/10835/5 "2018-05-11T06:40:13Z")

</div>

I’m not sure that benchmarking in a function will make a big difference in this specific case (although you should never benchmark in the global scope, and you should _definitely_ check out BenchmarkTools, I never benchmark with anything else).

What will make a difference is using StaticArrays, you’ll get 5-10x speedup, probably, at least if you avoid the splatting and the loop. You can unroll the loop using `@generated` functions.

But, question: If you always know the number composed rotations at compile time, can you just avoid the whole `compose` function, and just type out the multiplication inline in your code, that is, just put `D1 * D2 * D3 * D4` in there?

---

<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:** [May 11, 2018, 6:50am UTC](https://discourse.julialang.org/t/improve-the-performance-of-multiplication-of-an-arbitrary-number-of-matrices/10835/6 "2018-05-11T06:50:17Z")

</div>

BTW, for reasons that are unclear to me, I get significantly improved performance if, in addition to using StaticArrays, I add parens, like this:

```julia
((D1 * D2) * D3) * D4

```

---

<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:** [May 11, 2018, 7:00am UTC](https://discourse.julialang.org/t/improve-the-performance-of-multiplication-of-an-arbitrary-number-of-matrices/10835/7 "2018-05-11T07:00:49Z")

</div>

```julia
using StaticArrays
using BenchmarkTools

function composeN(D1, Ds...)
   for D in Ds
       D1 = D1 * D
   end
   return D1
end

compose4(D1, D2, D3, D4) = D1 * D2 * D3 * D4
compose4_(D1, D2, D3, D4) = ((D1 * D2) * D3) * D4
@inline compose4i(D1, D2, D3, D4) = ((D1 * D2) * D3) * D4

D1 = rand(3, 3)
D2 = rand(3, 3)
D3 = rand(3, 3)
D4 = rand(3, 3)

S1 = @SMatrix rand(3, 3)
S2 = @SMatrix rand(3, 3)
S3 = @SMatrix rand(3, 3)
S4 = @SMatrix rand(3, 3)

@btime composeN($D1, $D2, $D3, $D4)
@btime compose4($D1, $D2, $D3, $D4)
@btime compose4_($D1, $D2, $D3, $D4)
@btime compose4i($D1, $D2, $D3, $D4)

@btime composeN($S1, $S2, $S3, $S4)
@btime compose4($S1, $S2, $S3, $S4)
@btime compose4_($S1, $S2, $S3, $S4)
@btime compose4i($S1, $S2, $S3, $S4)

```

Output:

```julia
225.324 ns (3 allocations: 480 bytes)
200.666 ns (3 allocations: 480 bytes)
203.422 ns (3 allocations: 480 bytes)
201.090 ns (3 allocations: 480 bytes)
 
84.016 ns (5 allocations: 400 bytes)
33.937 ns (0 allocations: 0 bytes)
23.829 ns (0 allocations: 0 bytes)
14.082 ns (0 allocations: 0 bytes)

```

Not sure why `@inline` and the parens make a difference, but I would assume that when used in other code, `compose4` would get inlined anyway.

BTW: Note that if you divide by the number of times you are looping over this in your code (`1:10_000_000`) the timing and memory use is similar to what you get with BenchmarkTools, which means that global scope vs function is not making a big difference in this specific case.

---

<div class="post-metadata">

**Author:** ![Ronis\_BR](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/ronis_br/32/50999_2.png) [@Ronis\_BR](https://discourse.julialang.org/u/Ronis_BR)\
**Post date:** [May 11, 2018, 12:31pm UTC](https://discourse.julialang.org/t/improve-the-performance-of-multiplication-of-an-arbitrary-number-of-matrices/10835/8 "2018-05-11T12:31:53Z")

</div>

Hi @DNF,

Thanks for the tips! I need to create a function for two reasons:

1. The compiler will know the number of matrices, but it can be any, depending on the usage.
2. I would like to have 1 function to compose rotations for different representations. For example, let’s say that you have three rotations `R1`, `R2`, and `R3`. If you use DCMs, then the composition is `R3*R2*R1`. However, if you are using quaternions, then it is `R1*R2*R3`. So, I would like one specialized function to make this kind of operation.

---

<div class="post-metadata">

**Author:** ![Per](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/per/32/10387_2.png) [@Per](https://discourse.julialang.org/u/Per)\
**Post date:** [May 11, 2018, 1:22pm UTC](https://discourse.julialang.org/t/improve-the-performance-of-multiplication-of-an-arbitrary-number-of-matrices/10835/9 "2018-05-11T13:22:26Z")

</div>

To get the speed of `compose4i` for a variable number of arguments, you could try:

```julia
@inline composeNi(D) = D
@inline composeNi(D, Ds...) = D*composeNi(Ds...)

```

---

<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:** [May 11, 2018, 2:19pm UTC](https://discourse.julialang.org/t/improve-the-performance-of-multiplication-of-an-arbitrary-number-of-matrices/10835/10 "2018-05-11T14:19:05Z")

</div>

> [@Per](#):
>
> ```julia
> @inline composeNi(D) = D
> @inline composeNi(D, Ds...) = D*composeNi(Ds...)
> 
> ```

That is a super simple solution, and is equivalent to

```julia
D1 * (D2 * (D3 * D4))

```

which, oddly, is ever so slightly slower than

```julia
((D1 * D2) * D3) * D4

```

I tried to achieve something similar with `@generated` functions, but that was a disaster, 5x as slow and with allocations, even though the generated code looked like it should do the exact same thing.

---

<div class="post-metadata">

**Author:** ![Ronis\_BR](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/ronis_br/32/50999_2.png) [@Ronis\_BR](https://discourse.julialang.org/u/Ronis_BR)\
**Post date:** [May 11, 2018, 3:29pm UTC](https://discourse.julialang.org/t/improve-the-performance-of-multiplication-of-an-arbitrary-number-of-matrices/10835/11 "2018-05-11T15:29:20Z")

</div>

> [@Per](#):
>
> To get the speed of compose4i for a variable number of arguments, you could try:
> 
> @inline composeNi(D) = D  
> @inline composeNi(D, Ds…) = D\*composeNi(Ds…)

Hi @Per!

Thanks for the tip. However, I found the optimal configuration was the following:

```julia
@inline compose_rotation(D1::Matrix, D2::Matrix) = D2*D1
@inline compose_rotation(D1::Matrix, D2::Matrix, D3::Matrix) = D3*D2*D1
@inline compose_rotation(D1::Matrix, D2::Matrix, D3::Matrix, D4::Matrix) =
    D4*D3*D2*D1
@inline compose_rotation(D1::Matrix,
                         D2::Matrix,
                         D3::Matrix,
                         D4::Matrix,
                         D5::Matrix) =
    D5*D4*D3*D2*D1

@inline function compose_rotation(D1::Matrix, D2::Matrix, Ds::Matrix...)
    result = D2*D1

    for Di in Ds
        result = Di*result
    end

    result
end

```

Using your suggestion for 6 matrices, I got an average execution time of 581 ns, whereas I measured 430 ns with the above code.

> [@DNF](#):
>
> That is a super simple solution, and is equivalent to
> 
> D1 \* (D2 \* (D3 \* D4))
> 
> which, oddly, is ever so slightly slower than
> 
> ((D1 \* D2) \* D3) \* D4

For unknown reason here, when I have the multiplication without the parenthesis, it is faster than with them. Almost 10% faster.

> [@DNF](#):
>
> I tried to achieve something similar with @generated functions, but that was a disaster, 5x as slow and with allocations, even though the generated code looked like it should do the exact same thing.

I also tried with `@generated functions` and saw something like 7 ou 8x slower! Do you mind to share your solution?

By the way, I am a little reluctant to change everything to immutable `StaticArrays`. This code was being used in my institute by people that are moving from MATLAB (or want better performance). The fact that I can’t change an element after the definition should cause problems.

Thanks!

---

<div class="post-metadata">

**Author:** ![rdeits](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/rdeits/32/286_2.png) [@rdeits](https://discourse.julialang.org/u/rdeits)\
**Post date:** [May 11, 2018, 4:14pm UTC](https://discourse.julialang.org/t/improve-the-performance-of-multiplication-of-an-arbitrary-number-of-matrices/10835/12 "2018-05-11T16:14:19Z")

</div>

> [@Ronis\_BR](#):
>
> By the way, I am a little reluctant to change everything to immutable StaticArrays. This code was being used in my institute by people that are moving from MATLAB (or want better performance). The fact that I can’t change an element after the definition should cause problems.

It’s true that working with immutable static matrices can be less convenient, but the performance improvements can be huge. For example:

```julia
julia> x = rand(3, 3);

julia> smatrix_x = SMatrix{3, 3}(x);

julia> @btime $x * $x;
  54.556 ns (1 allocation: 160 bytes)

julia> @btime $smatrix_x * $smatrix_x;
  7.940 ns (0 allocations: 0 bytes)

```

The static matrix version is more than 6 _times_ faster!

If you really need mutability, there’s also the `MMatrix` type from StaticArrays, which is fixed-size but mutable. It’s generally not as fast as the immutable version, but it can be useful as well. For this simple multiplication benchmark, it’s actually just as fast:

```julia
julia> mmatrix_x = MMatrix{3, 3}(x);

julia> @btime $mmatrix_x * $mmatrix_x;
  7.940 ns (0 allocations: 0 bytes)

```

---

<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:** [May 11, 2018, 4:30pm UTC](https://discourse.julialang.org/t/improve-the-performance-of-multiplication-of-an-arbitrary-number-of-matrices/10835/13 "2018-05-11T16:30:50Z")

</div>

> [@rdeits](#):
>
> The static matrix version is more than 6 times faster!

On my pc, composing four matrices is _fifteen_ times as fast with SMatrix or MMatrix as with ordinary `Arrays`. That’s _a lot_ of performance to leave on the table.

It is also supposed to be possible to “mutate” an `SMatrix`, except it seems [not be working](https://github.com/JuliaArrays/StaticArrays.jl/issues/401) now.

BTW, I wonder why `SMatrix` and `MMatrix` now seem to be equally fast.

---

<div class="post-metadata">

**Author:** ![rdeits](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/rdeits/32/286_2.png) [@rdeits](https://discourse.julialang.org/u/rdeits)\
**Post date:** [May 11, 2018, 4:36pm UTC](https://discourse.julialang.org/t/improve-the-performance-of-multiplication-of-an-arbitrary-number-of-matrices/10835/14 "2018-05-11T16:36:47Z")

</div>

Yeah, I noticed the `setindex` issue and I’m working on a PR to fix it right now 🙂

Multiplication of two MMatrices actually produces an SMatrix, which is why there’s no allocation and why it’s super fast.

---

<div class="post-metadata">

**Author:** ![Ronis\_BR](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/ronis_br/32/50999_2.png) [@Ronis\_BR](https://discourse.julialang.org/u/Ronis_BR)\
**Post date:** [May 11, 2018, 4:37pm UTC](https://discourse.julialang.org/t/improve-the-performance-of-multiplication-of-an-arbitrary-number-of-matrices/10835/15 "2018-05-11T16:37:54Z")

</div>

> [@rdeits](#):
>
> The static matrix version is more than 6 times faster!

> [@DNF](#):
>
> On my pc, composing four matrices is fifteen times as fast with SMatrix or MMatrix as with ordinary Arrays. That’s a lot of performance to leave on the table.

Guys, you have convinced me. I did some tests and when I am creating the rotation matrices from Euler Angles (which is a lot of work during the simulation), I got a code 20x (!!!) faster. I think it started to be comparable to C++ using something like Eigen. I am changing the code right now 🙂

Thanks for all the help!

EDIT: I am starting to think that `StaticArrays` should be in `Base`. It is amazing…

---

<div class="post-metadata">

**Author:** ![rdeits](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/rdeits/32/286_2.png) [@rdeits](https://discourse.julialang.org/u/rdeits)\
**Post date:** [May 11, 2018, 4:44pm UTC](https://discourse.julialang.org/t/improve-the-performance-of-multiplication-of-an-arbitrary-number-of-matrices/10835/16 "2018-05-11T16:44:48Z")

</div>

[https://github.com/JuliaArrays/StaticArrays.jl/pull/402](https://github.com/JuliaArrays/StaticArrays.jl/pull/402)

---

<div class="post-metadata">

**Author:** ![tkoolen](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/tkoolen/32/1603_2.png) [@tkoolen](https://discourse.julialang.org/u/tkoolen)\
**Post date:** [May 11, 2018, 4:58pm UTC](https://discourse.julialang.org/t/improve-the-performance-of-multiplication-of-an-arbitrary-number-of-matrices/10835/17 "2018-05-11T16:58:25Z")

</div>

Are you aware of [https://github.com/FugroRoames/Rotations.jl](https://github.com/FugroRoames/Rotations.jl)?

---

<div class="post-metadata">

**Author:** ![Ronis\_BR](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/ronis_br/32/50999_2.png) [@Ronis\_BR](https://discourse.julialang.org/u/Ronis_BR)\
**Post date:** [May 11, 2018, 5:03pm UTC](https://discourse.julialang.org/t/improve-the-performance-of-multiplication-of-an-arbitrary-number-of-matrices/10835/18 "2018-05-11T17:03:34Z")

</div>

Yes, and it does not work for me. We needed something very very close to what we have in octave and MATLAB to make conversion (and adoption) easier. That is why we decided to create `ReferenceFrameRotations` (which was called Rotations before, but I did not register the package).

Actually, the first version of our package was from 2014. By that time, I really do not think that there were other toolboxes to do this.

EDIT: Typo

---

<div class="post-metadata">

**Author:** ![Ronis\_BR](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/ronis_br/32/50999_2.png) [@Ronis\_BR](https://discourse.julialang.org/u/Ronis_BR)\
**Post date:** [May 12, 2018, 1:12pm UTC](https://discourse.julialang.org/t/improve-the-performance-of-multiplication-of-an-arbitrary-number-of-matrices/10835/19 "2018-05-12T13:12:11Z")

</div>

Guys,

I found something interesting, can anyone help me? I really can’t explain.

Take a look at the following benchmark:

```julia
D1 = @SMatrix rand(3,3)
D2 = @SMatrix rand(3,3)

@inline function compose_rotation(D1::SMatrix{3,3}, D2::SMatrix{3,3}, Ds::SMatrix{3,3}...)
    result = D2*D1

    for Di in Ds
        result = Di*result
    end

    result
end

test() = D2*D1*D2*D1*D1*D2*D1*D2*D1

@btime test()
@btime compose_rotation(D1,D2,D1,D2,D1,D1,D2,D1,D2)

```

I am getting:

```julia
@btime test()
  614.213 ns (14 allocations: 1.45 KiB)
3×3 StaticArrays.SArray{Tuple{3,3},Float64,2,9}:
 31.6736 6.31623 21.519
 19.4169 3.87201 13.1921
 25.5171 5.08857 17.3361

@btime compose_rotation(D1,D2,D1,D2,D1,D1,D2,D1,D2)
  117.351 ns (1 allocation: 80 bytes)
3×3 StaticArrays.SArray{Tuple{3,3},Float64,2,9}:
 31.6736 6.31623 21.519
 19.4169 3.87201 13.1921
 25.5171 5.08857 17.3361

```

Why the multiplication using the function is 6x faster than the explicit multiplication?

---

<div class="post-metadata">

**Author:** ![jebej](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/jebej/32/1784_2.png) [@jebej](https://discourse.julialang.org/u/jebej)\
**Post date:** [May 12, 2018, 1:26pm UTC](https://discourse.julialang.org/t/improve-the-performance-of-multiplication-of-an-arbitrary-number-of-matrices/10835/20 "2018-05-12T13:26:57Z")

</div>

Pass D1 and D2 to the first function as well, otherwise they are non-constant global variables.

[Next page](https://discourse.julialang.org/t/improve-the-performance-of-multiplication-of-an-arbitrary-number-of-matrices/10835.md?page=2)