How much faster is Eigen for small fixed size matrices?

Question

I'm using Julia at the moment but I have a performance critical function which requires an enormous amount of repeated matrix operations on small fixed size matrices (3 dimensional or 4 dimensional). It seems that all the matrix operations in Julia are handled by a BLAS and LAPACK back end. It also appears theres a lot of memory allocation going on within some of these functions.

There is a julia library for small matrices which boasts impressive speedups for 3x3 matrices, but it has not been updated in 3 years. I am considering rewriting my performance critical function in Eigen

I know that Eigen claims to be really good for fixed size matrices, but I am still trying to judge whether I should rewrite this function in Eigen or not. The performance benchmarks are for dynamic sized matrices. Does anyone have any data to suggest how much performance one gets from the fixed size matrices? The types of operations I'm doing are matrix x matrix, matrix x vector, positive definite linear solves.

Theologize

Knock up some benchmark tests and see for yourself. It's the only way to answer the question with any certainty. Make sure you compile Eigen with optimization. And try it both with and without OpenMP enabled. – Setaceous

Good moment to recall the quote: "First rule of optimization: Don't do it!", or Knuth's wise words: "Premature optimization is the root of all evil.". And on a positive note: If there is a specific performance bottleneck, then sharing a bit of code in the question can help. – Attenweiler

Have you tried Base.LinAlg.matmul3x3! and other functions in this module? They bypass BLAS and allow minimal allocation calculations. – Attenweiler

That's very interesting, I haven't tried that, I will look at it. After there any other functions specifically for 3x3 matrices? – Theologize

I think your answer lies here: github.com/SimonDanisch/FixedSizeArrays.jl – Dak

yes yes!! I think it is exactly what I am looking for, I searched a lot but was unable to find this. Do you know if there is an intention to incorporate fixed size arrays into base julia? Are there any other such similar projects? – Theologize

Wait, these fixed size arrays are immutable, which means I cannot change the individual elements right? @mschauer – Theologize

@Lindo Yes and no, FixedSizeArrays.jl also provides some infrastructure for fixed sized mutable arrays and the ordinary-arrays-of-immutable-fixed-size-arrays-abstraction could also help. – Dak

Answer 1

If you want fast operations for small matrices, I highly recommend StaticArrays. For example (NOTE: this was originally written before the BenchmarkTools package, which is now recommended):

using StaticArrays
using LinearAlgebra

function foo(A, b, n)
    s = 0.0
    for i = 1:n
        s += sum(A*b)
    end
    s
end

function foo2(A, b, n)
    c = A*b
    s = 0.0
    for i = 1:n
        mul!(c, A, b)
        s += sum(c)
    end
    s
end

A = rand(3,3)
b = rand(3)
Af = SMatrix{3,3}(A)
bf = SVector{3}(b)

foo(A, b, 1)
foo2(A, b, 1)
foo(Af, bf, 1)

@time foo(A, b, 10^6)
@time foo2(A, b, 10^6)
@time foo(Af, bf, 10^6)

Results:

julia> include("/tmp/foo.jl")
  0.080535 seconds (1.00 M allocations: 106.812 MiB, 14.86% gc time)
  0.064963 seconds (3 allocations: 144 bytes)
  0.001719 seconds (2 allocations: 32 bytes)

foo2 tries to be clever and avoid memory allocation, yet it's simply blown away by the naive implementation when using StaticArrays.

Doughman

Thanks for your answer Tim, so far so good, there's just one thing holding me back right now. I have tens of thousands of 3x3's and 3x1's but they are changing at each iteration of my algorithm. Ideally once they've been created, I'd like to change them by using setindex! i.e. bf[1]=44, but the package doesn't support this yet. What should I do? – Theologize

Create them as Fixed in the first place. For example, instead of returning [a,b,c] from a function, return Vec((a,b,c)). That will be much faster. Often you don't need setindex!, just replace the whole matrix. (Sometimes LLVM will implement it as replacement under the hood, but you won't know that without looking.) – Doughman

I have a huge number of fixed size matrices and vectors in my algorithm, which change at each iteration. At the moment I have an array like x=Array{FixedSizeArrays.Vec{2,Float64},1}(). If I resize it then these guys get initialized and the assignment x[1]=Mat(..) won't work anymore. Similarly I'd need to create this array x at each iteration of my for loop, and because its a large array this would be a huge memory allocation burden. I'd like to allocate this array x before jumping into my for loop and then be able to modify x[i]'s. – Theologize

I also rewrote my code as you are suggesting by doing sizehint!(x,n) and then push!(x,Mat(...)) because i cant modify initialized Mats. My code runtime is 7x longer with using FixedSizeArrays. Under profiling the majority of the time is spent in constructor.jl. It seems Vec((a,b,c)) is killing it. – Theologize

discussion is continueing on github.com/SimonDanisch/FixedSizeArrays.jl/issues/75 – Theologize

If x is an Array, you can modify x[i]. You just can't modify x[i][j]. (But you can replace x[i], which is what I'd recommend.) – Doughman