Treating antipodal quaternions as equal breaks some operations

[bzinberg]

Currently, antipodal quaternions are treated as equal, i.e., q == -q. I can see the motivation for this -- they induce the same rotation -- but this leads to several issues:

• AFAIK, accessing the quaternion components is part of the public API. So it's confusing that two UnitQuaternions can be == but have different values of the field: q == -q but q.w != (-q).w.

• The exponential map is broken: exp(Rotations.UnitQuaternion(0, π, 0, 0)) gives incorrect result #128

• zero(UnitQuaternion) is not a valid unit quaternion. As described in Skew-symmetric matrices as generators for rotations #30 (comment)), there is an undocumented semantics of "unnormalized unit quaternion" which is used for some operations such as exp and log, but is not ultimately coherent. Better would be to have a separate type for infinitesimal unit quaternions, similar to what's being done in Skew-symmetric matrices as generators for rotations #30. That type has a zero; UnitQuaternion should not have a zero.

[hyrodium]

So it's confusing that two UnitQuaternions can be == but have different values of the field: q == -q but q.w != (-q).w.

I think the == operator here is for "equals as a matrix" because UnitQuaternion <: Rotation <: StaticMatrix{3,3,<:Real}. So, it's ok to have different values of the field: q == -q but q.w != (-q).w.

As a related example, the following code satisfies a == b && typeof(a) == typeof(b) && f(a) ≠ f(b):

julia> a, b = 0.0, -0.0
(0.0, -0.0)

julia> typeof(a), typeof(b)
(Float64, Float64)

julia> 1/a ≠ 1/b
true

If we'd like to deal with the problem correctly, I think it's better to add types for UnitQuaternion, MRP, (and etc.) and the types should not be subtype of StaticMatrix{3,3,<:Real}. (Maybe, it's ok to be a subtype of StaticMatrix{2,2,<:Complex} because unit quaternions and SU(2) are isomorphism.)

[bzinberg]

If we'd like to deal with the problem correctly, I think it's better to add types for UnitQuaternion, MRP, (and etc.) and the types should not be subtype of StaticMatrix{3,3,<:Real}.

I agree, I think subtyping is not appropriate here.

It seems sensible to me for RotMatrix{3, T} to be a subtype of StaticMatrix{3, 3, T} -- a RotMatrix "is" a matrix -- but the other Rotation types seem to fit the description "can be converted to a RotMatrix," rather than "is a matrix." WDYT?

[hyrodium]

but the other Rotation types seem to fit the description "can be converted to a RotMatrix," rather than "is a matrix." WDYT?

Indeed, but I think the relation between RotMatrix{3} <: Rotation and UnitQuaterion <: Rotation (and etc.) is similar to the relation between Vector <: AbstractVector and StepRange <: AbstractVector.
They have different internal implementation, and the latter uses less parameters (memory) and is faster.

However, I don't think it is necessary to go through RotMatrix when converting from RotX (and etc.) to UnitQuaternion.
For example, if we define the conversion like this, then UnitQuaternion(RotX(2π)) can keep its information and continuity.

function UnitQuaternion(r::RotX)
    θ = r.theta
    UnitQuaternion(cos(θ/2),sin(θ/2),0,0)
end

[bzinberg]

I'm totally on board with implementing efficient conversions where possible, such as UnitQuaternion(::RotX) as you showed above. To me, that is independent of the inheritance question. Inheritance decisions among "public" (in Julia parlance I guess we'd say exported, though I don't think of those two as being the same) types i.m.o. should be based solely on interfaces, not on implementations. (It might make sense to rely on inherited implementations internally, among "private" types, but i.m.o. even this is not always worth it.)

So, I think the decision to have StepRange be a subtype of AbstractVector is primarily one of whether it fits the interface of what an AbstractArray is supposed to support. Julia generally doesn't have very complete descriptions of its interfaces, even the most basic ones, which I'm guessing is partly intentional ("let's not become another Java") and partly something that will get resolved as Julia becomes more mature as a language.

So, making some or all of the concrete Rotation types subtypes of StaticMatrix{3, 3} is a judgment call, and above I shared some thoughts on what criteria I think should go into the decision: a deliberate consideration of what interface a StaticMatrix is supposed to satisfy (it's lightly documented here but probably there's room for our judgment too) and whether a user (not developer!) of a given concrete type such as AngleAxis would expect that type to satisfy that interface. For example, as a user, it would not be immediately obvious to me what it means to get the ith entry of an AngleAxis -- maybe it gets the angle or one of the coordinates of the axis? -- whereas it would be obvious to me what it means to get the ith entry of a RotMatrix.

[bzinberg]

Also, while I agree that the reasoning you give here is self-consistent:

I think the == operator here is for "equals as a matrix" because UnitQuaternion <: Rotation <: StaticMatrix{3,3,<:Real}. So, it's ok to have different values of the field: q == -q but q.w != (-q).w.

it is a pretty big red flag to me that the rest of the world already has a notion of what it means for unit quaternions to be equal, and the above described notion of equality differs substantially from that one, not just because of rounding errors but conceptually. One way to improve this i.m.o. would be to rename UnitQuaternion to something like QuatRotation. That way it's clear that quaternions are involved somewhere, but there isn't a strong signal that this type should directly represent the same concept as "unit quaternion" as a mathematical object (which it doesn't).

I don't think the above is a fire or anything, but I do think it would be beneficial to make a clarifying change such as the above, on some appropriate timeframe. And we can definitely give a generous deprecation timeline if we do decide to make breaking changes.

[c42f]

One way to improve this i.m.o. would be to rename UnitQuaternion to something like QuatRotation. That way it's clear that quaternions are involved somewhere

This is the main problem you're having with this package — a matter of terminology.

The API for this package has long taken the stance that "rotations are matrices which happen to be efficiently parameterized".

So the previous implementation (called Rotations.Quat) didn't intend to claim to be an quaternion algebraically.

A while back we were lucky to have many improvements merged in from @bjack205's package DifferentialRotations.jl.

At that point Quat became UnitQuaternion, and a few conceptual inconsistencies seem to have crept in at that time. For example, I think conj(::Quat) wasn't defined, but conj(::UnitQuaternion) now means quaternion conjugation which is quite inconsistent with the algebraic interpretation as a matrix. This should be fixed one way or another.

Overall I think the idea that "Rotations are just specially-parameterized matrices" could be questioned, but it's worth keeping the history of the package in mind — it's always had a practical focus on simple 3D geometry and optimization. If wholesale API changes made it more practical for those uses while improving the conceptual clarity, they could be worthwhile.

But if the perceived problems with UnitQuaternion could be solved with some documentation and a few careful renamings/deprecations that would be a lot better for the packages which depend on this one.

[bzinberg]

Thanks for this input @c42f, all of this sounds great.

if the perceived problems with UnitQuaternion could be solved with some documentation and a few careful renamings/deprecations that would be a lot better for the packages which depend on this one.

I think we should do two things:

• Documentation and a few careful renamings/deprecations
• Make a new type for infinitesimal rotations. (If UnitQuaternions are treated opaquely as matrices and accessing their components is not part of the public API, then the question of "which exponential map" goes away.)

Infinitesimal rotations are not rotations (skew-symmetric matrices are not orthogonal matrices), and rotations don't have a zero. And the thing returned by log(::Rotation) could have a better semantics if this type existed. The new type could have an implicit conversion to StaticMatrix{3, 3}, and exp(::Rotation) could be deprecated rather than removed.

One use case I have in mind by the way, is machine learning on orientation-valued data where the loss function includes geometric operations on the orientation. There, unit quaternions have some advantages over other representations, particularly regarding singularities such as gimbal lock. I don't see any problems with the current API for that -- the user can just keep track of the quaternion components themself and feed them into Rotations.jl when they need to do geometric operations -- but just wanted to bring that use case into the mix.

[andyferris]

Yeah, that rotations are just matrices is kinda fundamental here. It sounds like the current conj on quaternions is a bug. If you want to do actual quaternion algebra, being able to go between any 3d rotation and some kind of actual quaternion object, and back again, would be good (not via convert itself, though, which has it's own semantics). Obviously the conversion would be cheap (probably free) to and from UnitQuaternion.

I always hoped to add AntiSymmetric and AntiHermition (or skew, if you prefer) to LinearAlgebra, follow through that into StaticArrays, and define log and exp here. But even if I had the time to make the effort, it's a long burn. We could just do all here, for now.

[hyrodium]

One way to improve this i.m.o. would be to rename UnitQuaternion to something like QuatRotation.

I totally agree with that. Or maybe VersorRotation?

I just found out the word "versor" today. (same as unit quaternion)
https://en.wikipedia.org/wiki/Versor

[andyferris]

I had never heard of versor before - thanks.

It’s an apt name (I love the word!), but probably less well-known to other users so might hurt in discoverability…

[c42f]

Personally I'd go with QuatRotation over VersorRotation as it's unambiguous and more well known.

There, unit quaternions have some advantages over other representations, particularly regarding singularities such as gimbal lock. I don't see any problems with the current API for that -- the user can just keep track of the quaternion components themself and feed them into Rotations.jl when they need to do geometric operations -- but just wanted to bring that use case into the mix.

I think this makes sense — in this case, it's the derivative of the constructor with respect to the parameters which we care about? So calling rotation constructors as part of the model evaluation seems important.

Make a new type for infinitesimal rotations.

I've long thought that this package should have integration with ChainRulesCore.jl. I assume this would tie in naturally there?