Uncertainty annihilation in quaternions using Rotations.jl

[togo59]

using Measurements, Rotations

I create a quaternion with uncertainty:

julia> q1 = Quat(1.0 ± 0.01,0.0 ± 0.01,0.0 ± 0.001,0.0 ± 0.002)
3×3 Quat{Measurement{Float64}}(1.0±0.0, 0.0±0.01, 0.0±0.001, 0.0±0.002):
 1.0±0.0    0.0±0.004  0.0±0.002
 0.0±0.004  1.0±0.0    0.0±0.02 
 0.0±0.002  0.0±0.02   1.0±0.0 

.. get the rotation matrix

julia> matrix_immutable = SMatrix{3,3}(q1)
3×3 SArray{Tuple{3,3},Measurement{Float64},2,9}:
 1.0±0.0    0.0±0.004  0.0±0.002
 0.0±0.004  1.0±0.0    0.0±0.02 
 0.0±0.002  0.0±0.02   1.0±0.0 

.. create a new quaternion from the above rotation matrix

julia> q2 = Quat(matrix_immutable)
3×3 Quat{Measurement{Float64}}(1.0±0.0, 0.0±0.01, 0.0±0.001, 0.0±0.002):
 1.0±0.0    0.0±0.004  0.0±0.002
 0.0±0.004  1.0±0.0    0.0±0.02 
 0.0±0.002  0.0±0.02   1.0±0.0 

get q1 * inv(q2)

julia> q1 / q2
3×3 Quat{Measurement{Float64}}(1.0±0.0, 0.0±0.0, 0.0±0.0, 0.0±0.0):
 1.0±0.0  0.0±0.0  0.0±0.0
 0.0±0.0  1.0±0.0  0.0±0.0
 0.0±0.0  0.0±0.0  1.0±0.0

get inv(q1) * q2

julia> q1 \ q2
3×3 Quat{Measurement{Float64}}(1.0±0.0, 0.0±0.0, 0.0±0.0, 0.0±0.0):
 1.0±0.0  0.0±0.0  0.0±0.0
 0.0±0.0  1.0±0.0  0.0±0.0
 0.0±0.0  0.0±0.0  1.0±0.0

We would obviously expect unity matrices without uncertainty but I am surprised that the uncertainty is 0.0. For instance:

julia> a = 3 ± 1
3.0 ± 1.0

julia> b = 6 ± 2
6.0 ± 2.0

julia> a/b
0.5 ± 0.23570226039551584

[giordano]

I'm not sure what you expect here. The parameter of the type is Measurement, which is a number with uncertainty. The uncertainty of all elements is zero, which should be correct, right? Did you expect that the parameter automatically changed to Float64?

[togo59]

No. I don't understand why dividing one uncertain quaternion type (full of Float64s) by another annihilates the uncertainty whereas dividing one Float64 by another does not do this. The latter is expected. Perhaps it's a question of whether Rotations.jl can actually be used with Measurements.lj

[tkoolen]

I think a better comparison is:

julia> using Measurements

julia> a = 3 ± 1
3.0 ± 1.0

julia> a / a
1.0 ± 0.0

From the readme: "Functional correlation between variables is correctly handled, so x-x ≈ zero(x), x/x ≈ one(x), tan(x) ≈ sin(x)/cos(x), cis(x) ≈ exp(im*x), etc...".

[togo59]

Yes, because a and a are identically equal. I was under the impression that, as with floating point errors, that the propagation of uncertainty with (say) division would give rise to a result with greater uncertainty not less. Even if their values and uncertainties were numerically equal. E.g.

julia> using Measurements

julia> a = 1 ± 0.5
1.0 ± 0.5

julia> b= 1 ± 0.50000001
1.0 ± 0.50000001

julia> a/b
1.0 ± 0.7071067882576154

I've obviously misunderstood this.

[tkoolen]

Note that Measurements has a tracking mechanism; it doesn't matter that the values and derivatives are identically equal. This mechanism enforces the notion of independent vs. dependent variables, and is also in play with the Rotations example in the original post. So:

julia> q1 = Quat(1.0 ± 0.01,0.0 ± 0.01,0.0 ± 0.001,0.0 ± 0.002);

julia> q2 = Quat(1.0 ± 0.01,0.0 ± 0.01,0.0 ± 0.001,0.0 ± 0.002);

julia> q1 / q2
3×3 Quat{Measurement{Float64}}(1.0±0.0, 0.0±0.0141421, 0.0±0.00141421, 0.0±0.00282843):
 1.0±0.0         0.0±0.00565685  0.0±0.00282843
 0.0±0.00565685  1.0±0.0         0.0±0.0282843
 0.0±0.00282843  0.0±0.0282843   1.0±0.0

[togo59]

That's subtle and a potential gotcha (for me). My code produced q1===q2 (true), whereas your code produced q1===q2 (false). I thought that the way I had constructed my q2 would have guaranteed q1===q2 (false). I'll need to add checks to make sure my results don't get thwarted by the cleverness of the tracking system. Thank you for pointing this out!

[giordano]

When dealing with floating point numbers you never want to test for equality (==), let alone egality (===).

In addition to that, any result of mathematical operations involving Measurement objects returns a "derived quantity", which cannot be egal to any "direct measurement".

[togo59]

Of course! But note that when I created what I thought was a new variable called q2 out of q1 by way of a constructor that made a new matrix along the way (but actually didn't) , the result turned out to be the same object, ie yes q1===q2! Thus to be sure when I create a new object out of some operation, the first thing I will do, before any uncertainty calculation is check whether q1===q2, otherwise I have gone round in a circle and the result of a division will be the annihilation of the uncertainties, so yes I do need to use egality, as I demonstrated at the top of this thread.

And the whole point of using Measurements is the very fact that all engineering variables are uncertain not just the sensor values. I am building a model-based state-estimator and controller for a very real robotic system, weighing several tonnes and over 10m in length. A noob in Julia I may be, but I cut my programming teeth a very long time ago using FORTRAN IV; all my work is on real engineering systems. Noobs aren't necessarily pathologically stupid and inept. (Although, in my receding years, the latter frailties are beginning to take hold.) Fancy thinking I was actually creating a new quaternion!

[tkoolen]

Just to further clarify what's going on: in the usage section of the readme,

Measurement objects can be created with the two following constructors:

measurement(value, uncertainty)
value ± uncertainty

[ ...]

Every time you use one of the constructors above, you define a new independent measurement. Instead, when you perform mathematical operations involving Measurement objects you create a quantity that is not independent, but rather depends on really independent measurements.

Calling SMatrix{3,3}(q1) of course doesn't result in these constructors being called since StaticArrays and Rotations are Measurements-unaware, and as such this doesn't create new independent variables. You have to be very explicit when you want to create new independent variables. This is a very good thing, in my opinion.

[andyferris]

OK... this tagging system is really cool. :)

[giordano]

Noobs aren't necessarily pathologically stupid and inept.

@togo59 I never thought so and I'm sorry if I made you think otherwise. My point was that in Julia testing for egality should be done with much care, because you aren't just testing that two objects are superficially equal, but they're fundamentally the same thing, but that means to rely on the particular implementation of the data structure. I think that this shouldn't be done unless you know very well the data structure involved (and know it won't change in the future), because it may have, like in this case, some subtleties that can lead to unexpected results.