Do we really need normalized-flag in Quaternion?

[hyrodium]

The problems

The following are what I am concerned about.

• Someone just needs quaternion operations, without the flag.
• Performance disadvantage
• The flag is not well-documented, and some operations don’t work properly with the flag.

Someone just needs quaternion operations, without the flag.

See comments in Makie.jl for example.

Performance disadvantage

The following MyQuaternion is a quaternion without the flag. The benchmark shows -5.49% => improvement (5.00% tolerance).

julia> using Quaternions, BenchmarkTools

julia> struct MyQuaternion{T<:Real} <: Number
           s::T
           v1::T
           v2::T
           v3::T
       end

julia> Base.:*(q::MyQuaternion, w::MyQuaternion) = MyQuaternion(q.s * w.s - q.v1 * w.v1 - q.v2 * w.v2 - q.v3 * w.v3,
                                                      q.s * w.v1 + q.v1 * w.s + q.v2 * w.v3 - q.v3 * w.v2,
                                                      q.s * w.v2 - q.v1 * w.v3 + q.v2 * w.s + q.v3 * w.v1,
                                                      q.s * w.v3 + q.v1 * w.v2 - q.v2 * w.v1 + q.v3 * w.s)

julia> mq = MyQuaternion(1.0, 2.0, 3.0, 4.0)
MyQuaternion{Float64}(1.0, 2.0, 3.0, 4.0)

julia> q = Quaternion(1.0, 2.0, 3.0, 4.0)
Quaternion{Float64}(1.0, 2.0, 3.0, 4.0, false)

julia> a = @benchmark $mq*$mq
BenchmarkTools.Trial: 10000 samples with 1000 evaluations.
 Range (min … max):  2.793 ns … 31.009 ns  ┊ GC (min … max): 0.00% … 0.00%
 Time  (median):     4.260 ns              ┊ GC (median):    0.00%
 Time  (mean ± σ):   4.229 ns ±  0.486 ns  ┊ GC (mean ± σ):  0.00% ± 0.00%

                  ▁           ▃▃ ▆ ▄       ▆ █ █ ▅       ▂ ▁ ▂
  ▅▁▇▁▆▁▁▁▁▁▁▁▅▁█▁█▁▇▁▁▁▁▁▃▁▇▁██▁█▁█▁▁▁▁▁█▁█▁█▁█▁█▁▁▁▁▁█▁█▁█ █
  2.79 ns      Histogram: log(frequency) by time     4.75 ns <

 Memory estimate: 0 bytes, allocs estimate: 0.

julia> b = @benchmark $q*$q
BenchmarkTools.Trial: 10000 samples with 1000 evaluations.
 Range (min … max):  3.213 ns … 18.089 ns  ┊ GC (min … max): 0.00% … 0.00%
 Time  (median):     4.330 ns              ┊ GC (median):    0.00%
 Time  (mean ± σ):   4.474 ns ±  0.401 ns  ┊ GC (mean ± σ):  0.00% ± 0.00%

                                 ▅ █           ▂              
  ▂▁▂▁▂▁▂▁▁▁▁▁▁▁▂▁▂▁▂▁▂▁▁▁▁▁▂▁▄▄▁█▁█▁▅▁▁▁▁▁▃▁▇▁█▁▇▁▃▁▁▁▁▁▂▁▂ ▃
  3.21 ns        Histogram: frequency by time        5.17 ns <

 Memory estimate: 0 bytes, allocs estimate: 0.

julia> judge(mean(a),mean(b))
BenchmarkTools.TrialJudgement: 
  time:   -5.49% => improvement (5.00% tolerance)
  memory: +0.00% => invariant (1.00% tolerance)

The flag is not well-documented, and some operations don’t work properly with the flag.

See #36, #39, #51 and #59.

How to fix these problems

• If you need unit quaternions for rotations, I think Rotations.QuatRotation would be helpful.
  • Rotations.jl package have a dependency on Quaternions.jl, and conversion Quaterion <-> QuatRotation is supported.
• If you do not prefer the Rotations package, adding a new type Quaternions.UnitQuaternon seems better, I guess.
  • QuatRotaions is a rotation parametrization with Quaternion, but it is regarded as AbstractMatrix. i.e. the antipodal quaternions are equal. (Treating antipodal quaternions as equal breaks some operations Rotations.jl#171)

[sethaxen]

norm indicates that the quaternion/octonion is structurally normalized and allows the user to specify this constraint so that the resulting operations can be faster when possible. So a benchmark involving some unit quaternion-specific code would be I think more informative.

I'm surprised that simply the presence of the norm flag causes a decrease in performance. Do you have any idea why that might be?

I'd add to your concerns that there are several places in the code where type-stability is not guaranteed due to the use of norm. Of course, this could be repaired and tested without changing the user experience.

The choice of Rotations.jl to represent unit quaternions as a subtype of AbstractMatrix can have unintended consequences. In my experience, this plays poorly with ChainRules-based automatic differentiation engines, where in the pullbacks, matmuls are used instead of faster quaternion operations. So I found being able to use a Number type for unit quaternions to be safer, and then I could convert to a matrix representation when I needed one. One could make UnitQuaternion its own type; this makes the information about structural normalization known at compile time instead of runtime. But I expect this would increase the code complexity quite a bit.

What you seem to be proposing, dropping norm and potentially all unit quaternion/octonion-specific code, would be quite a breaking change. I'm pretty reluctant to make these kinds of changes without substantial agreement from the community.

[mikmoore]

What operations are actually made faster by the .norm field? So far the only functions that appear to take advantage of this are inv and the normalize family of functions, but I'm not sure these are really places where unitness is hard to determine or otherwise saves a lot. Checking numerically whether a quaternion is unit is very cheap

# should actually use a much tighter tolerance that depends on the type
#   5-10 ULPs should be plenty
@btime 0.99999<abs2(q)<1.00001 setup=(q=quatrand());
  2.100 ns (0 allocations: 0 bytes)

and could be even faster (and more accurate) if we used muladd to compute abs2(::Quaternion) using FMA instructions.

Meanwhile, the .norm flag grows the comfy 32-byte format Quaterion{Float64} into a 33-byte format that doesn't pack nicely in arrays (I'm guessing it pads to 36 bytes but haven't checked, but it probably messes up cache alignment in any case) and might negatively affect SIMD and other stuff.

Many places fail to set the .norm field accurately in the first place. For example Quaternion(1.0)/1.0 (equivalently, sign(Quaternion(1.0))) incorrectly sets .norm=false and normalize(Quaterniion(NaN)) claims to be unit (one could argue that it is and the NaNness has just been propagated, but this has other effects like isfinite incorrectly returning true). There's a lot of code dedicated to unitness already, but it is wrong or leads to type instability in a few places and is a pain to maintain in any case. All of this would "just work" if we determined unitness dynamically. The user probably knows what results should be unit anyway and doesn't need a stray byte to tell them so.

Meanwhile, unitness isn't something we actually work to preserve even when we do have it.

julia> normalize(Quaternion(1,1,0,0))^1000 |> abs
0.999999999999889 # but we claim this is a unit quaternion

Chain enough operations and the .norm flags can lead to results that are inaccurate enough to cause problems. I'm not suggesting we jam every unit operation through normalize, but as far as I can tell we only use unitness to short-circuit some operations (sometimes incorrectly, as in isfinite(Quaternion(NaN))) rather than actually take unitness seriously.

---

Fixing this doesn't have to be breaking, I think. We can overload getproperty to make the .norm "field" still accessible, we just re-route it to a call to a separate isunit function to decide the unit-ness numerically.
Quaternion(...,isnorm::Bool) can still exist (deprecated) and just ignore the flag. If we really wanted (not that I do want to), we could even use the isnorm=true case to do the user the favor of calling sign to be sure it's normalized or it could check that it truly is (within tolerance) unit or throw an error.

[sethaxen]

@mikmoore these arguments for dropping .norm are quite convincing.

Fixing this doesn't have to be breaking, I think. We can overload getproperty to make the .norm "field" still accessible, we just re-route it to a call to a separate isunit function to decide the unit-ness numerically.
Quaternion(...,isnorm::Bool) can still exist (deprecated) and just ignore the flag. If we really wanted (not that I do want to), we could even use the isnorm=true case to do the user the favor of calling sign to be sure it's normalized or it could check that it truly is (within tolerance) unit or throw an error.

In short, I agree we could add isunit now and deprecate the constructors with norm. I also agree we could reroute getproperty(::Quaternion, :norm) to isunit with a deprecation warning. The ColPrac guide gives no guidance on whether changing an (undocumented) field to a property should be considered breaking, but folks on Slack didn't seem to think it should be. So yes, this could all be non-breaking.

[sethaxen]

I started work on this in #75. @mikmoore, I'd appreciate it if you had any suggestions for what a good type-dependant approximate check for normalization would be.

[mikmoore]

Good question. I'm no expert on this, but I did some rough calculations and estimated that abs2(normalize(q)) should differ from 1 by at most something like 4eps(T) for a floating type T and so abs(normalize(q)) should be limited within 2eps(T) of error.

So a freshly normalized quaternion should safely satisfy abs(abs(q) - 1) <= 4eps(T). However, we don't want normalization to break after multiplying just a few unit quaternions either, and that could suffer some error propagation. I also reckoned that quaternion multiplication could give a non-unit result by a similar amount to abs or maybe a little higher (probably less than a factor of 2).

I don't really trust myself to have done the above analysis correctly, so I took a more empirical approach.

A numerical evaluation over 100e6 randn normalized quaternions never saw abs off by more than 1eps() or abs2 by more than 3eps(). E.g., (maximum(_->abs2(normalize(nquatrand())),1:100e6) -1). A test of 1e9 products of 2 random normalized quaternions returned a maximum deviation of 3eps(). 1e6 tests of products of 1000 normalized quaternions yielded a maximum deviation around 70eps() and 1e6 products of 10000 yielded around 180eps() for both Float32 and Float64 representations. These errors would grow much faster if every new argument was not normalized, as in power-by-squaring situations for example.

So for freshly-normalized quaternions, 10eps(T) should be plenty of margin for determining unitness. Products of thousands of normalized quaternions should still yield an error well-below, say 500eps(T). This might sound like a lot, but 500eps(Float64) is only 1.11e-13. and 500eps(Float32) is 5.96e-5. At some point, a person performing this many operations will have to manually renormalize their result, but they probably know this and will hopefully be responsible enough to do it more often than every few-thousand quaternion flops.

So I would propose a tolerance in the range of 50eps(T) to 1000eps(T) is reasonable, depending on how much margin we want to give people to go a bit without renormalizing their result.

That said, if we don't use isunit (or whatever we name it) for any computational shortcuts within the module itself (and I'd argue we shouldn't unless we can show a significant performance benefit), then the tolerance isn't really our concern. We could set a tight default tolerance (eg 20eps(T)) and provide a isunit(::Quaternion,tol) method where a tolerance can be manually specified by the user. Then isunit would basically serve as an analog of Base.isapprox.

Obviously, for Integer (and probably Rational) quaternions the tolerance should be zero. Since I think we've covered all the Base real types at this point, I guess the safe thing is to use zero tolerance for everything that isn't an AbstractFloat (or at least that doesn't have an eps() method).