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.

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.