How to optimize the speed of a numerical library in Haskell

Question

I have released a small numerical library for solving delay differential equations: http://github.com/masterdezign/dde

The main technical constraints:

1. Make the library flexible with respect to the number of dynamical variables (x(t), y(t),...), i.e. a State of the dynamical system for a given moment of time.
2. Easily integrate with libraries that exploit Data.Vector.Storable (e.g. hmatrix). Therefore, as input/output I extensively employ Data.Vector.Storable.

As a result, unlike in this solution: How do I optimize numerical integration performance in Haskell (with example), newtype State = State { _state :: V.Vector Double } is used instead of data State = State {-# UNPACK #-} !Double {-# UNPACK #-} !Double. However, the library runs twice slower now.

The question: is there any way to bring both the speed of data State = State {-# UNPACK #-}... and the flexibility of newtype State = State { _state :: V.Vector Double } for unspecified number of variables? Shall I consider Template Haskell to create data UNPACK-like structures during the compilation time?

Answer 1

I wouldn't use any particular vector implementation. Variable-length types like Data.Vector are a bad choice not only because the extra length information is a considerable overhead when the dimension of the space is low, also because you lose any type-system guarantee that the dimensions match.

Instead, you should make everything parametric on the choice of vector space. I.e, you effectively make the dimension a compile-time variable, plus you allow vector types with some meaningfully-named subvariables.

import Data.VectorSpace
import Data.AdditiveGroup

newtype Stepper1 state = Stepper1 {
  _stepper
    ::  Double
    -> RHS state   -- parameterised in a similar way
    -> state
    -> (Double, Double)
    -> (Double, Double)
    -> state
  }

rk₄ :: VectorSpace v => Stepper1 v
rk₄ = Stepper1 _rk4
 where _rk4 hStep rhs' y₀ ... = y₀ ^+^ (h/6)*^(k₁ ^+^ 2*^k₂ ^+^ 2*^k₃ ^+^ k₄)
         where k₁ = rhs' (y₀, ...)
               k₂ = rhs' (y₀ ^+^ (h/2)*^k₁, ...)
               k₃ = rhs' (y₀ ^+^ (h/2)*^k₂, ...)
               k₄ = rhs' (y₀ ^+^ h*^k₃, ...)

Then, what particular implementation it will be can be chosen by the user. For a 2-dimensional vector, a standard choice is V2 from the linear package; it is re-exported with VectorSpace instances from free-vector-spaces. But for testing, you could also just use plain old tuples, which have a VectorSpace instance right in vector-space. Of course, wrapping around the types from hmatrix would also be possible, but this would really not be good for performance – better just convert the final results, if necessary.

To get the best performance, you may need to use some {-# INLINE #-} pragmas. Bang patterns OTOH don't usually bring much performance advantage – most important is that the types are strict, and unboxed. It's certainly no use to pre-emptively put a bang before every variable definition – those won't have any effect anyway, because a CAF is only evaluated when it's used.

I would be excited to hear about the performance you get in the end! If it's notably worse than your original State {-# UNPACK #-} !Double {-# UNPACK #-} !Double, that's something we should investigate.