Spire is a Scala library that implements number types and numeric abstractions. With close to 1500 stars on github (as of writing), it seems to be the library around which the community has coalesced for those applications. However there are issues with the implementation of specific datastructures.

This article focuses on spire’s support for forward mode algorithmic differentiation, i.e. its implementation of Jets. There are three basic problems with the implementation - performance, mutability and numerical errors.

The source code referenced here can be found in the companion github repository. Also, a quick note on notation - this article uses residuals to refer to set of the first order sensitivities being computed.

Performance

A Jet is an interesting datastructure to implement because minimising latency requires maximising throughput of the operations over the residuals. In the general case, this means that the operations should use Single Instruction Multiple Data (SIMD) operations, a process often called vectorization. On the JVM, this is tricky in the best of times because programmers do not have explicit control over vectorization (unlike C++).

Spire’s Jet is parameterized over the type of the real, i.e. Double, spire.math.Complex etc. (that’s why the italicised real) . In principle, means that the same implementation can be used to compute sensitivies over multiple fields, such as real and complex; in practice this means that the the operations for the most common use case - Double - are painfully slow. They can be upto eight times slower than theoretical limits.

Performance is a deep issue that has no easy solution. The SimianQuant implementation of Jets runs within a few percent of the theoretical limit (while also providing a human readable API), but is built upon a massive infrastructure of metaprogramming and type algebra that would be impractical to develop for a simple library like spire.

Mutability

A spire Jet is a mutable object because it exposes a shallow copy of a mutable member variable (i.e. infinitesimal, the array of residuals). Moreover, this member is used for equality and hashcode computations which makes sets and maps unstable (in the general case). To illustrate:

val a = Jet(1.2) + Jet.h[Double](0)
val st = mutable.Set(a)
st.contains(a) // will be true

a.infinitesimal(1) = 11 // this is effectively a void setter
st.contains(a) // will be false

While not a dealbreaker in itself, this does break referential transparency, goes against the recommended Scala design patterns, and the emerging consensus on using immutable datastructures for atomic types.

A possible solution could be to add an extra function to return the residual at an index and return a deep copy if all of the residuals are required (this is the approach taken by the SimianQuant implementation).

Numerical Errors

Spire’s implementations of division and floor are incorrect. While division is a minor error, floor is an egregious one. A possible solution is to use a bridging library, like SimianQuant’s mathbridge, that fixes the numerical issues.

Division

Division is implemented as multiplication by the multiplicative inverse. While true algebraically, this is not true numerically and can lead to weird situations, such as when x / x is not scalar unity. To illustrate:

val a = Jet(0.9244501466787859, Array(0.21762703730059696))
val b = a / a
val c = Jet(1)
b == c // will be false

To be fair, for practical problems this issue is sufficiently minor that it can be ignored. That’s why mathbridge delegates to the default implementation instead of proving a custom one.

Floor

Since the floor of a function is piecewise constant, the derivative, where is is defined is zero. Spire implements floor as en elementwise floor of the residuals. This is obviously, and completely, wrong.

Harshad

In doing what we ought, we deserve no praise, for it is our duty