Taking Quaternions Seriously

[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 https://github.com/JuliaGeometry/Rotations.jl, 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. https://en.wikipedia.org/wiki/Split-quaternion

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.

[JeffreySarnoff]

Given q = one(Quaternion{Float64}), it seems to me that
eltype(q) === Float64 makes the most sense. Given the second half of the help text:

Determine the type of the elements generated by iterating a collection of the given type. For dictionary types, this will be a Pair{KeyType,ValType}. The definition eltype(x) = eltype(typeof(x)) is provided for convenience so that instances can be passed instead of types. However the form that accepts a type argument should be defined for new types.

Base.eltype(::Type{Quaternion{T}) where {T<:Number} = T makes sense imo. Also Base.eltype(::Type{QuaternionAsVector{T}}) = T.

[JeffreySarnoff]

@Orbots

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

My understanding has been that a pure quaternion is a quaternion. Implementing real(q) and imag(q) makes sense. Would you have imag(q) return a vector or a tuple? Do you agree that real(q) should return a scalar?

[chakravala]

In Grassmann you would get a bivector with imag of a quaternion, so a 2-form (neither vector nor tuple).

Of course, real would return a scalar, although that scalar might be wrapped in a 0-form, to preserve vector space the scalar exists in.

[Orbots]

@JeffreySarnoff

Would you have imag(q) return a vector or a tuple? Do you agree that real(q) should return a scalar?

Good questions. For the application I have in mind I actually would want both to end up as quaternions, since I’m feeding everything into a sparse linear system and the uniformity of dimensions helps book-keeping. But that’s related to my specific use case/implementation. Constructors with 1 and 3 airty could be used to easily convert.

[goretkin]

eltype(Quaternion{Float64}) === Float64 is only okay if and only if collect(q::Quaternion{Float64}) == [q.w, q.x, q.y, q.z] (modulo the field order and field access), and I believe also iff 1 .+ q == [1 + q.w, 1+q.x, 1+q.y, 1+q.z] (to say something about broadcasting).

I don’t think it should be done, and that the Quaternion type (or any type that implements the quaternion interface) should behave like Complex and all the types defined in ColorTypes.jl. That is, they should be treated like scalars for the purpose of iteration, indexing, and broadcasting.

Note that Base defines iterate(x::Number) = (x, nothing) and getindex(x::Number) = x (which is questionable, but I think outside the scope of this conversation), and I think it would basically be incorrect for a <:Number to define anything else.

[chakravala]

@goretkin

That is, they should be treated like scalars for the purpose of iteration, indexing, and broadcasting.

Exactly in agreement with my experiences. This is why I have valuetype instead of eltype.

However, I am willing to change it, if a convincing argument can be made for something else.

The situation should be similar to Complex.

[JeffreySarnoff]

I have become convinced (thank you @goretkin). Whether we call it valuetype or innertype or ?? … it should be distinct from eltype.

[elsuizo]

Hi, in some literature they define the Quaternion in term of a scalar part and a 3D vector

using StaticArrays
using LinearAlgebra

import Base.*

struct Quaternion{T<:Real} <: Number
   q₀::T
   q::SVector{3, T}
end

(*)(a::Quaternion{T}, b::Quaternion{T}) where T<:Real = Quaternion(a.q₀ * b.q₀ - dot(a.q, b.q), a.q₀ * b.q + b.q₀ * a.q + cross(a.q, b.q))

and they define the product in terms of dot and cross products

[yashi]

So, does anyone working on this?

I’ve been trying to use Makie.Quaternion with ReferenceFrameRotations.Quaternion and I am beaten.