Taking Quaternions Seriously

Specific Domains Numerics
21

I am doing something a little bit different. Instead of storing the quaternion elements inside a SVector, I created this structure:

struct Quaternion{T} <: AbstractVector{T}
    q0::T
    q1::T
    q2::T
    q3::T
end

Then I defined all the possible operations. It is really fast.

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

EDIT: However, it is very focused on satellite attitude representation. Can you explain with more details what you are trying to do? Now I am thinking that maybe this book will not be helpful depending on the case.

4 Likes
22

Thank you @Ronis_BR, it seems that chapter 2 is the core chapter about rotations, correct?

I don’t have very specific features in mind regarding these rotations yet, just investigating the possible approaches in other systems and fields of science.

3 Likes
23

Good! Well, you can also see my class notes:

https://www.ronanarraes.com/disciplinas/2020-01-cmc-202-4-movimento-de-um-solido/

Or if you understand Portuguese (which I know you do :slight_smile:), there is this nice document written by researchers of INPE:

http://mtc-m16d.sid.inpe.br/col/sid.inpe.br/mtc-m19/2012/01.26.19.13/doc/publicacao.pdf

2 Likes
24

Lovely :heart: :slight_smile:

1 Like
25

An SVector storage should be exactly equivalent, but with much less code.

5 Likes
26

It seems that there are different conceptions of quaternions and different implementations, appropriate for different purposes. Rather than trying to find one Version everyone is happy with, wouldn’t an approach like the Tables.jl interface work, where one meta package defines all operations on quaternions that everyone should support, and the implementation is left to each specific package? At least that way the different implementations might still mostly work across packages

8 Likes
27

IMHO this is the best thing to do. There are things that is not a consensus. For example, the real part is the first or fourth element? In this case, this package can abstract it.

2 Likes
28

100% agree with @jules.
My storage requirement for OpenGL + Makie pretty much exclude the performance/correctness requirements of other packages. But there shouldn’t be any reason, that Makie doesn’t work with any other Quaternion implementation out of the box.

5 Likes
29

It’s always been my desire to have Makie accept Grassmann.Chain{V,1} elements as replacements for the GeometryBasics.Point type, but converting the elements hasn’t been a big issue.

A quaternion is a combination of scalar and bivector, if you exclude the scalar, it is pure bivector. The pure bivector component is Grassmann.Chain{V,2} type, which is yet again the same type used for vectors, but with a different type parameter for the grade. So it would be called a 2-vector as opposed to a 1-vector.

However, once you add the scalar to a bivector, the element is no longer a bivector and becomes instead a general Grassmann.MultiVector, which has 2^3 elements in it. Currently, you need to rely on this representation for a full quaternion.

Alternatively, there is a MultiGrade type I have prototype for, with which I plan to make the scalar + bivector representation more efficient (this specialization task has not been completed yet). Although it’s not finished yet, this additional representation should be more efficient than the current usage of a full MultiVector type for quaternions.

So in Grassmann.jl alone there are already going to be multiple representations of quaternions, depending on whether it includes the scalar value or not, and with the scalar value added to it there are going to be multiple choices of representation for quaternions.

This wouldn’t quite work with my package the way you describe it. Although the AbstractTensors.jl does do something kind of similar to what you describe for my packages, it does not define the operations themselves, it only defines common generic methods.

2 Likes
30

Geometric Algebra doesn’t use complex numbers ( or need them ). By definition it’s the Clifford algebras with real number fields.

For anyone who hasn’t seen it, this is what quaternion basis looks like in GA

julia> using Grassmann
[ Info: Precompiling Grassmann [4df31cd9-4c27-5bea-88d0-e6a7146666d8]

julia> basis"+++"
(⟨+++⟩, v, v₁, v₂, v₃, v₁₂, v₁₃, v₂₃, v₁₂₃)

julia> 𝑖 = v12
v₁₂

julia> 𝑗 = v23
v₂₃

julia> 𝑘 = v13
v₁₃

julia> 𝑖*𝑖
-1v

julia> 𝑖*𝑗
1v₁₃

julia> 𝑘*𝑘
-1v

julia> 𝑗*𝑘
1v₁₂

julia> 𝑖*𝑗*𝑘
-1v

julia> q = 1/sqrt(2) - 1/sqrt(2)𝑖
0.7071067811865475 + 0.7071067811865475v₁₂

julia> q*1v2*~q
0.0 + 0.9999999999999998v₁

I think I got the signs right. If not, flip some stuff around :slight_smile:

3 Likes
31

As a convenience feature, you can use this command to get the correctly oriented quaternion basis

julia> i,j,k = hyperplanes(ℝ3) # using Grassmann
3-element Array{Simplex{⟨+++⟩,2,B,Int64} where B,1}:
 -1v₂₃
  1v₁₃
 -1v₁₂
3 Likes
32

I think I did not explain what I was after. I should clarify that I am not an expert on geometric algebra, but I am captivated by its generality and uniformity.

Geometric algebra subsumes complex numbers, quaternions, and trickier concepts like the distinction between a vector and a so-called “pseudovector”.

I believe it’s useful and necessary to regard of the state of affairs for doing spatial mathematics as a “legacy system”. People do spatial math using quaternions, cross products, triple products, and so on. There is an existing language, and from a certain perspective GA is a better language since the structure of all of these operations and objects exist within non-projective 2^3-dimensional geometric algebra (I guess notated G_{3,0,0}). There are uniform and systematic names for complex numbers, dual numbers, quaternions, etc.

But what about migration? A person can have never heard of geometric algebra, but still build a physics simulator using the language they already know. Pedagogically, I think it’s probably the right strategy to teach someone about complex numbers before you teach them about geometric algebra. The former is just a special case of the latter, and special cases matter.

So what’s the migration path? I have a function from a library that is not GA-aware and it operates with Base.Complex. There is another library that is “GA-aware” and because of that, more general. How do I make them interoperate?

I respect, and might even personally agree with

This is going to be my default algebra for everything, and I will rarely if ever use other ones.

but I don’t think it addresses “taking quaternions seriously”. I think Jiahao’s talk on Taking Vector Transposes Seriously was truly amazing for me to understand the spirit behind “taking [math thing] seriously”. The effort is strongly informed by history, practice, notation, and other programming languages to develop the taxonomy of “type 1” and “type 2” systems. To further the analogy, there are systems which as far as I can tell completely subsume the discussion on transposes in linear algebra. But it’s important to make this special case work well.

To be more concrete, if geometric algebra is the solution, I want something like the following to work. It’s likely that it shouldn’t work as I wrote it, but instead with the introduction of some shim type that wraps a certain kind of Multivector in a type so that it behaves like Complex.

(On Grassmann a79bff6d, the head of the default branch at time of writing)

julia> using Grassmann

julia> basis"-" # make something that behaves like complex numbers
(⟨-⟩, v, v₁)

julia> im′ = v1
v₁

julia> 2 + 3im′
2 + 3v₁

julia> (2 + 3im′) * (5 + 7im′)
-11 + 29v₁

julia> real(2 + 3im′)
2 + 3v₁

julia> imag(2 + 3im′)
ERROR: BoundsError: attempt to access (0, 1)
  at index [3]
Stacktrace:
 [1] getindex(::Tuple, ::Int64) at ./tuple.jl:24
 [2] getindex at /Users/goretkin/.julia/packages/StaticArrays/l7lu2/src/SVector.jl:39 [inlined]
 [3] imag(::MultiVector{⟨-⟩,Int64,2}) at /private/tmp/dev/Grassmann/src/parity.jl:246
 [4] top-level scope at REPL[15]:1
 [5] eval(::Module, ::Any) at ./boot.jl:331
 [6] eval_user_input(::Any, ::REPL.REPLBackend) at /Users/julia/buildbot/worker/package_macos64/build/usr/share/julia/stdlib/v1.4/REPL/src/REPL.jl:86
 [7] run_backend(::REPL.REPLBackend) at /Users/goretkin/.julia/packages/Revise/BqeJF/src/Revise.jl:1184
 [8] top-level scope at none:0

julia> (2 + 3im′) * (5 + 7im)
10 + 14im + (15 + 21im)v₁
5 Likes
33

Good question and relevant to Taking Quaternions Seriously.

Geometric Algebra has this meta-layer where you choose your algebra based on what you want to do. Usually you choose the smallest/simplest algebra you need. There are different ways to implement Quaternions in GA depending on what algebra you choose. In GA you would actually prefer to use a Rotor rather than a Quaternion ( multiply it all by the pseudoscalar ).

The best that can be done is a well thought out method interface as someone above mentioned. Complex has an interface I was able to quickly implement for your G(0,1,0) algebra of complex numbers.

I’m not a Grassmann expert or a GA expert, and I threw this together quickly, so don’t take it as the right way of doing it. But because both MultiVectors and Complex are at least in a duck-typing sense <:Number. I didn’t check, I believe they both are <:Number. Point is I didn’t have to check. I had a simple interface I could implement regardless of type hierarchy.

Now ideally if I implemented a Quaternion ( through methods ) using MultiVectors, I could then use fancy packages e.g. DifferentialEquations, JuMP with them and they would just work.

julia> Base.imag(M::MT) where {T,S, MT<:MultiVector{S,T,2}} = (-M⋅1.0v1)v1

julia> imag(2 + 3im′)
-0.0 + 3.0v₁

julia> ans*ans
-9.0v⃖

julia> (2 + 3im′)*(2 +3im)
4 + 6im + (6 + 9im)v₁

julia> Base.:*(M::MT, c::C) where {C<:Complex,S,T,MT<:MultiVector{S,T,2}} = M*(real(c)+imag(c)v1)

julia> (2 + 3im′)*(2 +3im)
-5 + 12v₁

julia> (2 +3im)*(2 +3im)
-5 + 12im

As another Quaternion use case: I’ve found use for Imaginary Quaternions ( no scalar part ) for conformal geometry processing using spin transforms. So quaternions are not just for rotating stuff in 3D.

< deleted some ramblings about Tensors and k-forms being possibly an elegant unification path >

1 Like
34

In that algebra generated by S"-", the appriopriate selectors would be scalar and vector. The imag definition is specifically defined in my paper, and is probably not appropriate here. Since it’s a 1 dimensional algebra, you can simply use the scalar and vector selector methods instead. For quaternions there is also a bivector selector method.

2 Likes
35

Offtopic: what kind of setup have you used to take those notes? I just realised you can even copy&paste text :wink:

2 Likes
36

I am using the app Notability in an iPad. It is awesome!

1 Like
37

Characterizing Quaternion Components

Hamilton identified a real number with a real quaternion:
a = [a, 0]. Nothing wrong here, but it invited Hamilton to go on and identify a pure quaternion (a quaternion with a null scalar part) with a vector: A = [0,A]. As the inventor of the vector he was entitled to call this object anything he wanted but the problem is that by this time people were already thinking about forces and such like objects very much as we think of vectors today and that the identification of Hamilton’s ‘vectors’ with what they had in mind created a great deal of confusion. [1]

Hamilton … professed that the quaternions make the study of vectors in three-space unnecessary since every vector can be considered as the vectorial part … of a quaternion … this interpretation is grossly incorrect since the vectorial part of a quaternion behaves with respect to coordinate transformations like a bivector or axial vector and not like an ordinary or polar vector. [2]

So let’s not further this terminological muddle.

[1] Simon L. Altmann (1989)
Hamilton, Rodrigues, and the Quaternion Scandal,
Mathematics Magazine, 62:5, 291-308, DOI: 10.1080/0025570X.1989.11977459
pp. 297-298

[2] M. Riesz, Clifford Numbers and Spinors. Lecture Series No 38,
The Institute for Fluid Dynamics and Applied Mathematics,
University of Maryland, Maryland, 1958.  p 21

Implementation Considerations

  • The Real abstract type presupposes commutative multiplication. Quaternion multiplication does not commute. Defining Quaternions as a subtype of Real lets them work incorrectly with other people’s code.

  • The Real abstract type does not supertype Complex numbers and, in general, does not supertype Dual numbers. Defining Quaternions with fields that are subtypes of Real excludes more general application of the type.

  • Performance is of substantial import. Any implementation should be as performant as one specialized to Float64 and Float32 fields.

Implementation Constructions

After reviewing the Julia quaternion implementations, the comments on this thread so far, and a number of papers, it seems to me that this is appropriate (field names are not the focus here):

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.

1 Like
38

That’s an interesting propositon, will need to ponder it, because it can make a difference in Grassmann.jl

I agree with the rest and I do plan to create something like a specialized specific Quaternion type in Grassmann.jl

I’m not proposing that my implementation should be the default for everyone, but it will be for my purposes.

PS, the pure quaternion component is a bivector by name and not a vector, just to repeat. Also, that paper you referenced by Riesz is a classic.

Also, since I don’t get paid, I can’t provide you with a date for when I will finalize my implementation variant, I just work on it for fun when I feel like it. The amount of work required to fully achieve what I envision with Grassmann.jl is non-trivial, and I should be charging many thousands of dollars for my overall effort. Yet, I’ve been making it available for free, even though the users on this forum want to prevent me from making money and stop me from being able to get a job or internship, I am still contributing this work for all of you anyways.

2 Likes
39

Yes, bivectors (thought it best not to alter the quote).

Insight that you have with regard to this developing implementation, apart from the nature of Grassmann.jl (we all understand that is done on your own time and its goals are wider than Quaternions), is likely helpful. As more detail becomes consensus, there may be opportunities for performance with which you are familiar.

1 Like
40

Please don’t do this. The most performant Quaternion is going to be realizable via SIMD. 4 scalars.

Additionally, it’s often ok to let your quaternions drift a certain amount before re-normalizing. How much? That’s up to the code using the Quaternions.

5 Likes