Taking Quaternions Seriously

There are at least 4 separate implementations of quaternions in Julia. Some of them are <:AbstractVector, some are <:AbstractMatrix (:neutral_face:), some are <:Number, some follow sijk convention, some follow ijks convention, etc. It would be no problem, but I’m constantly switching between them because different packages that I need to use together use different ones. - @jonniedie

It would be good to unify these with whichever choice is best all around. - @StefanKarpinski

Absatively

17 Likes

I’m sure the are specifics of each implementation that are good but I really like what’s in Rotations.jl.

3 Likes

Are there any advantages of using <: AbstractVector / <: AbstractMatrix over <: Number for quaternions?

3 Likes

Vectors are nice for differential equations

1 Like

My package ReferenceFrameRotations.jl also has quaternion algebra. The difference from Rotations.jl is that we treat all the rotations as passive (here we rotate reference systems, not vectors, which is what we usually do when working with satellite dynamics).

It is a structure that is also an AbstractVector. This is good because when you are propagating attitude, you compute a product between a matrix and a vector with the quaternion components. In this case, it is good that a quaternion can behave as a 4x1 vector.

3 Likes

Treating it as a number type can make it easily portable to new applications. A lot of problems can just run with a Vector{<:Number}, so you can expect a lot of Julia to just work. That said, the array of structs formulation might not be efficient for some calculations.

3 Likes

Numbers should be fine too (or any abstract vector space over the reals, for that matter — numbers are a vector space, and so are matrices … not all vector spaces are described by AbstractVector!).

If you don’t define a Quaternion as a Number type, many of the LinearAlgebra functions may not work for matrices of quaternions. And broadcasting should clearly treat a quaternion as a scalar, not as an AbstractVector (container). And so forth…

You’ll want to unroll/inline such a small matrix-vector product anyway (similar to StaticArrays), and hence you’ll probably want to define your own quaternion-operator type (similar to StaticMatrix) rather than using AbstractMatrix.

You could get the best of both worlds by defining it as

struct Quaternion{T<:Real} <: Number
    v::SVector{4,T}
end

and then you can define efficient operations on q::Quaternion by calling methods from StaticArrays.jl on q.v (and using q.v any time you need to iterate over components). For example

Base.+(q1::Quaternion, q2::Quaternion) = Quaternion(q1.v + q2.v)

and similarly defining your quaternion-operator type as an object encapsulating an SMatrix.

13 Likes

FWIW, I think someone did this once and reported back that the boundary value problem solvers just worked, even with ForwardDiff involved, so this style of implementation is what I’d suggest too since it’s fairly robust.

7 Likes

The kinematics of special relativity can be described using Quaternions with complex elements.

2 Likes

A unit complex number is a parameterization for SO(2) loosely the same way that unit quaternions [double] parameterizing SO(3).

But Complex is not defined to be a <:AbstractArray (not a AbstractVector and not a rotation matrix).

I don’t think the underlying “field” should be <:Real. Along the same lines, I don’t see a good reason for Complex{<:Real} either. For example, a dual quaternion could be represented by Quaternion{Dual{Float64}} (or isomorphically Dual{Quaternion{Float64}}).

2 Likes

Not that I disagree with your point, but note that in most implementations Dual <: Real (e.g. ForwardDiff.jl’s) for convenience. In my experience, Real in Julia often really means Abelian or similar.

I’ll just mention that my package Grassmann.jl supports quaternions and unifies them with vectors and other types of tensors. This package is my attempt ar unifying quaternions with linear algebra and tensor algebra. This is going to be my default algebra for everything, and I will rarely if ever use other ones.

AbstractTensors.TensorAlgebra <: Number

They are all subtype of Number, and it also supports Minkowski space and dual quaternions and even higher dimensions than just 3D.

I’ve been using it to develop a mesh interface for partial differential equations, similar to GeometryBasics but more flexible and generalized.

It also has faster performance than SMatrix for certain matrix operations and allows the construction of matrix like tensor elements too.

It’s still a work in progress and I plan on making more optimizations for quaternions specifically, but there are many possible things to focus on, so I first focused on the most general unification features and then specialized specific features as I go along.

16 Likes

Yeah, I think that arises in part from wanting to operate with other types and methods that force <:Real. I personally have yet to understand the benefit of this structure. Are there examples of two algorithms, one which assumes commutativity and one that doesn’t, and dispatch chooses the right one? I certainly could imagine that situation, and think it would be better served by a more general mechanism like using traits so that e.g. LinearAlgebra.Diagonal could be designated Abelian.

I could see even <:Number causing problems, since square matrices behave a lot like numbers. Indeed, complex, dual, [quaternion](https://en.wikipedia.org/wiki/Quaternion#Matrix_representations, split-complex numbers all have matrix representations. Perhaps the answer there is to wrap small matrices into a type that is <:Number.

This is probably outside the scope of “taking quaternions seriously”, but I think it would be very cool and useful to encode those matrix representations in Julia. In sloppy notation, matrep(Complex{Dual}) === kron(matrep(Complex), matrep(Dual)). It’s useful because you can fall back on many functions that are defined for square matrices (sqrt, exp, …). You also can then choose to use some operations, like \, which might do pivoting for better numerical accuracy (at the cost of performance).

2 Likes

Grassmann.jl already does that.

Again, this is something I am working on in Grassmann.jl but am still working on finalizing.

Most of the issues discussed here can be addressed with my packages, which have started with the most general possible foundation and are now working towards optimizing specialized cases.

1 Like

Can you give an example of, say, Complex{Float64} interoperating with a complex number as Grassmann.jl represents it?

Grassmann.jl is not designed to be automatically interoperable with the Julia’s Complex type currently. This is because in Grassmann.jl you can have a higher dimensional manifold with multiple imaginary units, much more than just 3 imaginary units (quaternions). So it would be ambiguous in general to know which imaginary basis unit it is supposed to convert to. Therfore, the user must make this choice personally.

Complex numbers are flat land, and live in a single plane, but in Grassmann.jl the choice of plane is ambiguous from a Julia Complex type, so you must personally choose how to interpret it.

Also, Julia’s Complex can be used as a coefficient type in the Grassmann algebra (as with Real), so automatic conversion would cause confusion.

The Grassmann algebra is designed to resolve this higher dimensional ambiguity by relying on combinatorics and “quantum computing.” But this requires an additional layer of context, which the Julia’s Complex lacks, so you personally choose how to interpret it in that case.

1 Like

Use the same convention as for complex numbers here. Complex numbers could be described as vectors, but that’s not useful in practice, hence they are numbers. Same with quaternions.

4 Likes

For what it’s worth, I think that while there may be some tantalizing historical reasons to think of Quaternions as being <: Vector (e.g. the word vector even comes from quaternions), it’s the least sensible of the supertype options Number, AbstractVector, AbstractMatrix.

Unlike vectors, a quaternion may be sensibly multiplied by another quaternion, and you can also divide by quaternions.

Algebraically, one can show that the reason people keep on finding a deep connection between quaternions and rotations is because the quaternions are actually a set of unit bivectors – oriented planes, rather than a set of unit vectors (oriented lines), further suggesting that Number or AbstractMatrix make more sense.

I’d definitely lean towards <: Number, but there are also distinct advantages to treating them as square matrices.

8 Likes

Julia’s Number type acts a lot like a matrix anyway.

julia> 1'
1

julia> 2[1]
2
7 Likes

@Ronis_BR do you have some reference that you could share? We are recently considering implementing some coordinate reference system rotations in GeoStats.jl and it would be really useful to understand the benefits of this approach.

1 Like