There are at least 4 separate implementations of quaternions in Julia. Some of them are <:AbstractVector, some are <:AbstractMatrix (
), 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
Iām sure the are specifics of each implementation that are good but I really like whatās in Rotations.jl.
Are there any advantages of using <: AbstractVector / <: AbstractMatrix over <: Number for quaternions?
Vectors are nice for differential equations
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.
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.
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.
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.
The kinematics of special relativity can be described using Quaternions with complex elements.
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}}).
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.
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](Quaternion - Wikipedia, 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).
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.
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.
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.
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.
Juliaās Number type acts a lot like a matrix anyway.
julia> 1'
1
julia> 2[1]
2
@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.