Skew-symmetric matrices as generators for rotations

[andyferris]

See #29.

I think it would be great to have a skew-symmetric matrix type that can be exponentiated into a Rotation, and vice-versa via log. These would form a Lie algebra to complement the Lie group of rotations.

Similarly to the rotations, these could have simplified parameterizations that are similar to AngleAxis and RodriguesVec.

[bzinberg]

It would be great to have the type of "Lie algebra element" be distinct from "Lie group element." Currently, log of a UnitQuaternion returns an "un-normalized unit quaternion," but there's no spec of what that means other than "for efficiency reasons we don't check unit length, so the caller should check that themself." Log of a unit quaternion does not belong in that type, because it is actually correct and necessary to allow non-unit length.

I think the conflating of Lie group elements with Lie algebra elements is the root cause of #126.

Also, there should be some type for which Base.zero returns the zero element of the Lie algebra.

[tkoolen]

One option could be to use RotationVec as the Lie algebra element, since it can be identified with the skew-symmetric matrix. That's kind of the approach I took with twists in RigidBodyDynamics.jl.

[bzinberg]

To fix #128, I think we may also need a distinct representation for the Lie algebra of unit quaternions. Related to that, I have some use cases [1] that rely on distinguishing between the quaternions q and -q. Maybe this package wants to specifically disclaim responsibility for maintaining that distinction; AFAICS it doesn't take a clear stance: q == -q but then we support exponential where that difference matters.

[1] I'm working with generative models of continuously changing quaternions, so I think I need to rely on UnitQuaternions remembering their components (i.e. not switching from q to -q, as long as the user doesn't round-trip convert them between other formats of course). Also, in SLERP interpolation, switching one of the endpoints from q to -q changes the path (I think this corresponds to traveling along a geodesic the short way vs. the long way, but not sure).

[andyferris]

I'd be tempted to build a family of representations of real 3x3 skew-symmetric matrices, starting with the obvious matrix-like one, but then also generators for RotationVec and Quaternion and so on. We can define exp and log to map between the group and algebra elements.

For unit quaternions our representation has an aditional Z(2) "gauge" (is it basically isomorphic to O(3) x Z(2)?). I could possibly imagine an additional group element which is basically an identity matrix but maps q to "-q" under multiplication (the quotes are because - is actually defined in terms of 3x3 matrices here...). I suppose there's an actual quaternion for that but I'm not sure if an explicit singleton or 1-bit struct would be good so each quaternion can be represented exactly as exp(quaternion_generator) * quaternion_guage. Would that be useful to anyone?

(We can't put that bit inside the quaternion generator because then addition becomes ambiguous).

[bzinberg]

Representing skew-symmetric matrices by their three free entries sounds great. For quaternions though, representation of the unit quaternion Lie algebra is even more straightforward than the so(3) Lie algebra: a tangent vector at the identity quaternion can be written as a pure-imaginary infinitesimal (a * ∂/∂x + b * ∂/∂y + c * ∂/∂z), and I think it is best to use this natural identification of the unit quaternion Lie algebra with R^3, rather than trying to build some other bijection (especially if it's discontinuous). Regarding the Lie groups, (if I understand the notation correctly) the unit quaternions are not isomorphic to SO(3) × (Z/2Z). The map from unit quaternions to SO(3) is a connected two-sheeted cover, so (very) schematically speaking it looks like this <https://en.wikipedia.org/wiki/Riemann_surface#/media/File:Riemann_sqrt.svg> rather than like this <https://en.wikipedia.org/wiki/Covering_space#/media/File:Covering_space_diagram.svg>. Any particular choice of how to locally decompose into a direct product would have singularities when one attempts to extend it to a global definition, thus losing the advantages quaternions have over other representations. So I think it's best to stick with representing unit quaternions directly as `(w, x, y, z)`.

(We can't put that bit inside the quaternion generator because then

addition becomes ambiguous). More ambiguous than it already is? The log function is already multi-valued. Actually I'm now realizing that even in the original formalism, a definition like q^t = exp(t * log(q)) depends on the branch cut when t isn't an integer, and I'm wondering what that means for the behavior of SLERP interpolation... [Edit: Ah, looks like changing the branch just causes SLERP to traverse the geodesic multiple times.]

…

On Fri, Jul 24, 2020 at 7:02 PM Andy Ferris ***@***.***> wrote: I'd be tempted to build a family of representations of real 3x3 skew-symmetric matrices, starting with the obvious matrix-like one, but then also generators for RotationVec and Quaternion and so on. We can define exp and log to map between the group and algebra elements. For unit quaternions our *representation* has an aditional Z(2) "gauge" (is it basically isomorphic to SU(2) x Z(2)?) I could possibly imagine an additional group element which is basically an identity matrix but maps q to "-q" under multiplication (the quotes are because - is actually defined in terms of 3x3 matrices here...). I suppose there's an actual quaternion for that but I'm not sure if an explicit singleton or 1-bit struct would be good so each quaternion can be represented exactly as exp(quaternion_generator) * quaternion_guage. Would that be useful to anyone? (We can't put that bit inside the quaternion generator because then addition becomes ambiguous). — You are receiving this because you commented. Reply to this email directly, view it on GitHub <#30 (comment)>, or unsubscribe <https://github.com/notifications/unsubscribe-auth/ABA7VV5GFHUEA5DOFDQHQXDR5IHH7ANCNFSM4DILZZ3Q> .

[hyrodium]

What would be an appropriate name for the abstract type to represent the Lie algebra (Skew-symmetric matrices)?
Maybe DiffRotation{d} <: StaticMatrix{d, d} or SkewSymmetric{d} <: StaticMatrix{d, d}?

[bzinberg]

I don't think DiffRotation (or InfinitesimalRotation or whatever we call it) should be subtype of StaticMatrix, because some operations that are naturally supported by matrices are not supported by elements of the Lie algebra (e.g., matrix multiplication: the product of skew-symmetric matrices is not necessarily skew-symmetric).

@hyrodium, are you planning to take a crack at implementing this?

[hyrodium]

some operations that are naturally supported by matrices are not supported by elements of the Lie algebra (e.g., matrix multiplication: the product of skew-symmetric matrices is not necessarily skew-symmetric).

Indeed, some operations are not algebraically closed in DiffRotation, but this property is same as Lie group Rotation <: StaticMatrix.
For example:

julia> using Rotations

julia> RotX(1) + RotY(2.3)
3×3 StaticArrays.SMatrix{3, 3, Float64, 9} with indices SOneTo(3)×SOneTo(3):
  0.333724  0.0        0.745705
  0.0       1.5403    -0.841471
 -0.745705  0.841471  -0.125974

Therefore, if the operation is not closed, I think we just need to return SMatrix.

are you planning to take a crack at implementing this?

Yes, I'm planning. But not sure the type namings and detailed implementation.

[bzinberg]

are you planning to take a crack at implementing this?

Yes, I'm planning. But not sure the type namings and detailed implementation.

I'd love to help!

[bzinberg]

Therefore, if the operation is not closed, I think we just need to return SMatrix.

This is a judgment call, but I think the implicit conversion to SMatrix introduces more potential for silent bugs (perhaps caused by the user needing to get clearer on the math they're trying to do) than it's worth -- I think it's better to have +(::Rotation, ::Rotation) return an error whose message includes "If you want to add matrices, you must explicitly cast to SMatrix." A lot of the Julia Base methods err on the side of safety in this way, and I think this is a good thing. We could even introduce a helper Rotations.to_matrix.

[bzinberg]

This also interacts with the issue described above in #30 (comment): UnitQuaternion being a subtype of Rotation is not quite right; in particular the conversion from UnitQuaternion to any other subtype of Rotation is a narrowing conversion. This leads to the strange situation that ::UnitQuaternion + ::RotMatrix{3}, if defined, has to narrow the first argument rather than widen the second argument.

And again, we might want to not implement infinitesimal unit quaternions in the first PR, instead sticking to just the other types which convert losslessly to elements of SO(3). One reason being that the exponential map is different for infinitesimal unit quaternions than it is for infinitesimal rotations. (The conversion is known and straightforward, but they are not the same return type.)

[hyrodium]

This is a judgment call, but I think the implicit conversion to SMatrix introduces more potential for silent bugs (perhaps caused by the user needing to get clearer on the math they're trying to do) than it's worth -- I think it's better to have +(::Rotation, ::Rotation) return an error whose message includes "If you want to add matrices, you must explicitly cast to SMatrix."

Hmm, that makes breaking changes, so I would like to keep Rotation <: AbstractMatrix for now.

A lot of the Julia Base methods err on the side of safety in this way

Do you have any examples? I think 1 / 2 isa Float64 is the simplest example that doesn't produce errors.

This leads to the strange situation that ::UnitQuaternion + ::RotMatrix{3}, if defined, has to narrow the first argument rather than widen the second argument.

I agree with that. But maybe, we can just return UnitQuaternion for ::UnitQuaternion + ::RotMatrix{3}.
(I'm also considering that UnitQuaternion(RotX(2π)) should be UnitQuaternion(-1,0,0,0) not UnitQuaternion(1,0,0,0). But this topic should be discussed in other issue.)

UnitQuaternion being a subtype of Rotation is not quite right; in particular the conversion from UnitQuaternion to any other subtype of Rotation is a narrowing conversion.

It's a small detail, but I think conversions from UnitQuaternion to MRP, RotX, RotXY, RotXYZ, AngleAxis and etc. are not narrowing.

And again, we might want to not implement infinitesimal unit quaternions in the first PR, instead sticking to just the other types which convert losslessly to elements of SO(3).

Okay, understood. I'll proceed that way:+1:

[hyrodium]

@andyferris
Do you have any thoughts about this?

[bzinberg]

Hmm, that makes breaking changes, so I would like to keep Rotation <: AbstractMatrix for now.

Agreed.

A lot of the Julia Base methods err on the side of safety in this way

Do you have any examples?

Main example that comes to mind (admittedly, it's the language itself, not Base methods 🙂) is not implicitly converting everything to a boolean:

julia> nothing || 5
ERROR: TypeError: non-boolean (Nothing) used in boolean context

julia> if [] 1 else 0 end
ERROR: TypeError: non-boolean (Vector{Any}) used in boolean context

Vaguely similar but not the same, checked casts:

julia> Int(3.3)
ERROR: InexactError: Int64(3.3)

It's a small detail, but I think conversions from UnitQuaternion to MRP, RotX, RotXY, RotXYZ, AngleAxis and etc. are not narrowing.

To me that is a big detail 🙂 I agree "narrowing" is not the right word, I'd say the conversion is co-widening. That is, a widening conversion is an injective function, whereas what we have here is a surjection, the 2-to-1 covering map. The non-injectivity of that map is why your above suggestion, "UnitQuaternion(RotX(2π)) should be UnitQuaternion(-1,0,0,0)" is currently impossible. Related is #128. But yeah, let's take that to a separate issue. (Update: filed #171.)

[hyrodium]

This leads to the strange situation that ::UnitQuaternion + ::RotMatrix{3}, if defined, has to narrow the first argument rather than widen the second argument.

I agree with that. But maybe, we can just return UnitQuaternion for ::UnitQuaternion + ::RotMatrix{3}.

Ah, I made a wrong reply here. I think it's ok to narrow UnitQuaternion because it is a subtype of AbstractMatrix. and the + result shuold be SMatrix.
I was thinking of ::UnitQuaternion * ::RotMatrix{3}, and at first, I thought narrowing the first unit quaternion can be evitable, but that was also wrong. Because we can't choose the mapping SO(3) → SU(2), it's not be able to widen RotMatrix{3} to UnitQuaternion properly.

[goretkin]

A lot of the Julia Base methods err on the side of safety in this way

Do you have any examples?

What about

julia> using Rotations

julia> Rotations.UnitQuaternion(1,0,0,0) + rand(3,3)
3×3 StaticArrays.SMatrix{3, 3, Float64, 9} with indices SOneTo(3)×SOneTo(3):
 1.72595   0.936594  0.125292
 0.197989  1.7558    0.827308
 0.258292  0.037621  1.90299

julia> Complex(1, 0) + rand(2, 2)
ERROR: MethodError: no method matching +(::Complex{Int64}, ::Matrix{Float64})
For element-wise addition, use broadcasting with dot syntax: scalar .+ array

A unit-magnitude complex (respectively quaternion) number has a representation as 2-by-2 (respectively 3-by-3) real orthogonal matrix.

Personally, I am in favor of a breaking release of this package that removes the behavior that conflates elements of the rotation group with their linear representation. Not to mention the double-cover issue.

[hyrodium]

A unit-magnitude complex (respectively quaternion) number has a representation as 2-by-2 (respectively 3-by-3) real orthogonal matrix.

It seems the example is not correct here, because Complex <: AbstractMatrix is false.
In this package, a 3D rotation is treated as a 3×3 matrix.

Personally, I am in favor of a breaking release of this package that removes the behavior that conflates elements of the rotation group with their linear representation. Not to mention the double-cover issue.

I think that will break the main concept of the package. Making a new package with a different name will be better, I guess.
Some of the problems will be fixed with #190.

[goretkin]

It seems the example is not correct here, because Complex <: AbstractMatrix is false.

That's exactly what I mean. A unit-complex/unit-quaternion is an efficient parameterization for a rotation in 2D/3D.

btw, this did introduce 2D Rotations type that <:AbstractMatrix, but parameterized by an angle, not a unit complex number.
#89

I think that will break the main concept of the package. Making a new package with a different name will be better, I guess.

It is a big part of Rotations.jl, true. I also value that the package provides many different methods for constructing rotations , and also types like RotX that encode information into the type system for the sake of multiple dispatch (in some cases? I have not checked recently.)