Taking Quaternions Seriously

[Orbots]

@JeffreySarnoff

struct Quaternion{S<:Number, T<:Number} <: Number
    s::S  # scalar value
    i::T  # imaginary i value
    j::T  # imaginary j value
    k::T  # imaginary k value
end

I think the idea is to not have a specific struct Quaternion but rather a set of methods.

In anycase. Assuming we are considering a reference implementation. Just wanted to point out that the type of the scalar part may be different than the type of the bivector part in some implementations.

If you really wanted a flag for normalization you could put it in S as well. Maybe a little too hacky?

[tuckermcclure]

@Ronis_BR

Sure! One very good book about reference systems, coordinate transformation, and also attitude dynamics is

I’ll second this book recommendation, and also add this paper by Shuster: A Survey of Attitude Representations (http://www.malcolmdshuster.com/Pub_1993h_J_Repsurv_scan.pdf). Glad he posted it on his site.

My own two cents about this topic:

• I seem to recall seeing some Julia package implement rotation operations by going to matrices and then back to the desired parameterization for the output. [Edit: At first, I thought I remembered which package, but I was wrong, so I removed the package name.] This loses more accuracy than dedicated operations for each parameterization. The difference would actually be the biggest source of error for, say, geographic surveying. The important part here doesn’t have anything to do with any particular package, but is simply that there are many, many methods for working with this stuff, according to different needs.
• Both aircraft and spacecraft folks tend to work with frame rotations, whereas (in my limited experience) robotics folks seem to like vector rotations (opposite conventions).
• Aircraft folks tend to put the scalar first, and spacecraft folks tend to put it last.
• As the paper points out, we might call this parameterization “Euler-Rodrigues parameters” to keep it separate from the more general idea of a quaternion. Further, typical math conventions for quaternions don’t actually tend to operate with the same order of operations that are desired for frame rotations, so functions specifically for an EulerRodriguesParameters type might be easier to use for working with frame rotations, and they could be totally separate from some more general Quaternion type. That might feel dirty since one’s a subset of the other in some sense, but on the other hand, dedicated implementations would almost certainly handle rotations more accurately and quickly.

Because there are so many competing conventions and methods for implementing the operations, I doubt we’ll unify this in a way that would satisfy broad audiences. Any time I see a “quaternion” on an interface, I tread carefully, double-check the data sheet, and then implement a bunch of tests to try to see what I actually have to do to work with it. I literally found a mixup of conventions in someone else’s code for this stuff yesterday.

[scheinerman]

I recently released SimpleQuaternions mostly as an exercise for myself, but thought it might be useful to others. Some thoughts:

• I see Quaternions as extensions of the Complex numbers. For that reason, they should be fully interoperable with Complex numbers and Quaternions should be a subtype of Number.
• Complex numbers use the symbol im for the square root of negative one. That is, a+bi is entered as a + b*im. My implementation exports jm and km for the quaternion units j and k, so the quaternion a+bi+cj+dk can be entered as a + b*im + c*jm + d*km.

[Ronis_BR]

@tuckermcclure

• I seem to recall seeing Rotations.jl implement rotation operations by going to matrices and then back to the desired parameterization for the output. (Please let me know if I’m mistaking this with another package.) This loses more accuracy than dedicated operations for each parameterization. The difference would actually be the biggest source of error for, say, geographic surveying. The important part here doesn’t have anything to do with Rotations.jl, but is simply that there are many, many methods for working with this stuff, according to different needs.
• Both aircraft and spacecraft folks tend to work with frame rotations, whereas (in my limited experience) robotics folks seem to like vector rotations (opposite conventions).

I completely agree. Those are the two reasons why I kept maintaining ReferenceFrameRotations.jl after Rotations.jl.

[bhawkins]

As others have mentioned, lots of varying conventions make this a difficult task. For example, there’s not consensus on whether i*j is +k or -k. I recently came across this paper (Sommer et al) that provides a nice summary.

[chakravala]

@bhawkins

For example, there’s not consensus on whether i*j is +k or -k.

Generally, the convention you want is i*j*k == -1.

[moble]

I agree strongly with @Orbots re: isnorm. In fact, in my experience, it’s best for numerical applications to ignore the possibility of a normalized quaternion; when normalization is relevant, it’s better to include the act of normalization explicitly. For example, when rotating a vector v by some quaternion R, it’s better to use R v R^{-1} than the usual R v \bar{R}. As you build up the computation of R, you’re freed from the irrelevant overhead of trying to maintain perfect normalization. Composing rotations by multiplying quaternions goes through just as before. You don’t get stiff ODEs if you need to integrate in time. And so on. (Though R v \bar{R} can still be useful for transformations that shouldn’t maintain the norm of the vector.)

Besides, if you really want a normalized type, the type system is where you should look. For example, we might have log(q::Quaternion) return Quaternion, implementing the general logarithm in the space of nonzero quaternions, while log(q::UnitQuaternion) returns a Bivector.

Of course, that line of reasoning leads rapidly to the deeper point, which is that we all need to start Taking Geometric Algebra Seriously. It would be best to design any quaternion package with an eye toward the future, so that types can be as interoperable and conventions as consistent as possible.

[moble]

@scheinerman

I recently released SimpleQuaternions …

Nice!

I’ll echo @chakravala in pointing out that the choice of plane when identifying complex numbers as a sub-algebra of the quaternions is arbitrary. Still, it’s entirely reasonable to just go ahead and make such a choice. Curiously, the conventional choice in the literature (at least the dark corners of it that I frequent) skews toward identifying the unit imaginary i \in \mathbb{C} with the generator of rotations about the z axis: k \in \mathbb{H}, or km in your code. Specifically, this pops up repeatedly in the 2-spinor formalism and the theory of (possibly spin-weighted) spherical harmonics. To put it another way, \mathbb{C} is regarded as the even subalgebra of the geometric algebra of the x-y plane, and i \in \mathbb{C} maps to x \wedge y.

[chakravala]

@moble

if you really want a normalized type, the type system is where you should look.

Indeed, I already know I would make this a parametric type parameter and not a struct value, but this was very obvious to me that I didn’t bother mentioning it yet.

@moble

we all need to start Taking Geometric Algebra Seriously.

Well, I have been doing that for a while now, but I know there seem to be many different reasons why people avoid taking my contribution seriously as of yet

@moble

the choice of plane when identifying complex numbers as a sub-algebra of the quaternions is arbitrary. Still, it’s entirely reasonable to just go ahead and make such a choice.

I’m not saying you can’t make a choice, but prefer leaving it as a choice for the user. My primary goal has been geometric algebra research, and I believe that geometric algebra researchers should be making those choices personally instead of locking the choices.

Also, by not making a choice, it keeps automatic conversions simpler, and allows the usage of complex numbers as coefficients instead of real numbers. If complex numbers automatically promote themselves into a specific algebra and plane, then it would start getting messy to handle complex coefficients. I don’t generally recommend complex coefficients, since real coefficients of a higher dimensional algebra would be the most natural choice, yet I want to keep complex coefficients available as another choice also (since they can also be useful)

[tkoolen]

Continuing the discussion from Taking Quaternions Seriously:

@tuckermcclure

I seem to recall seeing Rotations.jl implement rotation operations by going to matrices and then back to the desired parameterization for the output. (Please let me know if I’m mistaking this with another package.)

This is simply not true.

@JeffreySarnoff

struct Quaternion{T<:Number} <: Number
    s::T  # scalar value
    i::T  # imaginary i value
    j::T  # imaginary j value
    k::T  # imaginary k value
end

It may be performant to include an isnorm Bool field to flag realizations that have been normalized.

This is almost exactly what Quaternions.jl (old, practically unmaintained) did:

JuliaGeometry/Quaternions.jl/blob/01d45b4f345fb91945211d3c6387da4b5934e7d7/src/Quaternion.jl#L1-L7

To represent rotations, we’re really only interested unit quaternions, so in Rotations.jl the UnitQuaternion type basically asserts isnorm at compile time and tries to avoid renormalizing whenever that makes sense, avoiding performance costs but placing some burden on the user to ensure that numerical inaccuracies after several operations are kept in check. A type parameter that captures whether or not the quaternion is normalized (as mentioned) is also fine, but in Rotations.jl there has not really been a need for unnormalized quaternions.

@moble

For example, when rotating a vector vv by some quaternion RR , it’s better to use RvR−1R v R^{-1} than the usual Rv¯RR v \bar{R} .

This would have very severe performance consequences that are unacceptable in a lot of applications.

I think the fact that different users have different requirements in terms of performance, accuracy, memory layout, mathematical underpinnings, focus on rotation / interop with other rotation representations vs. support for unnormalized quaternions, as well as package dependencies means that it’s OK for there to be more than one implementation in the Julia package ecosystem. Of course, do consider contributing to an existing package before writing another implementation.

[moble]

@chakravala Sorry. I think there’s been some miscommunication. My first comment was really replying to JeffreySarnoff by way of Orbots; I didn’t mean to say anything about parametric types, or suggest that you’re not taking GA seriously. And in my second comment, I only meant to give credit to you for introducing that idea into this conversation. I agree with the goal of leaving choices to users as much as possible; I just meant to point out that such choices have actually been made before, and it’s good to be aware of existing conventions.

FWIW, I take your contribution quite seriously. In fact, Grassmann.jl is what first sparked my interest in Julia long ago because it’s such a nice piece of code, and I only regret that I don’t get more opportunity to use both Julia and your package.

[moble]

@tkoolen

@moble

For example, when rotating a vector v by some quaternion R , it’s better to use RvR^{−1} than the usual R v \bar{R} .

This would have very severe performance consequences that are unacceptable in a lot of applications.

I think it’s the opposite. Rather than repeatedly dividing by sqrt of the sum of squares of the components, as you manipulate the quaternions themselves, you just divide by the sum of squares (no sqrt!) once in this final operation. Of course, if you’re doing multiple rotations with the same quaternion, you should probably be converting to a matrix anyway, and again you just do this sum of squares once in that conversion.

[chakravala]

No worries, I recognized your user name, I check out all the people who interact with my github and immediate recognize you. Also, I don’t care whether you take it serious or not, I’m only saying that I understand there are various reasons why geometric algebra in general might not be take seriously by other people.

[tuckermcclure]

@tkoolen

This is simply not true.

My bad, @tkoolen! I’m truly sorry I got that wrong. I’m sure I saw the pattern I mentioned in a Julia package years ago and pulled it down and tested it at the time, and Rotations.jl is the only package I remember looking at. I tried to check by looking at GitHub - JuliaGeometry/Rotations.jl: Julia implementations for different rotation parameterizations, and it said, “At their heart, each rotation parameterization is a 3×3 unitary (orthogonal) matrix,” which caused me to believe I had remembered correctly. Looking at the source, however, I see I’m flat wrong. I removed the call-out in my original post to prevent confusing anyone else about Rotations.jl. No offense was intended! It was just intended to point out that there are multiple ways to go about working with these things, each with their own advantages.

[tkoolen]

@tuckermcclure

My bad, @tkoolen! I’m truly sorry I got that wrong.

No worries.

@moble

I think it’s the opposite. Rather than repeatedly dividing by sqrt of the sum of squares of the components, as you manipulate the quaternions themselves, you just divide by the sum of squares (no sqrt !) once in this final operation. Of course, if you’re doing multiple rotations with the same quaternion, you should probably be converting to a matrix anyway, and again you just do this sum of squares once in that conversion.

I’m not comparing normalization upon construction of R and then using R v \bar{R} to skipping normalization of R upon construction and using R v R^{-1}; I’m comparing both of those options to skipping normalization upon construction when performing operations that should preserve unit norm and then using R v \bar{R}. Even the computation of the squared norm, and especially the division needed for R v R^{-1} come at a performance cost that some users are not willing to pay.

[moble]

@tkoolen

I’m not comparing normalization upon construction of R and then using R v \bar{R} to skipping normalization of R upon construction and using R v R^{-1}; I’m comparing both of those options to skipping normalization upon construction when performing operations that should preserve unit norm and then using R v \bar{R}. Even the computation of the squared norm, and especially the division needed for R v R^{-1} come at a performance cost that some users are not willing to pay.

Oh, I see. Fair point. In my work, accuracy is paramount and quaternion operations are far from being the bottleneck, so I have a blind spot for these sorts of needs.

So we’re actually saying almost the same thing: operations on quaternions should usually ignore the norm. I’m just adding the point that this even works at the analytical level (including for things like slerp and squad) as long as you don’t hit 0, and use R v R^{-1} when applying the final rotation.

[JeffreySarnoff]

a first attempt to gather from all these most helpful contributions

There is a groundswell of support for a well-specified API to which different implementations can conform. Here is a first take on the least required. For your comments, improvements:

four conventions for brevity,

• const Q = Quaternion
• quaternion elements are as named fields: s, i, j, k
• T<:Number, so T could be Number, Real or Float64 …
• where {T} is omitted

The API autodefines sections are the default (fallback) definitions to be used where the API client does not provide specialized implementations.

imports

import Base: abs, abs2, conj, inv, +, -, *, /, \
import LinearAlgebra: normalize

Constructor

Q(s::T, i::T=zero(T), j::T=zero(T), k::T=zero(T))

API autodefines

Q(; s=0, i=0, j=0, k=0) = Q(promote(s,i,j,k)...)
Q(sijk::NTuple{4,T}) = Q(sijk...)
Q(s::T, ijk::NTuple{3,T}) = Q(s,ijk...)
Q(s, ijk::NTuple{3,T}) = Q(promote(s,ijk...))

Selectors

scalar(q::Q) = q.s
bivector(q::Q) = [q.i, q.j, q.k]
fieldvalues(q::Q) = (q.s, q.i, q.j, q.k) 
bivecvalues(q::Q) = (q.i, q.j, q.k)

API autodefines

abs2(q::Q) = sum(fieldvalues(q).^2)
abs(q::Q) = sqrt(abs2(q))
pure(q::Q{T}) = Q(zero{T}, bivecvalues(q))
conj(q::Q) = Q(scalar(q), (-).(bivecvalues(q)))

normalize(q::Q) = Q(fieldvalues(q) ./ abs(q))
inv(q::Q) = Q(fieldvalues(q) ./ abs2(q))

+(p::Q{T}, a::Real) = Q((+)(promote(scalar(p),a)...), bivecvalues(q))
+(a::Real, p::Q{T}) = Q((+)(promote(a,scalar(p))...), bivecvalues(q))
+(p::Q{T}, q::Q{T}) = Q((+)(fieldvalues(p) .+ fieldvalues(q))

-(p::Q{T}, a::Real) = Q((-)(promote(scalar(p),a)...), bivecvalues(q))
-(a::Real, p::Q{T}) = Q((+)(promote(a,scalar(p))...), bivecvalues(q))
-(p::Q{T}, q::Q{T}) = Q((+)(fieldvalues(p) .+ fieldvalues(q))

*(p::Q{T}, a::Real) = Q(fieldvalues(p) .* a)
*(a::Real, p::Q{T}) = Q(a .* fieldvalues(p))
*(p::Q{T}, q::Q{T}) =
    Q(p.s*q.s - p.i*q.i - p.j*q.j - p.k*q.k,
      p.s*q.i + p.i*q.s + p.j*q.k - p.k*q.j,
      p.s*q.j - p.i*q.k + p.j*q.s + p.k*q.i,
      p.s*q.k + p.i*q.j - p.j*q.i + p.k*q.s)

/(p::Q{T}, a::Real) = Q(fieldvalues(p) ./ a)
/(a::Real, p::Q{T}) = a * inv(p)
/(p::Q{T}, q::Q{T}) = (p * conj(q)) / abs2(q)

\(p::Q{T}, a::Real) = (a * p) / a^2
\(a::Real, p::Q{T}) = (conj(p) * a) / abs2(p)
\(p::Q{T}, q::Q{T}) = (conj(q) * p) / abs2(q)

To consider

• exp, log, trig
• rotation related ops
• how do we want to show Quaternions
• including definitional tests with the API

[chakravala]

@JeffreySarnoff

fieldvalues(q::Q) = (q.s, q.i, q.j, q.k)

in my API, this would be obtained with the already existing values method instead of fieldvalues.

in newest version of AbstractTensors there is a simplified variant of StaticVector implemented

Help with simplified StaticArrays implementation

Poblem: StaticArrays.jl is too large of a dependency and has too many extra features. To simplify my dependencies, I need a simplified variant of SVector, MVector, and SizedVector. Solution: In AbstractTensors.jl they are replaced with Values, Variables, and FixedVector. julia> using AbstractTensors # loads 10x faster than StaticArrays julia> a = Values(1,2,3) 3-element Values{3,Int64} with indices SOneTo(3): 1 2 3 My goal is to have an implementation that is robust as in StaticArrays b…

which loads much faster than StaticArrays and provides all the basic features needed for Values.

AbstractTensors.Values # simplified SVector variant
AbstractTensors.values # alternative to fieldvalues

my proposal is using simple Values and values for internal values of TensorAlgebra

[Orbots]

@JeffreySarnoff

pure(q::Q{T}) = Q(zero{T}, bivecvalues(q))

Usually this would be called an imaginary quaternion, I believe . Implementing Base.imag in place of pure would be good.

Any desire to support split quaternions here? They are very like quaternions but with a different metric. Split-quaternion - Wikipedia

Far less common that quaternions, but I thought I’d throw that out there.

In that case we’d want to support arbitrary ( or just 2 ) metrics.

[chakravala]

What about eltype of a Quaternion?

By default eltype(x) = eltype(typeof(x)) if they are a subtype of Number instead of AbstractArray.

At the moment, I am going with this default definition, since these are treated like Number and I am unsure if a change would introduce interoperability problems.

Instead, I currently use valuetype as an alternative.